1. Introduction

The progressively growing interest in the physics of fractional Schrödinger equation (FSE) [1–3], is motivated, on the one hand by their fundamental importance for the fractional field theory and the dynamics of fractional-spin particles [4], and, on the other hand, by the rich possibilities they offer for emulating propagation dynamics of beams in optics [5]. The last property allows one using optics to simulate the fractional quantum harmonic oscillator [5], what would be impossible in direct condensed matter experiments currently owing to the need to specially tailoring non-nearest neighbor hopping in the lattice [6].

In 2000, N. Laskin established the fractional quantum mechanics [1, 2]. He generalized the standard Schrödinger equation (SE) to the fractional Schrödinger equation (FSE) by replacing the second-order spatial derivative with a fractional one. The developed fractional path integral approach to quantum and statistical mechanics can easily be applied to d-dimensional consideration using the d-dimensional generalization of the fractional and the Lévy path integral measures [2]. The FSE describes the physics when the Brownian trajectories in Feynman path integrals are replaced by Lévy flights, which is crucial for understanding the phenomena involving in the fractional effects, including the fractional quantum Hall effect [7], the fractional Talbot effect [8], the fractional Josephson effect [9], and the fractional quantum oscillator [10].

However, the advance in this area is very slow due to the absence of effectively experimental observations. The theoretical difficult arises from the inherent nonlocal operator relevant to the fractional derivative sitting in the FSE [11–15]. Until 2015, by considering the similarity between the SE and the paraxial wave equation, Longhi introduced the FSE into optics and suggested a scheme to emulate the fractional quantum harmonic oscillator by an optical way [5].

This work opened new insights into the propagation dynamics of light beams in fractional dimensions. Even more recently, diffraction-free beams [16], chirped Gaussian beam propagation [17], $\mathcal{P}\mathcal{T}$ symmetry [18], “accessible solitons” (linear modes) [19, 20], and propagation management [21] were investigated in the FSE in linear regimes with or without a potential. Meanwhile, propagation of super-Gaussian beams [22] and gap solitons [23] were reported in the nonlinear fractional Schrödinger equation (NLFSE).

Surface waves are, by their very nature, guided waves propagating along the interface between two different media [24]. Nonlinear surface modes have attracted considerable interest in recent years due to the potential applications such as surface characterization, optical sensing, and switching [25]. They exhibit novel physics of nonlinear modes and have no analogs in schemes without interfaces. The presence of an interface drastically modifies the properties of nonlinear modes. It can enrich the families of solitons, change their existence domains, and alter their stability conditions. Various types of stable solitons at an interface such as, gap solitons [25], mixed-gap solitons [26], defect solitons [27], and disordered solitons [28] were predicted in different configurations.

Bearing in mind that surface solitons can exist in the conventional NLSE [25] and gap solitons are supported by the NLFSE [23], it is natural to ask: Can gap solitons be supported by the NLFSE with an interface? If yes, what is their dynamics and how does the Lévy index influence them? To answer these questions, we, for the first time, study the existence, stability, and propagation of spatial solitons in the NLFSE with an interface. We find that the Lévy index can be effectively used to suppress the instability of surface gap solitons, especially for solitons in higher bandgaps and solitons penetrating deeply in lattice regions.

2. Theoretical model

We start our analysis by considering beam propagation along the interface separating a defocusing uniform medium and a semi-infinite optical lattice imprinted in it. The balance between the fractional-order diffraction effects, nonlinearity, and modulation of refractive-index makes it possible to form surface solitons sitting on the interface. The dynamics of light beams is determined by the dimensionless NLFSE [13, 23]

For simplicity and comparison, we consider a particular form of potential R(x) = p[1−cos(Ωx)] at x ≥ 0 and R(x) = 0 at x < 0, where Ω ≡ 4 is the modulation frequency. The lattice depth p is assumed small compared with the unperturbed index and is of the order of nonlinear correction. A transition between the uniform media and the lattice can be realized by erasing part of the red-light imprinted lattice with an intense green-light wave [25]. The fractional-order diffraction effect may be implemented by the scheme proposed by Zhang et al. [16]. Equation (1) admits several conserved quantities, including the power or energy flow: $U={\displaystyle {\int }_{-\infty }^{\infty }{\left|q\left(x\right)\right|}^{2}dx}$.

Although the surface mode mainly resides at the interface, its wings or oscillatory tails unavoidably penetrate into the bulk of lattice. The existence region of solitons is consequently restricted in the finite gaps of lattice. Thus, before we discuss the properties of nonlinear modes, it is helpful to understand the Floquet-Bloch spectrum of the corresponding linear fractional system.

To obtain the dispersion relation, we search for solutions of the linear version of Eq. (1) in the form q(x, z) = w(x)exp(ikx + ibz), where k is the transverse Bloch wave number, b is the longitudinal wave number, and w(x) = w(x + T) is a periodic function with T being the lattice periodicity. According to the Bloch theorem, one only needs to consider the band-gap structure in the first Brillouin zone –π/T ≤ k ≤ π/T. Expanding w(x) and the potential R(x) in series of plane-waves, w(x) = ∑_n C_n exp(iK_nx) with K_n = 2πn/T and R(x) = ∑_m P_m exp(iK_mx) with $P_m=\frac{1}{T}{\displaystyle {\int }_0^{T}R\left(x\right)\exp \left(-i{K}_{m}x\right)dx}$ and plugging them into Eq. (1), one obtains: ${\displaystyle \sum _n\left[-b-\frac{1}{2}{\left|k+K_n\right|}^{\alpha }\right]C_n\exp \left[i\left(k+K_n\right)x\right]+{\displaystyle \sum _{m,n}{P}_m{C}_n\exp \left[i\left(k+K_n+K_m\right)x\right]=0}}$. Multiplying the above equation by exp[−i(k + K_q)x] and integrating over x ∈ (−∞, ∞), one ends up with an eigenvalue problem:

At a fixed lattice depth, the dependence of band-gap structure on the Lévy index is shown in Fig. 1(a). While the first finite gap shrinks with the decrease of α, the higher gaps expand. All finite gaps shift towards the direction of the semi-infinite gap. At α = 1.2, all bandgaps expand with the growth of lattice depth p [Fig. 1(b)]. It means the existence domain of nonlinear modes in bandgaps can be enlarged by increasing p.

 

Fig. 1 (a, b) Band-gap structure of periodic lattice. p = 4 in (a) and α = 1.2 in (b). Bands are marked with gray color; gaps are shown white. Bandgap spectrum at α = 2 in (c) and 1.2 in (d). p = 4 in (c, d) and Ω = 4 in all the panels.

Download Full Size | PDF

The detailed band-gap spectra of the periodic systems with α = 2 and 1.2 are displayed in Figs. 1(c) and 1(d). While the first gap between the lower edge of the first band and the upper edge of the second band in Fig. 1(d) [(2.619, 5.586)] is narrower than that in 1(c) [(0.282, 4.221)], the second gap of the system with α = 1.2 is obviously wider than that with α = 2. By comparing the corresponding band curves, one immediately finds that the curvature of dispersion curves decreases with the decrease of Lévy index α. When α → 1, the band curves become flatten, which indicates that a beam propagating in the present scheme is approximately diffraction-free, due to d²b/dk² → 0. Thus, the small Lévy index plays a role of suppressing the diffraction effect of the beam upon propagation. This property is very important for the stabilization of the evolution of nonlinear modes.

3. Numerical results and discussions

We now consider nonlinear waves sitting on an interface between a defocusing homogenous medium and a semi-infinite lattice imprinted in it. The stationary solution of Eq. (1) can be searched by assuming: q(x, z) = w(x)exp(ibz), where w(x) is the profile of surface gap soliton and b is a real propagation constant. Substitution of the expression into Eq. (1) results in a nonlinear fractional-order ordinary differential equation:

When both of forward and backward propagating waves in a lattice experience Bragg scattering, the nonlinear coupling between them results in the formation of gap solitons. Surface solitons in the first finite gap are shown in Figs. 2(a)–2(c). Gap solitons can exist only when b > 0 because the tails in the uniform medium exhibit an exponentially decaying behaviour. They are well localized although the nonlinearity is defocusing and the fractional diffraction still expands the beam (α > 1). It is the linear refractive-index modulation who compensates the above effects and affords the localization of solitons. Due to the asymmetry of the media beside the interface, the soliton profiles are also asymmetric. With the decrease of α, soliton peak increases and soliton becomes narrow [Fig. 2(a)]. Unlike the conventional surface gap solitons [25], in the presence of fractional diffraction effect, a dipole-like part appears in the second lattice channel and is in-phase with the main lobe in the first lattice site [Fig. 2(b)].

 

Fig. 2 Profiles of surface gap solitons in the first (a–c) and second (d) finite gaps. (a) b = 3.0. (b) b = 2.8. (c) A twisted soliton at b = 4. (d) b = 15, p = 12. p = 4 in (a–c) and α = 1.2 in (b–d).

Download Full Size | PDF

There exist another type of nonlinear localized modes constituting of two out-of-phase soliton units with almost equivalent amplitudes in the neighboring lattice sites [Fig. 2(c)]. The repulsive interaction between two out-of-phase components may lead to the instability of solitons. Solitons residing in the higher bandgap are also possible. In a deep lattice with p = 12, an example of solitons at b = 15 in the system with α = 1.2 is illustrated in Fig. 2(d). It consists of a main lobe in the first lattice site and an out-of-phase dipole in the second lattice site. One can expect that such solitons are stable, since the destabilization effect of the small-amplitude dipole is very weak.

To illustrate the role played by the Lévy index α, we show the width, peak value, power, and existence and stability regions of gap solitons in systems with α = 2 and 1.2 in Fig. 3. The integral form-factor of a soliton can be defined as $D={\displaystyle \int |x||w{|}^{2}dx/{\displaystyle \int |w{|}^{2}dx}}$, which is proportional to the effective width of solitons. When the propagation constant b is far from the band edges, solitons in systems with α = 2 and 1.2 are all well-localized in their respective gaps [see top plot in Fig. 3(a)]. At fixed b, solitons with a large α are broader than solitons with a small α [bottom plot in Fig. 3(a)]. Due to the nature of defocusing nonlinearity, the peak value of surface solitons decreases with the growth of b. Yet, in the common existence domain, the peak of solitons with a small α is larger than that of solitons with a large α [Fig. 3(b)].

 

Fig. 3 (a) Effective width of solitons versus propagation constant b (top plot) and Lévy index α (bottom). b = 3 in the bottom plot. (b) The peak value of solitons versus b for different α. Dependence of soliton power U versus b in the second gap at α = 2 (c) and 1.2 (d). The lower curves denote single-peaked solitons and the upper ones stand for two-peaked solitons. Solid: stable; dashed: unstable. p = 4 in (a, b) and 12 in (c, d).

Download Full Size | PDF

The stability is a fundamentally important issue in the study of spatial solitons. It can be analyzed by considering the perturbed stationary solution form as q (x, z) = w(x)exp(ibz) + u(x)exp[i(b + λ)z] + v^∗(x) exp[i(b – λ^*)z], where the perturbation components u, v could grow with a complex rate λ during propagation. The soliton is stable if the imaginary parts of perturbation eigenvalues equal zero. Substituting the perturbed solution into Eq. (1), we obtain the coupled fractional eigen-equations:

Surface solitons exist only when their power is beyond a threshold value. This means such nonlinear modes do not bifurcate from the Bloch modes at the band edges. It also indicates there is a lower cutoff of lattice depth p_cr below which no surface modes can be found. In the system with α = 2, surface solitons in the first gap are stable in almost their entire existence domain except for a very narrow region near the upper cutoff of propagation constant, where the power curve satisfies the condition dU/db > 0 [25]. The existence and stability properties of surface solitons in the system with a small α are similar. The difference comes from that the instability region becomes narrower, due to the weak fractional-order diffraction effect.

In the system with α = 2, while surface solitons with one peak in the second gap are still stable in almost their total existence domain, bound states with two in-phase peaks in the neighboring lattice sites exhibit an oscillatory instability when the propagation constant exceeds a certain value [Fig. 3(c)]. A fascinating feature in systems with a small α is that the instability of two-peaked solitons can be completely suppressed by the decrease of Lévy index α [Fig. 3(d)].

Though surface solitons have been intensively investigated in diverse optical schemes, gap solitons with more peaks in neighboring lattice channels near the interface have received little attention. We plot the multi-peaked surface gap solitons residing in the first gap in Fig. 4. Unlike the single-peaked solitons shown in Fig. 2, multi-peaked solitons penetrate into more lattice channels. The decaying tails adjacent to the main peak deviate from the interface. The effective width of solitons increases with the growth of the number of peaks. For solitons with a fixed number of peaks, the difference between the peak values of main peaks becomes undistinguishable with the increase of power U [Figs. 4(a) and 4(b)].

 

Fig. 4 Profiles of three-peaked solitons at b = 3 (a) and α = 1.2 (b). Profile of five-peaked soliton at b = 3.8 (c) and seven-peaked soliton at b = 3 (d) with α = 1.2. Power U versus b for multi-peaked surface gap solitons at α = 2 (e) and 1.2 (f). Solid: stable; dashed: unstable. All solitons reside in the first gap and p = 4 in all the panels.

Download Full Size | PDF

Generally to say, gap solitons with an arbitrary number of peaks can be supported by the present system. One may naturally associate such nonlinear modes with the “truncated” nonlinear Bloch waves localized in the gaps of the linear Bloch-wave spectrum, in Bose-Einstein condensates (BECs) loaded into optical lattices [30] and in optics [31, 32]. We stress that, however, though their appearances are similar, they belong to different families. The nonlinear gap waves in [30–32] originate from the band edges. The surface gap solitons here are purely nonlinear modes with a threshold power. They are not nonlinear waves bifurcating from band edges and cannot penetrate into bands. Another difference comes from that the peaks of “truncated” Bloch waves are of the same value and the peaks of surface solitons are slightly different [Figs. 4(c) and 4(d)].

We have stated that the fractional effect is helpful for the stabilization of surface solitons. For single-peaked solitons in the first gap, such advantage is not so obvious. At α = 2, surface solitons with several peaks are stable at moderate power. With the increase of the number of peaks n, the stability region shrinks, which indicates that, beyond a critical value of n, solitons will be completely unstable [Fig. 4(e)]. However, at α = 1.2, surface solitons with a different number of peaks can be stable provided that their power exceeds a critical value [Fig. 4(f)]. This vividly illustrates the role of the fractional effect. While the instability of high-power solitons is effectively suppressed by the decrease of α, solitons with a low power can also be stable.

At α = 2, the stability domains of three-, five-, and seven-peaked solitons occupy about 36.3%, 30.5%, and 19.3% of their individual existence domains [Fig. 4(e)]. The corresponding ratios in a system with α = 1.2 are 76.4%, 80.3%, and 86.8%, respectively [Fig. 4(f)]. Another issue we should emphasize is that, contrary to intuition, the stability region expands with the increase of the number of peaks. This is very important for the observation of stable surface solitons with many peaks. One can expect that, when α → 1, the stability region will expand and occupy the whole existence domain of solitons. The examples of detailed linear-stability spectra are shown in Fig. 5. The very small instability growth rate shown in Fig. 5(a) indicates that the corresponding soliton can propagate without obvious distortions for a very long distance, as shown in Fig. 6(c).

 

Fig. 5 Linear-stability spectra for the unstable three-peaked and stable seven-peaked solitons marked by open circles in Fig. 4(e) and Fig. 4(f), respectively. b = 0.5 in (a) and 4.7 in (d).

Download Full Size | PDF

 

Fig. 6 Propagation simulations of multi-peaked surface gap solitons. (a, b) Two-peaked solitons in the second gaps of lattices with p = 12 marked in Figs. 3(c) and 3(d). b = 3.8 in (a) and 14.9 in (b). (c–f) Multi-peaked solitons in the first gaps of lattices with p = 4 marked in Figs. 4(e) and 4(f). b = 0.5 in (c, d) and 4.7 in (e, f). α = 2 in (a, c) and 1.2 in other panels. x ∈ [−15, 15] and z = 6000 in (a–f).

Download Full Size | PDF

To verify the linear stability analysis results, we exhaustively perform propagation simulations on the stationary solutions by the modified beam propagation method with the input condition q(x, z = 0) = w(x)[1 + ρ(x)], here ρ(x) is a random noise with a Gaussian distribution whose variance ${\sigma }_{noise}^2=0.01$. Representative examples are shown in Fig. 5. In the system with a small α, the instability region of single-peaked solitons in the first gap is very narrow, similar to the cases in [25]. However, the unstable two-peaked solitons at α = 2 in the second gap [Fig. 6(a)] become stable when α decreases to 1.2 [Fig. 6(b)]. For a more direct comparison, we show some propagation examples of solitons with different peaks in the first gaps in Figs. 6(c)–6(f). Direct propagation simulations are in good agreement with the linear stability analysis results. Stable solitons propagate stably without any distortions even in the presence of strong initial perturbations. Unstable solitons can propagate robustly over a very long distance (thousands of diffraction lengths), greatly exceeding the present experimentally feasible sample lengths.

Finally, we briefly note the necessity of the study of spatial solitons supported by NLFSE. Although the fractional-order diffraction can suppress the diffraction of light beams in linear regimes [16, 17, 21], the investigation of beam propagation in nonlinear systems with a periodic refractive-index modulation is also an important issue. We make three comments here. First, linear modes can only exist in bands and solitons reside in bandgaps; Second, linear modes in periodic systems are Bloch waves expanding to infinity and solitons are well-localized modes; Third, the amplitudes of linear modes can be arbitrary and solitons usually have a threshold power.

4. Conclusions

To summarize, we studied the properties of surface solitons residing in the finite bandgaps of the corresponding periodic fractional-order system with a defocusing nonlinearity. The existence domain of different families of solitons can be shifted by the decrease of Lévy index. The most interesting finding is that the instability of solitons can be effectively suppressed by the small Lévy index. The stability region of multi-peaked surface solitons shrinks on the one hand by the decrease of Lévy index, and, on the other hand by the increase of the number of soliton peaks. We, thus, put forward the first example of nonlinear localized surface gap waves in fractional dimensions.