## Abstract

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 by 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 accessible soliton approximation suffers from systematic errors, when compared to the variational approach and the numerical solution. The errors increase with the dimension of the system. The variational highly nonlocal approximation provides more accurate results in any dimension and as such is more appropriate solution than the accessible soliton approximation.

© 2014 Optical Society of America

## 1. Introduction

Optical spatial solitons are self-localized wave packets that propagate in a nonlinear medium without changing their profile [1]. This is made possible by the balance between dispersion and nonlinearity or between diffraction and nonlinearity or between all three effects in the propagation of spatiotemporal localized optical fields. At most, the change in the profile may proceed in an oscillatory fashion, in which case one obtains soliton breathers. An important characteristic of many nonlinear media is their nonlocality, that is, the property that the size of the response of the medium is different from the size of the excitation [2]. High nonlocality is of special interest, because it is observed in many media. It means that the size of the response is much wider than the size of the excitation. For example, in nematic liquid crystals (NLCs) both experimental and theoretical investigations indicate that the nonlinearity is highly nonlocal [3–6].

Snyder and Mitchell introduced in 1997 a model of nonlinearity whose response is infinitely nonlocal [7]. They proposed an elegant theoretical model, intimately connected with the linear harmonic oscillator, that described complex soliton dynamics in simple terms, even though the model itself is actually linear. 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 very wide range of nonlocality, may lead to additional problems [8–10], because there exists no real physical medium without boundaries and the corresponding boundary conditions. Owing to the popularity of AS model and its widespread use, it is worthwhile to check how accurate it really is.

To include properly the dynamics of solitons within bounded media, as well the impact of the finite size of samples, we developed a variational approach to solitons in nonlinear media with long-range nonlocality. These include NLCs [11], materials with thermal nonlocality [12–16], photorefractive crystals [17], and Bose-Einstein condensates [18], among other. Starting from an appropriate ansatz, the variational approach (VA) delivers a stationary solution for the beam amplitude and width, as well as the period of small oscillations about the stationary state. It provides a natural explanation for oscillations seen when, e.g., noise is included into the nonlocal nonlinear models. Variational approach was already used quite successfully for the prediction of solitons in atomic BECs [19, 20].

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 calculation [11]. It may even destroy solitons after a prolonged propagation [21]. Although our VA results are corroborated by numerics and experiments, they still invited a heated exchange with other researchers [22, 23]. We have also investigated in detail the destructive influence of noise on the shape-invariant solitons in highly nonlocal NLCs in [21].

Having in mind 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 approaches 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 the multidimensional variational highly nonlocal approximation provides accurate results, very close to the numerical, while the beam parameters obtained by the AS approximation show systematic and predictable discrepancy with the numerical results. Except for our own preliminary analysis of *D* = 2 case in [24, 25], we could not locate any other publication that treats multidimensional VA and AS approximations on the same footing and under identical conditions, as done here. The closest in that respect was [26], where both VA and AS approximations were applied to a nonlocal nonlinear coupler, but in the opposite nonlocality regimes and for *D* = 1.

The paper is organized in the following manner. In Sec. 2 we introduce the model, Sec. 3 presents results obtained by the variational approach, and Sec. 4 discusses the accessible soliton approximation. In Sec. 5 we present numerical results and the comparison between the two approximate solutions. Finally, Sec. 6 brings conclusions.

## 2. The model

In this study we consider a (*D*+1)-dimensional highly nonlocal medium and adopt the following model of coupled dimensionless equations for the generation of scalar optical solitons [1,6,11]:

*D*-dimensional sphere (

*D*= 1, 2, 3). Here Δ is the

*D*-dimensional ”transverse” Laplacian and

*z*is the propagation coordinate. 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

*θ*. In the (3+1)-dimensional case one has the choice of considering the propagation of spatiotemporal solitons, in which case one of the 3 variables is the reduced time, or the temporal dynamics of three-dimensional wave packets according to the usual rules of quantum mechanics, in which case the marching variable

*z*stands for the regular time. Then, the model is relevant for the description of Bose-Einstein condensates [19,20], for example. In general, this is a fairly general model for nonlinear optical media with diffusive nonlocality, widely used in the literature [1,8,11,16].

In the local limit, the first term in Eq. (2) can be neglected and the model then 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 (high) nonlocality, we will omit in our analysis the *α* term from Eq. (2). In addition, we will consider VA and AS approximations to the fundamental soliton solutions only.

For radially-symmetric intensity distributions |*E*|^{2}, Eq. (2) can be solved using Green’s functions,

*G*(

_{D}*r*,

*ρ*) is the Green’s function of Eq. (2) in

*D*dimensions. An appropriate function can be written as follows [27] :

*d*is some characteristic (wide) transverse distance of interest. In the two-dimensional case (here and thereafter), the limit

*D*→ 2 should be taken, leading to

*δ*= max{exp(−

*d*

^{2}/

*T*

^{2})} ≪ 1 and max is taken over

*T*— some characteristic width of the optical field. The

*D*–dimensional power is given by

*Q*=

_{D}*ν*∫ |

*E*|

^{2}

*r*

^{D−1}

*dr*, where

*ν*= 2

*π*

^{D/2}/Γ(

*D*/2) is the size of the

*D*-dimensional unit sphere.

In an explicit form for different dimensions, Eq. (7) yields:

*r*= |

*x*|, in two dimensions $r=\sqrt{{x}^{2}+{y}^{2}}$, and in three $r=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}$ or $r=\sqrt{{x}^{2}+{y}^{2}+{t}^{2}}$, depending on the interpretation of the marching variable. In the equations above, $\text{erf}(z)=\left(2/\sqrt{\pi}\right){\int}_{0}^{z}\text{exp}(-{t}^{2})dt$ is the error function, $Ei=-{\int}_{-z}^{\infty}{t}^{-1}\text{exp}(-t)dt$ is the exponential integral function, and $\mathrm{\Gamma}(z)={\int}_{0}^{\infty}{t}^{z-1}\text{exp}(-t)dt$ is the Euler gamma function [28]. The form of

*θ*corresponds to a radially-symmetric solution of Eq. (2), with zero boundary conditions on a sphere 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 [11, 25] for the two-dimensional case.

Having found an approximate solution for the nonlocal response *θ* of the medium, it remains to do the same for the envelope *E* of the optical beam. This is done differently in the VA and AS approximations.

## 3. Variational approach

In the variational approach, to derive equations describing evolution of the electric field of the beam, expressed in an appropriate approximate form, one introduces a Lagrangian density, corresponding to Eqs. (1) and (2):

*A*is the amplitude,

*R*is the beam width,

*C*is the wavefront curvature along the transverse radial 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 of Eq. (7). 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}* =

*ν*∫

*ℒ*

_{D}r^{D−1}

*dr*is given by:

*z*and

*r*goes formally from 0 to ∞, but effectively from 0 to

*d*, owing to the boundary conditions imposed. Specifically for

*D*= 2:

*γ*is Euler’s constant.

In the process of optimization from the averaged Lagrangian, one obtains two algebraic equations whose solution is *Q _{D}* =

*P*and

_{D}*T*=

*R*, and four ordinary differential equations for the beam parameters:

*P*=

_{D}*ν*∫ |

*E*|

^{2}

*r*

^{D−1}

*dr*=

*π*

^{D/2}

*A*

^{2}

*R*is conserved. The system of Eqs. (19)–(21) describes the dynamics of the beam around a stationary state.

^{D}In the stationary state (*dR/dz* = *dC/dz* = *C* = 0), we find the equilibrium beam width *R*

*A*

*C*= 0) is given by the following relation:

*μ*=

*dψ/dz*in the stationary state:

Relations (23) – (26) completely define the VA approximate soliton solution in the highly nonlocal case. Now, one has to do the same for the AS approximation.

## 4. 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 and that the shape function of the field *E* is still a Gaussian given by Eq. (13). The only refractive index ”seen” by the beam is that confined near its propagation axis [7]. Expanding the solution of Eq. (7) in Taylor series up to the second-order terms, one finds

*θ*is

*can be obtained if Eqs. (28) and (6) are directly substituted into Eq. (2).*

_{D}The parameters of the trial function from Eq. (13) are given now by equations:

Equation (13), together with Eq. (28), exactly satisfies Eq. (1).The parameter *θ*_{Dmax} 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. (29). In this manner, Eq. (34) becomes

*R*in the AS approximation is and the amplitude

*A*:

That approach to AS approximation is based on the solution of Eq. (2) in the form of Eq. (28) when parameter Θ* _{D}* =

*A*

^{2}/(4

*D*) =

*P*/(4

*Dπ*

^{D/2}

*R*) is independent of

^{D}*z*. In that case, Eq. (34) for

*R*has an exact oscillatory solution

*C*immediately follows:

*R*and

_{*}*P*

_{*}are the width and the power of the AS solution, respectively. When

*P*=

*P*

_{*}, one obtains stationary AS; otherwise, the approximate solution oscillates, forming a breather. The quantity Λ = Λ

_{*}(

*P/P*

_{*})

^{2/D}represents the period of harmonic oscillations, while Λ

_{*}is the period of small oscillations of the perturbation to the beam width, Thus, in both approaches to AS approximation one obtains nice dependencies in closed form, but unfortunately of little real value, in view of the large discrepancy with the VA and the numerical solution of the full problem, to be addressed in the next section.

## 5. Numerical results

For numerical simulations of Eq. (1) we use finite difference time domain (FDTD) split-step method [30]. Equation (2) is solved by the tridiagonal matrix algorithm, which is also known as the Thomas algorithm [31]. In our simulations we applied the same mixed boundary conditions to both Eqs. (1) and (2):

*D*settings, except that in [6] they were written for an infinite sample. We need boundary conditions over a finite sample because it is there that the inadequacy of AS becomes apparent. The initial beam is the trial function given by Eq. (13), with parameters determined by the VA approximation, Eqs. (23) and (24). The numerical solutions slightly oscillate initially, but after a prolonged propagation settle into a steady state.

Figures (1,2,3) show analytical predictions of both approximations and numerical results for one, two, and three dimensions, respectively. One can see that the AS approximation is systematically off the numerical solution, becoming progressively worse as the dimension is increased. The beam width *R* of soliton solutions for AS approximation is 2^{D/(8−2D)} times smaller than the VA approximation at the 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*, the equilibrium amplitude

^{D}*A*in the AS approximation is 2

^{D2/(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) depicts an example of the soliton profile predicted by VA and AS, as well as the numerical output after a long propagation distance.

On the other hand, the values of quantities obtained by the VA approximation consistently follow the numerical values closely, regardless of the dimension. Essentially, in all the cases the VA predictions are confirmed by numerical results, which cannot be said for the AS approximation. The values of the beam parameters for AS approximation are systematically off the values for the VA approximation. 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 of course, the numerical solution.

Finally, we compare the stability of numerical to the VA solutions. The stability of AS is of no concern, because as a solution of the quantum-mechanical linear harmonic oscillator it is always stable. Even the stability of VA is beyond question, because it is clear that when perturbed, it can only oscillate about the equilibrium state. Still, it is advisable to check the stability of the numerical solution and to compare it with the VA approximation. This is demonstrated in Fig. 4. According to the Vakhitov-Kolokolov stability criterion, the necessary condition for the stability of solitons is that *dP/dμ* > 0, which is obviously satisfied.

## 6. Conclusion

In conclusion, we have discussed the differences between the VA and AS approximate solutions to the propagation 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. We also demonstrated that VA and numerical solutions are stable in the highly nonlocal regime.

## 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). 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.

## References and links

**1. **Y. S. Kivshar and G. Agrawal, *Optical solitons: from fibers to photonic crystals* (Academic press, 2003).

**2. **O. Bang, W. Krolikowski, J. Wyller, and J. J. Rasmussen, “Collapse arrest and soliton stabilization in nonlocal nonlinear media,” Phys. Rev. E **66**, 046619 (2002). [CrossRef]

**3. **J. Henninot, J. Blach, and M. Warenghem, “Experimental study of the nonlocality of spatial optical solitons excited in nematic liquid crystal,” J. Opt. A: Pure and Applied Opt. **9**, 20 (2007). [CrossRef]

**4. **X. Hutsebaut, C. Cambournac, M. Haelterman, J. Beeckman, and K. Neyts, “Measurement of the self-induced waveguide of a solitonlike optical beam in a nematic liquid crystal,” J. Opt. Soc. Am. B **22**, 1424–1431 (2005). [CrossRef]

**5. **J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, “Simulations and experiments on self-focusing conditions in nematic liquid-crystal planar cells,” Opt. Express **12**, 1011–1018 (2004). [CrossRef] [PubMed]

**6. **A. I. Yakimenko, Y. A. Zaliznyak, and Y. Kivshar, “Stable vortex solitons in nonlocal self-focusing nonlinear media,” Phys. Rev. E **71**, 065603 (2005). [CrossRef]

**7. **A. W. Snyder and D. J. Mitchell, “Accessible solitons,” Science **276**, 1538–1541 (1997). [CrossRef]

**8. **C. Conti, M. Peccianti, and G. Assanto, “Route to nonlocality and observation of accessible solitons,” Phys. Rev. Lett. **91**, 073901 (2003). [CrossRef] [PubMed]

**9. **C. Conti, M. Peccianti, and G. Assanto, “Observation of optical spatial solitons in a highly nonlocal medium,” Phys. Rev. Lett. **92**, 113902 (2004). [CrossRef] [PubMed]

**10. **J. Henninot, J. Blach, and M. Warenghem, “The investigation of an electrically stabilized optical spatial soliton induced in a nematic liquid crystal,” Journal of Optics A: Pure and Applied Optics **10**, 085104 (2008). [CrossRef]

**11. **N. B. Aleksić, M. S. Petrović, A. I. Strinić, and M. R. Belić, “Solitons in highly nonlocal nematic liquid crystals: Variational approach,” Phys. Rev. A **85**, 033826 (2012). [CrossRef]

**12. **D. Buccoliero, A. S. Desyatnikov, W. Krolikowski, and Y. S. Kivshar, “Boundary effects on the dynamics of higher-order optical spatial solitons in nonlocal thermal media,” J. Opt. A: Pure Appl. Opt. **11**, 094014 (2009). [CrossRef]

**13. **A. Litvak, “Self-focusing of powerful light beams by thermal effects,” Soviet Journal of Experimental and Theoretical Physics Letters **4**, 230 (1966).

**14. **F. Dabby and J. Whinnery, “Thermal self-focusing of laser beams in lead glasses,” Appl. Phys. Lett. **13**, 284–286 (1968). [CrossRef]

**15. **A. Litvak, V. Mironov, G. Fraiman, and A. Iunakovskii, “Thermal self-effect of wave beams in a plasma with a nonlocal nonlinearity,” Fizika Plazmy **1**, 60–71 (1975).

**16. **C. Rotschild, O. Cohen, O. Manela, M. Segev, and T. Carmon, “Solitons in nonlinear media with an infinite range of nonlocality: first observation of coherent elliptic solitons and of vortex-ring solitons,” Phys. Rev. Lett. **95**, 213904 (2005). [CrossRef] [PubMed]

**17. **M. R. Belić, D. Vujić, A. Stepken, F. Kaiser, G. F. Calvo, F. Agulló-López, and M. Carrascosa, “Isotropic versus anisotropic modeling of photorefractive solitons,” Phys. Rev. E **65**, 066610 (2002). [CrossRef]

**18. **L. Santos, G. V. Shlyapnikov, and M. Lewenstein, “Roton-maxon spectrum and stability of trapped dipolar bose-einstein condensates,” Phys. Rev. Lett. **90**, 250403 (2003). [CrossRef] [PubMed]

**19. **I. Tikhonenkov, B. Malomed, and A. Vardi, “Vortex solitons in dipolar bose-einstein condensates,” Phys. Rev. A **78**, 043614 (2008). [CrossRef]

**20. **I. Tikhonenkov, B. A. Malomed, and A. Vardi, “Anisotropic solitons in dipolar bose-einstein condensates,” Phys. Rev. Lett. **100**, 090406 (2008). [CrossRef] [PubMed]

**21. **M. S. Petrović, N. B. Aleksić, A. I. Strinić, and M. R. Belić, “Destruction of shape-invariant solitons in nematic liquid crystals by noise,” Phys. Rev. A **87**, 043825 (2013). [CrossRef]

**22. **G. Assanto and N. F. Smyth, “Comment on solitons in highly nonlocal nematic liquid crystals: Variational approach,” Phys. Rev. A **87**, 047801 (2013). [CrossRef]

**23. **N. B. Aleksić, M. S. Petrović, A. I. Strinić, and M. R. Belić, “Reply to comment on solitons in highly nonlocal nematic liquid crystals: Variational approach,” Phys. Rev. A **87**, 047802 (2013). [CrossRef]

**24. **B. Aleksić, N. Aleksić, M. S. Petrović, A. I. Strinić, and M. R. Belić, “Variational approach vs accessible soliton approximation in nonlocal nonlinear media,” arXiv preprint arXiv:1311.6840 (2013).

**25. **B. N. Aleksić, N. B. Aleksić, M. S. Petrović, A. I. Strinić, and M. R. Belić, “Variational approach versus accessible soliton approximation in nonlocal nonlinear media,” Phys. Scr. **T162**, 014003 (2014). [CrossRef]

**26. **X. Shi, B. A. Malomed, F. Ye, and X. Chen, “Symmetric and asymmetric solitons in a nonlocal nonlinear coupler,” Phys. Rev. A **85**, 053839 (2012). [CrossRef]

**27. **V. F. Zaitsev and A. D. Polyanin, *Handbook of Exact Solutions for Ordinary Differential Equations* (CRC Press, 2012).

**28. **M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Functions: With Formulars, Graphs, and Mathematical Tables*, vol. 55 (Dover Publications, 1964).

**29. **Q. Guo, B. Luo, F. Yi, S. Chi, and Y. Xie, “Large phase shift of nonlocal optical spatial solitons,” Phys. Rev. E **69**, 016602 (2004). [CrossRef]

**30. **B. Aleksić, N. Aleksić, V. Skarka, and M. Belić, “Using graphical processing units to solve the multidimensional ginzburg-landau equation,” Phys. Scr. **T149**, 014036 (2012). [CrossRef]

**31. **W. H. Press, *Numerical recipes 3rd edition: The art of scientific computing* (Cambridge University, 2007).