1. Introduction

In present years spatial solitons in Kerr-type nonlocal media have attracted a great amount of studies. The nonlocal bright spatial solitons have been experimentally observed in nematic liquid crystal[1, 2, 3] and in the lead glass[4]. The propagation and interaction properties of nonlocal bright/dark solitons are greatly different from that of local solitons. For example, the dependent functions of the beam power and phase constant on the beam width for nonlocal bright solitons are greatly different from those for local bright solitons[5, 6, 7]; There exist stationary higher order bright nonlocal soliton solutions[5, 6, 7 , 8, 9, 10, 11]; Bright nonlocal solitons with π phase difference attract rather repel each other[11, 12, 13, 14]; Two dark/gray nonlocal solitons can attract each other and form a bound state[15, 16, 17].

In this paper we use the spectral renormalization method[18, 19, 20] to compute dark/gray soliton solutions in Kerr-type nonlocal media and find it agrees well with the analytical results. Comparing to other numerical methods like the shooting method[21] and Newton relaxation method[22], it seems more convenient to use the spectral renormalization method to solve the nonlinear integro-differential equations and to solve the soliton profiles of multiple dimensions[18, 19, 20]. It was claimed in Ref. [17] that the gray soliton’s maximal velocity monotonically decreases with the characteristic nonlocal length. Contrarily, as will be shown in our paper by using the same model of Ref. [17], such maximal velocity does not vary with the characteristic nonlocal length when the latter is less than a certain critical value. Nevertheless this maximal velocity does decrease with the characteristic nonlocal length when the latter larger than this critical value. Such discrepancy will be discussed in this paper.

2. Nonlocal dark and gray solitons

The propagation of a (1+1) dimensional optical beam in Kerr-type nonlocal self-defocusing media is governed by the nonlocal nonlinear Schödinger equation(NNLSE)[11]

where u(x,z) is the complex amplitude envelop of the light beam, ∣u(x,z)∣² is the light intensity, x and z are transverse and longitude coordinates respectively, R((x), (∫R(x)dx = 1) is the real symmetric nonlocal response function, and n(x,z) = - ∫R(x - ξ) ∣u(ξ,z)² dξ is the light-induced perturbed refractive index. Note that not stated otherwise all integrals in this paper will extend over the whole x axis. Itis easy to prove that under the scaling transformation x̄ = x/η ξ,̄] = ξ/η,z̄ = z/η ²,ū = η, u, R̄ = ξR, the form of Eq.(1) keeps invariant. Owing to the scaling transformation, it is adequate to consider dark/gray soliton solutions with an unit background intensity ${\mid u\left(x,z\right)\mid }^2\stackrel{x\to \pm \infty }{\to }$1. When R(x) = δ(x), equation (1) will reduce to the NLSE $i\frac{\partial u}{\partial z}+\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}-{\mid u\mid }^{2}u=0$, which has dark/gray soliton solutions[23] u(x,z) = ψ(x)e ^iβz+iϕ(x), where ψ(x) = [1 - β ²sech²(Bx)]^1/2, $\varphi \left(x\right)=\sqrt{1-B^{2}x}+\arctan [B\tanh \left(\mathrm{Bx}\right)/\sqrt{1-B^2}]$ and $\beta =-\frac{1}{2}\left(3-B^2\right)$ is the propagation constant, B is the grayness. It is easy to check that $\frac{\mathrm{d\varphi }\left(x\right)}{\mathrm{dx}}\stackrel{x\to \pm \infty }{\to }\mu =\sqrt{1-B^2}$ and $\beta =-\left(\frac{{\mu }^2}{2}+1\right)$, where -μ is the gray soliton’s transverse velocity relative to its background intensity.

In this paper we numerically compute the dark/gray soliton solution of Eq. (1) that takes an ansatz u(x,z) = ψ(x)e ^iβz+iϕ(x), where ψ(x),ϕ(x) are real functions and β is a real constant, and $\psi \left(x\right)\stackrel{x\to +\infty }{\to }1,\frac{\mathrm{d\varphi }\left(x\right)}{\mathrm{dx}}\stackrel{x\to +\infty }{\to }\mu$. Inserting this ansatz into Eq. (1) and using asymptotic behaviors $\psi \left(x\right)\stackrel{x\to +\infty }{\to }\mathrm{}1$ and $\frac{\mathrm{d\varphi }\left(x\right)}{\mathrm{dx}}\stackrel{x\to +\infty }{\to }\mu ,$, we can find $\beta =-\left(\frac{{\mu }^2}{2}+1\right),\frac{\mathrm{d\varphi }}{\mathrm{dx}}=\frac{\mu }{{\psi }^2}$and

When μ ≠ 0, the term $\frac{{\mu }^2}{2{\psi }^3}$of Eq. (2) requires ψ(x) ≠ 0, which corresponds to even gray soliton solutions ψ(-x) = ψ(x). We study the case μ ≠ 0 firstly and leave the case μ = 0 to be discussed in the latter part of this section. By setting χ(x) = 1 - ψ(x), from Eq. (2) we have

where $N\left(\chi \right)\equiv \frac{{\mu }^2\left(3{\chi }^3-3\chi -{\chi }^3\right)}{2{\left(1-\chi \right)}^3}-\frac{{\mu }^2\chi }{2}-\left(1-\chi \right)\mathrm{\int R}\left(x-\xi \right)\left[{\chi }^2\left(\xi \right)-2\chi \left(\xi \right)\right]\mathrm{d\xi }$. Equation (3) is a nonlinear integro-differential equation to be solved under conditions χ′ (0) = 0 and $\chi \left(x\right)\stackrel{x\to +\infty }{\to }0$. Now we can use the spectral renormalization method presented in Ref. [18] to compute χ(x). To illustrate how to use this method to compute χ(x), we outline the main steps here. Define the Fourier transform by y0342;(k) = ∫y(x) exp(ikx)dx. Applying the Fourier transform to Eq. (3), we arrive at $\frac{k^2}{2}\tilde{\chi }+\tilde{N}\left(\chi \right)=0$. Let χ(x) hθ(x), where h is a constant to be determined, we have $\frac{k^2}{2}\tilde{\theta }+\frac{1}{h}\tilde{N}\left(\mathrm{h\theta }\right)=0$, or, alternatively $\tilde{\theta }={퓓}_h\left[\theta \right]\equiv \left[r\tilde{\theta }-\frac{1}{h}\tilde{N}\left(\mathrm{h\theta }\right)\right]/\left[r+\frac{k^2}{2}\right]$, where r is an arbitrary positive constant introduced to avoid a vanishing denominator[18]. Projecting $\frac{k^2}{2}\tilde{\theta }+\frac{1}{h}\tilde{N}\left(\mathrm{h\theta }\right)=0,$ onto θ, we obtain ${퓐}_{\theta }\left(h\right)\equiv \int {\tilde{\theta }}^{*}\left(k\right)\left[\frac{k^2}{2}\tilde{\theta }+\frac{1}{h}\tilde{N}\left(\mathrm{h\theta }\right)\right]\mathrm{dk}=0,$. We use 𝓐_θ(h) = 0 and θ͂ = 𝓓_h[θ] to iteratively compute h and θ. For an arbitrary initial function θ ₁(x), e.g. a Gaussian function, by solving equation 𝓐 _θ1(h ₁) = 0 we get the root h ₁. Then using θ ₁ and h ₁ we get another function θ͂₂ = 𝓓 h ₁[θ ₁]. For m ≥ 1, we get the iteration scheme 𝓐 _θ1(h_m) = 0, θ _m+1 = 𝓓_hm[θ_m]. Performing the iteration, we will find $\mid h_{m+1}-h_m\mid \stackrel{m\to \infty }{\to }0$ and $\parallel {\theta }_{m+1}-{\theta }_m\parallel \stackrel{m\to \infty }{\to }0$. To some accuracy we get a numerical solution ψ(x) = 1-h_mθ_m(x) and $\varphi \left(x\right)={\int }_0^x\frac{\mu }{{\psi }^2\left(\xi \right)}\mathrm{d\xi }$.

As two nonlocal examples we consider the exponentially decaying nonlocal response[16, 17, 25] $R\left(x\right)=\frac{1}{2w}\exp \left(-\frac{\mid x\mid }{w}\right)$ and the Gaussian nonlocal response[9, 12, 25]$R\left(x\right)=\frac{1}{w\sqrt{\pi }}\exp \left(-\frac{x^2}{w^2}\right)$, where w is the characteristic nonlocal length. For the exponentially decaying nonlocal response, the perturbative refractive index $n\left(x\right)=-\int \frac{1}{2w}\exp \left(-\frac{\mid x-\xi \mid }{w}\right){\mid u\left(\xi \right)\mid }^2\mathrm{d\xi }$ is the solution of differential equation $n-w^2\frac{d^{2}n}{d{x}^2}=-{\mid u\mid }^2$. For such two nonlocal responses some gray soliton solutions are shown in Fig. 1(a) and Fig. 1(b). From the numerical simulations shown in Fig. 1(c) and Fig. 1(d), we can find these numerical soliton solutions obtained by the spectral renormalization method can describe the nonlocal gray solitons very well. Based on a reduced weakly nonlocal model $i\frac{\partial u}{\partial z}+\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}-\left[{\mid u\mid }^2+\gamma \frac{{\partial }^2{\mid u\mid }^2}{\partial x^2}\right]u=0,$, where $\gamma =\frac{1}{2}\int R\left(x\right)x^2\mathrm{dx}$, Ref. [24] claimed that for an unit background intensity, when γ > 1/4 (or equivalently when w > 1/2 for the exponentially decaying nonlocal response or when w > 1 for the Gaussian nonlocal response), there is no dark/gray soliton and the plane-wave solution become modulationally unstable. However Ref. [25] indicated that for Eq. (1) the plane-wave solutions are always stable for such two nonlocal responses with any value of w, which is greatly different from that of the reduced model. Here as shown in Fig. 1 and Fig. 4 in our paper, for an unit background intensity Eq. (1) has gray/dark soliton solutions even when w > 1/2 for the exponentially decaying nonlocal response or when w > 1 for the Gaussian nonlocal response. In fact Eq. (1) has gray/dark soliton solutions with any background intensity for any value of w.

 

Fig. 1. (a), (b) Gray soliton solutions computed with the exponentially decaying or Gaussian nonlocal response respectively when μ = 1/3 and w = 0 (black line) or w = 1 (red line), or w = 3 (green line), or w = 5 (blue line). (c), (d) Numerical simulations computed with the exponentially decaying or Gaussian nonlocal response respectively when w = 5.

Download Full Size | PDF

 

Fig. 2. (a) The gray soliton (black line) computed with an exponentially decaying nonlocal response when μ = 0.43,w = 5, and its fitting functions (red line) 1 - 0.137exp(-0.0753∣x∣)cos(0.596∣x∣ + 0.055); (b) The gray soliton (black line) computed with a Gaussian nonlocal response when μ = 0.4, w = 5, and its fitting function (red line) 1 -0.238exp(-0.0733∣x∣)cos(0.713∣x∣ + 0.107).

Download Full Size | PDF

Now we investigate the gray soliton’s asymptotic behavior under x → ±∞. Since $\psi \left(x\right)\stackrel{x\mathrm{\to \pm \infty }}{\to }$ and χ(x) = 1 - ψ(x), we have ∣χ(x)∣ ≪ 1 when x → ±∞. In such a case, to the leading order Eq. (3) can be linearized to

 

Fig. 3. Areas below red lines are the existing parameter-spaces for gray soliton in a Kerr-type nonlinear media with (a) an exponentially decaying or (b) a Gaussian nonlocal response. Areas on the left sides or right sides of blue lines correspond to gray solitons with exponentially decaying or exponentially decaying oscillatory tails respectively.

Download Full Size | PDF

Since Eq. (4) is linear for χ(x), the superposition theorem applies. For the exponentially decaying nonlocal response, Eq. (4) can be expressed as two coupled equations $-\frac{1}{2}\chi ″-2{\mu }^2\chi +2f\left(x\right)=0$ and f(x)-w ² f″(x) = χ(x). By using an ansatz χ(x) = exp(λx), from such two coupled equations we have ${\lambda }^2=4\left(\frac{1}{1-{\lambda }^2{w}^2}-{\mu }^2\right)$, which has roots λ = ±(λ ₁ + λ ₂), where ${\lambda }_1=\sqrt{\left[1-4w^2{\mu }^2+4w\sqrt{1-{\mu }^2}\right]}/\left(4w^2\right)$, and ${\lambda }_2=\sqrt{\left[1-4w^2{\mu }^2+4w\sqrt{1-{\mu }^2}\right]}/\left(4w^2\right)$. For example when w = 5, μ, = 0.43, we have λ ₁ = 0.0753, λ ₂ = 0.596 which are used by the fitting function in Fig. 2(a). We can find the numerical solutions agree well with the analytical results. When $w\le \frac{1}{4},0\le {\mu }^2\le 1$, or when $\frac{1}{4}<w\le \frac{1}{2},\frac{4w-1}{4w^2}\le {\mu }^2\le 1$, the λ is a real number and the gray soliton has an exponentially decaying tail; When $w>\frac{1}{4},0\le {\mu }^2\le \frac{4w-1}{4w^2}$, the λ is a complex number and the gray soliton has an exponentially decaying oscillatory tail. In conclusion gray soltions can exist when $w\le \frac{1}{2},0\le \mid \mu \mid \le 1$ and when $w\ge \frac{1}{2},0\le \mid \mu \mid \le \sqrt{\frac{4w-1}{4w^2}}$. Such results are shown in Fig. 3(a). The red lines in Fig. 3 show the upper boundary of μ, below which the gray soliton can exist. There is a critical value of characteristic nonlocal length $w_c=\frac{1}{2}$, for w < w_c the maximal velocity (the upper bound of ∣μ∣) μ_max = 1 is independent of the characteristic nonlocal length w; for w > w_c, the ${\mu }_{max}=\sqrt{\frac{4w-1}{4w^2}}$ decreases with w. Contrarily, based on the same nonlocal model, it was claimed in Ref. [17] and shown in Fig. 4(a) of Ref. [17] that the μ_max monotonically decreases with w. There was no critical value for the characteristic nonlocal length in Ref. [17]. A possible reason for such a disagreement between Ref. [17] and our work could be that the essentially numerical results of Ref. [17] may be obtained by using a discrete nonlocal length larger than such a critical value. So the results of Ref. [17] can not show the real relationship of μ_max to w when w is less than such a critical value w_c.

Similar results are obtained by using a Gaussian nonlocal response. Using an exponentially decaying oscillatory ansatz χ(x) = exp(-λx)cos(κx), we have $\int \frac{1}{w\sqrt{\pi }}{e}^{-{\left(x-\xi \right)}^2/w^2}\chi \left(\xi \right)\mathrm{d\xi }=\mathrm{A\chi }\left(x\right)+\mathrm{B\chi }\prime \left(x\right)$, where $A=e^{\frac{w^{2\left({\lambda }^2-{\kappa }^2\right)}}{4}}\left(\cos \frac{w^2\mathrm{\lambda \kappa }}{2}-\frac{\lambda }{\kappa }\sin \frac{w^2\mathrm{\lambda \kappa }}{2}\right)$ and $B=-\frac{1}{\kappa }{e}^{\frac{w^{2\left({\lambda }^2-{\kappa }^2\right)}}{4}}\sin \frac{w^2\mathrm{\lambda \kappa }}{2}$. Then Eq. (4) turns into $-\frac{1}{2}\chi ″+2\mathrm{B\chi }\prime +\left(2A-2{\mu }^2\right)\chi =0,$, the eigenvalue problem of which provides two coupled equations for λ and κ: λ = -2B and $\kappa =\sqrt{-2\left(2A-2{\mu }^2\right)-4B^2}$. For example when w = 5, μ = 0.4, from such two coupled equations we get λ = 0.0733, κ = 0.713 which are used by the fitting function in Fig. 2(b). After some calculation we find when $w\le \frac{1}{\sqrt{e}},0\le {\mu }^2\le 1$ or when $\frac{1}{\sqrt{e}}<w<1,\frac{1+1n\left(w^2\right)}{w^2}<{\mu }^2\le 1$ the gray solitons have exponentially decaying tails; when $w>\frac{1}{\sqrt{e}},0\le {\mu }^2\le \frac{1+1n\left(w^2\right)}{w^2}$ the gray solitons have exponentially decaying oscillatory tails. In conclusion gray soltions can exist when w ≤ 1,0 ≤ ∣μ∣ ≤ 1 and when $w\ge \mathrm{1,0}\mid \mu \mid \le \sqrt{\frac{1+1n\left(w^2\right)}{w^2}}$.There also is a critical value of characteristic nonlocal length w_c = 1, for w < w_c the maximal velocity μ_max = 1 is independent of w; for w > w_c the maximal velocity ${\mu }_{max}=\sqrt{\frac{1+1n\left(w^2\right)}{w^2}}$ decreases with w. Such results are shown in Fig. 3(b).

 

Fig. 4. (a) Dark solitons computed with a Gaussian nonlocal response when w = 0(black line), or w = 1(red line ), or w = 3(green line), or w = 5(blue line). (b) The numerical simulation of a dark soliton computed with a Gaussian nonlocal response when w = 5.

Download Full Size | PDF

Now we discuss how to compute the nonlocal dark soliton solutions. In the case λ = 0, equation (2) reduces to

For dark soliton solutions we have ψ(-x) = - ψ(x). Acting $\frac{d}{\mathrm{dx}}$ on Eq. (5) and letting $\chi \left(x\right)=\frac{\mathrm{d\psi }\left(x\right)}{\mathrm{dx}}$, we obtain -χ″/2 + N(χ) = 0, where N(χ) = -χ+ χ ∫ R(x - ξ)ψ ²(ξ)dξ + 2ψ ∫ R(x - ξ) ψ(ξ)χ(ξ)dξ. Then we can use the spectral renormalization method to numerically compute ψ(x). Some dark nonlocal soliton solutions computed with a Gaussian nonlocal response are shown in Fig. 4(a). From the numerical simulation shown in Fig. 4(b), we can find this numerical dark soliton solution can describe the dark soliton very well. Similar results are obtained by using an exponentially decaying nonlocal response.

3. Conclusion

We use the spectral renormalization method to compute dark/gray soliton solutions in Kerr-type nonlocal media and find it agrees well with the analytical results. The analytical relationship for the maximal transverse velocity of gray soliton to the characteristic nonlocal length is obtained. It is indicated there is a critical value of characteristic nonlocal length, when the characteristic nonlocal length is less than such a critical value, the gray soliton’s maximal transverse velocity does not vary with the characteristic nonlocal length, otherwise such maximal velocity decreases with the characteristic nonlocal length.