Nonlocality-controllable Kerr-nonlinearity in nonlocally nonlinear system with oscillatory responses

We discuss the Kerr nonlinearities of the nonlocally nonlinear system with oscillatory responses by the variational approach. The self-focusing and self-defocusing states are found to dramatically depend on the degree of nonlocality. When the degree of nonlocality goes across a critical value, the nonlinearity can transit to its opposite counterpart, that is, focusing to defocusing or defocusing to focusing. The critical degree of nonlocality for the nonlinearities transition is given analytically, and confirmed by numerical simulations. As a versatile mathematical tool, we also employ the variational approach to analytically address the stabilities of solitons, and obtain the range of the degree of nonlocality for the stable solitons, which is confirmed by the linear stability analysis and the numerical simulation.


Introduction
In nonlinear optics, most of the nonlinear effects originate from the nonlinear refraction index (NRI) [1], a phenomenon that refers to the intensity dependence on the refractive index, which is the so-called optical Kerr effect. The focusing-defocusing property is an important property relevant to the Kerr effect. If the NRI changes uniformly with the intensity, that is, the NRI increases (decreases) as the intensity increases (decreases), this will cause the beam to focus, a phenomenon known as the self-focusing. Conversely, if the NRI changes inversely with the intensity, this will cause the beam to defocus, a phenomenon known as the self-defocusing. Self-focusing and self-defocusing have been considered as the intrinsic properties of materials, that is, they only depend on the media themselves.
If the NRI at a particular point is not determined solely by the optical intensity at that point, but also depends on that in its vicinity, the nonlinearity is spatially nonlocal, which holds for the media such as nematic liquid crystals (NLC) [2][3][4][5], photorefractive media [6] and materials with thermal or diffusion nonlinearity [7,8]. Spatial nonlocality has a considerable impact on many nonlinear phenomena, such as the suppression of collapse [9], and the support of multipole solitons [10], spiralling elliptic solitons [11] and the chaoticons [12]. The ratio of the width of the response function to that of the optical beam determines the degree of nonlocality, from which four categories of the nonlocality are divided [9]: local, weakly nonlocal, generally nonlocal, and strongly nonlocal.
Recently, focusing and defocusing properties were discovered to have a dramatic relation with the generalized degree of nonlocality (GDN) in the (1 + 1)-dimensional nonlocally nonlinear system with a sine-oscillation response function [13]. By changing the GDN across the critical value, the nonlinearities can transit between the focusing and defocusing states. In this paper, we make use of the variational approach to analytically discuss the focusing-defocusing property of both the (1 + 1)-dimensional and the (1 + 2)-dimensional nonlocally nonlinear system with oscillatory responses, and analytically obtain the ranges of the GDN for the focusing-defocusing states and the stable solitons in such a system. Although the

Model
We consider the propagations of (1 + D)-dimensional optical beams in nonlocally nonlinear system described by the coupled equations for dimensionless complex optical field amplitude φ (r, z) and the NRI Δn (r, z) given by [13][14][15][16] i ∂φ ∂z where z is the longitudinal (propagation direction) coordinate, r is D-dimensional transverse coordinate vectors (D = 1 or 2), ∇ 2 D is the D-dimensional transverse Laplacian operator, and s(= ±1) is the sign of Kerr coefficient, and w m is the characteristic length of the nonlinear response. Equations (1) and (2) with s = −1 resemble the model describing the propagations of optical beams in the NLC, where the unique tiny difference is that the sign before the second term is minus in equation (2) [2]. When s = 1, equations (1) and (2) are slightly different in their two last terms (φ * rather than φ and φ 2 rather than |φ| 2 ) from the well-known second-harmonic-generation model in media with the quadratic nonlinearity, however in the stationary state the two models are the same [17]. Therefore, the model described by equations (1) and (2) can be deemed a reasonable extension of the two physical models that describe the reorientational nonlinearity in the NLC and the parametric process in the quadratic nonlinear material, and its (1 + 1)-dimensional case has been numerically discussed in our recent work [13]. In fact, we have found that equation (2) with s = −1 might describe the response of the NLC with negative dielectric anisotropy to the optical beam [16], and the detailed discussion on the its Kerr nonlinearities will be reported elsewhere.
When w m = 0, equation (2) is reduced to Δn = s|φ| 2 , then the system above becomes the standard nonlinear Schrödinger equation which exhibits the self-focusing nonlinearity for s = 1 and the self-defocusing nonlinearity for s = −1 [18,19]. When w m → ∞, equation (2) changes to w 2 m ∇ 2 D Δn = s|φ| 2 , which can describe the self-focusing nonlinearity in lead glass for s = −1 [7] and the self-defocusing nonlinearity for s = 1 as will be discussed in the following. Therefore, the nonlocal nonlinear system governed by equations (1) and (2) exhibits opposite nonlinearities in the two nonlocal limits of w m = 0 and w m → ∞ for a given s. If we define the generalized degree of nonlocality (GDN) by the ratio of the characteristic length of the nonlinear response w m to the beam width that is, σ = w m /w, then it is well expected that there must exist some certain critical GDN, across which the nonlinearity can be switched to its opposite. From equation (2), the NRI can be expressed by the following convolution where R D is the D-dimensional response function. When D = 1, the response function R 1 exhibits sine-oscillatory structure obtained in our recent paper [13] While, in the (1 + 2)-dimensional case, the response function is obtained as with Y 0 being the zeroth order Bessel function of the second kind, which also has an oscillatory profile shown in figure 1.

Variational expression of the dynamics for optical beams
Now we investigate the evolutions of a Gaussian beam at different GDNs in order to determine its sampled focusing and defocusing states, which can be analytically achieved by using the variational approach. The Lagrangian density of the system described by equation (1) plus equation (5) is [20,21] Here, we introduce a trial Gaussian beam with the input power P 0 = +∞ −∞ |φ (r, z) | 2 d D r = π D/2 A 2 w D , A and α are the amplitude and phase of the complex amplitude, c represents the phase-front curvature. The phase-front curvature c(z) indicates the expansion or the contraction of the optical beam when its sign is positive or negative [11]. All of A, α, c and w are all allowed to vary with propagation distance z. The Lagrangian can be determined by the substitution of the trial solution (9) into the Lagrangian density (8) where with m and n being integers. According to the standard procedure of variational approach [20], we obtain the evolution of the beam width w (detailed derivation is given in appendix A) where

GDN-dependent focusing and defocusing states
We assume the Gaussian beam is incident at its waist where its phase front is a plane, therefore we let c(z)| z=0 = 0 and assume w 0 = w(z)| z=0 = 1, the solution of the differential equation (11) can be obtained by the direct twice integration  where (13) implicitly gives the dependence of the beam width w(z, σ 0 ) on the propagation distance z and the GDN σ 0 , which in fact represents the detailed evolution of a Gaussian beam in the nonlocal nonlinear media under some given (initial) GDN σ 0 . It should be noted that there is not any assumption made on the specific forms of the response function in the above derivation process. Therefore, equation (13) exhibits a certain universality and then can describe evolutions of the Gaussian beams in nonlocal nonlinear media with any-profiled response function. For the response function without singularities, such as the Gaussian-formed response [22,23], equation (13) can give an approximately analytical expression of the beam width under the high GDN, after expanding the response function, for example, to the fourth order [21] and the second order [24]. However, for the response functions with complicated profiles, such as the oscillatory responses discussed in the paper, it is difficult to obtain an analytical expression from equation (13) under the general GDN.
As an example, the dependence of the beam width at the propagation distance z = 1 on the GDN σ 0 is shown in figure 2. It is known that [25], after a linear propagation of the distance of z, the beam width of a Gaussian beam is widen to be w(z) = √ 1 + z 2 , which has the value of √ 2 at z = 1. Compared with the linear propagations, the Gaussian beam will evolve in different manners if the nonlinear effects are introduced. The optical beam sampling the defocusing nonlinearity expands faster than in the linear case; and the focusing case follows an opposite trend. After nonlinear propagations in one Rayleigh distance for the Gaussian beam, its beam width w at z = 1 as the function of the GDN σ 0 is shown by the color curves in figure 2. It can be found that the defocusing exist in the smaller side of the σ 0 -coordinate when s = −1, and the defocusing will transit to the focusing when the GDN σ 0 increases above the critical value. The inverse happens for case that s = 1: focusing will transit to the defocusing when GDN σ 0 increases above the same critical value. The critical GDNs possess different values for the (1 + 1)-and (1 + 2)-dimensional cases, that is, σ c = 0.82 for the former and σ c = 0.63 for the latter, respectively. Besides, the optical power P 0 strengthens the nonlinear effects for both the focusing and the defocusing nonlinearities, which is shown in figure 2. For comparison, we numerically simulate the nonlocally nonlinear system described by equations (1) and (5) by using the split-step Fourier method [18] (the main procedure of numerical simulations can be found in appendix B). A well agreement between the numerical results and the variational ones (13) can be found in figure 2.

Critical GDNs vs propagation distances
We discussed the critical GDNs for the nonlinearities transition by considering output beam width at one Rayleigh distance above. However, if we investigate the evolutions of the beam at different propagation distances (the beam width is given implicitly by equation (13)), the corresponding critical GDNs will be found to be different. In other words, the critical GDNs should be the functions of propagation distance z.
The beam in focusing states narrows its width for strong nonlinearity and weakens expanding tendency for weak nonlinearity. While in defocusing state, the expanding rate of beam width is always quicker than that in the linear case. Then, it is well expected that the critical case occurs when the width at z of the optical beam in nonlinear state is just equal to √ 1 + z 2 , which is the value in linear propagating case. Therefore, by replacing w with √ 1 + z 2 in equation (13), we obtain the critical GDNs at any propagation distance Equation (14) gives an implicit function of the critical GDN σ c on the propagation distance z, the optical power P 0 , the sign of Kerr coefficient s, and the dimension D of the nonlinear system. For the given s, D and P 0 , the function σ c (z) can be obtained by numerically solving the integral equation equation (14) for different propagation distance z. The critical GDN as a function of the propagation distance up to z = 2 is shown in figure 3, Where the critical GDN σ c is found to increase as the propagation distance z. Besides, at a certain z, the critical GDN σ c also depends on the optical power P 0 . An increasing P 0 can increase the critical GDN σ c for s = 1, while decrease σ c for s = −1.
We can understand the z-dependent critical GDNs by considering such a beam evolution at the GDN σ 0 (= 0.85) for the case of s = −1, which is a little larger than σ c (= 0.82) obtained at z = 1. Although the focusing effect is acting in this situation, the beam still expands due to the dominant diffraction effect, where the nonlinearity only weakens the expanding rate. As the beam propagating, the expansion of its beam width drops the GDN so that σ| z=z 0 < σ c | z=1 at some certain propagation distance z 0 (> 1), and makes the beam fall into the defocusing state. If we focus on the beam width at z = z 0 for different σ 0 , the GDN σ 0 (= 0.85) should be in the defocusing range, which also implies that the critical GDN σ c | z=z 0 should be larger than 0.85. This aforementioned kind of initial focusing evolving into the subsequent defocusing for s = −1 when the initial GDN σ 0 is little bit larger than the critical value is demonstrated in figure 4. The same analysis can be applied into the case of s = 1 as follows. Taking an initial GDN σ 0 a little larger than the critical value obtained at z = 1, we can expect the initial nonlinearity is defocusing at least during one Rayleigh distance, which makes optical beam expand quickly and then drops the GDN. When the GDN decreases below the critical value at some certain propagation distance z 0 , that is, σ| z=z 0 < σ c | z=1 , then the beam becomes to sample the focusing nonlinearity. Therefore, by discussing the evolutions of beam width at z 0 , we will infer that focusing nonlinearity instead of defocusing nonlinearity exist at z 0 for the GDN σ 0 , and the critical GDN σ c | z=z 0 should be larger than the value of σ c | z=1 . The initial defocusing evolving into the subsequent focusing for s = 1 when the initial GDN σ 0 is little bit larger than σ c | z=1 is demonstrated in figure 5. Therefore, the critical GDNs should be z-dependent, and rise as z increases for both s = 1 and s = −1.

Physical explanation of focusing-defocusing property
As mentioned above, when the GDN goes across σ c , the nonlinearity exhibited by the nonlocally nonlinear system described by equation (1) plus equation (5) can transit into its opposite counterpart. In principle, whether the nonlinearity is focusing or defocusing can be determined by the relative distributions of optical beam and its induced NRI. For the focusing nonlinearity, the NRI changes uniformly with the intensity, that is, the NRI increases as the intensity increases, and vice versa; the defocusing nonlinearity is on the contrary. We present the variation of the NRI in large and small ranges of the GDNs in figure 6 by substituting a Gaussian beam into the convolution (5) in the (1 + 1)-dimensional case.
According to the commutative law of the convolution, the NRI can be expressed by Δn(x) = s ∞ −∞ R 1 (x )|φ(x − x )| 2 dx , which mathematically is the area enclosed by the product of the response function R 1 (x ) and the optical intensity |φ(x − x )| 2 with the x -axis. For the large GDN, the beam is much narrower than w m , and it is a good approximation to take the response function outside the convolution, which yields The profile of the NRI is the same as the response function shown in figures 6(a)-(c). Therefore, the NRI in the region occupied by the Gaussian beam exhibits the bell and inverted bell shapes for cases of s = −1 and s = 1, respectively. Correspondingly, the Gaussian beam will be in the focusing and defocusing states [13]. The case of small GDNs is quite different, for which the beam is much wider than w m . As shown in figures 6(d) and (e), when the Gaussian beam centered at the origin of coordinate moves to the neighboring point x = x, it will 'see' more negative parts and less positive parts of the response function, which results in the decrease of the NRI in the case of s = 1 when the beam center is away from the origin of coordinate. According to the even symmetry, the NRI induced by the Gaussian beam will exhibit the bell-shaped profile shown in figure 6(f) when s = 1, and the Gaussian beam will be in the focusing state. Conversely, when s = −1 the Gaussian beam induces the inverted bell-shaped NRI, and will be in the defocusing state. Therefore, as the GDN goes down across σ c , the defocusing nonlinearity will transit to the focusing nonlinearity when s = 1, while the opposite happens when s = −1.
The (1 + 2)-dimensional case can also be discussed in the same way and is found to exhibit the similar characteristic behavior with the (1 + 1)-dimensional case, which is not given here.

Variational results of the existence range for bright solitons
Although the existence range of GDNs for the stable solitons has been numerically obtained for the (1 + 1)-dimensional nonlocally nonlinear system governed by equations (1) and (5) in cases of s = −1 [13] and s = 1 [14], enormous amounts of numerical iterations and simulations were needed there, which is quite time consuming. We will take advantage of the variational approach to approximately obtain the GDN-range of stable solitons supported by the (1 + D)-dimensional nonlocally nonlinear system, which can supply a guidance to numerical calculations and then save amounts of time for the numerical explorations. Furthermore, as will be shown in the following, the analytical (variational) results and numerical results can support and complement each other.
Integrating equation (11) with w 0 = 1 being assumed and N D given by equation (12). Soliton solutions correspond to the extremum points of the equivalent potential V D (w) [20,23]. Then by letting that ∂ w V D = 0 the critical power can be obtained as It should be noted that in equation (16) the critical power P c must be positive, which requires that N D > 0 for s = 1 and N D < 0 for s = −1.
As mentioned above, the first order derivative of the equivalent potential V D (w) can give its extreme point where solitons exist at P 0 = P c . While, the second order derivative of V D (w) can determine the stability of solitons. If the extreme point of V D (w) happens to be a minimum, i.e.  the soliton will be stable. Depending on D, Φ D in equation (17) has two different forms Specially for the case of D = 2, N 2 and Φ 2 can exhibit the following compact expressions (detailed derivation is given in appendix appendix C) with Ei being the exponential integral function [26]. From equation (17), the variational results of the GDN-range for stable solitons are obtained as shown in table 1. The numerical results are also given for comparison, the details of which will be given in the following. To confirm the variational predication, we can address the linear stability analysis of the nonlinear system (equation (1) plus equation (5)) [27] when D = 1. As for the case of D = 2, it is more difficult to perform the linear stability analysis because of the limitation of the computer resources. However, we can directly make the numerical simulations of iterated soliton profiles obtained by the perturbation-iteration method [28,29] to test their stabilities.
In the (1 + 1)-dimensional case of the nonlinear system, we search for perturbed solutions in the form φ (x, z) = [u (x) + μ (x, z) + iν(x, z)] exp (iλz), where perturbations μ and ν can grow with a complex rate δ on propagation. Linearization of equation (1) around u(x) yields: where ΔN = 2s ∞ −∞ R(x − x )u x μ x dx and Δn is given by equation (5). We numerically solve the system of equations (19) and (20), the result of which is summarized in table 2. The linear stability analysis indicates that solitons are stable when 1.05 σ 0 < +∞ in the case of s = −1 and when 0.05 σ 0 0.78 in the case of s = 1, the former case of which agrees well with the variational predication. As for the latter case of s = 1, although our numerical analysis in the GDN range of 0 σ 0 0.05 cannot be given due to the limitation of the computer resources (mainly because in this weakly nonlocal case it needs massive data points caused by the rapid oscillation of the response function 'seen' by the optical beam), our variational result can still work and give the correct prediction that the solitons are stable when 0 σ 0 0.05. Especially for the local case of σ 0 = 0, the standard nonlinear Schrödinger equation (3) has the stable sech-form bright soliton [25]. Then, it is well expected that as the GDN σ 0 decreases, the soliton profile will gradually approach to the sech-profile, which is shown in figure 7. Therefore, we can summarize the σ 0 spectrum of stable solitons for the (1 + 1)-dimensional nonlocally nonlinear system as 0 σ 0 0.78 when s = 1 and 1.05 < σ 0 < ∞ when s = −1.   In the (1 + 2)-dimensional case of the nonlinear system, the numerical simulations show that the bright solitons can propagate stably when 0.33 σ 0 < 0.56 in the case of s = 1 shown in figure 8. In fact, it is expected that solitons should be unstable under low GDN. Especially, it is well known that the (1 + 2)-dimensional solitons [30] (the self-trapping beams [31]) will collapse [32] in the local (σ 0 = 0) bulk Kerr media. In the case of s = −1, bright solitons propagate stably when σ 0 > 1.11 shown in figure 9. Obviously, all the numerical simulations again confirm our variational predication of the GDN range of the stable bright solitons.
Of course, the variational approach can be further applied into the prediction of the GDN existence range for solitons with complicated profiles like the Hermite-Gaussian-like ones, which is not the focus in the paper and will be reported elsewhere.

Conclusion
By using the variational approach, we discuss the focusing-defocusing property of a kind of nonlocally nonlinear system with oscillatory responses. In such a nonlocal nonlinear system, the focusing and defocusing states are found to depend closely on the generalized degree of nonlocality (GDN) of the system. The nonlinearities can transit between focusing and defocusing states when the GDN of the system goes cross a critical value. The critical GDN is obtained by the variational approach, which is the function of the propagation distances. The optical beam, when s = −1, samples self-defocusing and self-focusing when the GDN below and above the critical value, and the opposite happens when s = 1. The GDN-range of bright solitons is also obtained by the variational approach, which is consistent with the linear stability analysis and the numerical simulation. The quantitatively relation between the (focusing or defocusing) nonlinearities and the GDN, obtained by the variational approach, supply a theoretical guidance to the tailoring and engineering of nonlinearities exhibited by the nonlocally nonlinear system.
where the exponential operatorD can be evaluated in the Fourier domain In order to improve the accuracy, the split-step Fourier method can be further optimized to the following procedure of the symmetric form