The fractional nonlinear PT dimer

We examine a fractional Discrete Nonlinear Schrodinger dimer, where the usual first-order derivative 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 the Laplace 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 the direct explicit Grunwald algorithm. In general, the main effect of the fractional derivative is the onset of a monotonically decreasing time envelope for the amplitude of the oscillatory exchange. In the presence of PT symmetry, the dynamics shows damped oscillations for small gain/loss in both sites, while at higher gain/loss parameter values, the amplitudes of both sites grows unbounded. In the presence of nonlinearity, selftrapping is still possible although the trapped fraction decreases as the nonlinearity is increased past threshold, in marked contrast with the standard case.


I. INTRODUCTION
The topic of fractional calculus has experienced a rekindled interest in recent times. Essentially, it extends the notion of a derivative or an integral of integer order, to one of a fractional order, (d n /dx n ) → (d α /dx α ) for real α. The subject has a long history, dating back to letters exchanged between Leibnitz and L'Hopital, and later contributions by Euler, Laplace, Riemann, Liouville, and Caputo to name some [1][2][3][4][5]. The starting point was the computation of d α x k /dx α , where α is a non-integer number: For instance, (d 1/2 /dx 1/2 )x = (2/ √ π) √ x, and dx/dx = (d 1/2 /dx 1/2 )(d 1/2 /dx 1/2 )x = (2 √ π)(Γ(3)/Γ(1))x 0 = 1, as expected. From Eq.(1) the fractional derivative of an analytic function f (x) = k a k x k can be computed by deriving term by term. This basic procedure is not exempt from ambiguities. For instance, (d α /dx α ) 1 = (d α x 0 /dx α ) = (1/Γ(1 − α))x −α = 0, according to Eq.(1). However, one could have also taken (d α−1 /dx α−1 )(d/dx) 1 = 0. For the case of a fractional integral, a more rigorous starting point is Cauchy's formula for the integral of a function. From the definition we apply the Laplace transform L to both sides of Eq.(2) After n integrations, one obtains Extension to fractional α is direct: After noting that the RHS of Eq.(5) is the product of two Laplace transforms we have, after using the convolution theorem From this definition, it is possible to define the fractional derivative of a function f (x) as where, m − 1 < α < m. Eq. (7) is known as the Riemann-Liouville form. An alternative, closely related form, is the Caputo formula [5]: which has some advantages over the Riemann-Liouvuille form for differential equations with initial values. The various technical matters that arise in fractional calculus have prompted a whole line of research that has extended to current times. Long regarded as a mathematical curiosity, it has now regained interest due to its potential applications to complex problems in several fields: fluid mechanics [6,7], fractional kinetics and anomalous diffusion [8][9][10], strange kinetics [11], fractional quantum mechanics [12,13], Levy processes in quantum mechanics [14], plasmas [15], electrical propagation in cardiac tissue [16] and biological invasions [17]. In general, fractional calculus constitutes a natural formalism for the description of memory and non-locality effects found in various complex systems. Experimental realizations are not straightforward given the nonlocal character of the coupling, however some optical setups have been 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 continuum Schrödinger equation, or the paraxial wave equation. The effective discreteness comes from expanding the solution sought in terms of (continuous) modes that can be labelled unambiguously. The simplest of such examples is the bonding, anti-bonding electronic mode that one finds for a two-sites (dimer) molecule after diagonalizing the two-site Schrödinger equation in the tight-binding approach. Something similar happens in optics, where the paraxial equation is formally equivalent to the Schrödnger equation. In that case, for two optical waveguides, the total electric field is expanded in terms of the electromagnetic modes in each guide which interact through the evanescent field between the 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 the Schrödinger equation for N units [20,21]. Of course, at the end one has to collect all the discrete amplitudes and multiply them by the corresponding continuous mode profiles and superpose them, to obtain the final field. The simplest 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][23][24], and more recently, they have been considered in model of correlated disorder [26] and in magnetic metamaterial modeling [25].
In this work we consider the discrete Schrödinger equation for a dimer system, where the standard time derivative is replaced by a fractional one. The dimer considered is rather general and contains asymmetry, PT symmetry and nonlinearity ( Fig.1). Our main interest is in ascertaining the effect of the fractional derivative on the excitation exchange 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 where α is the fractional order of the (Caputo) derivatives with 0 < α < 1. Quantities C 1,2 are probability amplitudes in a quantum context, or electric field amplitudes, in an optical setting. Parameter V is the coupling term and χ is the nonlinearity parameter. Let us first consider the case of a general linear (χ = 0) dimer, and assume C 1 (0) = 1, C 2 (0) = 0. We will solve this case in closed form by the use of Laplace transforms: For 0 < α < 1, the Laplace transform of the Caputo fractional derivative of order α is given by After applying the Laplace transform L to both sides of Eq.(9) we have Solving for L(C 1 ) and L(C 2 ) gives: and Using the inverse Laplace formula[27] we obtain: where E γ α,β (z) is defined as where (γ) n = Γ(γ + n)/Γ(γ), and α, β, γ ∈ C, Re(α) > 0, Re(β) > 0, z ∈ C. Figure 2 shows examples of the time evolution of the square of the dimer amplitudes, for several site energy parameters, and fractional derivative orders. In general we observe that, as soon as α differs from unity, the dynamics is either bounded or unbounded, depending on the values of the site energy parameters. For the bounded cases, there is some oscillation initially, with a decreasing envelope towards zero.

A. The linear PT dimer
A particularly interesting case of Eq. (9), is the fractional PT -symmetric dimer. For systems that are invariant under the combined operations of parity (P) and time reversal (T ), it was shown that they display a real eigenvalue spectrum, even though the underlying Hamiltonian is not hermitian [29,30]. In these systems there is a balance between gain and loss, leading to a bounded dynamics. However, as the gain/loss parameter exceeds a certain value the system undergoes a spontaneous symmetry breaking, where two or more eigenvalues become complex. At that point, the system loses its balance and its dynamics becomes unbounded. According to the general theory, for our system to be PT symmetric, the real part of the site energies in Eq.(9) must be even in space while the imaginary part must be odd: Re( 1 ) = Re( 2 ) and Im( 1 ) = −Im( 2 ). For simplicity we take the real parts of 1 , 2 as zero and thus, 1 = − 2 ≡ i , where is the gain/loss parameter. This leave us with the equations: Dimer amplitudes for the linear case (χ = 0) and several site energy parameters 1, 2 and different fractional derivative orders α. Solid(dashed) line denotes |C1| 2 (|C2| 2 ). whose exact solutions can be extracted from the general solution, Eqs.(15),(16) as where, E α,β (z) = E 1 α,β (z) is known as the generalized Mittag-Leffler function 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 time evolutions for |C 1 | 2 , |C 2 | 2 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, converging to zero at long times. The asymptotic behavior of C 1 (t), C 2 (t) depends on the behavior of the Mittag-Leffler functions E α,β (z) at large values of |z|. After writing z = |z| exp(iφ), we have [28] where Q = z 1/α = exp((1/α) log(|z|) + i φ) and |φ/α| ≤ π (φ = 2 − V 2 ). This implies, exp(Q) = exp(|z| 1/α cos((1/α)φ)) × exp(i |z| 1/α sin((1/α)φ)) Thus, bounded behavior in time will occur for (π/2) < |φ/α| < π, while unbounded behavior occurs for 0 < |φ/α| < π/2. In our case, φ = arg( 2 − V 2 ) = 0, π, implying that C 1 (t) and C 2 (t) will increase (decrease) asymptotically in time if ( 2 − V 2 ) is positive (negative). This behavior is sketched in Fig.4.   . Asymptotic stability for the amplitudes C1(t), C2(t) for the fractional PT dimer system (19). Here, U ≡ unbounded, B ≡ bounded, and z = ( 2 − V 2 )t 2α . The dots denote the position of our two phases, phase( 2 − V 2 ) = 0 for 2 − V 2 > 0, and phase( 2 − V 2 ) = π for 2 − V 2 < 0.

B. The nonlinear PT dimer
We now explore a PT dimer in the presence of nonlinearity, and subject to fractional evolution equations In the absence of PT symmetry ( = 0), and for order α = 1, eqs. (23) have been explored before in the literature [31][32][33]. For initial conditions C 1 (0) = 1, C 2 (0) = 0, it was shown that they lead to the phenomenon of a seltrapping transition: The existence of a critical nonlinearity parameter χ c /V = 4 below which, the long-time average of the square of the amplitudes, |C 1,2 | 2 = (1/T ) decreases towards zero. The trapped fraction at the initial site, |C 1 | 2 , increases abruptly as the critical nonlinearity is crossed.
For a fractional order derivative (0 < α < 1), where we take the Caputo version of the fractional derivative, and in the presence of PT symmetry, we resort to the Grunwald algorithm [34] to compute the time evolution of C 1 (t), C 2 (t) for initial conditions C 1 (0) = 1, C 2 (0) = 0. This approach is based on finite differences, and in our case leads to the following difference equations: where, X ≡ C 1 , Y ≡ C 2 , and Numerical results are shown in Fig.5. In panel (a) we show the behavior of |C 1 (t)| 2 in the linear limit (χ = 0), for a fixed α and several different gain/loss parameter values. We see that the effect of increasing is to augment the amplitude and decrease the frequency of the oscillation. When approaches V , the amplitude grows unbounded and the oscillation stops. In panel (b) we take χ = 0 as before, with a fixed gain/loss and for several order α values. As we noticed before, the presence of 0 < α < 1 induces a decreasing oscillation behavior in |C 1 (t)| 2 . If we reduce now the value of α, we see a further decrease of the oscillation amplitude, with little effect on the frequency. We now move to the nonlinear case. In panels (c, d) we show the behavior of |C 1 (t)| 2 in time for fixed α, parameters and for several χ values. Roughly speaking, what we observe here is that there seems to exist a special nonlinearity value below which the curves decrease steadily to zero at long times, and above which they approach a constant nonzero value in time. To help understand this, we show in Figure 5e |C 1 | 2 = (1/T ) T 0 |C 1 | 2 dt (T 1), the time-averaged fraction remaining at the initial site vs the nonlinearity strength. As soon as the nonlinearity increases from zero there is a finite amount of trapping at the initial site that increases monotonically with nonlinearity. As the nonlinearity parameter reaches a critical value χ c whose precise value depends on α, there is an abrupt increase in |C 1 | 2 signaling a seltrapping transition, like in the standard (α = 1) nonlinear dimer. What is interesting though, is that if we continue increasing the nonlinearity past the critical point, the trapped fraction begins to decrease instead of increasing towards unity as in the standard case. This fragility of the trapping could perhaps be a manifestation of the tendency of α to decrease the 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
We have examined the excitation dynamics in a nonlinear PT dimer when the evolution equations are ruled by a fractional-order time derivative, instead of the usual first-order. In the absence of gain/loss and nonlinearity, the effect of the fractional derivative alone is to induce a damping of the oscillatory exchange between the two sites. When PT symmetry is added, two oscillation regimes are present, a damped oscillatory regime for small gain/loss parameter, and a monotonic growth for large gain/loss parameters. However, there is no stable oscillatory regime, i.e., when there is complete balance between gain and losses, unless α = 1, the standard case.
Finally, when nonlinearity is added to the picture, we observe selftrapping at long times at the initial site, that increases steadily as the nonlinearity reaches a critical value of approximately, χ ∼ 4V . Above this threshold, the trapped fraction at the initial site decreases monotonically as nonlinearity is increased further. This is in marked contrast with the standard case where selftrapping keeps on increasing towards unity with increasing nonlinearity strength.
In spite of these interesting differences with the standard integer derivative case, some general features remain more or less similar to its standard counterpart: (1) There is exchange between the sites, (2) In the presence of PT symmetry there are two clearly marked regimes where the amplitudes decrease to zero or diverge to infinity (3) There is still a selftraping transition, whose finer details depend on the value of the fractional derivative order as well as on the strength of the gain/loss parameter. The persistence of this general phenomenology in the face of a different mathematical rate of change, suggests that this phenomenology is robust against "mathematical perturbations". This kind of robustness has been also observed for the fractional DNLS equation, where the Laplacian is replaced by a fractional one [35,36]. Thus, we expect that concepts of power exchange between sites, selftrapping transitions and existence of nonlinear excitations (discrete solitons) to be concepts of wide physical validity for nonlinear discrete systems.