TThe fractional nonlinear PT dimer Mario I. Molina
Departamento de F´ısica, Facultad de Ciencias, Universidad de Chile, Casilla 653, Santiago, Chile (Dated: February 5, 2021)We examine a fractional Discrete Nonlinear Schrodinger dimer, where the usual first-order deriva-tive of the time evolution is replaced by a non integer-order derivative. The dimer is nonlinear (Kerr)and PT -symmetric, and we examine the exchange dynamics between both sites. By means of theLaplace transformation technique, the linear PT dimer is solved in closed form in terms of Mittag-Leffler functions, while for the nonlinear regime, we resort to numerical computations using thedirect explicit Grunwald algorithm. In general, the main effect of the fractional derivative is theonset of a monotonically decreasing time envelope for the amplitude of the oscillatory exchange. Inthe presence of PT symmetry, the dynamics shows damped oscillations for small gain/loss in bothsites, while at higher gain/loss parameter values, the amplitudes of both sites grows unbounded. Inthe presence of nonlinearity, selftrapping is still possible although the trapped fraction decreases asthe nonlinearity is increased past threshold, in marked contrast with the standard case. a r X i v : . [ n li n . PS ] F e b I. INTRODUCTION
The topic of fractional calculus has experienced a rekindled interest in recent times. Essentially, it extends the notionof a derivative or an integral of integer order, to one of a fractional order, ( d n /dx n ) → ( d α /dx α ) for real α . Thesubject has a long history, dating back to letters exchanged between Leibnitz and L’Hopital, and later contributionsby Euler, Laplace, Riemann, Liouville, and Caputo to name some[1–5]. The starting point was the computation of d α x k /dx α , where α is a non-integer number: d n x k dx n = Γ( k + 1)Γ( k − n + 1) x k − n → d α x k dx α = Γ( k + 1)Γ( k − α + 1) x k − α . (1)For instance, ( d / /dx / ) x = (2 / √ π ) √ x , and dx/dx = ( d / /dx / )( d / /dx / ) x = (2 √ π )(Γ(3) / Γ(1)) x = 1, asexpected. From Eq.(1) the fractional derivative of an analytic function f ( x ) = (cid:80) k a k x k can be computed by derivingterm by term. This basic procedure is not exempt from ambiguities. For instance, ( d α /dx α ) 1 = ( d α x /dx α ) =(1 / Γ(1 − α )) x − α (cid:54) = 0, according to Eq.(1). However, one could have also taken ( d α − /dx α − )( d/dx ) 1 = 0. For thecase of a fractional integral, a more rigorous starting point is Cauchy’s formula for the integral of a function. Fromthe definition I x f ( x ) = (cid:90) x f ( s ) ds, (2)we apply the Laplace transform L to both sides of Eq.(2) L { I x f ( x ) } = (1 /s ) L{ f ( x ) } . (3)After n integrations, one obtains L { I nx f ( x ) } = (1 /s n ) L{ f ( x ) } . (4)Extension to fractional α is direct: L { I αx f ( x ) } = (1 /s α ) L{ f ( x ) } . (5)After noting that the RHS of Eq.(5) is the product of two Laplace transforms we have, after using the convolutiontheorem I αx f ( x ) = 1Γ( α ) (cid:90) x f ( s )( x − s ) − α ds. (6)From this definition, it is possible to define the fractional derivative of a function f ( x ) as d α f ( x ) dx α = (cid:18) d m dx m (cid:19) I m − αx f ( x ) == d m dx m (cid:20) m − α ) (cid:90) x ( x − s ) m − α − f ( s ) ds (cid:21) , (7)where, m − < α < m . Eq.(7) is known as the Riemann-Liouville form. An alternative, closely related form, is theCaputo formula[5]: d α f ( x ) dx α = I m − αx (cid:18) d m dx m (cid:19) f ( x ) == 1Γ( m − α ) (cid:90) x ( x − s ) m − α − f ( m ) ( s ) ds, (8)which has some advantages over the Riemann-Liouvuille form for differential equations with initial values. The varioustechnical matters that arise in fractional calculus have prompted a whole line of research that has extended to currenttimes. Long regarded as a mathematical curiosity, it has now regained interest due to its potential applications tocomplex problems in several fields: fluid mechanics[6, 7], fractional kinetics and anomalous diffusion[8–10], strangekinetics[11], fractional quantum mechanics[12, 13], Levy processes in quantum mechanics[14], plasmas[15], electricalpropagation in cardiac tissue[16] and biological invasions[17]. In general, fractional calculus constitutes a naturalformalism for the description of memory and non-locality effects found in various complex systems. Experimental V e , c e , c Figure 1. Nonlinear anisotropic fractional dimer where the excitation is on site 1 initially ( t = 0). In the PT case, (cid:15) = − (cid:15) = iγ realizations are not straightforward given the nonlocal character of the coupling, however some optical setups havebeen suggested that could measure the effect of fractionality on new beam solutions and new optical devices[18, 19]On the other hand, when dealing with effectively discrete, interacting units, as one encounters in atomic physics(interacting atoms), or in optics ( coupled optical fibers), it is common to deal with discrete versions of the continuumSchr¨odinger equation, or the paraxial wave equation. The effective discreteness comes from expanding the solutionsought in terms of (continuous) modes that can be labelled unambiguously. The simplest of such examples is thebonding, anti-bonding electronic mode that one finds for a two-sites (dimer) molecule after diagonalizing the two-siteSchr¨odinger equation in the tight-binding approach. Something similar happens in optics, where the paraxial equationis formally equivalent to the Schr¨odnger equation. In that case, for two optical waveguides, the total electric field isexpanded in terms of the electromagnetic modes in each guide which interact through the evanescent field betweenthe two guides giving rise to a transversal dynamics for the optical power. The procedure can be extended to N interacting units, either atoms or waveguides, where the relevant dynamics is given by a discretized version of theSchr¨odinger equation for N units[20, 21]. Of course, at the end one has to collect all the discrete amplitudes andmultiply them by the corresponding continuous mode profiles and superpose them, to obtain the final field. Thesimplest case N = 2 is termed a dimer and oftentimes constitute a basic starting point when studying an interacting,discrete system. Ensembles of interacting dimers have been studied before in classical and quantum statistics[22–24], and more recently, they have been considered in model of correlated disorder[26] and in magnetic metamaterialmodeling[25].In this work we consider the discrete Schr¨odinger equation for a dimer system, where the standard time derivativeis replaced by a fractional one. The dimer considered is rather general and contains asymmetry, PT symmetry andnonlinearity (Fig.1). Our main interest is in ascertaining the effect of the fractional derivative on the excitationexchange between the sites, its stability and selftrapping behavior, for several cases of interest. II. THE FRACTIONAL DIMER
Let us consider the fractional evolution equations for a general nonlinear dimer i d α C ( t ) dt α + (cid:15) C ( t ) + V C ( t ) + χ | C ( t ) | C ( t ) = 0 i d α C ( t ) dt α + (cid:15) C ( t ) + V C ( t ) + χ | C ( t ) | C ( t ) = 0 (9)where α is the fractional order of the (Caputo) derivatives with 0 < α <
1. Quantities C , are probability amplitudesin a quantum context, or electric field amplitudes, in an optical setting. Parameter V is the coupling term and χ isthe nonlinearity parameter.Let us first consider the case of a general linear ( χ = 0) dimer, and assume C (0) = 1 , C (0) = 0. We will solve thiscase in closed form by the use of Laplace transforms: For 0 < α <
1, the Laplace transform of the Caputo fractionalderivative of order α is given by L{ f ( t ) } = s α L{ f ( t ) } − s α − f (0 + ) . (10)After applying the Laplace transform L to both sides of Eq.(9) we have i ( s α L ( C ) − s α − ) + (cid:15) L ( C ) + V L ( C ) = 0 i s α L ( C ) + (cid:15) L ( C ) + V L ( C ) = 0 . (11)Solving for L ( C ) and L ( C ) gives: L ( C ) = − i ( (cid:15) + is α ) s α − s α − i ( (cid:15) + (cid:15) ) s α + V − (cid:15) (cid:15) (12)and L ( C ) = i V s α − s α − i ( (cid:15) + (cid:15) ) s α + V − (cid:15) (cid:15) . (13)Using the inverse Laplace formula[27] L − (cid:26) s ρ − s α + as β + b (cid:27) = t α − ρ ∞ (cid:88) r =0 ( − a ) r t ( α − β ) r E r +1 α,α +( α − β ) r − ρ +1 ( − bt α ) , (14)we obtain: C ( t ) = ∞ (cid:88) r =0 ( i ( (cid:15) + (cid:15) )) r t αr E r +12 α,αr +1 (( (cid:15) (cid:15) − V ) t α ) − i(cid:15) t α ∞ (cid:88) r =0 ( i ( (cid:15) + (cid:15) )) r t αr E r +12 α,α (1+ r )+1 (( (cid:15) (cid:15) − V ) t α ) (15) C ( t ) = iV t α ∞ (cid:88) r =0 ( i ( (cid:15) + (cid:15) )) r t αr E r +12 α,α (1+ r )+1 (( (cid:15) (cid:15) − V ) t α ) (16)where E γα,β ( z ) is defined as E γα,β ( z ) = ∞ (cid:88) k =0 ( γ ) k z k k ! Γ( αk + β ) (17)where ( γ ) n = Γ( γ + n ) / Γ( γ ), and α, β, γ ∈ C , Re( α ) > , Re( β ) > , z ∈ C . Figure 2 shows examples of the timeevolution of the square of the dimer amplitudes, for several site energy parameters, and fractional derivative orders. Ingeneral we observe that, as soon as α differs from unity, the dynamics is either bounded or unbounded, depending onthe values of the site energy parameters. For the bounded cases, there is some oscillation initially, with a decreasingenvelope towards zero. A. The linear PT dimer A particularly interesting case of Eq.(9), is the fractional PT -symmetric dimer. For systems that are invariantunder the combined operations of parity ( P ) and time reversal ( T ), it was shown that they display a real eigenvaluespectrum, even though the underlying Hamiltonian is not hermitian[29, 30]. In these systems there is a balancebetween gain and loss, leading to a bounded dynamics. However, as the gain/loss parameter exceeds a certain valuethe system undergoes a spontaneous symmetry breaking, where two or more eigenvalues become complex. At thatpoint, the system loses its balance and its dynamics becomes unbounded. According to the general theory, for oursystem to be PT symmetric, the real part of the site energies in Eq.(9) must be even in space while the imaginarypart must be odd: Re( (cid:15) ) = Re( (cid:15) ) and Im( (cid:15) ) = − Im( (cid:15) ). For simplicity we take the real parts of (cid:15) , (cid:15) as zero andthus, (cid:15) = − (cid:15) ≡ i (cid:15) , where (cid:15) is the gain/loss parameter. This leave us with the equations: i d α C ( t ) dt α + i(cid:15)C ( t ) + V C ( t ) = 0 i d α C ( t ) dt α − i(cid:15)C ( t ) + V C ( t ) = 0 (18) | C , | C | C , | C | C , | C ✏ = 2 , ✏ = , ↵ = 1 ✏ = 2 , ✏ = ↵ = 0 . ✏ = 2 , ✏ = ↵ = 0 . ✏ = 2 , ✏ = ↵ = 0 . ✏ = 0 . , ✏ = 0 . ↵ = 0 . ✏ = 0 . , ✏ = 0 . ↵ = 0 . ✏ = 0 ✏ = 0 ↵ = 0 . ✏ = 1 . ✏ = 2 ↵ = 0 . ✏ = 1 . ✏ = 2 ↵ = 0 . Figure 2. Dimer amplitudes for the linear case ( χ = 0) and several site energy parameters (cid:15) , (cid:15) and different fractionalderivative orders α . Solid(dashed) line denotes | C | ( | C | ). whose exact solutions can be extracted from the general solution, Eqs.(15),(16) as C ( t ) = − (cid:15)t α E α,α +1 (( (cid:15) − V ) t α ) + E α, (( (cid:15) − V ) t α ) C ( t ) = iV t α E α,α +1 (( (cid:15) − V ) t α ) , (19)where, E α,β ( z ) = E α,β ( z ) is known as the generalized Mittag-Leffler function E α,β ( z ) = (cid:88) k z k Γ( αk + β ) . (20)It is the natural extension of the exponential function and plays the same rol for fractional differential equations,as the exponential function does for the standard integer differential equations. Figure 3 shows examples of timeevolutions for | C | , | C | for several fractional orders and several gain/loss parameter values. As we can see, as α decreases from 1, both amplitudes decrease too, with oscillations that become monotonically decreasing, convergingto zero at long times.The asymptotic behavior of C ( t ) , C ( t ) depends on the behavior of the Mittag-Leffler functions E α,β ( z ) at largevalues of | z | . After writing z = | z | exp( iφ ), we have[28] E α,β ( z ) ≈ (1 /α ) Q − β exp( Q ) (21)where Q = z /α = exp((1 /α ) log( | z | ) + i φ ) and | φ/α | ≤ π ( φ = (cid:15) − V ). This implies,exp( Q ) = exp( | z | /α cos((1 /α ) φ )) × exp( i | z | /α sin((1 /α ) φ )) (22)Thus, bounded behavior in time will occur for ( π/ < | φ/α | < π , while unbounded behavior occurs for 0 < | φ/α | <π/
2. In our case, φ = arg( (cid:15) − V ) = 0 , π , implying that C ( t ) and C ( t ) will increase (decrease) asymptotically intime if ( (cid:15) − V ) is positive (negative). This behavior is sketched in Fig.4. | C , | C | C , | C | C , | C ( a ) ( b ) ( c )( d ) ( f )( g ) ( h ) ( i )( e ) Figure 3. Dimer amplitudes | C ( t ) | (continuous line) and | C ( t ) | (dashed line) for the linear PT case ( χ = 0 , (cid:15) (cid:54) = 0) forseveral fractional derivative orders and various gain/loss parameters. (a) α = 1 , (cid:15) = 0, (b) α = 0 . , (cid:15) = 0, (c) α = 0 . , (cid:15) = 0,(d) α = 0 . , (cid:15) = 0, (e) α = 1 , (cid:15) = 0 .
5, (f) α = 0 . , (cid:15) = 0 .
5, (g) α = 0 . , (cid:15) = 0 .
5, (h) α = 0 . , (cid:15) = 0 .
5, (i) α = 0 . , (cid:15) = 0 . C ( t ) , C ( t ) for the fractional PT dimer system (19). Here, U ≡ unbounded,B ≡ bounded, and z = ( (cid:15) − V ) t α . The dots denote the position of our two phases, phase( (cid:15) − V ) = 0 for (cid:15) − V >
0, andphase( (cid:15) − V ) = π for (cid:15) − V < B. The nonlinear PT dimer We now explore a PT dimer in the presence of nonlinearity, and subject to fractional evolution equations i d α C ( t ) dt α + i(cid:15)C ( t ) + V C ( t ) + χ | C ( t ) | C ( t ) = 0 i d α C ( t ) dt α − i(cid:15)C ( t ) + V C ( t ) + χ | C ( t ) | C ( t ) = 0 (23)In the absence of PT symmetry ( (cid:15) = 0), and for order α = 1, eqs.(23) have been explored before in the literature[31–33]. For initial conditions C (0) = 1 , C (0) = 0, it was shown that they lead to the phenomenon of a seltrappingtransition: The existence of a critical nonlinearity parameter χ c /V = 4 below which, the long-time average of thesquare of the amplitudes, (cid:104)| C , | (cid:105) = (1 /T ) (cid:82) T | C , | dt (with T (cid:29)
1) is the same: (cid:104)| C | (cid:105) = (cid:104)| C | (cid:105) = 1 /
2. Atnonlinearity values above χ c /V , (cid:104)| C | (cid:105) increases past 1 / χ/V values, while (cid:104)| C | (cid:105) | C | ✏ = 0 . = 0 ↵ = 0 . ↵ = 0 . ↵ = 0 . ↵ = 0 . ↵ = 0 . | C | = 0 ✏ = 0 . ↵ = 0 . ✏ = 0 . | C | = 1 | C | = 5 = 6 = 11 = 17 ↵ = 0 . ✏ = 0 . ( a )( c ) ( b )( d ) < | C > ( e ) ↵ = 0 . ✏ = 0 . ↵ = 0 . ✏ = 0 . Figure 5. Dimer amplitude at initial site | C ( t ) | for the nonlinear PT case for several fractional derivative orders α , variousgain/loss parameters (cid:15) and different nonlinearities χ . Lower right: Example of a time-averaged, trapped fraction at initial siteas a function of the nonlinearity parameter ( (cid:15) = 0 . , α = 0 . decreases towards zero. The trapped fraction at the initial site, (cid:104)| C | (cid:105) , increases abruptly as the critical nonlinearityis crossed.For a fractional order derivative (0 < α < PT symmetry, we resort to the Grunwald algorithm[34] to compute the time evolution of C ( t ) , C ( t )for initial conditions C (0) = 1 , C (0) = 0. This approach is based on finite differences, and in our case leads to thefollowing difference equations: X n +1 = n +1 (cid:88) ν =1 Φ αν X n +1 − ν + ih ( Y n + i(cid:15)X n + χ | X n | X n ) + r αn +1 X Y n +1 = n +1 (cid:88) ν =1 Φ αν Y n +1 − ν + ih ( X n − i(cid:15)Y n + χ | Y n | Y n ) + r αn +1 Y (24)where, X ≡ C , Y ≡ C , and Φ αν = ( − ν − (cid:18) αν (cid:19) r αν = ν − α Γ(1 − α ) . (25)Numerical results are shown in Fig.5. In panel (a) we show the behavior of | C ( t ) | in the linear limit ( χ = 0), fora fixed α and several different gain/loss parameter values. We see that the effect of increasing (cid:15) is to augment theamplitude and decrease the frequency of the oscillation. When (cid:15) approaches V , the amplitude grows unbounded andthe oscillation stops. In panel (b) we take χ = 0 as before, with a fixed gain/loss (cid:15) and for several order α values. Aswe noticed before, the presence of 0 < α < | C ( t ) | . If we reduce now thevalue of α , we see a further decrease of the oscillation amplitude, with little effect on the frequency. We now move tothe nonlinear case. In panels (c, d) we show the behavior of | C ( t ) | in time for fixed α, (cid:15) parameters and for several χ values. Roughly speaking, what we observe here is that there seems to exist a special nonlinearity value below whichthe curves decrease steadily to zero at long times, and above which they approach a constant nonzero value in time. Tohelp understand this, we show in Figure 5e (cid:104)| C | (cid:105) = (1 /T ) (cid:82) T | C | dt ( T (cid:29) χ c whose precise value depends on α , there is an abrupt increase in (cid:104)| C | (cid:105) signaling a seltrappingtransition, like in the standard ( α = 1) nonlinear dimer. What is interesting though, is that if we continue increasingthe nonlinearity past the critical point, the trapped fraction begins to decrease instead of increasing towards unity asin the standard case. This fragility of the trapping could perhaps be a manifestation of the tendency of α to decreasethe amplitude of oscillations in general. Thus, what we are seeing here is the interplay of two opposing tendencies:Trapping by nonlinearity and amplitude decay by α . III. CONCLUSIONS