On a class of derivative Nonlinear Schrödinger-type equations in two spatial dimensions
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  in the context of nonlinear fiber optics. In contrast to the usual nonlinear Schrödinger equation, this new model incorporates the additional effects of self-steepening and partial off-axis variations of the group velocity of the laser pulse. We prove global in-time existence of the corresponding solution for various choices of parameters, extending earlier results of . In addition, we present a series of careful numerical simulations concerning the (in-)stability of nonlinear ground states and the possibility of finite-time blow-up.
Key words and phrases:nonlinear Schrödinger equation, derivative nonlinearity, ground states, finite-time blow-up, self-steepening, spectral resolution, Runge-Kutta algorithm
2000 Mathematics Subject Classification:65M70, 65L05, 35Q55.
- 1 Introduction
- 2 Ground states
- 3 Numerical method for the time-evolution
- 4 Global well-posedness with full off-axis variation
- 5 (In-)stability properties of ground states with full off-axis variation
- 6 Well-posedness results for the case with partial off-axis variation
- 7 Numerical studies for the case with partial off-axis variation
This work is devoted to the analysis and numerical simulations for the following class of nonlinear dispersive equations in two spatial dimensions:
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
Indeed, we shall mainly be concerned with (1.1) rewritten in its evolutionary form:
Here and in the following, , for any , is the non-local 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 well-known Benjamin-Bona-Mahoney equation for uni-directional 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  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.
The NLS is a canonical model for slowly modulated, self-focusing 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 -super-critical for . It is well known that in these regimes, solutions to (1.4) may not exist for all , due to the possibility of finite-time blow-up. 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 .
In the -critical case, there is a sharp dichotomy characterizing the possibility of this blow-up: 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 self-similar blow-up 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 self-steepening of the laser pulse in the direction . Secondly, the operator describes off-axis variations of the group velocity of the beam. The case is thereby referred to as full off-axis dependence, whereas for the model incorporates only a partial off-axis variation. Both of these effects become more pronounced for high beam intensities (see ) and both are expected to have a significant influence on the possibility of finite-time blow-up. 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 super-critical regimes for equation (1.3). At least formally, though, equation (1.3) admits the following conservation law
generalizing the usual mass conservation. In the case of full-off axis dependence, (1.5) yields an a-priori bound on the -norm of , ruling out the possibility of finite-time blow-up. However, the situation is more complicated in the case with only a partial off-axis variation.
The latter was studied analytically in the recent work , but only for the much simpler case without self-steepening, i.e., only for . It was rigorously shown that in this case, even a partial off-axis variation (mediated by with ) can arrest the blow-up 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 finite-time blow-up as soon as .
Furthermore, the current work aims to extend the analysis of  to situations with additional self-steepening, 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 well-posedness versus finite-time blow-up 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 one-dimensional derivative NLS is known to be completely integrable. Nevertheless, there has only very recently been a breakthrough in the proof of global in-time existence for this case, see [15, 16]. In contrast to that,  gives strong numerical indications for a self-similar finite-time blow-up in derivative NLS with . The blow-up thereby seems to be a result of the self-steepening effect in the density , which generically undergoes a time-evolution 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 in-time 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 two-dimensional derivative NLS given by (1.3). Except for its physical significance, this class of models also has the advantage that the inclusion of (partial) off-axis variations via are expected to have a strong regularizing effect on the solution, and thus allow for several stable situations without blow-up.
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 time-integration 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 off-axis 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 self-steepening. 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 time-periodic solutions to (1.1) given by
The function then solves the following nonlinear partial differential equation:
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 non-negative solution to (2.1), cf. [13, 35]. Ground states will be an important benchmark for our numerical simulations below.
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 self-steepening and off-axis variations. Indeed, in one spatial dimension, equation (1.1) simplifies to
Seeking the ground state solution in the form (2.1) thus yields the following ordinary differential equation:
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
while the phase is given a-posteriori through
where . In view of (2.6), this implies that the phase function is given by
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 self-steepening , the phase . Thus, and we find
For , this is the well-known ground state solution to (1.4) in one spatial dimension, cf. [13, 35]. We notice that the effect of adding the off-axis 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.
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
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.
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
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 . 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 two-dimensional 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 non-vanishing and , as follows:
Step 1: In the case without self-steepening , 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 off-axis variation, the solution is elongated in the -direction. In the case of full off-axis 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 self-steeping , 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 non-vanishing . In Fig. 4 we show on the left the ground state for , , and , when the action of is orthogonal to the self-steepening. When compared to the case with , the solution is seen to be elongated in the -direction. Next, we simulate when acts parallel to the self-steepening, 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 time-evolution
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
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 semi-classical NLS is considered. Driscoll’s composite Runge-Kutta (RK) method  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 in-time 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.
The evolutionary form of our model (1.3) is in many aspects similar to the well-known Davey-Stewartson (DS) system, which is a non-local NLS type equation in two spatial dimensions, cf. [9, 35]. In [21, 23, 25], the possibility of self-similar blow-up 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 time-steps for times . In this case, we know that the exact time-dependent 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 time-evolution and the one for the ground state which in itself is obtained with an accuracy of order . Thus, the time-evolution 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 wave-numbers 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 time-integration 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 time-integration 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 in-time 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 (self-similar) blow-up in finite time. This behavior can be reproduced in our simulations.
To do so, we take initial data of the form
and work on the numerical domain given by (2.9) with . We will use time-steps 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.
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
we again use time steps for . As can be seen in Fig. 6 on the right, there is numerical indication for finite-time blow-up. 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 self-similar blow-up 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.
We want to point out that there are certainly more sophisticated methods available to numerically study self-similar blow-up, 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 self-similar blow-up in NLS and KdV type models. As a result, all our numerical findings concerning finite-time blow-up 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 blow-up, 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. Time-dependent change of variables in the case with self-steepening
In the case of self-steepening, the ability to produce an accurate numerical time-integration 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
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 Runge-Kutta 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, two-dimensional derivative NLS of the following form
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 time-steps for . We also apply a Krasny filter , 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 self-focusing behavior within . Moreover, its is no longer positive due to phase modulations.
Surprisingly, however, there is no indication of a finite-time blow-up, 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 finite-time blow-up. 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 finite-time blow-up, 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 .
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 finite-time blow-up, see .
4. Global well-posedness with full off-axis variation
In this section we will analyze the Cauchy problem corresponding to (1.3) in the case of full off-axis 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 in-time existence result, we rewrite (1.3), using Duhamel’s formula
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 off-axis variation, does not allow for any Strichartz estimates, see . 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 off-axis variations).
Let , and . Then for any and any , there exists a unique global in-time solution to (1.3), depending continuously on the initial data. Moreover,
Let . We aim to show that is a contraction on the ball
To this end, let us shortly denote
where for , we write
Now, let and recall that is an isometry on . Using Minkowski’s inequality, followed by Cauchy-Schwarz yields
To bound the two terms on the right hand side, we first note that
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 in-time solution . Standard arguments (see, e.g., ) 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
For , this conservation law yields a uniform bound on the -norm of , since
We consequently can re-apply 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 time-reversible 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 , 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. ): Let for . We first rewrite Duhamel’s formula (4.1), using the continuity of the semigroup to propagate backwards in time
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 Cauchy-Schwarz we find that this quantity is indeed finite, since
Denoting for simplicity , we find after a lengthy computation (see  for more details) that the integral
We can express using the integral formulation (4.5) and write
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). ∎
This proof can easily be extended to the analogous model problem in arbitrary space dimensions , provided
for (no self-steepening),
for (with self-steepening).
For this is the usual -subcritical restriction.
5. (In-)stability properties of ground states with full off-axis variation
In this section, we shall perform numerical simulations to study the (in-)stability properties of nonlinear ground states in the case with self-steepening and of full-off axis variation . In view of Theorem 4.1, we know that there cannot be any strong instability of , i.e., instability due to finite-time blow-up. 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.