Variational approach to multidimensional solitons in highly nonlocal nonlinear media

We apply the variational approach to solitons in highly nonlocal nonlinear media in $D=1,2,3$ dimensions. We compare results obtained by the variational approach with those obtained in the accessible soliton approximation, by considering the same system of equations in the same spatial region and under the same boundary conditions. To assess the accuracy of these approximations, we also compare them with the numerical solution of the equations. We discover that the variational highly nonlocal approximation provides more accurate results and as such is more appropriate solution than the accessible soliton approximation.


I. INTRODUCTION
Optical spatial solitons are self-localized wave packets that propagate in a nonlinear medium without changing their structure [1]. This is made possible by the robust balance between dispersion and nonlinearity or between diffraction and nonlinearity or between all three effects in the propagation of spatiotemporal localized optical fields. An important characteristic of real nonlinear media is their nonlocality, that is, the fact that characteristic size of the response of the medium is wider than the size of the excitation itself. Strong nonlocality is of special interest, because it is observed in many media. For example, in nematic liquid crystals (NLCs) both experimental and theoretical studies indicated that the nonlinearity is highly nonlocal -meaning that the size of response is much wider than the size of excitation [2][3][4].
In 1997, Snyder and Mitchell introduced a model of nonlinearity whose response is highly nonlocal [5] -in fact, infinitely nonlocal. They proposed an elegant theoretical model, intimately connected with the linear harmonic oscillator, that described complex soliton dynamics in simple terms, even in two and three dimensions. Because of the simplicity of the theory, they coined the term "accessible solitons" (ASs) for these optical spatial solitary waves. However, straightforward application of the AS theory, even in nonlinear media with almost infinite range of nonlocality, inevitably led to additional problems [6][7][8], because there exists no real physical medium without boundaries and without noise.
To include interactions between solitons within boundaries, as well the impact of the finite size of the sample, we developed a variational approach (VA) to solitons in nonlinear media with long-range nonlocality, such as NLCs [9], materials with thermal nonlocality [10], photorefractive crystals [11], and Bose-Einstein condensates [12]. Starting from an appropriate ansatz, this approach delivers a stationary solution for the beam amplitude and width, as well as the period of small oscillations about the stationary state. It provides for natural explanation of oscillations seen when, e.g., noise is included into the nonlocal nonlinear models. The noise is inevitable in any real physical system and causes a regular oscillation of soliton parameters with the period well predicted by our VA calculus [9]. It may even destroy solitons [13]. Although our VA results were corroborated by numerics and experiments, they still attracted a heated exchange with other researchers [14,15]. We have investigated in detail the destructive influence of noise on the shapeinvariant solitons in a highly nonlocal NLCs in [13].
Owing to great impact and practical relevance of the AS model, in this paper we systematically compare it with the VA approximation in highly nonlocal nonlinear media in D = 1, 2, 3 dimensions. We consider both systems using the same general equations, in the same spatial region, and under the same boundary conditions. To check the accuracy of both approximations, we compare them with the numerical solution of the equations. We find that multidimensional variational highly nonlocal approximation provides very accurate results, while the beam parameters obtained by using the AS approximation show systematic and predictable discrepancy with numerical results.
The paper is organized in the following manner. In Sec. 2 we introduce the model, Sec. 3 presents results obtained by the variational approach, Sec. 4 discusses the accessible soliton approximation, Sec. 5 presents numerical results, and Sec. 6 brings conclusions.

II. THE MODEL
For the study of VA and AS approximations to the fundamental soliton solutions in a (D+1)-dimensional highly nonlocal medium we adopt the following model of coupled dimensionless equations [1,9]: with zero boundary conditions on a D-dimensional sphere (D = 1, 2, 3). Here ∆ is the D-dimensional Laplacian. The system of equations of interest consists of the nonlinear Schrödinger equation for the propagation of the optical field E and the diffusion equation for the nonlocal response of the medium θ. This is a fairly general model for nonlinear optical media with diffusive nonlocality, widely used in the literature [1,6,9]. In the local limit, the first term in Eq. (2) can be neglected and the model reduces to the Schrödinger equation with Kerr nonlinearity. In the opposite limit, the second term in Eq.
(2) can be neglected and the highly nonlocal model is reached. Since we are interested in the strong nonlocality, we will omit in our analysis the α term from Eq.
(2). For radially-symmetric intensity distributions |E| 2 , Eq. (2) can be solved using Green's functions, where G D (r, ρ) is the Green's function in D dimensions, defined as follows: (4) Here d is some characteristic transverse distance of interest. In the two-dimensional case (here and thereafter), the limit D → 2 should be taken, leading to For the Gaussian-shaped beams with the field intensity: the solution of Eq. (3) can be written as where δ = max{exp(−d 2 /T 2 )} ≪ 1 and T is the characteristic width of the optical field. The D−dimensional power is given by In an explicit form for different dimensions, Eq. (7) yields: for one, two, and three dimensions, respectively. In the equations above, erf( is the Euler gamma function, and Γ(a, z) = ∞ z t a−1 exp(−t)dt is the incomplete gamma function [16]. The form of θ corresponds to a radially-symmetric solution of Eq. (2), with zero boundary conditions on a circle of radius d ≫ T , and in the limit of a thick transverse width. We can take these expressions as an approximate solution on a D−dimensional cube, as demonstrated in [9,17] for the two-dimensional case.
In the AS approximation, for the shape of the nonlocal response θ(r) of the medium one uses a parabolic function of the transverse distance. Expanding the solution of Eq. (7) in Taylor series up to the second terms, one finds where the parabolicity coefficient is and the maximum value of θ is In two dimensions this becomes: which agrees with the value determined in [9,17]. Here γ is Euler's constant.

III. VARIATIONAL APPROACH
In the variational approach, to derive equations describing evolution of the field beam expressed in an appropriate approximate form, one introduces a Lagrangian density, corresponding to Eq. (1): Thus, the problem is reformulated into a variational problem whose solution is equivalent to Eq. (1). To obtain evolution equation for an approximate field in the highly nonlocal region, an ansatz is introduced in the form of a Gaussian beam for the field [9]: in which A is the amplitude, R is the beam width, C is the wave front curvature along the transverse coordinate, and ψ is the phase shift. Variational optimization of these beam parameters will lead to the most appropriate VA solution of the problem. Likewise, a trial function for the nonlocal response of the medium is introduced, in the form Eq. (7) which is characterized by the power Q D = P D and the thickness T = R of the beam. Again, the form of θ corresponds to a radially-symmetric solution of Eq. (2), with zero boundary conditions on a D−dimensional sphere of radius d ≫ R (the limit of a thick cell δ ≪ 1). The averaged Lagrangian L D = ν L D r D−1 dr is given by: where the prime in Eq. (18) denotes the derivative with respect to z and Specifically for D = 2: In the process of optimization from the averaged Lagrangian, one obtains four ordinary differential equations (ODEs) for the beam parameters: According to Eq. (21), the beam power P D = ν |E| 2 r D−1 dr = π D/2 A 2 R D is conserved. The system of Eqs. (21-23) describes the dynamics of the beam around a stationary state.
In the stationary state (dR/dz = dC/dz = C = 0), we find the equilibrium beam width R: and the amplitude A: as functions of the beam power.
The period of small oscillations of the perturbation around the equilibrium position (C = 0) is given by the following relation: From relations (19,24,25) we also find that the propagation constant µ = dψ/dz, in the stationary state, can be written as: For the two dimensional case we have: Note that the integral quantity , which is proportional to the power, is also conserved.
Relations (25 -28) completely define the VA approximate soliton solution in the highly nonlocal case. It remains to do the same for the AS approximation.

IV. ACCESSIBLE SOLITON APPROXIMATION
In the AS approximation, the basic assumptions are that the shape of the nonlocal response of the medium is a parabolic function of the transverse distance, Eq. (11), and that the shape function of the field E is still a Gaussian given by Eq. (17). The only refractive index "seen" by the beam is that confined near its propagation axis [5].
The parameters of the trial function Eq. (17) are given now by equations: Equation (17), together with Eq. (11), exactly satisfies Eq. (1). The parameter θ D max is only a phase shift and as such quite arbitrary in the AS approximation. On the other hand, the value of Θ D is much more important; it is determined by Eq. (12). In this manner, Eq. (32) becomes: The equilibrium width R in the AS approximation is: and the amplitude A: The period of small oscillations of the width perturbation around the equilibrium can be obtained from (34): .
(37) It remains to compare these expressions with the ones obtained in the VA approach. Another useful approximation in AS approximation is based on the solution of Eq. (2) in the form of Eq. (11) when parameter Θ D = A 2 / (4D) = P D / 4Dπ D/2 R D is independent of z. In contrast to Eq. (34), in which the condition Θ D may depend on z, Eq. (32) for R has an exact oscillatory solution [18] R = R * cos 2 (πz/Λ) + P * P D sin 2 (πz/Λ), from which the solutions for C and ψ immediately follow: Here R * and P * are the width and the power of the AS solution, respectively. When P D = P * , one obtains stationary AS; otherwise, the approximate solution oscillates. The quantity Λ = Λ AS P * /P D represents the period of harmonic oscillations around the equilibrium (soliton) state, while Λ AS is the period of small oscillations of the width perturbation: Thus, one obtains nice dependencies in closed form, but unfortunately of little practical value, in view of the large discrepancy with the VA and the numerical solution of the full problem, to be presented in the next section.

V. NUMERICAL RESULTS
For numerical simulations of Eq. (1) we use finite difference time domain (FDTD) split-step method [19]. Equation (2) is solved by the tridiagonal matrix algorithm (TDMA). In our simulations we applied the same boundary conditions for both Eqs. (1,2): and ∂θ ∂r r=0 = 0, θ r=d = 0 (43) The initial beam is the trial function given by Eq. (17), with parameters determined by VA approximation, Eqs.
(25,26). Figures (1,2,3) show analytical predictions of both approximations and numerical results for one, two, and three dimensions, respectively. The beam width R of soliton solutions for AS approximation is 2 D/(8−2D) times smaller than the VA approximation at a same power, see panel (a) in Figs. (1,2,3). The corresponding stationary amplitudes in both approximations are presented in panel (b) of Figs.(1,2,3). Because P D = π D/2 A 2 R D , the equilibrium amplitude A in the AS approximation is 2 D 2 /(16−4D) times greater than the one in the VA approximation at the same power. The period in the AS approximation is 2 D/(4−D) times less than that in the VA approximation at the same power, see panel (c) in Figs. (1,2,3). Panel (d) in Figs. (1,2,3) shows an example of the soliton shape predicted by VA and AS, as well as the numerical output after a long propagation distance.
Thus, the values of the beam parameters for AS are systematically off the values for the VA approximation. In all the cases the VA predictions are confirmed by numerical results. This confirms that, while AS approximation may represent a valuable aid for an easy understanding of highly nonlocal solitons, it is a poor approximation to the more exact approaches, such as the VA approximation, and to the numerical solution.

VI. CONCLUSION
In conclusion, we have discussed the differences between the VA and AS approximate solutions to the prop-   Fig. (1), but for D = 3. The propagation distance is now z = 10 4 and P = 60. agation of solitons in highly nonlocal nonlinear media. The AS model provides a radical simplification of the problem and allows for an elegant description of solitons, but has a limited practical relevance, mainly because of the competition between the nonlocality and the finite size of the sample, which leads to systematic errors. The VA solution is not so simple, but it works very well in the limited region of large nonlocality. We have found that the AS approximation can differ up to eight times, when compared to the more realistic VA approximation and the numerical solution.

ACKNOWLEDGMENTS
This publication was made possible by NPRP Grants # 5 -674 -1 -114 and # 6-021-1-005 from the Qatar National Research Fund (a member of the Qatar Foundation). The statements made herein are solely the responsibility of the authors. Work at the Institute of Physics Belgrade was supported by the Ministry of Science of the Republic of Serbia under the projects No. OI 171033, No. 171006 and III 46016.