# An perturbation-iteration method for multi-peak solitons in nonlocal nonlinear media

###### Abstract

An perturbation-iteration method is developed for the computation of the Hermite-Gaussian-like solitons with arbitrary peak numbers in nonlocal nonlinear media. This method is based on the perturbed model of the Schrödinger equation for the harmonic oscillator, in which the minimum perturbation is obtained by the iteration. This method takes a few tens of iteration loops to achieve enough high accuracy, and the initial condition is fixed to the Hermite-Gaussian function. The method we developed might also be extended to the numerical integration of the Schrödinger equations in any type of potentials.

## I Introduction

Nonlocal spatial optical solitons have been a hot topic in the current research on nonlinear optics in the past few decades. In the Snyder and Mitchell’s workA. W. Snyder-Science-97 (), the nonlocal nonlinear Schrödinger equation (NNLSE), which models the optical beam propagation in the nonlocal nonlinear mediaAssanto-book (); Qi Guo-book-01 (), can be simplified to a linear equation, so-called Snyder-Mitchell (SM) model, in the strongly nonlocal case, and then an exact Gaussian-shaped stationary solution called an accessible soliton has been obtained. This pioneer work raised much attention, and many theoretical A. W. Snyder-Science-97 (); Assanto-book (); Qi Guo-book-01 (); A.W. Snyder-josab-97 (); A.W. Snyder-ijstqe-00 (); C. Conti-prl-03 () and experimentalAssanto-book (); Qi Guo-book-01 (); C. Conti-prl-04 (); M. Peccianti-ol-02 (); M. Peccianti-apl-02 (); M. Peccianti-pre-03 () works have been focused on the propagation of nonlocal spatial solitons.

It is well known that the complex-form solitons with multi-peak structure cannot be self-guided in the local nonlinear media due to the natural repulsion existing between lobes of opposite phase. It has been found by Christodoulides et al. in 1998 that multi-peak solitons are possible in saturable nonlinear mediaChristod-prl-98 (). Meanwhile, it is also found in nonlocal nonlinear media that the repulsion can be overcome by the nonlocalityA. V. Mamaev-pra-97 (); X. Hutsebaut-oc-04 (); A. I. Yakimenko-pre-06 (), and the multi-peak solitons have been observed both numerically D. W. McLaughlin-pd-96 (); D. Deng-josab-07 (); D. Buccoliero-prl-07 ()and experimentallyX. Hutsebaut-oc-04 (); C. Rotschild-ol-06 (), and their stability has also been discussedZ. Xu-ol-05 (). In the strongly nonlocal case, Multi-peak solitons can also be analytically obtained based on SM modelD. Deng-josab-07 (); Dongmei-jpb-2008 (); Dongmei-ol-2007 (); Dongmei-ol-2009 ().

Generally, most of the solitons defy analytical expressions and have to be computed numerically. So far, a number of numerical methods have been developed. Examples include the shooting methodJ. Yang-pre-02 (), the Petviashvili-type methodV.I. Petviashvili-pp-76 (); M.J. Ablowitz-ol-05 (); T.I. Lakoba-jcp-07 (), the imaginary time evolution method and its accelerated versionJ.J. Garcia-Ripoll-sjsc-01 (); J. Yang-sam-08 (), the Newton’s methodJ.P. Boyd-book-01 (); J.P. Boyd-jcp-02 (); William H. Press-book-92 (), the squared-operator iteration method and its modified versionJ. Yang-sam-07 (), Newton-conjugate-gradient (Newton-CG) methodsJ. Yang-jcp-09 (), etc. Generally, the first three methods can only converge to the ground statesT.I. Lakoba-jcp-07 (); J. Yang-sam-08 (); D.E. Pelinovsky-sjna-04 (); J. Yang-prl-99 (), i.e., the single-peak solitons. The Newton’s method has been widely used for iterating both single- and multi-peak solitons. However, it requires that the initial condition should be close enough to the exact solutionWilliam H. Press-book-92 (). The squared-operator iteration method and its modified version, based on the idea of time-evolving a “squared” operator equation, can converge to any solitons, including multi-peak solitonsJ. Yang-sam-07 (). However, these methods are quite slow, especially when the propagation gets near the edge of the continuous spectrum (detailed discussions can be found in Ref.J. Yang-jcp-09 () ). The Newton-CG method is based on Newton iteration, coupled with conjugate-gradient iterations to solve the resulting linear Newton-correction equation. It can be applied to compute both single- and multi-peak solitons and is faster than the other leading methods as concluded by the author. Although this method can tolerate wider ranges of initial conditions than that of the Newton method, it still require that the initial condition is reasonably close to the exact solutionJ. Yang-jcp-09 ().

In this paper, we develop a method for the computation of the Hermite-Gauss-like solitons with arbitrary peak numbers for arbitrary degree of nonlocality (as long as the nonlocality can support the soliton). The idea is inspired by the work of Ouyang et. al.Shigen Ouyang-pre-06 (), where the NNlSE is approximated to a linear equation of a perturbed harmonic oscillator, and the approximate analytical soliton solutions with the peak numbers up to three have been obtained in the second-order approximation using the standard procedure of the perturbation theory in quantum mechanicsW.Greiner-book-01 (). In the method we developed, named as perturbation-iteration (PI) method, we start from the perturbed model of the harmonic oscillator, determine the “minimum” perturbation by means of the weighted least-squares (WLS) methodWilliam H. Press-book-92 (), use the formal expression of infinite-order perturbation expansionsA.Messiah-book-1965 () to numerically calculate the eigenfunctions and eigenvalues of the perturbed model, and iterate this perturbation to obtained the multi-peak solitons with enough high accuracy. The initial condition of this method is fixed to the Hermite-Gaussian functions, since it iterates the perturbation rather than iterating the profile directly.

## Ii Methods

We start from the (1+1)-D dimensionless NNLSE which readsQi Guo-book-01 ():

(1) |

with the complex amplitude envelop of the optical beam, and the transverse and longitude coordinates respectively, and the nonlocal response function of the media. One can look for the multi-peak-soliton solutions of Eq. (1) with the form

(2) |

where is a coefficient related to the power of soliton, and are respectively the transverse distribution and propagation constant of soliton, and is the order of soliton. Note that , and are all real. corresponds to the first-order (fundamental) soliton with one peak, corresponds to the second-order soliton with two peaks, and so on. Then Eq. (1) becomes a stationary equation which reads:

(3) |

It has been discussed in Ref. Shigen Ouyang-pre-06 () that Eq. (3) can be treated as a perturbed model of the harmonic oscillator, and the approximate analytical soliton solutions up to has been obtained there by means of the perturbation theory in Quantum MechanicsW.Greiner-book-01 (). In order to apply the perturbation theory for searching the numerical soliton solution with arbitrary order , we start from a perturbed model of the harmonic oscillator with the following form

(4) |

where and are respectively the eigenfunctions and the eigenvalues of Eq. (4), is the perturbation compared to the term (the potential of the harmonic oscillator), and is expressed as with a constant. The eigenfunction series and the eigenvalues of the unperturbed equation for the harmonic oscillator, which is

(5) |

are W.Greiner-book-01 () the Hermite-Gaussian function (HGF) and respectively, where are the -order Hermite polynomials and are normalization coefficients such that . We emphasize here the meaning of the function : mathematically, is the difference between Eq. (3) and Eq. (5) at the zero-order approximation about the perturbation such that and ; physically, is the difference between the soliton-induced nonlinear refraction and the parabolic refraction .

For the given , and can be obtained using the standard procedure of the perturbation theoryW.Greiner-book-01 (). The formal expressions of the infinite-order perturbation expansions of and can be respectively expressed asA.Messiah-book-1965 ():

(6) |

(7) |

and

(8) |

with , and , where the summation indices and in Eqs.(6)-(8) are respectively the order of the expansion and the order of the HGF. Note that although the convergency of the expansions (6) and (7) is lacking in mathematical rigor, our numerical computations (presented in the next section) indicate that the expansions are in physical precision.

In the following, we are going to search a under the expected degree of the nonlocality , with the characteristic length of the media and the root-mean-square(RMS) width of the soliton (expressed as ). Since is unknown, we can use the RMS width of its zeroth-order approximation [] instead, which is W.Greiner-book-01 (), and then give an apriori value of which reads

(9) |

The reason is that mainly depend on according to the expansion (7), and the higher-order corrections in the expansion only slightly deviate the RMS width from , as will be discussed below. Therefore, by setting , we determine a “minimum” by means of the WLS methodWilliam H. Press-book-92 (), of which the procedure is as follows. For a given transverse distribution of the optical beam , the deviation, which is sampled by the optical beam, between the nonlinear refraction and the parabolic refraction is

(10) |

By letting the first derivative of with respect to and that to equal zero respectively, we have

(11) |

where , , , , , and . By easily solving this equation set, and can be obtained as

(12) |

Then, the perturbation which makes the deviation reaches its minimum is obtained through the relationships of Eq. (12). As a result, we can calculate and via Eqs. (6) and (7) by truncating the summations by the maxima of and (marked as and , respectively) when the error of the numerical solution can hardly decreases.

We can repeat the above procedure to achieve higher accuracy by the iteration. First, we let as an input wave function [ is the real wavefunction in the - loop, and algorithmize this procedure according to Eqs. (6), (7), and (12) through the following iterating process: within the - loop, we have

(13) |

(14) |

(15) |

(16) |

(17) |

At last, the soliton solutions of Eq. (1) are obtained as with , where is the number of iteration. Note that inside the loop is used to compute the error within the loop. In principle, for every -order of the expansion(17) can be individually chosen. For the practical convenience, however, is set to the same value for every -order in our computation below. .

For the obtained soliton solution with the form given by Eq. (2), the corresponding residual is

(18) | ||||

with

(19) |

Although the residual is the function of , we consider its average value within the region quadruple larger than the width of the soliton, since the data in the domain far from the soliton are meaningless for discussing the error. Obviously, consists of and that we defined as the relative residual, and is more reasonable than to be used to discuss the error of the numerical results, since of the solitons with different are quite different. The iteration process are terminated when the criterion is sufficiently small.

## Iii Examples

We show in Fig. 1 the computation of -peak solitons () for Gaussian and exponential-decay response functions which respectively readsQi Guo-book-01 (); Qi Guo-book-02 () and , where is the characteristic length of the media. In our computation, for convenience, we let to make the RMS widths of equal to , and then [see Eq. (9) above)], which is chosen to be in this example. The computational domain is taken as , discretized by 8192() points (The below computations are the same with these parameters). Figure 1(a) are the amplitude distributions of the solitons for the two responses. In order to chose proper and , we present the diagrams versus for different in Figs. 1(b) and (c) respectively for and . Here we terminate the iterations when is sufficiently small. For both two cases, it is clear that rapidly declines when increases and becomes nearly unchangeable. As keeps on increasing, also declines and becomes unchangeable. The results above indicate the best accuracy can be achieved when and are both large enough. In detail, for the -peak solitons, and are found for to reach the best accuracy of and and for to reach . However, for other achievable accuracies one required, it is difficult to obtain the proper values of and in a fast way. One can start for small values of and (typically, and ), fix the value of and increase , until reaches the required value; if cannot reach the required value, then increase and repeat the procedure above. By setting and for the best accuracy mentioned above, we further calculate the diagrams versus the number of iterations shown in Fig. 1(d) for the cases of the two responses. It is shown that the rapidly declines with the iteration, and becomes nearly unchangeable when the number of the iteration is larger than and respectively for the cases of and .

Our method can be used to calculate multi-peak solitons with an arbitrary order of soliton . For the given , the computer times for the solitons with different are shown in Figs. 2(a) and (b) for and , respectively. In addition, the needed values of the pair (, ) and the iteration number for every soliton are marked in the figures. For comparison, the results taken by Newton-CG methodJ. Yang-jcp-09 () are also presented. It is shown that our method takes even shorter times than the Newton-CG method for , but much longer times for , since the case of takes much larger values of , and than that of . Even so, our method is quite accessible, since the cost computer times are in the order of seconds, and the initial condition is fixed to , rather than careful choosing of an initial profile. Note that the Newton-CG method requires that the initial condition is reasonably close to the exact solutionJ. Yang-jcp-09 (). In our realization of the Newton-CG method, therefore, we choose (with the initial amplitude) as the initial condition, and should be carefully chose to ensure the convergence of the iteration. Moreover, the sensitivity of also increases as the peak number increases, according to our computation.

We also present the RMS widths of the solitons with different as a function of the characteristic length , i.e, the degree of the nonlocality, for and , which are shown in Figs. 3(a) and (b), respectively. The request values of are the same with those in Fig. 2. As expected above, the higher-order corrections in the expansion (7) only slightly deviate the RMS widths from . As increases, for the case of tends to , while that for the case of tends to a value deviated from . This asymptotic behavior can be understood as followsQi Guo-book-01 (); Qi Guo-book-02 (): as increases, namely, the degree of the nonlocality increases, Eq. (3) with tends to the SM model, i.e., the model for harmonic oscillator described by Eq. (5) A. W. Snyder-Science-97 (), and then tends to ; for the case of , in the strongly nonlocal limit, i.e., , Eq. (3) does not tend to Eq. (5) due to the discontinuity in the derivative of at the origin (detailed discussions can be found in Qi Guo-book-01 (); Qi Guo-book-02 ()), therefore this deviation is in expectation. It is worth mentioning that the best of achievable is not sensitive to , i.e., the degree of the nonlocality, according to our computation. Figs. 3(c) and (d) are the power of the solitons (shown in logarithmic coordinate) with different as the function of for and , respectively. As and increase, the power of the solitons increases rapidly, which is the reason why is more proper than to discuss the error of the solitons, as mentioned above.

Finally, we show in Fig. 4 the computation of the 2-peak soliton with the sine-oscillation response sine1 (); sine2 (), which could hardly be obtained by other known methods as far as we have tried. The best accuracy of is achieved when and . The physical properties for such multi-peak solitons with will be further discussed otherwhere.

## Iv Conclusion

In this work, we develop a method for the numerical computations of the Hermite-Gauss-like solitons with arbitrary peak numbers in nonlocal nonlinear media. We start from the perturbed model of the harmonic oscillator, determine the “minimum” perturbation by means of the WLS method, and use the formal expression of infinite-order perturbation expansions to numerically calculate the eigenfunctions and eigenvalues of the perturbed model. Finally we iterate the procedure above to achieve a sufficient high accuracy. This method takes a few tens of iteration loops, and the initial condition is fixed to the Hermite-Gaussian function. It is worth mentioning that our method is not restricted to the NNLSE, but might also be able to be extended to the Schrödinger equations in any type of potentials.

## V Funding Information

National Natural Science Foundation of China (Grant No. 11474109).

## References

- (1) A. W. Snyder , D. J. Mitchell, “Accessible solitons,” Science, 276, 1538-1541 (1997).
- (2) G. Assanto, Nematicons: Spatial Optical Solitons in Nematic Liquid Crystals (New Jersey: John Wiley & Sons, 2013).
- (3) Qi Guo, Daquan Lu, and Dongmei Deng, Nonlocal spatial optical solitons, Chapter 4, in Advances in Nonlinear Optics, edited by X. Chen, Qi Guo, W. She, H. Zeng, and G. Zhang, (De Gruyter, Berlin, Feb. 2015), p.p.227-305.
- (4) A.W. Snyder and Y. Kivshar, “Fright spatial solitons in non-Kerr media: stationary beams and dynamical evolution ,” J. Opt. Soc. Am. B 11, 3025-3031 (1997).
- (5) A.W. Snyder, “Guiding light into the millennium ,” IEEE J. Sel. Top. Quantum Electron. 6, 1408-1412 (2000).
- (6) C. Conti, M. Peccianti, G. Assanto, “ Route to nonlocality and observation of accessible solitons,” Phys. Rev. Lett. 91, 073901 (2003).
- (7) C. Conti, M. Peccianti, G. Assanto, “Observation of optical spatial solitons in highly nonlocal medium,” Phys. Rev. Lett. 92, 113902 (2004).
- (8) M. Peccianti, C. Conti, G. Assanto,“All-optical switching and logic gating with spatial solitons in liquid crystals ,” Appl. Phys. Lett. 81, 3335-3337 (2002).
- (9) M. Peccianti, K. A. Brzdakiewicz, and G. Assanto,“Nonlocal spatial soliton interactions in nematic liquid crystals, ” Opt. Lett. 27, 1460-1462 (2002).
- (10) M. Peccianti, C. Conti, G. Assanto,“Optical modulational instability in a nonlocal medium, ” Phys. Rev. E 68, 025602 (2003).
- (11) D. N. Christodoulides, T. H. Coskun, M. Mitchell, Z. G. Chen, M. Segev, “Theory of incoherent dark solitons,” Phys. Rev. Lett. 80, 2310-2313 (1998).
- (12) A. V. Mamaev, A. A. Zozulya, V. K. Mezentsev, D. Z. Anderson, and M. Saffman,“ Bound dipole solitary solutions in anisotropic nonlocal self-focusing media ,” Phys. Rev. A 56, 1110-1113 (1997).
- (13) X. Hutsebaut, C. Cambournac, M. Haelterman, A. Adamski, and K. Neyts,“Single-component higher-order mode solitons in liquid crystals,” Opt. Commun. 233, 211-217(2004).
- (14) A. I. Yakimenko, V. M. Lashkin, and O. O. Prikhodko, “Dynamics of two-dimensional coherent structures in nonlocal nonlinear media,” Phys. Rev. E 73, 066605 (2006).
- (15) D. W. McLaughlin, D. J. Muraki, M. J. Shelley, “Self-focussed optical structures in a nematic liquid crystal,” Physica D 97, 471-497 (1996).
- (16) D. Deng, X. Zhao, Q. Guo , S. Lan, “Hermite-Gaussian breathers and solitons in strongly nonlocal nonlinear media,” J. Opt. Soc. Am. B 24,2537-2544 (2007).
- (17) D. Buccoliero, A. S. Desyatnikov, W. Krolikowski, and Y. S. Kivshar, “Laguerre and Hermite soliton clusters in nonlocal nonlinear media,” Phys. Rev. Lett. 98, 053901 (2007).
- (18) C. Rotschild, M. Segev, Z. Y. Xu, Y. V. Kartashov, L. Torner, and O. Cohen, “Two-dimensional multipole solitons in nonlocal nonlinear media,” Opt. Lett. 31, 3312-3314 (2006).
- (19) Z. Xu, Y. V. Kartashov, and L. Torner, “Upper threshold for stability of multipole-mode solitons in nonlocal nonlinear media,” Opt. Lett. 30, 3171-3173 (2005).
- (20) Dongmei Deng, Qi Guo, and Wei Hu, “Hermite-Laguerre-Gaussian beams in strongly nonlocal nonlinear media,” J. Phys. B: At. Mol. Opt., 41, 225402 (2008).
- (21) Dongmei Deng, and Qi Guo, “Ince-Gaussian solitons in strongly nonlocal nonlinear media,” Opt. Lett., 32, 3206-3208 (2007).
- (22) Dongmei Deng, Qi Guo, and Wei Hu, “Comeplex-variable-function-Gaussian solitions,” Opt. Lett, 34, 43-45 (2009).
- (23) J. Yang,“Internal oscillations and instability characteristics of (2+1)-dimensional solitons in a saturable nonlinear medium,” Phys. Rev. E 66, 026601 (2002).
- (24) V.I. Petviashvili, “On the equation of a nonuniform soliton,” Plasma Phys. 2, 469-472 (1976).
- (25) M.J. Ablowitz, Z.H. Musslimani,“Spectral renormalization method for computing self-localized solutions to nonlinear systems,” Opt. Lett. 30, 2140-2142 (2005).
- (26) T.I. Lakoba, J. Yang,“A generalized Petviashvili iteration method for scalar and vector Hamiltonian equations with arbitrary form of nonlinearity,” J. Comput. Phys. 226, 1668-1692 (2007).
- (27) J.J. Garcia-Ripoll, V.M. Perez-Garcia, “Optimizing Schrodinger functionals using Sobolev gradients: Applications to quantum mechanics and nonlinear optics,” SIAM J. Sci. Comput. 23, 1315-1333 (2001).
- (28) J. Yang, T.I. Lakoba,“Accelerated imaginary-time evolution methods for the computation of solitary waves,” Stud. Appl. Math. 120, 265-292 (2008).
- (29) J.P. Boyd, Chebyshev and Fourier Spectral Methods (Dover Publications, second edition, 2001).
- (30) J.P. Boyd, “Deleted residuals, the QR-factored Newton iteration, and other methods for formally overdetermined determinate discretizations of nonlinear eigenproblems for solitary, cnoidal, and shock waves,” J. Comput. Phys. 179, 216-237(2002).
- (31) William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T. Vetterling, Numerical Recipes in C: the Art OF Scientific Computing (Cambridge University Press,second edition,1992)
- (32) J. Yang, T. I. Lakoba, “Universally-convergent squared-operator iteration methods for solitary waves in general nonlinear wave equations,” Stud. Appl. Math. 118, 153-197 (2007).
- (33) Jianke Yang, “Newton-conjugate-gradient methods for solitary wave computations,” J. Comput. Phys. 228, 7007-7024 (2009).
- (34) D.E. Pelinovsky, Y.A. Stepanyants, “Convergence of Petviashvili’s iteration method for numerical approximation of stationary solutions of nonlinear wave equations,” SIAM J. Numer. Anal. 42, 1110-1127 (2004).
- (35) J. Yang, B.A. Malomed, D.J. Kaup, “Embedded Solitons in Second-Harmonic-Generating Systems,” Phys. Rev. Lett. 83, 1958 (1999).
- (36) Shigen Ouyang, Qi Guo, and Wei Hu, “Perturbative analysis of generally nonlocal spatial optical solitons,” Phys. Rev. E 74, 036622 (2006).
- (37) W.Greiner, Quantum Mechanics An Introduction (Springer-Verlag, 4th edition,Berlin,2001).
- (38) A. Messiah, Quantum Mechanics (John Wiley & Sons, Inc, New York,1965), Vol.2, p.688.
- (39) In the Newton-CG method, the residuals diverges for the sine oscillation response as far as we can try.
- (40) Qi Guo, Wei Hu, Dongme Deng, Daquan Lu, and Shigen Ouyang, Features of strongly nonlocal spatial solitons, Chapter 2, in Nematicons: Spatial Optical Solitons in Nematic Liquid Crystals, edited by Gaetano Assanto (John Wiley & Sons, New Jersey, 2013), p.p.37-69.
- (41) B. K. Esbensen, M. Bache, W. Krolikowski, and O. Bang, “Quadratic solitons for negative effective second-harmonic diffraction as nonlocal solitons with periodic nonlocal response function,” Phys. Rev. A 86, 023849 (2012).
- (42) J. Wang, Y. Li, Q. Guo, and W. Hu, “Stabilization of nonlocal solitons by boundary conditions,” Opt. Lett. 39, 405-408 (2014).