Supercritical blowup in coupled parity-time-symmetric nonlinear Schrödinger equations

# Supercritical blowup in coupled parity-time-symmetric nonlinear Schrödinger equations

João-Paulo Dias, Mário Figueira, Vladimir V. Konotop, and Dmitry A. Zezyulin
Avenida Professor Gama Pinto 2, Lisboa 1649-003, Portugal
Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa,
Avenida Professor Gama Pinto 2, Lisboa 1649-003, Portugal
Campo Grande, Ed. C8, Lisboa 1749-016, Portugal
###### Abstract

We prove finite time supercritical blowup in a parity-time-symmetric system of the two coupled nonlinear Schrödinger (NLS) equations. One of the equations contains gain and the other one contains dissipation such that strengths of the gain and dissipation are equal. We address two cases: in the first model all nonlinear coefficients (i.e. the ones describing self-action and non-linear coupling) correspond to attractive (focusing) nonlinearities, and in the second case the NLS equation with gain has attractive nonlinearity while the NLS equation with dissipation has repulsive (defocusing) nonlinearity and the nonlinear coupling is repulsive, as well. The proofs are based on the virial technique arguments. Several particular cases are also illustrated numerically.

Keywords: blowup, parity-time symmetry, nonlinear Schrödinger equation.

## 1 Introduction and preliminary observations

In this paper we consider the Cauchy problem for two coupled nonlinear Schrödinger (NLS) equations

 iut=−Δu+iγu+κv−(g1|u|2+g|v|2)u, (1a) ivt=−Δv−iγv+κu−(g|u|2+g2|v|2)v, (1b) where x∈RN, t≥0, Δ is the N-dimensional Laplacian, and u(x,0)=u0(x),v(x,0)=v0(x),u0,v0∈H1(RN),|x|u0(x),|x|v0(x)∈L2(RN). (1c)

We consider , and assume that and [the cases and (or) can be considered analogously after simple redefinition of the functions and ]. Our main interest is in establishing sufficient conditions for the initial conditions and to blow up in finite time.

The system (1a)–(1b) is referred to as parity-time () symmetric because of its linear counterpart (), which under certain conditions allows for stable propagating of the linear waves, which justifies the particular relevance of the model. The concept of symmetry can be formulated also in the nonlinear case referring to the formal property as follows: if functions and solve the pair of equations in (1a)–(1b), then the new functions and also solve the same equations, provided that (hereafter an overbar stands for the complex conjugation).

The concept of symmetry was originally employed in quantum mechanics to construct non-Hermitian potentials with pure real spectra [1, 2]. Soon after that, the ideas of symmetry were transferred to the optics [3]. Using the mathematical analogy between the Schrödinger equations in quantum mechanics and paraxial approximation in optics, it was shown that a medium with mutually balanced gain and losses allows for stable propagation of linear waves. Further, it was suggested in [4] that a simple implementation of symmetry in optical systems can be achieved using two waveguides. This geometry was proven to be very suitable for experimental observation of the effects related to the symmetry [5, 6]. Optical applications which naturally admit inclusion of the nonlinearity have stimulated further studies of the nonlinear -symmetric system [7]. In the two-waveguide geometry described by the coupled NLS equations with gain and dissipation, there have been considered bright [8, 9] and dark [10] solitons, breathers [11], and rogue waves [12].

Interplay between the symmetry and nonlinearity was also considered in the context of Bose–Einstein condensates (BECs) [13, 14], where the coupled NLS (alias Gross–Pitaevskii) equations appear as a natural model for two BECs split in an external double-well trap. Alternatively, the symmetry in BECs can be implemented using spinor components of a single atomic specie in two different ground states, one of which is pumped by an external laser and another one has (either natural or induced) loss of atoms [15]. Since a three-dimensional realization of a BEC appears in the most of experimental settings, it is of interest to consider dynamics of two coupled three-dimensional, and more generally -dimensional, NLS equations (1) one of which accounts for gain [equation (1a)] and the other one has dissipation [equation (1b)].

Turning now to the blowup phenomenon, in coupled conservative NLS equations it was previously discussed in Refs. [16, 17] without linear coupling and in Ref. [18] with linear coupling (see also [19] and references there in for the critical collapse). On the other hand, blowup in a single NLS equation with homogeneous and inhomogeneous linear damping was considered in [20, 21] and [22], respectively.

As mentioned above, particular interest in analyzing the blowup in the -symmetric coupled NLS equations (1a)–(1b) stems from the fact that depending on the relation among parameters the system can display features either of a conservative medium, where linear waves [i.e. solutions of the Eqs. (1) with ] can propagate stably or of a dissipative one where a linear mode is unstable (i.e. it either decays or infinitely grows with time). The former case is referred to as the unbroken -symmetric phase [1], and corresponds to

 κ>γ>0. (2)

If , then the symmetry is said to be broken. The particular case corresponds to the exceptional point [23] of the underlying linear operator.

Thanks to the possibility of stable propagation of linear modes, the case of unbroken phase is of special physical interest. It also allows one to make some general preliminary conclusions about blowing up and global existence of solutions. Indeed, if (2) is satisfied one can perform the transformation

 (uv)=(eiα/2−e−iα/2e−iα/2eiα/2)(UV), (3)

where the constant phase is determined by the relation

 eiα=−κ√κ2−γ2−iγ. (4)

We notice that (2) implies that . The newly introduced functions and solve the system

 iUt=−ΔU+ωU−F1,iVt=−ΔV−ωV−F2 (5)

with ,

 F1=G|U|2U+G+|V|2U+2Q|U|2V+M|V|2V+PU2¯V+(G−2g)V2¯U, F2=¯G|V|2V+G+|U|2V+2¯Q|V|2U+M|U|2U+¯PV2¯U+(¯G−2g)U2¯V,

and the other parameters being defined by

 G±=12(g1±g2),G=g+G+−iG−tanα, M=−G−cosα,Q=−M+igsinα,P=Mcos(2α)−2iG+sinα

and the initial conditions

 U(x,0)=U0(x),V(x,0)=V0(x),U0,V0∈H1(RN). (6)

Notice that the performed transformation (3) does not depend on the dimensionality of the space .

If we have and the corresponding system (5) has particular solutions which can be found in the form with solving the equation

 iUt=−ΔU+ωU−(g+g1)|U|2U. (7)

(analogously one can consider a particular solution with ). This leads to a number of conclusions (which follow from the well known results on the NLS equation, see. e.g. [24]) as follows.

First, if and , then there exist solutions blowing up in finite time. This collapse is characterized by the conserved squared norms and simultaneous blowup of both the fields and . Notice that hereafter we use the abbreviated notation for the standard -norm: .

Second, if then there exist global (dispersive) solutions.

Third, subject to sufficiently small initial conditions and (the critical case) one can find solutions existing globally. Moreover, in the case a solution with either or with smooth localized initial data exists globally, since it is described by the linear Schrödinger equation.

Below in this paper we concentrate on corresponding to the super-critical case.

## 2 Main results

In order to formulate our main results it is convenient to define the following integral Stokes components for the NLS equation:

 S0=∥u∥22+∥v∥22,S1=2Re∫u¯vdx, S2=2Im∫u¯vdx,S3=∥u∥22−∥v∥22.

Hereafter for the sake of brevity we use the notation . Obviously, is a conserved quantity for , but becomes time-dependent for . Indeed, it is straightforward to compute

 12dS0(t)dt=γS3(t), (8)

 S0(t)≤S0(0)e2γt. (9)

We also define the energy

 (10)

which is another conserved quantity in the conservative case (). For arbitrary the energy changes with time according to

 dEdt=2γ∫(|∇u|2−|∇v|2−g1|u|4+g2|v|4)dx. (11)

Next, we define the mean squared width of the solution

 X(t)=∫|x|2(|u|2+|v|2)dx, (12)

and its derivative :

 Y(t)≡dX(t)dt=4Im∫(ux⋅∇¯u+vx⋅∇¯v)dx+2γ∫|x|2(|u|2−|v|2)dx.

For the second derivative of one can compute

 d2X(t)dt2=4NE(t)+4∫(2−N)(|∇u|2+|∇v|2)dx+4γ2X+ +16γIm∫(¯ux⋅∇u−¯vx⋅∇v)dx+8γκIm∫|x|2v¯udx− −8κNRe∫u¯vdx. (13)

Now we introduce functions

 F(t) = X(0)+Y(0)t+8NN+2E(0)t2+4κγ2S0(0)(e2γt−2γt−1), (14) M(t) = supτ∈[0,t]F(τ)+1>0, (15) G(t) = M(t)(c1t22+exp(c3γtc2)−1), (16)

where is defined by -norms of the initial conditions and , as well as constants

 c1 = 4γκ+4γ25N+6N−2, c2 = ⎧⎪⎨⎪⎩45min{1,g1,g2},ifg≥0,45min{1,g1+g√g1√g2,g2+g√g2√g1},if−√g1g2

Our main result consists in the following Theorem which provides sufficient conditions for the finite-time blowup of solution for the problem (1) in the supercritical case.

###### Theorem 1

Let and

 g1,2>0,g>−√g1g2. (17)

Assume that the initial conditions and of the Cauchy problem (1) are chosen such that there exists for which the following two conditions hold:

 F(T0)+1<0, (18) G(T0)<1. (19)

Then the solution of the problem (1) does not exist in the interval .

Before presenting the proof of the Theorem 1, it is useful to ensure that the conditions (18) and (19) are consistent and that the initial conditions and satisfying (18)–(19) do exist. We illustrate the consistency of the conditions (18) and (19) formulating the two lemmas as follows. Lemma 1 establishes that the blow up occurs if the initial energy is negative (with large enough absolute value), while Lemma 2 shows that the blow up conditions can be satisfied if is negative (with large enough absolute value).

###### Lemma 1

Let ,

 C0=|Y(0)|2γ+4κγ2S0(0),~M(t)=1+X(0)+C0(e2γt−1),

and

 T0,max=1βln(1+β2(1+X(0))(β2+c1)), (20) T0,min=1βln(1+β2~M(T0,max)(β2+c1)). (21)

Then there exists such that the blow-up conditions (18)–(19) are satisfied at provided that

 E(0)<−(N+2)~M(T0,max)8NT20,min. (22)

Proof. Let us introduce

and define as the smallest solution of the equation (which implies that . Using that we find that . Introduce . For we have . With simple transformations we obtain that and therefore for all one has

 M(t)≤supτ∈[0,t]~F(τ)+1≤supτ∈[0,t]~M(τ)=~M(T0,max), (24)

which implies that .

In order to satisfy the condition (18) we require which is satisfied automatically if

 E(0)<−(N+2)~M(T0)8NT20. (25)

If the condition (22) holds then the requirement (25) also holds. Therefore, both conditions of the blow-up are satisfied at .

###### Lemma 2

Redefine the constant as follows:

 C0=4N|E(0)|(N+2)γ2+4κγ2S0(0), (26)

and keep others definitions the same as in Lemma 1. The conditions (18)–(19) are satisfied if

 Y(0)<8κγS0(0)−~M(T0,max)T0,min. (27)

Proof. The proof is almost identical to that for Lemma 1 except for definition of the function . Now it is defined as . Then for any .

## 3 Proof of Theorem 1

We start with the following estimate which is obtained from (11) and from the definition of the constant :

 E(t)≤2γ∫t0(∥∇u∥22+∥∇v∥22+g1∥u∥44+g2∥v∥44)dτ+E(0) ≤γc3N+216N∫t0(∥∇u∥22+∥∇v∥22+∥u∥44+∥v∥44)dτ+E(0). (28)

Rearranging the terms two in the r.h.s. of (2) as follows

 4NE(t)+4∫(2−N)(|∇u|2+|∇v|2)dx=1N+2{16NE(t)− 8(N−2)∫(|∇u|2+|∇v|2)dx−2N(N−2)[g1∫|u|4dx+g2∫|v|4dx]− 4Ng(N−2)∫|u|2|v|2dx+8N(N−2)κRe∫u¯vdx}, (29)

and taking into account that

 16γIm∫(¯ux⋅∇u−¯vx⋅∇v)dx≤4(N−2)N+2(∥∇u∥22+∥∇v∥22)+16γ2(N+2)N−2X, (30)

and

 8γκIm∫|x|2v¯udx≤4γκX, (31)

from (2) we obtain the following estimate:

 d2Xdt2+2(N−2)N+2(2∥∇u∥22+2∥∇v∥22+Ng1∥u∥44+Ng2∥v∥44+2Ng∫|u|2|v|2dx)≤ ≤c1X+16NN+2E(t)+16κS0(t), (32)

which implies

 d2Xdt2+c2(∥∇u∥22+∥∇v∥22+∥u∥44+∥v∥44)≤c1X+16NN+2E(t)+16κS0(t). (33)

Then using (9) and (3) we arrive at

 d2Xdt2+c2(∥∇u∥22+∥∇v∥22+∥u∥44+∥v∥44)≤c1X+16NN+2E(0)+ +16κS0(0)e2γt+c3γ∫t0(∥∇u∥22+∥∇v∥22+∥u∥44+∥v∥44)dτ (34)

Next, introducing

 ρ(t)=∫t0∫σ0(∥∇u∥22+∥∇v∥22+∥u∥44+∥v∥44)dτdσ, (35)

and using function defined in (14), we rewrite (3) in the form

 d2dt2(X(t)+c2ρ(t))≤c1X(t)+c3γdρ(t)dt+d2F(t)dt2. (36)

From this inequality we obtain

 X(t)+c2ρ(t)≤F(t)+c1∫t0∫σ0X(τ)dτdσ+c3γ∫t0ρ(τ)dτ (37)

In order to complete the proof, we use the arguments of reductio ad absurdum. Let the conditions (18)–(19) hold, but a solution of the Cauchy problem (1) nevertheless exists for all . Then we can define

 T1=sup{t∈[0,T0]:X(s)≤M(T0)\,\,% for any \,\,s∈[0,t]}. (38)

From (37) and (18) for all one has

 X(t)+c2ρ(t)≤F(t)+c1M(T0)T202+c3γ∫t0ρ(τ)dτ ≤M(T0)−1+c1M(T0)T202+c3γ∫t0ρ(τ)dτ

Since , by Gronwall’s inequality

 ρ(t)≤M(T0)c2exp(c3c2γt). (40)

Using this estimate back in right hand side of (3) and using function defined by (16), we obtain

 X(t)≤F(t)+G(T0). (41)

Therefore, using (19) we conclude that and thus . Hence

 X(T0)≤F(T0)+1<0, (42)

which is impossible because .

## 4 “Early-collapse”

Now we consider the case, when the NLS equation with gain, i.e. Eq. (1a), contains attractive nonlinearity, while the nonlinearity of the system with loss i.e. Eq. (1b), as well as nonlinear coupling, are either repulsive or zero, i.e. we set

 g1>0,g2≤0,g≤0. (43)

For this set of parameters, one can use the idea of the “early-collapse” suggested in [25]. The method is based on the fact that even if the energy is a growing function, its growth can be controlled and it is hence possible to choose the initial conditions and such that the blowup occurs at sufficiently early times of the evolution.

In order to obtain the growth rate of the energy, we use (11) which subject to (43) allows one to obtain

 dEdt≤2γ∫(|∇u|2+|∇v|2+κ(¯uv+u¯v)−g12|u|4−g22|v|4−g|u|2|v|2)dx+ +2γκ(∥u∥22+∥v∥22)≤2γE(t)+2κγS0(t). (44)

Thus using (9) and the definition of the energy (10) we derive

 E(t)≤(E(0)+2κγS0(0)t)e2γt=:Emax(t). (45)

Next, we derive the following estimate

 16γIm∫(¯ux⋅∇u−¯vx⋅∇v)dx≤16γ2N−2X+4(N−2)(∥∇u∥2+∥∇v∥2), (46)

and from the identity (2) we obtain

 d2Xdt2≤c24X+4NEmax(t)+4κNS0(t),c4=2γ√κγ+N+2N−2. (47)

Now one can obtain the upper bound for :

 X(t)≤Z(t)ec4t, (48)

where

 Z(t):=X(0)+∫t0e−2c4s[Y(0)−c4X(0)+4N∫s0ec4σ[Emax(σ)+κS0(0)e2γσ]dσ]ds, (49)

and the function was defined in (45).

The obtained result can be reformulated as the following Theorem.

###### Theorem 2

Let , the coefficients , and satisfy (43), and initial conditions and in the Cauchy problem (1) are such that the function defined by (49) has a real positive zero . Then the solution of the problem (1) does not exist in the interval .

Since the occurrence of the blowup is now reduced to the study of the zeros of and this function depends on several parameters of the problem, we limit further consideration of this section by the numerical analysis. To this end now, as well as in all other numerical examples below, we limit the analysis to the Gaussian initial conditions (as the most typical for experimental settings):

 u0(x)=AπN/4aN/2exp(−|x|22a2),v0(x)=BπN/4bN/2exp(−|x|22b2), (50)

where , , , and are positive constants.

In Fig. 1 we show plots of the function for several particular choices of the parameters in the case , where we observe that for a proper choice of the initial conditions and the system parameters zeros of the function indeed exist. The time where becomes zero gives the respective (with is the minimal root obtained numerically). In Fig. 1(a) we show the behavior of for different values , which corresponds to different initial conditions for the equation (1b) with dissipation, provided the input of the NLS with gain, i.e. is fixed. Increase of the initial amplitude of the pulse subjected dissipation [i.e ] results in larger values and eventually leads to nonexistence of the zeros of . In the cases where does not exist, the question about the finite time blowup remains open.

In Fig. 1(b) we address different values of the gain/loss coefficient , and in Fig. 1(c) we illustrate different values of the coupling coefficient at fixed initial conditions. We observe increase of for stronger gain and dissipation [panel (b)], as well as for stronger coupling [panel (c)].

## 5 A remark on the Manakov model

The sufficient blowup conditions formulated in Sec. 2 allow for a simplification in the case of equal nonlinear coefficients: (below we refer to this model as Manakov’s one, after the original work [26]). Without loss of generality now we can set (which can be achieved by the renormalization , ).

An interesting property of the Manakov model is the existence of (at least) two integrals of motion. Indeed, straightforward computations show that the quantities

 S1(t)=S1(0)andS(t):=κS0(t)−γS2(t)=κS0(0)−γS2(0) (51)

do not depend on time. Then, from (8) one obtains

 d2S0dt2+4ω2S0=4κS, (52)

where as defined in (5). Thus in the case of unbroken the -symmetric phase (2) the quantity undergoes oscillatory dynamics:

 S0(t)=κSω2+S01cos(2ωt)+S02sin(2ωt), (53) S01=S0(0)(1−κω2)+S2(0)γκω2,S02=S3(0)γω, (54)

and hence there exists an upper bound for [c.f. (9)]:

 S0(t)≤S0,max:=κSω2+√S201+S202. (55)

In order to obtain sufficient conditions of the blowup in the Manakov case, we again employ Eq. (2) and Eq. (3) and, using that the quantity is conserved, we can replace the estimate (3) by the following one:

 d2Xdt2+2(N−2)N+2(2∥∇u∥22+2∥∇v∥22+Ng1∥u∥44+Ng2∥v∥44+2Ng∫|u|2|v|2dx)≤ ≤c1X+16NN+2(E(t)−κS1(0)). (56)

One can repeat all the subsequent steps of the proof of Theorem 1 using the newly obtained estimate (5). Therefore, in the Manakov case Theorem 1 can be reformulated with functions , and replaced by , and , respectively, where the new functions are defined as

 ^F(t) = X(0)+Y(0)t+8NN+2(E(0)−κS1(0))t2, (57) ^M(t) = supτ∈[0,t]^F(τ)+1, (58) ^G(t) = ^M(t)(c1t22+exp(48NγtN+2)−1). (59)

## 6 Numerical illustrations

The analytical results obtained above give sufficient conditions for the finite time blowup, but do not describe the blowup dynamics and its dependence on the parameters of the model, i.e. , , , and . In order to understand better the effect of those parameters on the phenomenon of blowup, now we resort to numerical analysis of the Cauchy problem (1). We concentrate on the case which is the most interesting one from the physical point of view. We also set and consider how the blowup depends on the nonlinear coupling and on the gain/loss coefficient . As we mentioned above, without loss of generality we can set (which is achieved by the renormalization , , ). Using another evident renormalization, without loss of generality we can set . Therefore we concentrate our attention on the system

 iut=−Δu+iγu+v−(|u|2+g|v|2)u, (60a) ivt=−Δv−iγv+u−(g|u|2+|v|2)v, (60b)

where and .

We first notice that the blowup conditions (18)–(19) can be easily satisfied by a proper choice of the initial pulses and . Several examples are presented in Fig. 2 which shows spatial profiles of initial conditions and chosen in the form of the Gaussian beams (50) and the behavior of the associated functions and defined by (14) and (16). Panels Fig. 2 (a–b) correspond to a situation when the blowup conditions are not satisfied as functions and grow monotonously. However, increasing powers and of the input beams [or one of the beams, see Fig. 2 (c–d)] or decreasing the characteristic widths and of the beams [Fig. 2(e–f)], one can easily satisfy the blowup conditions.

We also performed numerical simulations of the three-dimensional system (60) in the spherical symmetric case, i.e. assuming that the functions and depend only on the radius in the spherical coordinates and do not depend on the polar and azimuthal angles. In this case it is convenient to introduce new functions and and to reformulate the problem as follows

 ipt = −prr+iγp+q−r−2(|p|2+g|q|2)p, (61a) iqt = −qrr−iγq+p−r−2(g|p|2+|q|2)q. (61b)

We also limited the numerical study to a finite interval subject to the zero boundary conditions . The interval width is taken sufficiently large such that increase of practically does not affect the results presented below.

We solved system (61) using a semi-implicit finite-difference scheme proposed in [27]. The numerical simulations were carried on up to the time when the ratios and were of order , and then the simulations were interrupted. We used the spatial step of order and the adaptively decreasing temporal step which was of order in the beginning of simulations (i.e. at ) and of order at the time of the termination of the simulations. We also performed several additional runs with smaller spatial steps and checked convergence of the numerical solution.

In the numerical simulations we have observed different scenarios of the dynamics: dispersion of the initial pulses, growth of the solution and its derivative at the origin in one of the components (either with gain or with dissipation), and the simultaneous growth in both the components. The dispersion corresponds to the spreading of the initial pulses (which eventually occupy the entire computational domain ), while the observed growth is a presumable manifestation of the blowup. Our numerical results allow us to conjecture that depending on the initial conditions and parameters of the model the blowup can occur either in both or in only one component. In the latter case this can be a component either with gain or with damping.

Summarizing results of dynamical simulation of system, we first of all conjecture that the obtained above blowup conditions (18)–(19) are not necessary for the blowup to occur: it is easy to find initial conditions that do not satisfy (18)–(19) but nevertheless grow in numerical simulations. Moreover, our numerics indicate that the growth can occur even for initial conditions with and , i.e. in the case when conditions (18)–(19) obviously can not be satisfied for any .

Several (the most interesting) examples of the blowup-like dynamics obtained from numerical solution of system (61) are presented in Fig. 3. The case of the unbroken symmetry is illustrated in Fig. 3(a) from where one observes that positive and negative values of result in different scenarios of the blowup-like dynamics. Namely, for (this is the attractive inter-species nonlinearity) we observe that the the intensity growth occurs in both the components. However for and the same initial conditions display growth only in the dissipative component while in the gain component the intensity at the origin decays. Notice also that the maximal squared amplitude in the component with gain for and (not shown on the plot) was of the order of at the moment when the simulations were terminated. It is also interesting to notice that for and the only component “responsible” for the blowup is the one with dissipation [i.e. ] and that the growth with (the dotted curves) occurs earlier than that with (the dashed curves) in spite of the additional “defocusing” that large negative values of induce. We finally notice that the case is not covered by the condition of the Theorem 1, which is another illustration for the fact that the obtained blowup conditions (18)–(19) are sufficient but not necessary.

Repeating the same simulations for [the broken symmetry, Fig. 3(b)] we do not observe any quantitative difference in the blowup scenarios with respect to the case . However it is interesting to notice that for the chosen initial conditions and for the three different values of addressed in Fig. 3(b) the numerical blowup for the broken symmetry occurs later than that with the unbroken symmetry in Fig. 3(a) (this is especially well-visible for and ). Ho