Self-accelerating beam dynamics in the space fractional Schrödinger equation

Self-accelerating beams are fascinating solutions of the Schrödinger equation. Thanks to their particular phase engineering, they can accelerate without the need for external potentials or applied forces. Finite-energy approximations of these beams have led to many applications, spanning from particle manipulation to robust in vivo imaging. The most studied and emblematic beam, the Airy beam, has been recently investigated in the context of the fractional Schrödinger equation. It was notably found that the packet acceleration decreased with a reduction in the fractional order. Here, I study the case of a general n th-order self-accelerating caustic beam in the fractional Schrödinger equation. Using a Madelung decomposition combined with the wavelet transform, I derive the analytical expression of the beam’s acceleration. I show that the nonaccelerating limit is reached for inﬁnite phase order or when the fractional order is reduced to 1. This work provides a quantitative description of self-accelerating caustic beams’ properties.


I. INTRODUCTION
The Schrödinger equation has been at the heart of quantum mechanics for nearly a century.However, it is only two decades ago that this fundamental equation of physics has been extended to fractional calculus [1], thanks to Laskin [2][3][4].He generalized Feynman's path integral formulation of quantum mechanics to Lévy flights, i.e. beyond Brownian motion or Wiener stochastic processes that are based on usual Gaussian statistics.This led to establish the basis of fractional quantum mechanics and to derive the space fractional Schrödinger equation (SFSE).Fractional derivatives in partial differential equations have a long history in accurately modeling a wide range of physical systems, where their integer-order counterparts failed [5].This includes systems experiencing anomalous diffusion, like fluids in heterogeneous porous media with long range spatial correlation decaying as a power law [6].It also concerns several linear or nonlinear viscoelastic and wave propagation problems (acoustic, water...) in lossy media, or even industrial controllers for vehicule suspensions [7].Following Laskin's work, various schemes for a physical implementation of the SFSE have been discussed [8][9][10].
Amongst the many solutions of the Schrödinger equation, the Airy beams discovered in 1979 by Berry and Balazs [11] share a particular place.These intriguing wave packet solutions notably possess the apparent property of accelerating without external potentials or applied forces.They were long considered to be a mathematical curiosity until their physical realisation, using approximations of Airy beams with a finite energy [12][13][14].They were subsequently observed in platforms other than optically-based ones, using e.g.electron beams [15] or surface plasmon polariton [16], and used in applications such as optical particle manipulation [17], light sheet mi-Fundamentally, accelerating beams can be understood in the broader context of the catastrophe theory, introduced by Thom [23,24].Their caustic properties arise from diffraction integrals which define the different sets of catastrophe in Thom's theory.Aside from the widely studied Airy beam that connects from the fold catastrophe, only little attention was brought to caustic beams of higher order, such as Pearcey and Swallowtail beams [25][26][27].These beams share similar accelerating properties but with a different acceleration parameter.Recently, Airy beams physics has been investigated in the context of the SFSE [28][29][30].For a free particle, it has been qualitatively shown that the acceleration of the Airy packet decreases as the fractional index is reduced.These last developments invite for a more general and quantitative description of self-accelerating beams in the SFSE.
In this paper, I study the properties of general onedimensional self-accelerating beams in the context of the SFSE.The wave packet acceleration along with the peaks' trajectory is obtained analytically using a simple method based on the phase decomposition of the wave function in momentum space [31].This method, combined with the use of the wavelet transform (WT), allows to understand self-accelerating beams as the result of a self-interference of the wave function, and to monitor the trajectory of their individual modes.This effect was notably found, following a similar approach, in excitonpolaritons [32], atomic condensates [33] or to be at the origin of the formation of nonlinear X-waves [34].
The paper is organised as follows.In Section II, I introduce the formalism of the SFSE and I derive a general wave packet solution.In Section III, I introduce self-accelerating beams as derived from catastrophe theory.In Section IV, I present a detailed analysis of the Airy beam as solution of the SFSE.In Section V, I extend these results to the case of a general self-accelerating beam within the SFSE, deriving the general wave packet acceleration.Finally, Section VI concludes the paper.
The prefactor D α is the quantum diffusion constant with dimension [energy] 1−α [lenght] α [time] −α .With α = 2 and D α = 1/2m, m being the mass of the particle, one recovers the conventional one-dimensional Schrödinger equation.Using the definition of the Riesz derivative [35], the SFSE can be more conveniently written in momentum space as with = 1.In Eq. ( 1), the fractional derivative in real space is intrinsically a non-local operator [35].However, regarded from a solid-state physics point of view, it simply corresponds to a modification of the dispersion relation in momentum space, here defined as 36].It spans from a parabola (α = 2) to a symmetric linear dispersion (α → 1), much as the one experienced by conduction electrons on a Dirac cone [37] or by Bogoliubov excitations in an interacting Bose-Einstein condensate [38].
For a given initial condition ψ(k, 0), the solution of Eq. ( 2) is obtained by simple integration: In his seminal article [3], Laskin solved Eq. ( 2) for a symmetric Lévy stable distribution, or "Lévy packet", that can be expressed analytically through its characteristic function, i.e. its Fourier transform, as where a µ is a constant scaling the amplitude of the distribution's (Pareto) tails.Analytical expressions in position space for Eq. ( 4) only exist for the case µ = 2 where L 2 (x) is a Gaussian distribution, and for µ = 1 where L 1 (x) is a Lorenztian distribution.With such an initial condition, the solution of Eq. (2) for a Lévy packet reads As I discuss in the following, this choice of the initial condition does not influence the mode propagation, which here only depends on the properties of the dispersion, but only the mode population through the spread of the wave function in momentum space.8) and derived from the dispersion relation only.Supplemental Movie S1 provides an animation of the Gaussian wave packet dynamics with its WT [39].
I first consider the simple case of a Gaussian wave packet, i.e. with µ = 2 in Eq. ( 4), as initial condition: ).The solution of Eq. ( 2) is now To analyse and understand the individual mode properties of the wave packet solution given in Eq. ( 6), I employ the method developed in Ref. [31].It consists in studying the phase properties of the wave function ψ(k, t), since only a phase dynamics takes place in momentum space.The wave function is decomposed as , with an amplitude √ N and a phase term φ.This is analogous to a Madelung transformation in position space, notably performed in hydrodynamics studies [40] and where the gradient of the (real space) phase corresponds to the fluid velocity.In the present case, the initial condition (either Gaussian or Lévy) does not have a complex phase, so the whole phase dynamics arises as a consequence of the dispersion relation The gradient of this phase, with respect to k, gives the wave packet's mode displacement [31] This quantity has a simple interpretation, that a particle in a given mode k will propagate at a distance d (α) after a time t.This has been equivalently introduced by Laskin as the expectation value of the space position for a single mode k 0 of the dispersion relation [3]: using the k-dependent group velocity dispersion Like for the phase in Eq. ( 7), the mode displacement in Eq. ( 8) is independent of the initial condition.The choice of ψ(k, 0) only governs the population of the modes, but not their propagation.A Gaussian distribution has a well-defined variance σ 2 k so the populated modes are then well-located in a certain range of momentum, whereas a more general Lévy stable distribution has no defined variance (for µ < 2) due to the presence of Pareto tails.The dynamics of |ψ(x, t)| 2 for a Gaussian wave packet in the SFSE of order α = 1.5 is shown in Fig. 1(a).Unlike the conventional case of α = 2, the packet splits into two entities at long times, which is a direct consequence of the modified dispersion.
To better visualise and understand the effect of the mode displacement d (α) (k, t), I use the wavelet transform (WT) which permits an alternative representation of the wave function in both position (x) and momentum (k).The WT is defined as [41] (11) where I use Gabor wavelet family which are Gaussian functions with an internal frequency w G .The choice of the Gabor wavelet is arbitrary but here natural since Gaussians are simple wave packet solutions of the SE.One could however use other wavelet families (Morlet, Mexican hat etc.) to obtain equivalent results.I apply the WT to the Gaussian packet previously calculated and show its corresponding wavelet energy densities |W(x, k)| 2 at different times of the evolution in Fig. 1(b)-(d).The density is well-fitted by the mode displacement d (α) (k, t), here shown as dashed green lines.This representation permits an efficient way to visualize, at a given instant of time, both the population of the modes and the distance they have propagated.As I will show in Section IV, the mode displacement can be strongly affected by the presence of a complex phase in the initial condition, leading to interesting transient effects on the wave packet dynamics.

III. SELF ACCELERATING BEAMS
The first self-accelerating beam, the Airy beam, was discovered by Berry and Balazs as a solution of the free one-dimensional Schrödinger equation, using an Airy function as initial condition [11].However, these beams are not physical solutions as they do not possess a squareintegrable wave function, i.e they would require an infinite energy to maintain their properties as they propagate.Physical approximations of Airy beams have been experimentally realised by including an exponential cutoff on the beam's tail [12], in order to ensure its squareintegrability.The initial condition thus becomes where a controls the exponential cut-off of the wave function and b the width of the peaks of the Airy function.
Other examples of the family of caustic beams that emerge in the catastrophe theory can be constructed from canonical diffraction integrals of codimension K [23,24] with the associated potential functions These integrals here depend on one active variable u and on a certain number of control parameters r n .The first example with K = 1 corresponds to the fold catastrophe, and it is related to the (rescaled) Airy function as The second example with K = 2 corresponds to the cusp catastrophe and defines the Pearcey function An interesting application of the fold and cusp catastrophes has been recently found as quantum caustics in spin chains systems [42].The third integral with K = 3 corresponds to the swallowtail catastrophe e i(u 5 +u 3 z+u 2 y+ux) du = Sw(x, y, z) .
(18) One could also construct the four-control parameters integral ξ 4 , from the Butterfly catastrophe.Thom's theory shows that it only exists seven fundamental types of catastrophe, the three remaining ones, hyperbolic, elliptic and parabolic umbilic catastrophes, are obtained with two active variables [23].While multi-dimensional wave packets can be difficult to engineer experimentally, for obvious reasons, onedimensional versions of high-order catastrophe can be more easily realised.For instance, one-dimensional versions of a finite-energy Pearcey beam and its dual selfaccelerating properties have been recently studied [27].The square-integrability of the beam was notably ensured by a symmetrical cut-off on the initial condition, which becomes Pe(x, 0) exp(−ax 2 ).One can then similarly construct the next order of one-dimensional accelerating beams, that is a 1D-Swallotail beam Sw(x, 0, 0), or any beam of order n arising from the following integral The integral defined in Eq. ( 19) can be notoriously difficult to compute, due to its highly oscillatory integrand.Different methods and algorithms based on sets of differential equations or using contour integration techniques have been developed to evaluate such cupsoid canonical integrals [43][44][45].
The function ξ n (x) is here defined for an integer-order n, but it worth mentioning that such a function could also be defined through a fractional order in its integral representation, similarly to so-called fractional Airy beams [46].

IV. AIRY BEAMS IN THE FRACTIONAL SCHR ÖDINGER EQUATION
In this Section, I solve the SFSE for a finite-energy Airy beam as initial condition, see Eq. ( 13), which can be written in momentum space as The solution of Eq. ( 2) is now (21) The real-space solution ψ To understand the wave packet dynamics, I follow the same method as in Section II, first by taking the complex phase of the solution in Eq. ( 21) from which one obtains the mode displacement ) I also compute the WT for the Airy beams with the different fractional orders and show their corresponding wavelet energy densities at three selected times of the evolution in Fig. 2(d)-(l).The mode displacement d (α) Ai previously calculated is plotted on top of the wavelet energy density as orange and purple dashed lines, with again an excellent agreement.
As shown in Ref. [31], the presence of fringes in the Airy beam's density |ψ The displacement of the extremum mode d Ai (k ext ) is found by substituting back Eq. ( 24) into Eq.( 23), which gives The point Ai (k ext ), k ext } is marked as a red dot in the WT panels of Fig. 2, and a green line indicates its trajectory in the x-k phase space.At long times, the self-interference vanishes as the point P I drifts away to the region of large momenta, where the wavelet energy density is zero, leaving no signal to participate into the self-interference.This process is faster as α is large, as seen from Fig. 2(d)-(l) or from Eq. (25).Physically, this means that the fringes in the beam's density "survive" a longer time as α is reduced to 1.

Since d (α)
Ai (k ext ) measures the distance traveled by the extremum mode, it corresponds to the beam's wave front trajectory in real space.Trajectories for the three cases of Airy beams in Fig. 2(a)-(c) and computed from Eq. ( 25) are reported as dashed blue lines.For the case of α = 2 in Eq. ( 25), one recovers the usual t 2 parabolic acceleration of the Airy beam.
The last discussion of this Section will concern the splitting of the wave packet at long times, observed either with the Gaussian packet (see Fig. 1) or the Airy beam (see Fig. 2) and that has been previously discussed [28,29,47] .It is worth looking at the effect of lowering the fractional order on the far-field density, representing the wave function in both momentum k and energy E. The density |ψ for the three different Airy cases I have considered so far.In this phase space, the density follows the dispersion relation E (α) (k), here plotted as a red line.One can appreciate its deviation with the parabolic dispersion of the Schrödinger equation, shown as a dashed blue line.
The splitting of the beams at long times is only the consequence of the linearisation of the dispersion relation, no matter what the initial condition is.Indeed, even when it contains a complex phase, as for the Airy beam, at long times the effect of the dispersion is dominant on the mode propagation.In Eq. ( 23), the only time-dependent term is the one that derives from the dispersion, on the LHS.As the dispersion becomes linear (E (α→1) = |k|), only two group velocities remain accessible, that are v ± = ±αD α .So if the initial wave function ψ(k, 0) is distributed around k = 0, at long times the packet inevitably splits into two parts, travelling with the same velocity but in opposite directions.In the wavelet picture, this would translate as two vertical lines for the mode displacement on the x-k phase space, that is also the signature of non-diffusing wave packets [31].

V. SELF ACCELERATING BEAMS IN THE FRACTIONAL SCHR ÖDINGER EQUATION
Let's now generalise the results obtained in the previous Section for the Airy beam, to a general 1D n th -order self-accelerating caustic beam, as defined by the integral in Eq. (19).The method employed in this paper to calculate the mode displacement (see Eqs. ( 8) and ( 23)), which leads to the expression of the acceleration, requires the knowledge of the initial condition in momentum space.Fortunately, the integral ξ n (x, 0) possesses a handily Fourier transform.Indeed, one can find that as shown in Appendix A.
It is instructive to first examine the symmetry of the function in Eq. ( 26) in order to deduce some of the beams' properties.When n is even (n = 4, 6, 8...), ξn (k) possesses both an even real and an even imaginary part, which leads to ξ n (x) having also an even real and an even imaginary part.This follows that the density |ξ n (x)| 2 is an even function when n is an even number.Such wave functions, like the Pearcey beam (n = 4), thus possess a symmetrical double-accelerating wave front [31].On the other hand, when n is odd (n = 3, 5, 7...), ξn (k) possesses an even real part and an odd imaginary part, which leads ξ n (x) to be a purely real function with an even real part and an odd real part.The density |ξ n (x)| 2 thus has no particular symmetry when n is an odd number.This is why the wave functions of the Airy or the Swallowtail beams (n = 3 and 5) only possess a single accelerating wave front.
With ξn (k) as initial condition [48], the solution of the SFSE (see Eq. ( 2)) now reads with some amplitude √ N 0 and the associated phase The corresponding mode displacement, obtained by differentiating the previous phase with respect to k, is The extremum mode is obtained by solving n (k, t) derived from Eq. ( 29).The point PI , around which the self-interference occurs, is shown as a red dot, and the green line stands for its trajectory in the x-k phase space.Both packets accelerate as t 4/3 .Supplemental Movie S3 provides an animation of the Pearcey and Swallowtail beams dynamics with their WT [39].
Substituting back Eq.(30) into Eq.( 29), one obtains the mode displacement of the extremum mode which corresponds to the wave front acceleration.From the previous formula, one finds that the wave packet accelerates as t β , with β = (n − 1)/(n − α).The acceleration parameter β as function of the fractional order α is plotted in Fig. 3 for different beam orders n.The values at α = 2 and marked by a dot correspond to the acceleration obtained for the different beams in the conventional Schrödinger equation.One thus retrieves the previously derived t (n−1)/(n−2) acceleration law for general 1D caustic beams [31].
From Eq. ( 31), one finds that the non-accelerating limit (β → 1) is reached either when α → 1 or when n → ∞.It also follows from Eq. ( 31) that it exists different ways to obtain a given acceleration with different types of caustic beams.A given acceleration, for example β = 4/3, can be realised in three ways, (1) with an Airy beam (n = 3) and a fractional order α = 1.5, (2) with a Percey beam (n = 4) and a fractional order α = 1.75, or (3) with a Swallowtail beam (n = 5) and a fractional order α = 2.This Airy beam configuration was already shown in Fig. 2(b).Cases of the Pearcey and Swallowtail beams showing the same acceleration are presented in Fig. 4. One can note the symmetrical double-accelerating wave front in the Pearcey beam.
effect of vanishing self-interference is particularly well visible in the supplemental movie S3 that timeanimates Fig. 4. Indeed, the Pearcey beam possesses here a higher spread in momentum space than the Swallowtail beam, due to their different cut-offs.The selfinterference, and thus the fringes in the real space density, are preserved for longer times in the Pearcey beam as the point P I with coordinates {d

VI. CONCLUSION
In conclusion, I have studied general one-dimensional self-accelerating caustic beams in the SFSE.I showed that the previously observed effect of the wave packet splitting at long times is a natural consequence of the linearisation of the dispersion relation, which occurs when the fractional order is reduced down to 1. Using a method of analysis based on a combination of spectral techniques, I have calculated the analytical expression for the mode propagation, from which the wave packet acceleration is derived.In both the limits of an infinite caustic order or a fractional order of 1, the beams loose their accelerating property.This method of analysis could be employed in the study of other accelerating beams, or beams with a different phase engineering, and with other types of dispersion relations, such as with coupled or periodic systems.

FIG. 1 .
FIG. 1. Diffusion of a Gaussian wave packet in the SFSE of order α = 1.5, and with initial width σ k = 1/1.5.(a) Wave function density |ψ(x, t)| 2 .(b)-(d) Wavelet energy density |W(x, k)| 2 computed at selected times.The dashed green lines stand for the mode displacement d (α) (k, t) computed from Eq. (8) and derived from the dispersion relation only.Supplemental Movie S1 provides an animation of the Gaussian wave packet dynamics with its WT[39].

FIG. 2 .
FIG. 2. Airy beam dynamics for different fractional orders α.Each row corresponds to a different value of α. (a)-(c) Wave function density |ψ (x, t)| 2 .The dashed blue line shows the wave packet front trajectory, computed from the extremum mode displacement in Eq. (25).(d)-(l) Corresponding wavelet energy density |W(x, k)| 2 at selected times.The dashed purple and orange lines are the mode displacement d (k, t) computed from Eq. (23).The point PI , around which the self-interference occurs, is shown as a red dot, and the green line stands for its trajectory in the x-k phase space.(m)-(o) Far-field density |ψ (k, E)| 2 .The field's density follows the dispersion relation E (α) (k), shown as a solid red line.The parabolic dispersion E (α=2) (k) is plotted as a dashed blue line for comparison.Parameters are as follows: a = 0.01, b = 0.3, Dα = 0.5.Supplemental Movie S2 provides an animation of the different Airy beam dynamics with their WT [39].
t) is simply obtained by inverse Fourier transform.Typical dynamics for |ψ(α) Ai (x, t)| 2 are shown in Fig. 2(a)-(c) for three different fractional orders α.As the fractional order is reduced, the wave front acceleration diminishes, and at long times, the packet undergoes the same reshaping as the one observed in the Gaussian case, see Fig. 1.
(x, t)| 2 arises from a transient self-interference of the wave function.The mode displacement d (α) Ai is a multivalued function in the x-k phase space.A self-interference thus occurs when the wavelet energy density spreads around the point where d (α) Ai becomes multivalued.Indeed, for a given position x, the wave function can contain two different modes k which overlap in real space and lead to the self-interference.The self-interference occurs around a point P I with coordinates {d (α) Ai (k ext ), k ext }, where k ext is the momentum of the extremum mode of the branch, i.e. the mode with the largest displacement, or also max[d (α) Ai ].The extremum mode is found solving ∂ k d (α) Ai (k, t) = 0 for k, which gives

FIG. 3 .
FIG. 3. Exponent of the wave front acceleration t β as function of the fractional order α for self-accelerating beams of different orders n.The three cases with n = 3, 4 and 5 correspond to an Airy, a Pearcey and a Swallotail beam, respectively.These beams can reach the same acceleration for a different fractional order α, for example β = 4/3, see the dashed grey line.The dots at α = 2 indicate the acceleration obtained with the conventional Schrödinger equation.The non-accelerating limit (β → 1) is reached either when α → 1 or n → ∞.

FIG. 4 .
FIG. 4. High-order accelerating beams dynamics.The first column corresponds to the case of a Pearcey beam with n = 4 and α = 1.75, and the second column to a Swallowtail beam with n = 5 and α = 2. (a)-(b) Wave function density |ψ (α) n (x, t)| 2 .The dashed blue lines indicate the trajectory of the extremum mode, derived from Eq. (31).(c)-(d) Wavelet energy density |W(x, k)| 2 computed at t = 150.The dashed orange and purple lines stand for the mode displacement d (α)

n
(k ext ), k ext } drifts to large momenta.