Some case example exact solutions for quadratically nonlinear optical media with \mathcal{PT}-symmetric potentials

# Some case example exact solutions for quadratically nonlinear optical media with PT-symmetric potentials

Y.N. Truong Vu Department of Mathematics and Statistics, Amherst College, Amherst, MA, USA    J. D’Ambroise Department of Mathematics and Statistics, Amherst College, Amherst, MA, USA    F.Kh. Abdullaev Department of Physics, Faculty of Sciences, IIUM, Jln. Indera Mahkota, Sultan Ahmad Shah, 25200, Kuantan, Malaysia    P.G. Kevrekidis Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA USA
###### Abstract

In the present paper we consider an optical system with a -type nonlinearity and unspecified -symmetric potential functions. Considering this as an inverse problem and positing a family of exact solutions in terms of cnoidal functions, we solve for the resulting potential functions in a way that ensures the potentials obey the requirements of -symmetry. We then focus on case examples of soliton and periodic solutions for which we present a stability analysis as a function of their amplitude parameters. Finally, we numerically explore the nonlinear dynamics of the associated waveforms to identify the outcome of the relevant dynamical instabilities of localized and extended states.

###### pacs:
42.65.Ky, 42.65.Sf, 42.65.Tg, 42.81.Dp

## I Introduction

Recently a great deal of attention has been devoted to the investigation of quantum and classical systems with -symmetric non-Hermitian Hamiltonians. It was shown in the seminal paper Bender () that such type of Hamiltonians can have real eigenvalues. Due to the analogy between the Schrodinger equation and the paraxial wave equation in optics, this result has applications in optics; this is one among numerous other areas studied over the past decade. For optical beams -symmetry imposes the condition on the complex refractive index : even in space for the real part of the refractive index and odd in space for the imaginary part . Recently effects of -symmetry have been observed in optical experiments Gao (). However, optics is certainly not the sole area where -symmetric applications have recently emerged. More specifically, a mechanical system realizing -symmetry has been proposed and realized in pt_mech (), while a major thrust of efforts has focused on the context of electronic circuits; see e.g. the original realization of tsampikos_recent () and the more recent review of this activity in tsampikos_review (). Additionally, we note that further intriguing realizations of -symmetry have also emerged e.g. in the realm of whispering-gallery microcavities pt_whisper (). While many of these experimental realizations have been chiefly explored at the level of linear dynamics, the intrinsic nonlinearity of optical systems Gao () and the potential nonlinearity also of electrical ones (e.g. in the form of of a -symmetric dimer of Van-der-Pol oscillators bend_kott ()) have prompted a considerable amount of work at the interface of nonlinearity and -symmetry.

Nonlinearity leads to new effects in the -symmetric systems, such as solitons (and vortices) in continuous Musslimani08 (); achilleos () and discrete nonlinear optical media with -symmetric potentials Dmitriev (), gap solitons in media with -symmetric periodic potentials Aceves (), non-reciprocity, instabilities and nonlinear -transitions in -symmetric (nonlinear) couplers Ramezani (); Sukhorukov (); likev (); pickton (); barasflach (); barasflach2 (), as well as the smoothing of the spectral singularity in transmission arxive (), among many others. Important applications to nonlinear plasmonic systems and metamaterials are under recent investigation plasm (); metha (). Additional developments are connected with nonlinear -symmetric lattices Zez (); KevPel (); tyugin1 (); pelinovsky () and -symmetry management DM (); Longhi (); Horne () etc.

In nonlinear wave equations with -symmetric terms, solitonic solutions can exist as was shown recently e.g. in the works of Musslimani08 (); Abdullaev10 (); Abdullaev11 (); Konotop11 (); Malomed11 (). In the case of an NLS system with inhomogeneous in space loss/gain parameters such solutions were found for linear -potentials in the work Musslimani08 (); Abdullaev10 (); Salerno13 (); Kevrekidis13 (), and for nonlinear -potentials in the works Abdullaev10 (); Abdullaev11 (), for a nonlocal NLS equation in nonlocal () and for a cubic- quintic model in QC (). Naturally, it is also of interest to find exact solutions for solitons in other physically important systems. Recent numerical simulations of the system with -symmetric potential showed the existence of stable bright solitons Moreira (), as well as of gap solitons in periodic -symmetric potentials MKM (). Discrete -symmetric systems with quadratic nonlinearity were also explored at the level of oligomer systems kaidima ().

In the present work we will study the system with -symmetric potentials, describing wave processes in quadratically nonlinear media with spatially distributed gain/loss parameters. Such systems are, in principle, of interest for nonlinear optics and potentially even for atomic-molecular Bose-Einstein condensates trapped in complex potentials (see, for example, a realization of imaginary potential in Oberthaler ()). In this paper, we find exact solitary and periodic solutions of the system. We begin in Section II by outlining the mathematical model. In Section III we derive exact solutions in terms of the Jacobi elliptic cnoidal function, and we present special case solutions for which we later examine stability properties. We essentially follow an inverse-function approach, somewhat reminiscent of flach () to obtain such exact solutions for suitably tailored -symmetric potentials in the presence of the quadratic nonlinearity. In Section IV we present the stability analysis of the obtained solutions as a function of the solution parameters (such as their amplitudes), and we show the results of propagation of the solutions in the (analogous to “time”) variable . Finally, in Section V we make our concluding statements and present a number of possibilities for future work.

## Ii The model

Let us consider the system describing the first harmonic (FH) and second harmonic (SH) propagation in quadratically nonlinear media with -symmetric potentials as follows:

 iuz + d1uxx+V1(x)u+iW1(x)u=u∗v (1) ivz + d2vxx+κv+V2(x)v+iW2(x)v=u2

where are even functions of corresponding to real parts of the refraction index, and are odd functions of pertaining to imaginary parts thereof. describe the inhomogeneous in space gain/loss. Seeking standing waves in the form: and we obtain the system

 ωU + d1Uxx+V1(x)U+iW1(x)U=U∗V (2) σV + d2Vxx+V2(x)V+iW2(x)V=U2

where . It is useful to introduce the amplitude-phase decomposition in the form . Solitonic solutions for have been reported e.g. in KS (); WD94(2) (); WD94 (); MST94 (); Hayata (), cnoidal wave solutions in Parker (), and solitons in the conservative 2D system with a potential were considered recently in SM12 (). For a review of solitary wave dynamics in quadratic systems see e.g. buryak ().

Assuming that are real, and that is either real or purely imaginary, we obtain the system

 ρ∗1ρ2 = ωρ1+d1ρ1,xx−d1ρ1(θx)2+V1(x)ρ1 (3) ρ21 = σρ2+d2ρ2,xx−4d2ρ2(θx)2+V2(x)ρ2 Wj(x)ρ2j = −jdj(ρ2jθx)x

for .

## Iii Solutions in terms of the cnoidal function

We begin by writing and for . This gives the system

 F∗1(y)F2(y) = ωF1(y)+d1r2Γ1(y)−d1F1(y)G2(y)+V1(x)F1(y) (4) F21(y) = σF2(y)+d2r2Γ2(y)−4d2F2(y)G2(y)+V2(x)F2(y) Wj(x)Fj(y) = jrdjdn(rx,k)sn(rx,k)(2F′j(y)G(y)+G′(y)Fj(y))

for

 Γj(y) = y(2k2−1−2k2y2)F′j(y)+(1−y2)(k2y2+1−k2)F′′j(y) (5)

with and where the primes denote differentiation with respect to . Notice that by writing (4) in terms of both and we avoid restrictions on the domain which would be applicable if we composed with an inverse function.

To find exact solutions, we apply the reverse engineering approach Abdullaev10 (); BK (); this type of technique was applied much earlier in order to obtain exact traveling wave solutions in dynamical lattices flach (). Our general strategy is to first specify the form of the functions and then use (4) to solve for appropriate potentials . Thus, we rewrite the above equations as:

 V1(x) = F∗1(y)F1(y)F2(y)−d1r2Γ1(y)F1(y)+d1G2(y)−ω (6) V2(x) = F21(y)−d2r2Γ2(y)F2(y)+4d2G2(y)−σ Wj(x) = jrdjdn(rx,k)sn(rx,k)(2F′j(y)G(y)Fj(y)+G′(y))

for and . In each of the following subsections we make specific choices of and in such a way that the resulting potentials , in (6) do not contain terms with denominators (that may lead to singularities) and they also obey the requirements of -symmetry ( even functions of and odd functions of ). In other words, we require at least that the conditions

 F1(y) ∣ Γ1(y) (7) F2(y) ∣ (F21(y)−d2r2Γ2(y)) (8) Fj(y) ∣ F′j(y)G(y) (9)

for are met for any choices of that we specify.

### iii.1 Polynomial Functions

Consider in the form of generalized polynomials in

 F1(y)=i{0,1}k1∑n=s1Cnyn,F2(y)=k2∑m=s2Dmym (10)

with coefficients , integer indexing bounds with , and . Notice that we restrict our attention here to polynomials with at least two terms. The case of a monomial type solution will be included in the next section where we consider a more general class of power functions. Recall that, in the derivation of (3), is required to be either real or purely imaginary. To show this in (10) we have included an optional multiple of in the definition of . In the following subsections we outline the process of solving for . We separate into two cases which are convenient based on the resulting maximal power of the polynomial conditions (7)-(8).

#### iii.1.1 Cnoidal parameter k≠0

To proceed in solving for using (6) and assuming (10) we must satisfy condition (7). That is, we must have that is a factor of the polynomial . For , will have maximal power so that (7) amounts to the condition

 k1+2∑n=s1−2((n+2)(n+1)(1−k2)Cn+2+n2(2k2−1)Cn−(n−2)(n−1)k2Cn−2)yn=(α1y2+β1y+γ1)k1∑n=s1Cnyn (11)

for some with . For convenience we use the convention that for any index . Equating the coefficients of (11) then gives

 (n2(2k2−1)−γ1)Cn=(α1+(n−2)(n−1)k2)Cn−2+β1Cn−1−(n+2)(n+1)(1−k2)Cn+2 (12) α1=−k1(k1+1)k2,s1(s1−1)(1−k2)=0,(s1+1)s1(1−k2)Cs1+1=0

where the first equation is a recursion relation that holds for and the latter three equations are obtained from equating the coefficients of the , terms in (11), respectively. Note that the latter equations have made use of the conditions and for .

From the first equation of the latter three in (12), we now have that the coefficient is determined by the highest degree chosen for . The latter two equations in (12) then give us a starting point for finding more specific solutions. That is, we can look for solutions with in terms of and these latter two equations are satisfied. Alternatively, we can look for solutions with in which case the latter equations of (12) give that either so that the first term in the polynomial is required to be a constant, or the first term is required to be of degree with coefficient . We will proceed here to outline the general solution for cnoidal parameter . Later towards the end of this section we will focus primarily on for the special case where are quadratic polynomials.

The recursive equations in (12) give us conditions for the constants . Later we will choose and then use conditions (12) to solve for the coefficients of and also . Now that (7) is satisfied by imposing (12), we have via (6) as

 V1(x)=±F2(y)−d1r2(α1y2+β1y+γ1)+d1G2(y)−ω (13)

with as usual. The plus sign in (13) applies to real (using in (10)) and the minus sign applies to purely complex (using in (10)). Since the cnoidal function is an even function of , is an even function of so that this potential is compatible with -symmetry. Notice that may not be a polynomial if is not a polynomial. In Section IV, we will specify choices for the function and we will choose so that as .

To solve for we must have that statisfies condition (8). One way to proceed is to require that

 F21(y)=F2(y)P(y) (14)

where is a polynomial in . Then, similar to the case, we may also impose that so that for some constants with . Using a similar procedure as in the case, now equations (12) must hold after performing the replacements , and in the subscripts . Using this -version of equation (12), now the coefficient is determined by the highest degree of the polynomial . Since here the -version of the latter two equations in (12) either requires us to take so that the first term in the polynomial is required to be a constant, or alternately the first term is required to be of degree with coefficient . Also, the constants are required to satisfy the same recursive equations in (12) but with appropriate -version described above.

The most obvious choice in order to satisfy both (14) and the -version of (12) is to take as scalar multiples of each other. In other words,

 F2(y)=i{0,1}AF1(y) \ and \ P(y)=F1(y)i{0,1}A (15)

for some . Since are required to be real-valued the multiple of in front is included only if it’s included in the definition of in (10). Now we have the real-valued potential function

 V2(x)=F1(y)i{0,1}A−d2r2(α2y2+β2y+γ2)+4d2G2(y)−σ (16)

where later in specific examples we will choose to be such that as .

Next we want to determine an appropriate form for the function with that will satisfy (9). Since have been chosen to be scalar multiples of each other, if (9) holds for then it holds for . So, we take

 G(y)=T(y)F1(y) (17)

for a function and this gives via (6)

 Wj(x) = jrdjdn(rx,k)sn(rx,k)(2F′1(y)T(y)+G′(y)) (18)

for and . Since is an odd function of and is even, are odd functions of as required by -symmetry as long as the quantity is an even function of . This is reasonable since cn is an even function of .

Now we have a complete solution of (4) given by in (10), in (13) and (16), and in (18), all under the conditions seen in (12), (15), (17). To be more explicit, let us focus on details in the case where is a quadratic function and so that . Consider of the form for . Then and so that the latter two equations in (12) are satisfied. Now the remaining equations in (12) give , and the conditions

 C1=β1C0,C2=−6C0+β1C14,C3=0=(α1+2)C1+β1C29, (19)

that we may solve for the four constants .

In the case that , (19) gives and so that combining with (15) we have

 F1(y)=C0(1−32y2),F2(y)=AF1(y). (20)

Proceeding with as in (17) for any function we have by (13), (16) and (18) that

 V1(x) = AC0(1−32y2)+6d1r2y2+d1G2(y)−ω (21) V2(x) = C0A(1−32y2)+6d2r2y2+4d2G2(y)−σ Wj(x) = jrdjtanh(rx)sech(rx)(−6C0yT(y)+G′(y)).

Equations (20)-(21) give us a solution that we will refer to as the dark-dark soliton case. In Section IV we show the dark soliton shape, analyze the stability of the dark-dark soliton, and show plots over the propagation variable for a specific choice of the function and other parameters.

We also consider here the case of , for which (19) gives two possibilities for the coefficients of the polynomial . Then combined with (15) we have

 F1(y)=C0(1±√22y+4y2),F2(y)=AF1(y). (22)

Letting be as in (17) for some function we then obtain

 V1(x) = AC0(1±√22y+4y2)−d1r2(±√22y−6y2)+d1G2(y)−ω (23) V2(x) = C0A(1±√22y+4y2)−d2r2(±√22y−6y2)+4d2G2(y)−σ Wj(x) = jrdjtanh(rx)sech(rx)(2C0(8y±√22)T(y)+G′(y))

for . These solutions are quite interesting in their own right, as the one with the sign corresponds to an antidark-antidark soliton setting of a pair of bright solitary waves on top of a non-vanishing background. On the other hand, the solution with the sign is especially structurally complex, resembling a conglomeration of multiple –more specifically of 4– dark solitons.

#### iii.1.2 Cnoidal parameter k=0 and y=cos(rx)

For of the polynomial form (10) now consider the case of , or . In solving for the polynomial condition analogous to (11) has maximal power . This roughly makes sense because when we differentiate a or the result is a function of the same overall power (in contrast to derivatives of sech and tanh, for example). The analogue of (11) in this case is

 k1∑n=s1−2((n+2)(n+1)Cn+2−n2Cn)yn=γ1k1∑n=s1Cnyn (24)

for . Equating coefficients we get

 Cn = (n+2)(n+1)Cn+2n2+γ1 \ for \ n∈{s1,…,k1−1} γ1 = −k21,s1(s1−1)=0,(s1+1)s1Cs1+1=0 (25)

where the latter three equations came from equating the coefficients of the , terms in (24), respectively. As before, the first equation of the latter three in (III.1.2) shows that the coefficient is determined in terms of the maximal power of the polynomial . The latter two equations in (III.1.2) then show that either so must have a constant term, or alternatively the first term is required to be of degree with coefficient . The remaining recursive equations in (III.1.2) then give us conditions for unknowns . Choosing one of these coefficients will lead us to find the others. In solving for we have similar conditions to (III.1.2) for the constants and where in (III.1.2) one should replace , and in the subscripts . We proceed in a similar way as in the case above, assuming the forms of as seen in (14), (15), (17) and finally we have

 V1(x) = ±F2(y)−d1r2γ1+d1G2(y)−ω (26) V2(x) = F1(y)i{0,1}A−d2r2γ2+4d2G2(y)−σ Wj(x) = jrdjsin(rx)(2F′1(y)T(y)+G′(y))

for .

Let us focus on the details of the quadratic case where for and , . (III.1.2) then gives and . Then, we have

 F1=C0(1−2y2),F2(y)=AF1(y). (27)

We also have by the analogue of (III.1.2) (see description above). Proceeding with as in (17) for some function we have

 V1(x) = AC0(1−2y2)+4d1r2+d1G2(y)−ω (28) V2(x) = C0A(1−2y2)+4d2r2+4d2G2(y)−σ Wj(x) = jrdjsin(rx)(−8C0yT(y)+G′(y))

for . Equations (27)-(28) give us a solution that we call the quadratic oscillatory case. In Section IV we show the shape of the solution, analyze the stability, and explore its dynamics over the evolution variable ().

### iii.2 Power Functions

Next we take and to be power functions

 F1(y)=i{0,1}Cyp1,F2(y)=Dyp2 (29)

for and . This special case considerably simplifies the relevant compatibility conditions. In particular, substituting (29) into (6) and examining conditions (7)-(8) we require that either and , or .

#### iii.2.1 Cnoidal parameter k=1 and y=sech(rx)

In the case of we can have non-integer and ; this is in contrast to the polynomial case. For this case, we apply (6) and find

 V1(x) = ±Dyp2−d1r2p1(1−2y2)−d1r2p1(p1−1)(1−y2)+d1G2(y)−ω (30) V2(x) = ±C2Dy2p1−p2−d2r2p2(1−2y2)−d2r2p2(p2−1)(1−y2)+4d2G2(y)−σ Wj(x) = jrdjtanh(rx)(2pjG(y)+yG′(y))

for . As long as are even functions in , which is easy to choose since is an even function of , then and are odd and compatible with the symmetry criterion. We refer to the solutions given in (29) and (30) as the bright-bright soliton case. In section IV we show the wave’s shape, analyze its stability and explore its direct numerical evolution. Note that the case corresponds to solitonic solutions found by Karamzin-Sukhorukov in KS (), and to solitonic solutions found by Menyuk et al. MST94 ().

#### iii.2.2 Powers p1=p2=1

In the case , we have

 V1(x) = ±Dy−d1r2(2k2−1−2k2y2)+d1G2(y)−ω (31) V2(x) = ±C2Dy−d2r2(2k2−1−2k2y2)+4d2G2(y)−σ Wj(x) = jrdjdc(rx,k)sn(rx,k)(2G(y)+yG′(y))

where . Similar to the previous section, we only need to choose a function so that is an even function of in order to satisfy the symmetry criterion. Equations (29) and (31) give us a solution we call the linear oscillatory case. More details are included in Section IV.

### iii.3 Other solutions

Here, we introduce a possibility which is distinct from the previous ones as follows. We introduce and in the form

 F1=iAyp(1−y2)1/2,F2=Byq (32)

for i.e., a non-polynomial form. By (6) and examining conditions (7)-(8) we find that we must have and . That is, we require that either with , or , or . We will focus on the former two cases. As for in (6) and condition (9), we find that must be in the form

 G(y)=Cya(1−y2)b (33)

where .

#### iii.3.1 Cnoidal parameter k=1 and y=sech(rx)

In the case, we have the solutions

 ω =−d1r2p2, V1(x)=−Byq+d1r2y2(p+1)(p+2)+d1G2(y) (34) σ =−d2r2q2, V2(x)=−A2By2p−q(1−y2)+d2r2y2q(q+1)+4d2G2(y)

with and in the form of (33). Since for this family of solutions the form of is specified, it is immediately clear which choices of will give as . In contrast to previous sections, those choices have been made in (34). Also, we have by (6)

 W1 = Crd1sech(rx)tanh(rx)(ya−1(1−y2)b(a+2p)−2ya+1(1−y2)b−1(b+1)) W2 = 2Crd2tanh(rx)(ya(1−y2)b(a+2q)−2bya+2(1−y2)b−1). (35)

Equations (32)-(III.3.1) give us a solution that bears a bright soliton coupled with a dark-in-bright soliton. The latter involves a pair of bright solitary waves coupled in a bound state anti-symmetric (i.e., they bear a phase difference of ) configuration; another example of this form has been previously reported e.g. in njp (). More details on the propagation of this solution and its stability are included in Section IV.

#### iii.3.2 Powers p=q=1

In the case where , we obtain

 ω =−d1r2(5k2−4), V1(x)=−By+6d1r2k2y2+d1G2(y) (36) σ =−d2r2(2k2−1), V2(x)=−A2By(1−y2)+2d2r2k2y2+4d2G2(y).

Again here since is known we have made choices of reflected in (36) so that as . With as in (33) we have

 W1 = rCd1sn(rx,k)dn(rx,k)(ya−1(1−y2)b(a+2)−2ya+1(1−y2)b−1(b+1)) W2 = 2rCd2dc(rx,k)sn(rx,k)(ya(1−y2)b(a+2)−2bya+2(1−y2)b−1). (37)

We will refer to the solution in (32)-(33) and (