Tripole-mode and quadrupole-mode solitons in (1 + 1)-dimensional nonlinear media with a spatial exponential-decay nonlocality

The approximate analytical expressions of tripole-mode and quadrupole-mode solitons in (1 + 1)-dimensional nematic liquid crystals are obtained by applying the variational approach. It is found that the soliton powers for the two types of solitons are not equal with the same parameters, which is much different from their counterparts in the Snyder-Mitchell model (an ideal and typical strongly nolocal nonlinear model). The numerical simulations show that for the strongly nonlocal case, by expanding the response function to the second order, the approximate soliton solutions are in good agreement with the numerical results. Furthermore, by expanding the respond function to the higher orders, the accuracy and the validity range of the approximate soliton solutions increase. If the response function is expanded to the tenth order, the approximate solutions are still valid for the general nonlocal case.

solitary wave in a nematic liquid crystal 33 , based on the previous work for the NLS equation performed by Kath et al. 34 . The work in refs 33 and 34 does more than find the steady solitary wave, it also finds the evolution to this solitary wave from an initial condition 33,34 . Malomed presented a general review of these variational methods in nonlinear fiber optics and related fields 35 . Recently, Aleksić et al. analytically investigated the fundamental solitons based on the variational approach in (2 + 1)-dimensional NLCs 36 . Especially, MacNeil et al. obtained exact solutions of the nematicon equations in (1 + 1) and (2 + 1) dimensions for fixed parameter values, and they also got approximate solutions based on the variational approach method 37 . Furthermore, Panayotaros and Marchant addressed the existence of a solitary wave solution of the nematic equations mathematically 38 . It has been proven that in the media with an exponential-decay nonlocal response, the soliton bound states are stable if the solitons contain fewer than five-poles 21 . Namely the fundamental, dipole-mode, tripole-mode and quadrupole-mode solitons can all propagate stably in such media. The fundamental and dipole-mode solitons in nonlinear media with an exponential-decay nonlocal response have been investigated analytically by various mathematical methods, such as the classical Lie-group method 39 , the perturbative analysis method 40 , and the variational approach method 36,37,41,42 . But so far, no one gives analytical expressions of tripole-mode and quadrupole-mode solitons in NLCs. Presenting an analytical solution is helpful for one to get a good understanding of the dynamics of nonlocal solitons. In this paper, based on the variational approach, we study the tripole-mode and quadrupole-mode solitons in nonlinear media with an exponential-decay nonlocal response. The approximate expressions of such solitons are obtained and the characteristics of them are investigated in detail. The approximate results are confirmed by the numerical ones which are obtained using the iterative numerical technique based on the NNLSE directly. Since a surface soliton in nonlocal nonlinear media can be regarded as a half of a bulk soliton with an antisymmetric amplitude distribution 43,44 , the results on quadrupole-mode solitons here may also be helpful for the investigation of the surface dipole nonlocal solitons.

Results
Variational method for NNLSE and tripole-mode soliton solutions. First of all, let us roughly recall the derivation of the dimensionless NNLSE for the (1 + 1)-dimensional NLCs based on refs 1, 31 and 45. Considering only one transversal dimension 45 , an external quasistatic electric field E LF is applied in the transversal direction to control the initial tilt angle of the NLCs. The evolution of an optical beam Q(X, Z) in the paraxial approximation and optically induced reorientation angle perturbation Ψ X Z ( , ) of the liquid crystal molecules can be described as follows 1, 31 where k 0 and k are, respectively, the wave numbers in vacuum and the NLCs. Δε HF and Δε LF are, respectively, the anisotropy of the liquid crystal at the optical frequency and the anisotropy of the liquid crystal at the frequency of the quasistatic electric field. ε 0 is the vacuum permittivity. K is the relevant elastic constant taken equal for splay, bend, and twist. Introducing the normalization: x = X/W 0 , z = Z/Z R , q = Q/Q 0 , and ψ = Ψ Ψ / 0 , where W 0 is the full width at half maximum of t he amplitude of t he opt ic a l b e am, , one can obtain the following dimensionless equations where σ denotes the degree of nonlocality which is expressed as 31 In Eq. (4), the term ∂ 2 ψ/∂z 2 has been canceled since , which is a part of the coefficient multiplying the derivative ψ zz in the non-dimensional director equation 1,3 . Based on Eq. (4), combined with Fourier transformation and the convolution theorem, one can find where I = |q| 2 is the intensity of the optical beam, and the normalized nonlocal response function R takes an exponential-decay function as follows Scientific RepoRts | 7: 122 | DOI:10.1038/s41598-017-00197-6 Several kinds of solitons in (1 + 1)-dimensional nonlinear media with a spatial exponential-decay nonlocality has been investigated in the past years 46,47 . In experiments 4,[9][10][11] , the typical values of the parameters are: W 0 = 10 μm, K = 10 −11 N, Δε HF = 0.64ε 0 , Δε LF = 15ε 0 , E LF = 10 4 V/m. For the visible wavelengths, ε 0 is 8.85 × 10 −12 in MKS units, and the diffraction length is about 1 mm. With the above parameters, the degree of nonlocality is about 12 , i.e. σ ≈ 12 2 , which belongs to the sub-strongly nonlocal case. If σ → w m and ψ → Δn, the dimensionless NNLSE, which governs the beam propagation in (1 + 1)-dimensional nonlinear media with an exponential-decay nonlocal response, can be rewritten phenomenologically as follows 17,31  where x and z denote, respectively, the normalized transversal and longitudinal coordinates; q denotes the complex amplitude of the optical beam; Δn denotes the nonlinear perturbation of the refractive index of the nonlocal medium; w m denotes the characteristic length of the nonlocal material response. If w m → 0, Eq. (8) is simplified to the well-known nonlinear Schrödinger equation in local nonlinear media; if w m~wR (where w R is the second-order moment width of multipole solitons), it represents the general nonlocal case; and if w m → ∞, it represents the strongly nonlocal case 48,49 . In experiments, the magnitude of w m can be controlled by changing the pretilt angle of NLCs via a bias voltage 1,3,17 . Based on the above derivation, Eq. (8), together with Eq. (9), denotes a NNLSE with an exponential-decay nonlocal response 31 , which can be restated as follows, where the Lagrangian density can be expressed as  (12), the reduced variational problem can be obtained as follows 0 where the average Lagrange is expressed as For the strongly nonlocal media (especially for the Snyder-Mitchell model), the higher-order Gaussian solitons (such as Hermite-Gaussian solitons, Laguerre-Gaussian solitons, Ince-Gausian solitons etc.) are the exact soliton solutions 12,15,50,51 . For (1 + 1)-dimensional nonlinear media with a spatial exponential-decay nonlocality, we consider that solitons come into being from the Hermite-Gaussian beams. In our previous research, we have proven that the first-order Hermite-Gaussian function can be used to describe the dipole solitons in such nonlocal media 52 . Hence, we take the second-and third-order Hermite-Gaussian functions as the trial functions of tripole and quadrupole solitons, respectively. The trial solution of tripole-mode solitons takes the following form where a is the amplitude, θ(z) is the phase, and w is the width of a Gaussian soliton. Because of the complexity of the intensity distribution, the second-order moment beam width is adopted to describe the width of multipole solitons. The second-order moment beam widths is defined as Thus the width of a tripole-mode soliton is = w w 10 R . For the sake of convenience in the following discussion, we introduce a nonlocal parameter α to define the degree of the material nonlocality, i.e., α = w m /w R . The larger the nonlocal parameter, the stronger the degree of nonlocality.
In theory, based on Eqs (11), (13), (15) and (16), one can obtain the expression of [ ]  . However, the integrals in the averaged Lagrangian based on the trial function could not be calculated explicitly due to the inability to find closed form integrals. Fortunately, for the strongly nonlocal case, we can calculate it by expanding the response function. If it is expanded to the second order, one can get Substituting Eqs (13), (16) and (18) into Eq. (15), the expression of [ ]  is obtained as follows . Based on the corresponding Euler-Lagrangian equations, one can get As we know, for the soliton case, θ′(z) = βz, where β is the propagation constant. Combining Eqs (20) and (21) It is evident that Eqs (22)(23)(24) are valid only when π < w ww 32 2 29 , a becomes an imaginary number, β < 0, and P < 0, which is impossible in physics. Figure 1 shows the propagation constant of tripole-mode solitons versus the soliton powers. In Fig. 1(a), the degree of nonlocality is 7, which belongs to the strongly or at least sub-strongly nonlocal case. It is found that the approximate result is in good agreement with the numerical one which is obtained directly based on Eqs (8) and (9) using the iterative numerical technique 53 . When w m is fixed at 10, the approximate result is also in good agreement with the numerical ones as shown in Fig. 1(b). Figure 1(b) also shows that the approximate result it is invalid when β < 3.68 as the variational solution (22)-(24) breaks down, as discussed after Eq. (24). The reason for the invalidity is that the response function is only expanded to second order [see Eq. (18)], which leads to the inaccuracy. By expanding the respond function to higher orders, one can improve the accuracy of the approximate solutions, which will be discussed in the following Section. In addition, one can find from Fig. 1 that the slope of the power versus propagation constant is positive, which implies that the soliton propagation is stable. It is also found that if w m takes a larger value, the valid region of β also becomes larger. The accuracy of approximate results is only dependent on the degree of nonlocality. The stronger the degree of nonlocality, the more accurate the approximate results. As a result, when the degree of nonlocality is still fixed at 7, the analytical result is always accurate independent of the soliton width [see Fig. 2(a)]. Nevertheless, if w m is fixed, the degree of nonlocality decreases with the increase of the soliton width, so the validity of approximate results is getting declined continuously [see Fig. 2(b)]. The variational approximation shows bistability in Fig. 2(b), while the numerical solution does not. This is an important point as it shows that the variational approximation can predict behaviour which is not actually present. Caution should then be exercised with variational approximations. Figure 3 shows the profiles of the tripole-mode soliton with different soliton powers. It is found that the approximate solutions agree well with the numerical solutions for β = 30 and 20. When β decreases to 7, the approximate solution becomes a little inaccurate, and when β decreases to 4, it becomes even worse. The degrees of nonlocality are 3.19, 2.78, 1.93, and 1.59 for β = 30, 20, 7, and 4, respectively. It is found that when the nonlocality degree is 1.59, i.e. β = 4, an obvious deviation appears between the variational solution and the numerical one. The reason is that the Taylor series truncation (18) is starting to break down for this low w m . So we can conclude that if the response fuction is expanded to the second order, the approximate solutions are valid only for the strongly nonlocal case.

Quadrupole-mode soliton solutions.
For the quadrupole-mode solitons in (1 + 1)-dimensional NLCs, we take the following ansatz solution   . Comparing with the SMM which is valid for the strongly nonlocal case, it is found that the soliton powers are different. The soliton powers for tripole-mode and quadrupole-mode solitons in (1 + 1)-dimensional NLCs are not equal, and it is larger for quadrupole-mode solitons than tripole-mode solitons with the same parameters. However the multipole solitons and the higher-order solitons are all the same in SMM [12][13][14][15][16] . Furthermore, for the the strongly nonlocal limit (i.e. , which indicates the nonlocality degree decreases. As a results, the approximate result becomes invalid gradually. Especially when β < 7.17, approximate results do not exist any more [see Fig. 4(b)]. Figure 4 shows that the slope of the power versus propagation constant is positive, similar to the case of tripole solitons, which implies a stable propagation of solitons. Figure 6 presents the profiles of the quadrupole-mode soliton with different soliton powers and propagation constants. It is found that the approximate solutions agree well with the numerical solutions for β = 50 and 30. When β decreases to 12, the approximate solution has a little inaccuracy, and when β decreases to 8, it becomes even worse. We also calculate the nonlocal parameter α, and find that it is equal to 3.02, 2.53, 1.85, and 1.61 for β = 50, 30, 12, and 8, respectively.
More accurate approximate solutions. In the above sections, the response function is only expanded to the second order. So with the degree of nonlocality decreasing, the approximate results become invalid gradually. In order to get a more exact approximate solution, one can expand the response function to the higher orders. The  higher orders the response function is expanded into, the more exact the approximate solutions are. Of course, the calculations are more complicated. As an example, we expand R(x) to the fourth order, i.e., . In the previous section, it is found that the approximate tripole-mode solution is inaccurate when β = 4, and even it is not obtained when β < 3.68. In order to show the improvement of the approximate solutions obtained by expanding the respond function to the fourth order, Fig. 7 shows the comparison between the approximate results and the numerical ones. It is found that the approximate solutions are in good agreement with the numerical ones when β = 4 and even β = 2. For the case of β = 2, the degree of nonlocality is about 1.23, which already belongs to the general nonlocality. When β = 1, the approximate solution has a little inaccuracy, and when β = 0.8, it becomes even worse. The degrees of nonlocality for β = 1 and β = 0.8 are, respectively, 0.91 and 0.76. Furthermore, if R(x) is expanded to the tenth order, i.e.,  and It should be noted that Eqs (35)(36)(37) are valid only when A + B > 0. Figure 8 shows the profiles of the tripole-mode soliton obtained by expanding the response function to the tenth order. When β = 0.5, 0.35 and 0.12, the corresponding degrees of nonlocality are, respectively, 0.753, 0.657   and 0.417, which all belong to the case of the general nonlocality. Therefore, the approximate solutions can be improved by expanding the response function to the higher orders.
Although it has been proven that in a nonlinear medium with an exponential-decay nonlocal response, the soliton bound states are stable if the solitons contain fewer than five-poles 21 . In order to confirm the validity of our results, we take the tripole-mode soliton as an example and simulate the propagation based on Eqs (8) and (9) directly. Here we take the split-step Fourier method 54 to simulate the soliton propagation. As expected, for the strongly nonlocal case, the approximate results obtained by expanding the response function to the second order are accurate, and the solitons can stably propagate for a long distance [see Fig. 9(b)]. With the decrease of the nonlocality degree, the approximate results become inaccurate gradually, and irregular oscillations occur during propagation [see Fig. 9(a)]. However, if the nonlocal response function is expanded to the tenth order, the accuracy and the validity range of the approximate solutions increase. The approximate soliton can still keep a relatively stable propagation [see Fig. 9(d)], but it can not even be obtained if the response function is only expanded to the second order. In Fig. 9(c), because the degree of nonlocality is weaker than that in Fig. 9(d), the irregular oscillations appear more obviously. The oscillations of Fig. 9(c,d) are just typical behaviour for NLS-type equations. For such equations an initial condition near a solitary wave will evolve to the solitary wave with the amplitude and width displaying decaying oscillations. All these oscillations shows the degree of accuracy of the variational solutions. Therefore, we can conclude that by expanding the respond function to the higher orders, the accuracy of the approximate soliton solutions is improved. Corresponding to the cases of Figs 9(a,b) and 10 illustrates the propagation of the tripole-mode solitons by expanding the response function to the tenth order, which shows a stable propagation of solitons. Especially, when β = 4, the approximate soliton obtained by expanding the response function to the tenth order is still valid [see Fig. 10(a)]. Contrarily, it is already invalid when the response function is expanded to the second order. As another example, Fig. 11 illustrates the stable propagation of a quadrupole-mode soliton, which confirms the validity of the approximate variational quadrupole-mode soliton solutions.
At last, for the completeness, we also give the solutions of quadrupole-mode solitons as follows when the response function is expanded to the tenth order.  Figure 9. Propagation of tripole-mode solitons in the presence of 1% white input noise. The profiles of Fig. 3(a,d) are employed as the input shapes of solitons in (a,b), and the profiles of Fig. 8(a,c) are employed as the input shapes of solitons in (c,d).

Discussion
By applying the variational approach, we obtain the approximate analytical expressions of tripole-mode and quadrupole-mode solitons in nonlinear media with an exponential-decay nonlocal response. It is found that with the same parameters, the soliton power of the quadrupole-mode solitons is larger than that of the tripole-mode solitons, which is much differnt from the SMM (In SMM, the soliton powers with different multipoles are the same 12,15 ). The numerical simulations are carried out to illustrate the accuracy of the approximate solutions. The results show that the accuracy of the approximate solutions is only related with the degree of nonlocality. For the strongly nonlocal case, if the response function is expanded to the second order, the approximate soliton solutions are in good agreement with the numerical ones. With the degree of nonlocality decreasing, the approximate solutions become invalid gradually. Furthermore, by expanding the respond function to the higher orders, one can improve the accuracy of the approximate solutions. The higher orders the response function is expanded to, the more exact the approximate solutions are. If the response function is expanded to the tenth order, the approximate solutions are still valid for the general nonlocal case. Since a surface soliton in nonlocal nonlinear media can be regarded as a half of a bulk soliton with an antisymmetric amplitude distribution 43,44 , the results on quadrupole-mode solitons here may also be helpful for the investigation of the surface dipole nonlocal solitons.