Spatiotemporal Airy Ince–Gaussian wave packets in strongly nonlocal nonlinear media

The self-accelerating Airy Ince–Gaussian (AiIG) and Airy helical Ince–Gaussian (AihIG) wave packets in strongly nonlocal nonlinear media (SNNM) are obtained by solving the strongly nonlocal nonlinear Schrödinger equation. For the first time, the propagation properties of three dimensional localized AiIG and AihIG breathers and solitons in the SNNM are demonstrated, these spatiotemporal wave packets maintain the self-accelerating and approximately non-dispersion properties in temporal dimension, periodically oscillating (breather state) or steady (soliton state) in spatial dimension. In particular, their numerical experiments of spatial intensity distribution, numerical simulations of spatiotemporal distribution, as well as the transverse energy flow and the angular momentum in SNNM are presented. Typical examples of the obtained solutions are based on the ratio between the input power and the critical power, the ellipticity and the strong nonlocality parameter. The comparisons of analytical solutions with numerical simulations and numerical experiments of the AiIG and AihIG optical solitons show that the numerical results agree well with the analytical solutions in the case of strong nonlocality.

Ince-Gaussian (IG) beams, which constitute the third complete family of transverse eigenmodes of stable resonators, have attracted extensive attention from research communities all over the world since it was introduced by Bandres and Gutiérrez-Vega 1,2 . The transverse structure of IG modes is described by the Ince polynomials, and the limiting cases of IG modes are Laguerre Gaussian (LG) modes and Hermite Gaussian (HG) modes when the ellipticity parameter of IG modes tends to zero or infinity, respectively. Meanwhile, Schwarz et al. have generated single high order IG modes with very high quality by slightly breaking the symmetry of the cavity of a diode pumped Nd:YVO 4 laser and its pump beam configuration 3,4 .
Soon afterwards, Deng et al. have discovered the IG solitons in strongly nonlocal nonlinear media (SNNM) [5][6][7] . Recent experimental results demonstrate that nematic liquid crystals 8 and lead glass 9 are strongly nonlocal nonlinear media. Propagation properties of Airy beams in the SNNM has been studied by Zhou et al. 10 . Then, the anomalous interaction of Airy beams in nonlocal nonlinear media have been reported by Shen et al. 11 , in which nonlocal nonlinearity also affects the interaction of out-of-phase bright solitons and dark solitons. In addition, owing to the SNNM, Zhu et al. have obtained that the Airy vortex beams follow a periodic trajectory when propagating through a SNNM 12 .
Over the past decades, tremendous efforts have been made by researchers to generate the three-dimensional (3D) localized optical solitons [13][14][15][16][17][18][19][20][21] , which also called spatiotemporal light bullets, are localized in both space and time. Among which, Airy wave packets were widely reported in free space [16][17][18] and quadratic index medium 20,21 due to the self-accelerating [22][23][24][25][26] , self healing 27 and no diffraction 28 features of Airy distribution. It is worth mentioning some pioneer works in Airy related wave packets. The spatiotemporal Airy-Bessel light bullets by combining an Airy pulse with a two-dimensional Bessel beam have been studied by Chong et al. 14 18 . Then, the chirped Airy Gaussian vortex wave packets in quadratic index medium 20 and Airy Gaussian and Airy Gaussian vortex light bullets in harmonic potential 21 have been reported. However, what will happen when the 3D localized Airy Ince-Gaussian (AiIG) wave packets propagate in the SNNM? It will be significant to investigate the AiIG and Airy helical Ince-Gaussian (AihIG) optical breathers and solitons and their propagation dynamics in the SNNM.

Methods
In the paraxial approximation and in the nonlocal nonlinear media, the propagation of the wave packet obeys the (1 + 3)D nonlocal nonlinear Schrödinger equation [5][6][7]29,30  where U(r, t, z) = U(x, y, t, z) is a paraxial wave packet, r = x y 2 2 + , x and y are the transverse coordinates, z is the longitudinal coordinate (the propagation distance in nonlocal nonlinear materials), k is the wave number in the media without nonlinearity, k 0 ′′ is the dispersion coefficient of the media evaluated at ω 0 . Δn = n 2 ∫N(r − r′)|U(r′, z)| 2 d 2 r′ is the spatial nonlinear perturbation of the refraction index, n 2 is the nonlinear index coefficient, r and r′ are two-dimensional transverse coordinates, N is the normalized symmetrical real spatial response function of the media. Without loss of generality, we assume the material response to be the Gaussian function where w m is the characteristic length of the material response function. Underlying the dimensionless coordinates (X, Y, T, Z) = (x/w 0 , y/w 0 , t/t 0 , z/L), w 0 is the spatial scaling parameter, t 0 is the temporal scaling parameter, = L kw 0 2 is the diffraction length. In the case of the strong nonlocality 29 , Eq. (1) can be simplified into the normalized dimensionless Snyder-Mitchell model where P 0 is the input power at Z = 0, P c = n 0 /(γn 2 L 2 ) is the critical power for the soliton propagation, γ is the material parameter relating to N. Eq. (2) is an equation in the case of the nonlinearity limit with the degrees of nonlocality approaching to infinity, the field can change the refractive index of the medium while propagating, which is creating a structure similar to a graded-index fiber. Schrödinger equation in ref. 20 . is from the quadratic index medium, while Eq. (2) is from strongly nonlocal nonlinear medium condition, they are two different questions. The solution of Eq. (2) can be obtained from the method of separation of variables as in the transverse plane perpendicular to Z, the elliptic coordinates [1][2][3][4][5] denote the radial and angular elliptic variables, respectively. Semifocal separation f(Z) = f 0 w(Z)/w 0 , f 0 denotes semifocal separation at waist plane Z = 0. ξ and J satisfy continuity in the whole space. φ G is the Gaussian solution 5,18 (2), we have the following equations where a and p are separation constants and ε = f w 2 / 0 2 0 2 is the ellipticity parameter of the IG wave packet. In Eq. (7), we deal with finite energy Airy wave packet for the physical reality, which can be expressed as A(T, is the decay factor. By using the Fourier-transform method, one can obtain the solution of Eq. (7) as The solution of Eq. (8) can be given as , where m p 0 ⩽ ⩽ for an even function, m p 1 ⩽ ⩽ for an odd function, and indices (p, m) have the same parity. Then, the spatial even IG, odd IG and helical IG (hIG) distribution can be expressed as ]. The even AiIG, odd AiIG and AihIG wave packets in SNNM can be expressed as , ,

Results and Discussion
Having found the analytical solutions of Eqs (16)(17)(18) in SNNM, we concentrate on the comparison of analytical solutions with numerical simulations and numerical experiments of the spatiotemporal AiIG and AihIG wave packets in the following. We discuss various examples of the localized wave packets for different parameters. The evolution properties of the AiIG and AihIG wave packets can be easily manipulated and modulated through adjusting the ratio between the input power and the critical power, the ellipticity and the strong nonlocality parameter. By comparing Fig. 1(a1) with Fig. 1(a2), there is very good qualitative agreement between analytical result and numerical simulation of the temporal Airy wave packets. Propagation of the spatial odd IG wave packets is shown in Fig. 2(a1-d1), for different choices of the parameter P 0 /P c . Without doubt, the wave is free to broaden when P 0 approaches to zero. IG breathers initially broaden because wave packet diffraction initially overcomes wave packet induced refraction as P 0 = 0.5P c in Fig. 2(b1). The comparison shows that the IG solitons propagate stably with the input power equaling the critical power in Fig. 2(c1), which means that diffraction is exactly balanced by nonlinearity. And the IG breathers initially narrow as P 0 = 1.5P c in Fig. 2(d1). The numerical simulation of the spatial propagation distributions is shown in Fig. 2(a2-d2), which agrees well with the related analytical result. The spatiotemporal AiIG solitons can be realised theoretically when the input power equals the critical power, with steady state in spatial and almost non-dispersion in temporal dimensions. On the other hand, we can achieve 3D localized AiIG breathers periodically oscillating when the balance between diffraction and nonlinearity is broken.
As revealed, the initial spatiotemporal AiIG and AihIG soliton patterns vary with the ellipticity ε are shown in Fig. 3. Under the condition of the input power equaling the critical power, we obtain an approach to generate approximately steady 3D AihIG solitons, which can maintain the shape after propagating many Rayleigh lengths by compare Fig. 4(a1-c1) with Fig. 3(a3-c3). Actually, the peaks along the T direction go ahead after propagation due to the self-accelerating character of Airy wave packets, which can also be found by comparing different propagation distances in Fig. 1(a1) and (a2). In addition, we utilize split step Fourier transform to show the numerical simulation, α = w 0 /w m denotes the degree of the material nonlocality. The less α is, the stronger the nonlocality is. The numerical results (Fig. 4(a2-c2)) agree well with the analytical solutions ( Fig. 4(a1-c1)) in the case of strong nonlocality with α = 0.01. However, at the same propagation distance, the numerical simulation results present unstable properties when the degree of nonlocality becomes weaker (α = 0.5), as shown in Fig. 4(a3-c3).
Likewise, one can easily obtain the self-accelerating AiIG breathers periodically oscillating when the balance between diffraction and nonlinearity is broken. The phenomenon shows the spatial periodically oscillating (breather state) or steady (soliton state), temporal self-accelerating and approximately non-dispersion properties of AiIG breathers and solitons reported in SNNM are different from other spatiotemporal Airy related wave packets diffracting in spatial dimension while propagating in free space 17,18 . Without loss of generality, the relationship of AiIG, Airy Laguerre Gaussian (AiLG), and Airy Hermite Gaussian (AiHG) breathers and solitons will be analysed as follow. When ε → 0, the transition from IG pm to LG nl occurs. Simultaneously, the indices of both modes are related as: l = m and n = (p − m)/2 in the limit. When ε → ∞, the transition from IG pm to HG l l x y occurs. In the limit, the indices are related as: for even AiIG breathers and solitons l x = m and l y = p − m, on the contrary, for odd AiIG breathers and solitons l x = m − 1 and l y = p − m + 1.
We formulate the spatial hIG wave packet analytically in SNNM and present related numerical experiment results of the hIG wave packet generation, corroborating the properties we describe. The intensity patterns of the initial input of hIG beams are shown in Fig. 5(a1-a3). The phase patterns of the initial input of hIG beams are shown in Fig. 5(b1-b3), where the locations of the optical vortexes are obtained. For numerical experimental generation, we launch an initial beam [see Fig. 5(a1-a3)] to reconstruct the computer-generated holograms [see gray insets display of Fig. 5(c1-c3)] of the desired beam profiles 18     To obtain the stability analysis of the hIG wave packets, in Fig. 6, we numerically simulate the intensity distribution, which is excited by an initial perturbation 31 . The initial condition is supposed to be IG h (X, Y, Z = 0) + ρψ(X, Y, Z = 0), where ψ(X, Y, 0) is a random complex function whose maximum amplitude is less than that of the hIG beam, and ρ denotes the perturbation parameter. The difference between the output intensity distribution of the hIG beams and the solitons is huge when ρ = 0.05 in Fig. 6(c1-c3); the difference becomes so small that one can hardly distinguish it when ρ = 0.01 in Fig. 6(b1-b3). Stable soliton propagation in a numerical experiment where the initial beam is perturbed by noise, which may not constitute a rigorous proof of the stability, but does provide strong support for the existence of observable nonlinear modes in laboratory experiments. The related numerical simulations of spatiotemporal AihIG wave packets with ρ = 0.01 and ρ = 0.05 are shown in Fig. 7, while the comparison with ρ = 0 can be found in Fig. 4(a2-c2).
The local energy flow is usually expressed in terms of the Poynting vector. The Poynting vector has a magnitude of energy per unit area (or per unit time), and a direction which represents the energy flow at any point in the field. The vector potential is a linear polarization state. Given an X polarized vector potential V = IG h (X, Y, Z)e ikZ e X , where e X is the unit vector along the X direction. The time-averaged Poynting vector can be expressed as

32
, where c 0 is the velocity of light in a vacuum. Figures 8(a1-c1) show that the transverse energy flow appears to rotate counterclockwise around the vortex (distinct circulation of current). The length and direction of the arrows represent the intensity and direction of the Poynting vector respectively. As in mechanics, the time-averaged angular momentum density for the electromagnetic field is the angular momentum per unit area (per unit time), obtained by forming the cross product of the position vector with the time-averaged momentum density 〈j〉 = r × 〈E × B〉 33 . The spatial angular momentum distributions in Fig. 8(a2-c2) are more coherent than the intensity distributions in Fig. 5(a1-a3), and the positions of optical vortex may not be obvious.

Conclusion
In summary, we have obtained a novel class of self-accelerating AiIG and AihIG optical breathers and solitons in SNNM from the method of separation of variables. The evolution properties of the AiIG and AihIG wave packets can be easily manipulated and modulated through adjusting P 0 /P c , the ellipticity ε and the strong nonlocality parameter α. The AiIG and AihIG optical solitons can be obtained when the input power equals the critical power, while the AiIG and AihIG breathers can be achieved with P 0 < P c or P 0 > P c . The comparisons of analytical solutions with numerical simulations and numerical experiments of the AiIG and AihIG optical solitons show that the numerical results agree well with the analytical solutions in the case of strong nonlocality. For the first time, the spatial periodically oscillating (breather state) or steady (soliton state), temporal self-accelerating and approximately non-dispersion properties of AiIG and AihIG breathers and solitons are reported in SNNM, which is quite different from the case in free space 18 . If the last term of Eq. (1) in left part disappears, it becomes a free space condition, and the propagation of the spatial part is similar shown in Fig. 2(a1), which will expand due to the diffraction, and never show breathers and solitons.
In general, the analytical solution described here is applicable to other optical breathers and solitons such as self-accelerating AiHG optical breathers and solitons, self-accelerating AiLG optical breathers and solitons, and related self-accelerating Airy elegant breathers in SNNM. We foresee potential applications in signal processing due to the stabilized propagation properties of these wave packets.