On a class of derivative Nonlinear Schrödingertype equations in two spatial dimensions
Abstract.
We present analytical results and numerical simulations for a class of nonlinear dispersive equations in two spatial dimensions. These equations are of (derivative) nonlinear Schrödinger type and have recently been obtained in [11] in the context of nonlinear fiber optics. In contrast to the usual nonlinear Schrödinger equation, this new model incorporates the additional effects of selfsteepening and partial offaxis variations of the group velocity of the laser pulse. We prove global intime existence of the corresponding solution for various choices of parameters, extending earlier results of [2]. In addition, we present a series of careful numerical simulations concerning the (in)stability of nonlinear ground states and the possibility of finitetime blowup.
Key words and phrases:
nonlinear Schrödinger equation, derivative nonlinearity, ground states, finitetime blowup, selfsteepening, spectral resolution, RungeKutta algorithm2000 Mathematics Subject Classification:
65M70, 65L05, 35Q55.Contents
 1 Introduction
 2 Ground states
 3 Numerical method for the timeevolution
 4 Global wellposedness with full offaxis variation
 5 (In)stability properties of ground states with full offaxis variation
 6 Wellposedness results for the case with partial offaxis variation
 7 Numerical studies for the case with partial offaxis variation
1. Introduction
This work is devoted to the analysis and numerical simulations for the following class of nonlinear dispersive equations in two spatial dimensions:
(1.1) 
where , is a given vector with , and is a parameter describing the strength of the nonlinearity. In addition, for , we denote by the following linear differential operator
(1.2) 
Indeed, we shall mainly be concerned with (1.1) rewritten in its evolutionary form:
(1.3) 
Here and in the following, , for any , is the nonlocal operator defined through multiplication in Fourier space using the symbol
where is the Fourier variable dual to . For this obviously yields a bounded operator . In addition, is seen to be uniformly elliptic provided .
The inclusion of implies that (1.1), or equivalently (1.3), shares a formal similarity with the wellknown BenjaminBonaMahoney equation for unidirectional shallow water waves [3, 4]. However, the physical context for (1.1) is rather different. Equations of the form (1.1) have recently been derived in [11] as an effective description for the propagation of high intensity laser beams. This was part of an effort to remedy some of the shortcomings of the classical (focusing) nonlinear Schrödinger equation (NLS), which is obtained from (1.1) when , i.e.
(1.4) 
The NLS is a canonical model for slowly modulated, selffocusing wave propagation in a weakly nonlinear dispersive medium. The choice of thereby corresponds to the physically most relevant case of a Kerr nonlinearity, cf. [13, 35]. Equation (1.4) is known to conserve, among other quantities, the total mass
A scaling consideration then indicates that (1.4) is critical for and supercritical for . It is well known that in these regimes, solutions to (1.4) may not exist for all , due to the possibility of finitetime blowup. The latter means that there exists a time , depending on the initial data , such that
In the physics literature this is referred to as optical collapse, see [13].
In the critical case, there is a sharp dichotomy characterizing the possibility of this blowup: Indeed, one can prove that the solution to (1.4) with exist for all , provided its total mass is below that of the nonlinear ground state, i.e., the least energy solution of the form
Solutions whose norm exceeds the norm of , however, will in general exhibit a selfsimilar blowup with a profile given by (up to symmetries), see [30, 31]. In turn, this also implies that stationary states of the form are strongly unstable. For more details on all this we refer the reader to [6, 13, 35] and references therein.
In comparison to (1.4), the new model (1.3) includes two additional physical effects. Firstly, there is an additional nonlinearity of derivative type which describes the possibility of selfsteepening of the laser pulse in the direction . Secondly, the operator describes offaxis variations of the group velocity of the beam. The case is thereby referred to as full offaxis dependence, whereas for the model incorporates only a partial offaxis variation. Both of these effects become more pronounced for high beam intensities (see [11]) and both are expected to have a significant influence on the possibility of finitetime blowup. In this context, it is important to note that (1.3) does not admit a simple scaling invariance analogous to (1.4). Hence, there is no clear indication of sub or supercritical regimes for equation (1.3). At least formally, though, equation (1.3) admits the following conservation law
(1.5) 
generalizing the usual mass conservation. In the case of fulloff axis dependence, (1.5) yields an apriori bound on the norm of , ruling out the possibility of finitetime blowup. However, the situation is more complicated in the case with only a partial offaxis variation.
The latter was studied analytically in the recent work [2], but only for the much simpler case without selfsteepening, i.e., only for . It was rigorously shown that in this case, even a partial offaxis variation (mediated by with ) can arrest the blowup for all . In particular, this allows for nonlinearities larger than the critical case, cf. Section 6 for more details. One motivation for the present work is to give numerical evidence for the fact that these results are indeed sharp, and that one can expect finitetime blowup as soon as .
Furthermore, the current work aims to extend the analysis of [2] to situations with additional selfsteepening, i.e., , and to provide further insight into the qualitative interplay between the latter effect and the one stemming from . From a mathematical point of view, the addition of a derivative nonlinearity makes the question of global wellposedness versus finitetime blowup much more involved. Derivative NLS and their corresponding ground states are usually studied in one spatial dimension only, see e.g. [1, 8, 12, 14, 27, 28, 36, 37] and references therein. For , the classical onedimensional derivative NLS is known to be completely integrable. Nevertheless, there has only very recently been a breakthrough in the proof of global intime existence for this case, see [15, 16]. In contrast to that, [28] gives strong numerical indications for a selfsimilar finitetime blowup in derivative NLS with . The blowup thereby seems to be a result of the selfsteepening effect in the density , which generically undergoes a timeevolution similar to a dispersive shock wave formation in Burgers’ equation. To our knowledge there is, however, no rigorous proof of this phenomenon is currently available.
In two and higher dimensions, even the local intime existence of solutions to derivative NLS type equations seems to be largely unknown, let alone any further qualitative properties of their solutions. In view of this, the present paper aims to shine some light on the specific variant of twodimensional derivative NLS given by (1.3). Except for its physical significance, this class of models also has the advantage that the inclusion of (partial) offaxis variations via are expected to have a strong regularizing effect on the solution, and thus allow for several stable situations without blowup.
The organization of our paper is then as follows:

Certain perturbations of such ground states will form the class of initial data considered in the numerical timeintegration of (1.3). The numerical algorithm used to perform the respective simulations is detailed in Section 3. In it, we also include several basic numerical tests which compare the new model (1.3) to the classical (derivative) NLS.

In the former case, the picture is much more complete, which allows us to perform a numerical study of the (in)stability properties of the corresponding nonlinear ground states in Section 5.

In the case with only partial offaxis variations, the problem of global existence is more complicated and one needs to distinguish between the cases where the action of is either parallel or orthogonal to the selfsteepening. Analytically, only the former case can be treated so far (see Section 6). Numerically, however, we shall present simulations for both of these cases in Section 7.
2. Ground states
In this section, we focus on timeperiodic solutions to (1.1) given by
(2.1) 
The function then solves the following nonlinear partial differential equation:
(2.2) 
subject to the requirement that as . A solution to (2.1) with minimal energy will be called a ground state. For the classical NLS, i.e. , the ground state is the unique, radial and nonnegative solution to (2.1), cf. [13, 35]. Ground states will be an important benchmark for our numerical simulations below.
Remark 2.1.
2.1. Explicit ground state solutions in 1D
In one spatial dimension, equation (2.2) allows for explicit formulas for the ground state, which will serve as a basic illustration for the combined effects of selfsteepening and offaxis variations. Indeed, in one spatial dimension, equation (1.1) simplifies to
(2.3) 
Seeking the ground state solution in the form (2.1) thus yields the following ordinary differential equation:
(2.4) 
To solve this equation, we shall use the polar representation for
and impose the requirement that Plugging this ansatz into (2.4) and isolating the real and imaginary part yields the following coupled system:
Multiplying the second equation by and integrating from to gives
Using the above, we infer that the amplitude solves
(2.5) 
while the phase is given aposteriori through
(2.6) 
After some lengthy computation, similar to what is done for the usual NLS, cf. [13], the solution to (2.5) can be written in the form
(2.7) 
where . In view of (2.6), this implies that the phase function is given by
(2.8) 
where we omitted a physically irrelevant constant in the phase (clearly, is only unique up to multiplication by a constant phase).
Note that in the case with no selfsteepening , the phase . Thus, and we find
For , this is the wellknown ground state solution to (1.4) in one spatial dimension, cf. [13, 35]. We notice that the effect of adding the offaxis dispersion () widens the profile, causing it to decay more slowly as as can be seen in Fig. 1 on the left. On the right of Fig. 1, it is shown that the maximum of the ground state decreases with but that the peak becomes more compressed.
Remark 2.2.
2.2. Numerical construction of ground states
In more than one spatial dimension, no explicit formula for is known. Instead, we shall numerically construct by following an approach similar to those in [22, 24]. Since we can expect to be rapidly decreasing, we use a Fourier spectral method and approximate
by a discrete Fourier transform which can be efficiently computed via the Fast Fourier Transform (FFT). In an abuse of notation, we shall in the following use the same symbols for the discrete and continuous Fourier transform. To apply FFTs, we will use a computational domain of the form
(2.9) 
and choose , sufficiently large so that the obtained Fourier coefficients of decrease to machine precision, roughly , which in practice is slightly larger due to unavoidable rounding errors.
Now, recall that for a solution of the form (2.1) to satisfy (1.1), the function needs to solve (2.2). In Fourier space, this equation takes the simple form
where
For , the solution can be chosen to be real, but this will no longer be true for . In the latter situation, we will decompose
and separate (2.2) into its real and imaginary part, yielding a coupled nonlinear system for . By using FFTs, this is equivalent to the following system for and :
Formally, the system can be written as where and solved via a Newton iteration. One thereby starts from an initial iterate and computes the th iterate via the well known formula
(2.10) 
where is the Jacobian of with respect to . Since our required numerical resolution makes it impossible to directly compute the action of the inverse Jacobian, we instead employ a Krylov subspace approach as in [34]. Numerical experiments show that when the initial iterate is sufficiently close to the final solution, we obtain the expected quadratic convergence of our scheme and reach a precision of order after only 4 to 8 iterations.
As a basic test case, we compute the ground state of the standard twodimensional focusing NLS with , using the initial iterate
on the computational domain (2.9) with . By choosing many Fourier modes, we have after seven iterations of (2.10) a residual smaller than . The obtained solution is given on the left of Fig. 2. As expected, the solution is radially symmetric.
The numerical ground state solution hereby obtained then be used as an initial iterate for the situation with nonvanishing and , as follows:
Step 1: In the case without selfsteepening , the iteration is straightforward even for relatively large values such as . It can be seen in the middle of Fig. 2, that the ground state for and is no longer radially symmetric. As an effect of the partial offaxis variation, the solution is elongated in the direction. In the case of full offaxis dependence, the ground state for the same value of can be seen in Fig. 2 on the right. The solution is again radially symmetric, but as expected less localized than the ground state of the classical standard NLS. This is consistent with the explicit formulas for found in the one dimensional case above.
Step 2: In the case with selfsteeping , smaller intermediate steps have to be used in the iterations: We increment , by first varying only in steps of , always using the last computed value for as an initial iterate for the slightly larger . The resulting solution can be seen in Fig. 3. Note that the imaginary part of is of the same order of magnitude as the real part.
Step 3: In order to combine both effects within the same model, we shall use the ground state obtained for and as an initial iterate for the case of nonvanishing . In Fig. 4 we show on the left the ground state for , , and , when the action of is orthogonal to the selfsteepening. When compared to the case with , the solution is seen to be elongated in the direction. Next, we simulate when acts parallel to the selfsteepening, that is when , , and . The result is shown in the middle of Fig. 4. In comparison to the former case, the imaginary part of the solution is essentially rotated clockwise by 90 degrees. The elongation effect in the direction is still visible but less pronounced.
Step 4: For the ground states become increasingly peaked, as is seen from the 1D picture in Figure 1. Hence, to construct ground states for higher nonlinear powers in 2D, we will consequently require more Fourier coefficients to effectively resolve these solutions. To this end, we work on the numerical domain (2.9) with and Fourier modes. We use the ground state obtained for as an initial iterate for the case and follow the same program as outlined above.
3. Numerical method for the timeevolution
3.1. A Fourier spectral method
In this section, we briefly describe the numerical algorithm used to integrate our model equation in its evolutionary form (1.3). After a Fourier transformation, this equation becomes
Approximating the above by a discrete Fourier transform (via FFT) on a computational domain given by (2.9), yields a finite dimensional system of ordinary differential equations, which formally reads
(3.1) 
Here is a linear, diagonal operator in Fourier space, and has a nonlinear and nonlocal dependence on . Since can be large, equation (3.1) belongs to a family of stiff ODEs, for which several efficient numerical schemes have been developed, cf. [17, 21] where the particular situation of semiclassical NLS is considered. Driscoll’s composite RungeKutta (RK) method [10] has proven to be particularly efficient and thus will also be applied in the present work. This method uses a stiffly stable third order RK method for the high wave numbers of and combines it with a standard explicit fourth order RK method for the low wave numbers of and the nonlinear part . Despite combining a third order and a fourth order method, this approach yields fourth order intime convergence in many applications. Moreover, it provides an explicit method with much larger time steps than allowed by the usual fourth order stability conditions in stiff regimes.
Remark 3.1.
The evolutionary form of our model (1.3) is in many aspects similar to the wellknown DaveyStewartson (DS) system, which is a nonlocal NLS type equation in two spatial dimensions, cf. [9, 35]. In [21, 23, 25], the possibility of selfsimilar blowup in DS is studied, using a numerical approach similar to ours.
As a first basic test of consistency, we apply our numerical code to the cubic NLS in 2D i.e. equation (1.4) with . As initial data we take the ground state , obtained numerically as outlined in Section 2 above. We use timesteps for times . In this case, we know that the exact timedependent solution is simply given by . Comparing this to the numerical solution obtained at yields an difference of the order of . This verifies both the code for the timeevolution and the one for the ground state which in itself is obtained with an accuracy of order . Thus, the timeevolution algorithm evolves the ground state with the same precision as with which it is known.
For general initial data , we shall control the accuracy of our code in two ways: On the one hand, the resolution in space is controlled via the decrease of the Fourier coefficients within (the finite approximation of) . The coefficients of the highest wavenumbers thereby indicate the order of magnitude of the numerical error made in approximating the function via a truncated Fourier series. On the other hand, the quality of the timeintegration is controlled via the conserved quantity defined in (1.5). Due to unavoidable numerical errors, the latter will numerically depend on time. For sufficient spatial resolution, the relative conservation of will overestimate the accuracy in the timeintegration by orders of magnitude.
3.2. Reproducing known results for the classical NLS
As already discussed in the introduction of this paper, the cubic NLS in two spatial dimensions is critical and its ground state solution is strongly unstable. Indeed, any perturbation of which lowers the norm of the solution below that of itself, is known to produce purely dispersive, global intime solutions which behave like the free time evolution for large . However, perturbations that increase the norm above that of are expected to generically produce a (selfsimilar) blowup in finite time. This behavior can be reproduced in our simulations.
To do so, we take initial data of the form
(3.2) 
and work on the numerical domain given by (2.9) with . We will use timesteps within . We can see on the right of Fig. 5 that the norm of the solution decreases monotonically, indicating purely dispersive behavior. The plotted absolute value of the solution at confirms this behavior. In addition, is conserved to better than , indicating that the problem is indeed well resolved in time.
Remark 3.2.
Note that we effectively run our simulations on , instead of . As a consequence, the periodicity will after some time induce radiation effects appearing on the opposite side of . The treatment of (large) times therefore requires a larger computational domain to suppress these unwanted effects.
Next, for initial data of the form
(3.3) 
we again use time steps for . As can be seen in Fig. 6 on the right, there is numerical indication for finitetime blowup. The code is stopped at when the relative error in the conservation of drops below . The solution for can be seen on the left of Fig. 6. This is in accordance with the selfsimilar blowup established by Merle and Raphaël, cf. [30, 31]. In particular, we note that the result does not change notably if a higher resolution in both and is used.
Remark 3.3.
We want to point out that there are certainly more sophisticated methods available to numerically study selfsimilar blowup, see for instance [24, 28, 35] for the case of NLS type models, as well as [19, 20] for the analogous problem in KdV type equations. However, these methods will not be useful for the present work, since, as noted before, the model (1.1) does not admit a simple scaling invariance, which is the underlying reason for selfsimilar blowup in NLS and KdV type models. As a result, all our numerical findings concerning finitetime blowup have to be taken with a grain of salt. An apparent divergence of certain norms of the solution or overflow errors produced by the code can indicate a blowup, but might also just indicate that one has run out of resolution. The results reported in this paper therefore need to be understood as being stated with respect to the given numerical resolution. However, we have checked that they remain stable under changes of the resolution within the accessible limits of the computers used to run the simulations.
3.3. Timedependent change of variables in the case with selfsteepening
In the case of selfsteepening, the ability to produce an accurate numerical timeintegration in the presence of a derivative nonlinearity () becomes slightly more complicated. The inclusion of such a nonlinearity can lead to localized initial data moving (relatively fast) in the direction chosen by . In turn, this might cause the numerical solution to “hit” the boundary of our computational domain .
To avoid this issue, we shall instead perform our numerical computations in a moving reference frame, chosen such that the maximum of remains fixed at the origin. More precisely, we consider the transformation
and denote . Under this transformation (1.3) becomes
(3.4) 
The quantity is then determined by the condition that the density has a maximum at for all . We get from (3.4) the following equation for :
Differentiating this equation with respect to and respectively, and setting yields the desired conditions for and .
Note that the computation of the additional derivatives appearing in this approach is expensive, since in practice it needs to be enforced in every step of the RungeKutta scheme. Hence, we shall restrict this approach solely to cases where the numerical results appear to be strongly affected by the boundary of . In addition, we may always choose a reference frame such that either one of the two components of is zero which consequently allows us to set either , or .
3.4. Basic numerical tests for a derivative NLS in 2D
As an example, we consider the case of a cubically nonlinear, twodimensional derivative NLS of the following form
(3.5) 
which is obtained from our general model (1.3) for , and . We take initial data given by (3.3). Here, is the ground state computed earlier for this particular choice of parameters, see Fig. 3. We work on the computational domain (2.9) with , using Fourier modes and timesteps for . We also apply a Krasny filter [26], which sets all Fourier coefficients smaller than equal to zero. The real and imaginary part of the solution at the final time can be seen in Fig. 7 below. Note that they are both much more localized and peaked when compared to Fig. 3, indicating a selffocusing behavior within . Moreover, its is no longer positive due to phase modulations.
Surprisingly, however, there is no indication of a finitetime blowup, in contrast to the analogous situation without derivative nonlinearity (recall Fig. 6 above). Indeed, the Fourier coefficients of at are seen in Fig. 8 to decrease to the order of the Krasny filter. In addition, the norm of the solution, plotted in the middle of the same figure, appears to exhibit a turning point shortly before . Finally, the velocity component plotted on the right in Fig. 8 seems to slowly converge to a some limiting value . The latter would indicate the appearance of a stable moving soliton, but it is difficult to decide such questions numerically. All of these numerical findings are obtained with conserved up to errors of the order .
It might seem extremely surprising that the addition of a derivative nonlinearity is able to suppress the appearance of finitetime blowup. Note however, that in all the examples above we have used only (a special case of) perturbed ground states as initial data. For more general initial data, the situation is radically different, as can be illustrated numerically in the following example: We solve (3.5) with purely Gaussian initial data of the form
on a numerical domain with , using Fourier coefficients and time steps for . This case appears to exhibit finitetime blowup, as is illustrated in Fig. 9. The conservation of the numerically computed quantity drops below at which indicates that plotting accuracy is no longer guaranteed. Consequently we ignore data taken for later times, but note that the code stops with an overflow error for .
Remark 3.4.
These numerical findings are consistent with analytical results for derivative NLS in one spatial dimension. For certain values of and certain velocities the corresponding ground states, or more general solitary wave solutions, are found to be orbitally stable, see [8, 14, 27]. However, for general initial data and large enough, one expects finitetime blowup, see [28].
4. Global wellposedness with full offaxis variation
In this section we will analyze the Cauchy problem corresponding to (1.3) in the case of full offaxis dependence, i.e. , so that
In this context, we expect the solution of (1.3) to be very well behaved due to the strong regularizing effect of the elliptic operator acting in both spatial directions.
To prove a global intime existence result, we rewrite (1.3), using Duhamel’s formula
(4.1) 
Here, and in the following, we denote by
the corresponding linear propagator, which is easily seen (via Plancherel’s theorem) to be an isometry on for any . It is known that in the case with full offaxis variation, does not allow for any Strichartz estimates, see [5]. However, the action of allows us to “gain” two derivatives and offset the action of the gradient term in the nonlinearity of (4.1). Using a fixed point argument, we can therefore prove the following result.
Theorem 4.1 (Full offaxis variations).
Let , and . Then for any and any , there exists a unique global intime solution to (1.3), depending continuously on the initial data. Moreover,
Proof.
Let . We aim to show that is a contraction on the ball
To this end, let us shortly denote
(4.2) 
where for , we write
Now, let and recall that is an isometry on . Using Minkowski’s inequality, followed by CauchySchwarz yields
To bound the two terms on the right hand side, we first note that
(4.3) 
In general spatial dimensions , the dual of Sobolev’s embedding then yields
provided , when . In addition, if we impose , then we have
again by Sobolev’s embedding. This allows us to also estimate
Together with a Hölder’s inequality in , we can consequently bound
By choosing sufficiently small, Banach’s fixed point theorem directly yields a unique local intime solution . Standard arguments (see, e.g., [33]) then allow us to extend this solution up to a maximal time of existence and we also continuous dependence on the initial data.
Next, we shall prove that
(4.4) 
For , this conservation law yields a uniform bound on the norm of , since
We consequently can reapply the fixed point argument as many times as we wish, thereby preserving the length of the maximal interval in each iteration, to yield . Since the equation is timereversible modulo complex conjugation, we obtain a global solution for all , provided (4.4) holds.
To prove (4.4), we adapt and (slightly) modify an elegant argument given in [32], which has the advantage that it does not require an approximation procedure via a sequence of sufficiently smooth solutions (as is classically done, see e.g. [6]): Let for . We first rewrite Duhamel’s formula (4.1), using the continuity of the semigroup to propagate backwards in time
(4.5) 
As is unitary in , we have . The latter can be expressed using the above identity:
We want to show that . In view of (4.2) we can rewrite
By CauchySchwarz we find that this quantity is indeed finite, since
Denoting for simplicity , we find after a lengthy computation (see [2] for more details) that the integral
We can express using the integral formulation (4.5) and write
(4.6) 
Next, we note that the particular form of our nonlinearity implies
Here, the first expression in the last line is obviously zero, whereas for the second term we compute
for solutions . In summary, the second term on the right hand side of (4.6) simply vanishes and we find
This finishes the proof of (4.4). ∎
Remark 4.2.
This proof can easily be extended to the analogous model problem in arbitrary space dimensions , provided

for (no selfsteepening),

for (with selfsteepening).
For this is the usual subcritical restriction.
5. (In)stability properties of ground states with full offaxis variation
In this section, we shall perform numerical simulations to study the (in)stability properties of nonlinear ground states in the case with selfsteepening and of fulloff axis variation . In view of Theorem 4.1, we know that there cannot be any strong instability of , i.e., instability due to finitetime blowup. Nevertheless, we shall see that there is a wealth of possible scenarios, depending on the precise choice of parameters, , , and on the way we perturb the initial data.