Unique Features of Nonlocally Nonlinear Systems with Oscillatory Responses

Abstract

We review the recent investigation of a new form of nonlocally nonlinear system with oscillatory responses. The system has various new features, such as the nonlocality-controllable transition of self-focusing and self-defocusing nonlinearities, a unique modulational instability and new forms of solitons. We also discuss the propagation of the optical beam in a nematic liquid crystal with negative dielectric anisotropy and demonstrate theoretically that propagation can be modelled by the system.

1. Introduction

The optical Kerr effect (OKE), a phenomenon that refers to the dependence of the refractive index on the optical intensity, is one of the most important effects in nonlinear optics [1,2]. The physical mechanism of OKE includes molecular reorientation, thermal nonlinearity, photorefractive effect, electronic contribution and electrostriction, etc. No matter what the physical mechanism is, the caused refractive index can always be phenomenologically represented as $n(\mathit{r},z,t)=n_0(\mathit{r},z)+n_{nl}(\mathit{r},z,t)$, with $n_0$ and $n_{nl}$ being, respectively, the linear refractive index and nonlinear refractive index (NRI). If NRI at a certain point in space is not determined solely by the optical intensity at that point but also depends on its vicinity, nonlinearity is spatially nonlocal [3,4,5]. Mathematically, NRI can be expressed as an integral of a kernel function (also called the response function) and an optical intensity [3]. The media of nonlocal nonlinearity are referred to as nonlocally nonlinear media [6,7,8,9], which ranges from nematic liquid crystal (NLC) with positive dielectric anisotropy [10,11,12,13,14], lead glass [15], nonlinear ion gas [16], photorefractive crystal [17], dipolar Bose–Einstein condensate [18] and so on. For the physical systems mentioned above, response functions are positive definite, for example, the exponential-decay response function for the (1 + 1)-dimensional case [19] and the zeroth-order modified Bessel function for the (1 + 2)-dimensional infinite case [10,20].

In fact, there is the other kind of response function without positive definiteness: the sine-oscillation function, brought out in the study of quadratic solitons by the formal equivalence in mathematics between quadratic and nonlocal solitons [21,22]. This kind of sine-oscillation function was also obtained in a system of coupled Gross–Pitaevskii–Poission equations [23], which govern the evolutions for matter–wave components and the microwave magnetic field in atomic Bose–Einstein condensates. Recently, we investigated the nonlocally nonlinear system with oscillatory responses [24,25,26,27,28,29,30,31,32] and found various new features such as the nonlocality-controllable transitions between focusing and defocusing nonlinearities, unique modulational instabilities (MI) and new kinds of solitons. The nonlocally nonlinear system with oscillatory responses can model the propagation of optical beams in NLC with negative dielectric anisotropy.

In this review article, we provide a brief overview of unique features of nonlocally nonlinear system with oscillatory responses, including its mathematical model, the treatment of the model and the nonlocality-controllable transition between focusing and defocusing nonlinearities in Section 2, the short-term and long-term MI properties in Section 3, different kind of soliton solutions in Section 4 and a special algorithm for finding soliton solutions in Section 5. In Section 6, we discuss the nonlinearity characteristic of NLC with negative dielectric anisotropy and show that it can be modeled by the nonlinear system with oscillatory responses.

2. Nonlocality-Controllable Kerr Nonlinearities

The model we considered is the dimensionless system in the form of [27,28]

$$i\frac{\partial u}{\partial z}+\frac{1}{2}{\nabla }_D^{2}u+\Delta nu=0,$$

(1)

$$w_m^2{\nabla }_D^2\Delta n+\Delta n-s{\left|u\right|}^2=0.$$

(2)

The equations are the coupled equations for a dimensionless optical beam $u(\mathit{r},z)$ and its (dimensionless) induced NRI $\Delta n(\mathit{r},z)$, where $\mathit{r}$ is D-dimensional transverse coordinate vectors ($D=1$ or 2), ${\nabla }_D^2$ is the D-dimensional transverse Laplacian operator, $s(=\pm 1)$ is the sign of Kerr coefficient and $w_m$ is the nonlinear characteristic length. We will show in Section 6 that the model, Equation (1) plus (2), can describe the evolutions of optical beams in the NLC with negative dielectric anisotropy. However, the model discussed extensively before is the following [3,5]:

$$w_m^2{\nabla }_D^2\Delta n-\Delta n+s{\left|u\right|}^2=0,$$

(3)

of which the second term has different signs from that of Equation (2). It is such a difference in the signs that makes the two nonlinear systems described by Equations (2) and (3) exhibit utterly different features.

Although the NRI described by Equation (2) in an infinite space can also be expressed by the following convolution:

$$\Delta n(\mathit{r},z)=s{R}_D\otimes {\left|u\right|}^2,$$

(4)

similarly to the case of Equation (3) [3,5], the D-dimensional response function $R_D$ is completely different. When $D=1$, the response takes the sin-oscillatory function as [27]

$$R_1=\frac{1}{2w_m}sin\left(\frac{\left|x\right|}{w_m}\right).$$

(5)

While, when $D=2$ [28], the response function is expressed by the following:

$$R_2=\frac{1}{4w_m^2}{Y}_0\left(\frac{r}{w_m}\right)$$

(6)

with $Y_0$ being the zeroth order Bessel function of the second kind. Substituting Equation (4) into Equation (1) yields a nonlocally nonlinear Schrödinger equation in the form of

$$i\frac{\partial u}{\partial z}+\frac{1}{2}{\nabla }_D^{2}u+su{R}_D\otimes {\left|u\right|}^2=0.$$

(7)

We can glimpse at the self-focusing and self-defocusing property of the nonlinear system described by Equation (7) from two limits where $w_m=0$ and $w_m\to \infty$. When $w_m=0$, Equation (2) is reduced to $\Delta n={s\left|u\right|}^2$. Clearly, the focusing nonlinearity occurs for $s=1$ and the defocusing nonlinearity does for $s=-1$. When $w_m\to \infty$, Equation (2) changes to $w_m^2{\nabla }_D^2\Delta n=s{\left|u\right|}^2$, which models the focusing nonlinearity in lead glass for $s=-1$ [15] and the defocusing nonlinearity for $s=1$. Therefore, the nonlinearities reverse if $w_m$ proceeds from 0 to ∞ for both cases of $s=\pm 1$. This indicates that the self-focusing and self-defocusing property of the nonlinear system described by Equation (7) depends on the degree of nonlocality, which is defined by $\sigma =w_m/w$ with w being the beam width [27]. In contrast, for the nonlinear system given by Equation (3), the self-focusing and self-defocusing property is different. In the case where $s=1$, Equation (3) is reduced to $\Delta n={\left|u\right|}^2$ when $w_m=0$, and to $w_m^2{\nabla }_D^2\Delta n=-{\left|u\right|}^2$ when $w_m\to \infty$. Under the two limits, the system (3) exhibits both self-focusing nonlinearity: that is, the local self-focusing nonlinearity when $w_m=0$ and the (thermally induced) nonlocal self-focusing nonlinearity in lead glass $w_m\to \infty$ [15]. Conversely, in the case where $s=-1$, the self-defocusing nonlinearity exists for the two limits of $w_m$. Therefore, no nonlinearities transition can take place in the nonlinear system given by Equation (3).

The dependence of the self-focusing and self-defocusing property on the degree of nonlocality can be obtained by both the variational approach and by numerical simulations. The Lagrangian density of the system described by Equation (7) is the following: [33,34]

$$\mathcal{L}=\frac{i}{2}\left(u^{\ast }\frac{\partial u}{\partial z}-u\frac{\partial u^{\ast }}{\partial z}\right)-\frac{1}{2}|{\nabla }_D{u|}^2+\frac{s}{2}{\left|u\right|}^2\left(R_D\otimes {\left|u\right|}^2\right).$$

(8)

We introduce a trial Gaussian beam:

$$u(\mathit{r},z)=A\left(z\right)exp\left[i\alpha \left(z\right)\right]exp\left[-\frac{r^2}{2w{\left(z\right)}^2}+ic\left(z\right)r^2\right],$$

(9)

where all meanings of parameters $A,\alpha$ and c can be found in Ref. [28]. According to the standard procedure of variational approach [33,34], we obtain the evolution of the beam width w [28]:

$$\frac{d^{2}w}{d{z}^2}=\frac{1}{w^3}-\frac{2s{P}_0}{{\pi }^D}{N}_D(w,w_m),$$

(10)

where $P_0={\int }_{-\infty }^{+\infty }{\left|u(\mathit{r},z)\right|}^2{d}^D\mathit{r}={\pi }^{D/2}{A}^2{w}^D$ is the input power, and the mathematical expression of $N_D$ can be found in Ref. [28]. By assuming ${w|}_{z=0}=w_0=1$ and ${dw/dz|}_{z=0}=0$, Equation (10) can be solved:

$$z={\int }_1^w\frac{d{w}_2}{\sqrt{1-1/w_2^2-(4s{P}_0/{\pi }^D){\int }_1^{w_2}{N}_D(w_1,{\sigma }_0)d{w}_1}},$$

(11)

where ${\sigma }_0=w_m/w_0{|}_{w_0=1}=w_m$. For linear propagations of the input Gaussian beam, the beam width $w_0=1$ is widened to be $w\left(z\right)=\sqrt{1+z^2}$ [35]. If the beam is in a focusing state, its width is $w\left(z\right)<\sqrt{1+z^2}$, while in a defocusing state it is $w\left(z\right)>\sqrt{1+z^2}$. Therefore, when nonlinearities transit between focusing and defocusing as the degree of nonlocality changes, by replacing w with $\sqrt{1+z^2}$ in Equation (11) transition point ${\sigma }_c$ can be obtained as

$$z={\int }_1^{\sqrt{1+z^2}}\frac{d{w}_2}{\sqrt{1-1/w_2^2-(4s{P}_0/{\pi }^D){\int }_1^{w_2}{N}_D(w_1,{\sigma }_c)d{w}_1}}.$$

(12)

Equation (12) provides an implicit function of the critical degree of nonlocality ${\sigma }_c$ on z, $P_0$, s and D. For a given $s,D$ and $P_0$, the function ${\sigma }_c\left(z\right)$ can be obtained by numerically solving the integral Equation (12) at different propagation distances z, which is shown in Figure 1.

At the specific distance $z=1$, the critical degree of nonlocality for the $(1+1)$- and $(1+2)$-dimensional cases can be obtained from Equation (12) in that ${\left.{\sigma }_c\right|}_{z=1}=0.82$ and ${\left.{\sigma }_c\right|}_{z=1}=0.63$, respectively. The dependence of beam width on the degree of nonlocality is given in Figure 2. For the case of $s=-1$, optical beams experience self-defocusing nonlinearity when ${\sigma }_0<{\sigma }_c$, and self-focusing nonlinearity when ${\sigma }_0>{\sigma }_c$. The case of $s=1$ is on the contrary. Figure 3 shows the evolutions of the optical beam for different degree of nonlocality in two cases of $s=1$ and $s=-1$, where the degree-of-nonlocality-dependence self-focusing and self-defocusing effects can be obviously observed.

3. Unique Modulational Instability

In this section, we briefly review the modulational instabilities (MI) of the (1 + 1)-dimensional nonlocally nonlinear system (7) with oscillatory response (5). The evolution of perturbations we investigated includes two processes: short-term evolution and long-term evolution. During the short-term evolution of MI, the perturbation is small enough so that the linear analysis on the NNLSE (7) is valid [31]. However, after the optical field evolves for long term, the increasing perturbation is not much lower than the amplitude of plane wave, and the linear approximation method is not applicable, and the nonlinear evolution of MI should be considered: that is, the long-term evolution of MI [32].

3.1. Short-Term Evolution of MI

According to the standard procedure [36], we add a random perturbation $\psi$ (${\left|\psi (x,z)\right|}^2\ll I_0$) on the plane wave:

$$u(x,z)=\left[I_0^{1/2}+\psi (x,z)\right]exp\left[i2\pi s{I}_0\tilde{R}\left(0\right)z\right],$$

(13)

where $\tilde{R}\left(0\right)=\tilde{R}\left(k_x\right){|}_{k_x=0}$ and $\tilde{R}\left(k_x\right)$ represents the Fourier transform of the response function (5). After calculations, the gain coefficient of perturbation is obtained as follows [31]:

$$\mathrm{g}=\mathrm{Re}\left[2|k_x|{\left(\frac{s{I}_0}{1-w_m^2{k}_x^2}-\frac{{k_x}^2}{4}\right)}^{1/2}\right],$$

(14)

which only exists when $s{I}_0/(1-w_m^2{k}_x^2)>k_x^2/4$. When $s=-1$, MI occurs when the following is the case:

$$1<w_m^2{k}_x^2<\frac{1}{2}{\left(1+16w_m^2{I}_0\right)}^{1/2}+\frac{1}{2},$$

(15)

which is shown Figure 4. When $s=1$ and the light intensity is small enough to meet the following condition:

$$0<I_0\le \frac{1}{16w_m^2},$$

(16)

there are two MI gain bands, as shown in Figure 5, at each side of the origin.

$$\left\{\begin{array}{c}\hfill 0<w_m^2{k}_x^2<\frac{1}{2}-\frac{1}{2}{\left(1-16w_m^2{I}_0\right)}^{1/2},\\ \\ \hfill 1>w_m^2{k}_x^2>\frac{1}{2}+\frac{1}{2}{\left(1-16w_m^2{I}_0\right)}^{1/2}.\end{array}\right.$$

(17)

Furthermore, after $I_0$ exceeds the critical value of $1/16w_m^2$, the MI gain bands combine into one as follows:

$$0<w_m^2{k}_x^2<1,$$

(18)

which is shown Figure 5. MI in the nonlocally nonlinear system with oscillatory responses is found to have two unique properties. First, MI exists both when the Kerr coefficient is positive and when it is negative. Second, the maximum gain points of MI do not shift with light intensity. The physical mechanism behind the properties of MI has been revealed by utilizing the theory of four-wave mixing in Ref. [31].

3.2. Long-Term Evolution of MI

To show the long-term evolution of MI, we take an infinite plane wave with perturbation [32]:

$${u(x,z)|}_{z=0}=1+{10}^{-4}cos\left(k_{x}x\right),$$

(19)

as the initial input of Equation (7) with $s=-1$ and $D=1$, where the first term, the second term and $k_x$ represent plane wave, perturbation and frequency of perturbation, respectively. In the following numerical simulation, we assume that $w_m=1$ and $I_0=1$. Then, MI appears in the range of $1<|k_x|<1.61$. $|k_x|=1$ is the singular point at which the gain coefficient g is infinite and $|k_x|=1.61$ is the cutoff point. The nonlinear evolutionary process of initial input (19) of different perturbation frequencies is shown in Figure 6. In the short-term evolution, the curve of intensity obtained through linear processing is consistent with the numerical simulation, which means that linear approximation is applicable. However, in long-term evolution, the results obtained by linear approximation gradually deviated from numerical simulations. Additionally, the perturbation of the frequency close to the singular point increases more quickly to a larger peak; by contrast, the perturbation of the frequency close to the cutoff point presents a slower growth of a lower peak. The reason why the numerical simulations deviate from the analytical result is that, during the nonlinear evolutions of perturbation higher-order, harmonic waves appear. In this case, the linear approximation of perturbation is not appropriate, and harmonic waves with higher frequencies should be taken into account.

We can obtain the stable solution of induced MI by considering the input modulated wave to be composed of harmonics with various frequencies:

$$u(x,z)=\sum _{n=-\infty }^{+\infty }{A}_n\left(z\right)exp\left(i{k}_{x}nx\right).$$

For stable solutions, the amplitudes $A_n$ of various harmonic do not vary with propagation distance z; then, we assume $A_n=a_{n}exp\left(i\beta z\right)$. Substitution into Equation (7) with $D=1$ and $s=-1$ yields the following:

$$-\beta a_n-\frac{1}{2}{a}_n{k}_x^2{n}^2+s\sum _{m=-\infty }^{+\infty }{a}_m{I}_{n-m}\left(z\right)\frac{1}{1-w_m^2{k}_x^2{(n-m)}^2}=0.$$

(20)

If only the first harmonic is taken into account, the analytical solution to Equation (20) can be obtained as

$$\begin{array}{ccc}\hfill u(x,z)& =& \left[\sqrt{\frac{\left(\frac{1}{1-w_m^2{k}_x^2}-\frac{k_x^2}{4}\right)(1-w_m^2{k}_x^2)}{2}}+\frac{\sqrt{\left(\frac{1}{1-w_m^2{k}_x^2}+\frac{k_x^2}{4}\right)(1-w_m^2{k}_x^2)}}{4}cos\left(k_{x}x\right)\right]\hfill \\ & & exp\left[-i\left(\frac{1}{1-w_m^2{k}_x^2}+1+\frac{k_x^2}{4}\right)z\right].\hfill \end{array}$$

(21)

By using $u(x,0)$ as the initial conditions, the evolution of stable solution is displayed in Figure 7. As shown in the figure, only the solutions when the modulation frequency is near the cutoff frequency are able to be stably propagated, while the ones when the modulation frequency departs from the cutoff frequency will fluctuate greatly during propagations; in this case, an increased number of harmonics should be considered. In this case, it is impossible to obtain the analytical solution, and the numerical method is required in order to solve Equation (20) [32].

4. Characteristics of Solitons

The MI of the nonlocally nonlinear system (7) always exists for cases of $s=\pm 1$ [31]. In consequence, due to the MI [36], the dark solitons [27,37] in self-defocusing sides of such a nonlocally nonlinear system are unstable [27], while the bright solitons existing in self-focusing sides can propagate stably. Concretely, in the case of $s=-1$, self-focusing and self-defocusing nonlinearities are, respectively, exhibited in high $({\sigma }_0>{\sigma }_c)$ and low $({\sigma }_0<{\sigma }_c)$ ranges of the degree of nonlocality, while in the case of $s=1$, the self-focusing and self-defocusing properties are on the contrary. In the following, we will review our works on the bright solitons in the self-focusing sides of such a system described by Equation (7) in cases of $s=\pm 1$. In addition to fundamental solitons [27,28], there are new kinds of solitons including the in-phase and out-of-phase bound-state solitons [29] and multi-peak (more than four peaks) solitons [30] in the nonlocally nonlinear system with oscillatory responses.

The variational approach is used to obtain the approximate solutions of bright solitons, including the fundamental solitons [28] and and Hermite–Gaussian-type multi-peak solitons [30]. The variational results have been confirmed by the numerical ones. We used the imaginary time evolution method to iterate the fundamental solitons [27,28] and the in-phase and out-of-phase bound-state solitons [29]. To iterate the multi-peak solitons with arbitrary peak numbers [30], we specifically developed a perturbation-iteration method [38].

4.1. Fundamental Solitons

Integrating Equation (10) once yields $\frac{1}{2}{\left(\frac{dw}{dz}\right)}^2+V_D\left(w\right)=0$, where $V_D\left(w\right)$ is the equivalent potential in the form of [28]

$$V_D\left(w\right)=-{\int }_1^w\left[\frac{1}{w_1^3}-\frac{2s{P}_0}{{\pi }^D}{N}_D(w_1,{\sigma }_0)\right]d{w}_1.$$

(22)

The soliton solutions correspond to the extremum points of the equivalent potential [28,33]; therefore, by allowing ${\partial }_w{V}_D=0$ the critical power of the solitons can be obtained as

$$P_c=\frac{{\pi }^D}{2s{N}_D(1,{\sigma }_0)}.$$

(23)

It should be noted that, in Equation (23), the critical power $P_c$ must be positive, which requires $N_D>0$ for $s=1$ and $N_D<0$ for $s=-1$. Furthermore, the solitons are stable at the minimum points of $V_D\left(w\right)$, then we have:

$${\left.\frac{d^2{V}_D}{d{w}^2}\right|}_{w=1,P_0=P_c}>0.$$

(24)

From Equation (24), the existence ranges of the degree of nonlocality ${\sigma }_0$ for the solitons can be determined, which depend on the values of D and the sign of s. When $D=1$, the solitons exist if ${\sigma }_0\in [0,0.76]$ for $s=1$ and $[1.05,+\infty )$ for $s=-1$. When $D=2$, the existence ranges of ${\sigma }_0$ are $[0.33,0.61]$ for $s=1$ and $[0.90,+\infty )$ for $s=-1$. These variational results agree well with the numerical ones, which are given in the following.

In order to numerically iterate the soliton solution, we assume $u(x,z)=u\left(x\right)exp\left(i\lambda z\right)$. Substitution of the solution into Equations (1) and (2) yields the following:

$$\begin{array}{ccc}\hfill \frac{1}{2}\frac{d^{2}u}{d{x}^2}-\lambda u+\Delta nu& =& 0,\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill w_m^2\frac{d^2\Delta n}{d{x}^2}+\Delta n-s{u}^2& =& 0.\hfill \end{array}$$

(26)

By the imaginary-time method, the solitons can be numerically found in higher ranges of the degree of nonlocality, ${\sigma }_0\in [1.05,+\infty ]$ for $D=1$ and ${\sigma }_0\in [1.11,+\infty ]$ for $D=2$ in the case of $s=-1$. These ranges of ${\sigma }_0$ are all within the self-focusing sides of such a nonlinear system and agree well with the variational results. The bright solitons have two abnormal properties. The first one is the negative propagation constants, while it is the opposite of the cases in local nonlinear media [39] and nonlocally nonlinear media with the positively defining attenuating response functions [3]. The other is the negative slope of the $P_c\left(\lambda \right)$ ($P_c$ and $\lambda$ being the critical power and the propagation constant), as shown in Figure 8. By means of linear stability analysis, all solitons are shown to be stable [27]. This means that the stability criterion obeys an anti-Vakhitov–Kolokolov stability criterion [40], which is also obeyed by the in-phase and out-of-phase bound-state solitons and the multi-peak solitons reviewed in the next two sections. However, in nonlocally nonlinear media with the positively defining attenuating response functions, stable solitons should comply with the Vakhitov–Kolokolov criterion [9]. We have noticed that the evolutions for matter–wave components and the microwave magnetic field in the atomic Bose–Einstein condensates can also be governed by the mathematical model of Equation (1) plus (2) with $s=-1$ in some specific conditions [23,41,42], and in such a system, bright solitons are obtained in strongly nonlocal cases.

Although what Equations (25) and (26) determine are the solitons of the nonlocally nonlinear system (7), the coupled equations with $s=1$ are formally the same as the model governing parametric spatial solitary waves in quadratic nonlinear materials [21,22,43]. The fundamental and the second harmonic waves are mathematically equivalent to the paraxial beam and its induced NRI [21]. Esbensen et al. employed the formal equivalence between two distinct physical systems to investigate the stabilities of solitons and found that the solitons are all unstable. However, it should be noted that equivalence is limited only to the solitary waves, while their dynamic properties, such as the stability, are completely different. In quadratic nonlinear materials, there is the energy and phase exchange between the fundamental and the second harmonic waves. Therefore, for the second harmonic wave, its z-dependant evolution should be taken into account, which is not the case for the NRI shown in Equation (2). Detailed discussions on this point can be found in Ref. [10]. Our group also made some improvements on the issue of quadratic solitons [24,25,26]. We found that the boundary confinement of media can support stable solitons [24]. Furthermore, the fundamental waves even can be multipole solitons [25], the existence of which depends closely on the sample size and the degree of nonlocality. If nonlocality is fixed and the sample size is varied, soliton width varies piecewise and approximately periodically [26].

4.2. In-Phase and Out-of-Phase Bound-State Solitons

In the case of $s=1$, the nonlinear system exhibits self-focusing nonlinearity in lower ranges of degree of nonlocality, that is, ${\sigma }_0<{\sigma }_c$. In these ranges, we obtained bright solitons with complicated structures, which are the in-phase bound-state solitons when ${\sigma }_0\in \left[0.28,0.78\right]$ and the out-of-phase bound-state solitons when ${\sigma }_0\in \left[0.38,0.78\right]$ [29]. The in-phase bound-state soliton, shown in Figure 9, owns the symmetrical profile and the nonzero central value. On the other hand, the out-of-phase bound-state soliton, shown in Figure 10, owns the antisymmetrical profile and the zero central value. When the degree of nonlocality decreases, the in-phase bound-state soliton approaches the sech profile shown in Figure 9, and the out-of-phase bound-state soliton tends toward the first-order Hermite–Gaussian profile shown in Figure 10. The above-mentioned abnormal properties in the previous section, that is, the negative propagation constants and the negative slope of the power-propagation constant (the so-called anti-Vakhitov–Kolokolov stability criterion), also exist for bound-state solitons, which is shown in Figure 11a. Furthermore, it can be found from Figure 11b that the power-propagation constant diagrams are the same for both the in-phase and the out-of-phase bound-state solitons. In other words, we can say that two forms of bound-state solitons form degenerate modes.

4.3. Multi-Peak Solitons

The Hermite–Gaussian soliton solutions with multipeaks in fact have been discussed in nonlocally nonlinear media with positively-definiting attenuating response functions [44,45]. However, the upper thresholds of the peak number of the stable solitons in nonlocally nonlinear media are different, which closely depend on the response function of media. For instance, in the nonlocal system with the Gaussian response, the multipeak solitons with any peak number [4] are stable [44]. In nematic liquid crystals with positive dielectric anisotropy, the response function is an exponential decay one for the $(1+1)$ case, and the peak-number of the solitons supported is only less than five [44]. In the nonlocally nonlinear system with oscillatory responses described by Equations (1) and (2), we found the upper thresholds of the peak number of the stable solitons are five and four in the cases of $s=-1$ and 1, respectively [30].

The variational approach can be applied to find solitons and discuss their stabilities. The Hermite–Gaussian solitons of the $(1+1)$-dimensional NNLSE (7) are of the following profiles:

$$u_n={\left(\frac{P_{nc}\sqrt{2n+1}}{2^{n}n!\sqrt{\pi }}\right)}^{1/2}{H}_n\left(\sqrt{2n+1}x\right)exp\left[-\frac{(2n+1)x^2}{2}\right].$$

(27)

The critical power is obtained as follows:

$$P_{nc}=\frac{(2n+1){(\sqrt{\pi }{2}^{n}n!)}^2}{s{N}^{\left(n\right)}(1,{\sigma }_0)},$$

(28)

where function $N^{\left(n\right)}$ is the $N_n$ given in Ref. [30]. The degree of nonlocality range within which HG solitons exist is obtained by the variational approach, which is summarized in

Table 1. The ranges of ${\sigma }_0$ within which the HG-type solitons exist.

		$\mathit{n}=0$	$\mathit{n}=1$	$\mathit{n}=2$	$\mathit{n}\ge 3$
$s=-1$	variational	$(1.05,+\infty )$	$(1.06,+\infty )$	$(1.06,+\infty )$	$(1.06,+\infty )$
	numerical	$(1.05,+\infty )$	$(1.05,+\infty )$	$(1.05,+\infty )$	$(1.05,+\infty )$
$s=+1$	variational	$(0,0.77)$	$(0.38,0.79)$	$(0.39,0.79)$	$(0.39,0.79)$
	numerical	$(0.05,0.78)$	$(0.38,0.78)$	$(0.39,0.78)$	$(0.40,0.78)$

Note: self-focusing nonlinearity is exhibited when σ₀ ∈ (σ_c, +∞) for s = −1, and [0, σ_c] for s = 1 (for example, ${{\sigma }_c|}_{z=1}$ = 0.82.)

, where the numerical results are also given for a comparison.

Multipeak solitons can be numerically obtained for two cases of $s=\pm 1$, which are shown in Figure 12 and Figure 13, respectively. Clearly, the multipeak solitons in the case of $s=-1$ exhibit Hermite–Gaussian structures and are in good agreement with numerical ones. While, the profiles of solitons in the case of $s=1$ deviate from the Hermite–Gaussian form, especially when ${\sigma }_0$ becomes large. Likely, for multipeak solitons, the propagation constants and the slope of the power-propagation constant are both negative, which is shown in Figure 14. By linear stability analyses, the ranges of the degree of nonlocality ${\sigma }_0$ within which the stable multi-peak solitons exist are summarized in

Table 2. The ranges of ${\sigma }_0$ within which the HG–type solitons are stable.

	$\mathit{n}=0$	$\mathit{n}=1$	$\mathit{n}=2$	$\mathit{n}=3$	$\mathit{n}=4$	$\mathit{n}\ge 5$
$s=-1$	$(1.05,+\infty )$	$(1.41,+\infty )$	$(1.10,+\infty )$	$(1.10,+\infty )$	$(1.10,1.61)$	no
$s=1$	$(0.05,0.78)$	$(0.38,0.78)$	$(0.42,0.45)$	$(0.48,0.49)$	no	no

.

5. Perturbation-Iteration Method

In the nonlocally nonlinear system with oscillatory responses governed by Equation (7), the Hermite–Gauss-like solitons with arbitrary peak numbers are hardly obtained by the usual numerical algorithms such as the Newton-conjugate-gradient method. Therefore, we developed a special numerical method by which we called the perturbation-iteration method to solve this problem [38]. The key idea behind the perturbation-iteration method is to treat the NNLSE as a perturbed model of the harmonic oscillator, and soliton solutions can be obtained by perturbation theory in quantum mechanics [46]. The procedure for the perturbation-iteration method is briefly reviewed. Substituting the following soliton solutions:

$$u_n(x,z)=\sqrt{A_n}{\psi }_n\left(x\right)exp(-i{\beta }_{n}z)(n=0,1,2,3\dots )$$

(29)

into the NNLSE (7) yields the following:

$$\left[-\frac{1}{2}\frac{d^2}{d{x}^2}-A_n{\int }_{-\infty }^{\infty }R(x-\xi ){\psi }_n^2\left(\xi \right)d\xi \right]{\psi }_n\left(x\right)={\beta }_n{\psi }_n\left(x\right),$$

(30)

where $A_n$ is a coefficient related to the power of soliton, ${\psi }_n\left(x\right)$ and ${\beta }_n$ are, respectively, the transverse distribution and propagation constant of soliton and n is the order of soliton. We can treat Equation (30) as a perturbed model of the harmonic oscillator:

$$\left[-\frac{1}{2}\frac{d^2}{d{x}^2}+\frac{1}{2{\mu }^4}{x}^2+f_n\left(x\right)\right]{\psi }_n\left(x\right)=E_n{\psi }_n\left(x\right),$$

(31)

where ${\psi }_n\left(x\right)$ and $E_n$ are, respectively, eigenfunctions and eigenvalues, $f_n\left(x\right)$ is the perturbation compared to the potential of the harmonic oscillator. Starting from the perturbed model of the harmonic oscillator, we determine the “minimum” perturbation, then use the formal expression of infinite-order perturbation expansions to numerically calculate the eigenfunctions and eigenvalues of the perturbed model and iterate this perturbation to obtain multipeak solitons with enough high accuracy.

This form of perturbation-iteration method developed for the $(1+1)$-dimensional NNLSE also has been extended to the $(1+2)$-case [47]. In fact, the method we developed might also be extended to the numerical integration of the Schrödinger equations in any type of potentials.

6. Optical Beams in NLC with Negative Dielectric Anisotropy

6.1. Evolution Equation for Optical Beams in NLC

The physical mechanism of nonlinearity in NLC is the optically induced molecular reorientation and nonlocality comes from the interactions between NLC molecules [3]. We consider the propagation of optical beams in the sample cell of NLC shown in Figure 15. NLC molecules exhibit anisotropy in both the low-frequency domain and the optical frequency domain, which is expressed by ${ϵ}_a^{rf}(={ϵ}_{\Vert }-{ϵ}_{\perp })$ and ${ϵ}_a^{op}(=n_{\Vert }^2-n_{\perp }^2)$, respectively. Anisotropy ${ϵ}_a^{rf}$ can be either positive [10,48] or negative [48,49], while ${ϵ}_a^{op}$ is always positive.

When an external low-frequency electrical field ${\overrightarrow{E}}_{rf}=E_{rf}{\overrightarrow{e}}_X$ and an optical field ${\overrightarrow{E}}_{op}=E_{op}{\overrightarrow{e}}_X$ are present, the molecular director $\overrightarrow{n}$ expressed by $\overrightarrow{n}\left(\theta \right)=(sin\theta ,0,cos\theta )$ ($\theta$ is the tilt angle between $\overrightarrow{n}$ and z) will be reorientated. NLC molecules approach the equilibrium state when the three torques are balanced, i.e., ${\overrightarrow{\Gamma }}_{op}+{\overrightarrow{\Gamma }}_{lf}+{\overrightarrow{\Gamma }}_{el}=0$, where ${\overrightarrow{\Gamma }}_{op}$, ${\overrightarrow{\Gamma }}_{lf}$ and ${\overrightarrow{\Gamma }}_{el}$ are the optical field-induced torque, low-frequency electrical field induced torque and the elastic torque, respectively. Using the one-constant approximation [48], we can have the following:

$${\overrightarrow{\Gamma }}_{el}=K{\nabla }^2\theta {\overrightarrow{e}}_Y,$$

where K is the NLC average elastic constant. We can also obtain the following [48]:

$$\begin{array}{ccc}\hfill {\overrightarrow{\Gamma }}_{lf}& =& {\overrightarrow{D}}_{rf}\times {\overrightarrow{E}}_{rf}={ϵ}_0{ϵ}_a^{rf}sin\theta cos\theta E_{rf}^2{\overrightarrow{e}}_Y,\hfill \\ \hfill {\overrightarrow{\Gamma }}_{op}& =& <{\overrightarrow{D}}_{op}\times {\overrightarrow{E}}_{op}>=\frac{1}{2}{ϵ}_0{ϵ}_a^{op}sin\theta cos\theta {\left|E_{op}\right|}^2{\overrightarrow{e}}_Y,\hfill \end{array}$$

where ${\overrightarrow{D}}_{rf}={ϵ}_0{ϵ}_{\perp }{\overrightarrow{E}}_{rf}+{ϵ}_0{ϵ}_a^{rf}(\overrightarrow{n}·{\overrightarrow{E}}_{rf})\overrightarrow{n}$ and ${\overrightarrow{D}}_{op}={ϵ}_0{n}_{\perp }^2{\overrightarrow{E}}_{op}+{ϵ}_0{ϵ}_a^{op}(\overrightarrow{n}·{\overrightarrow{E}}_{op})\overrightarrow{n}$ are the electric displacements, and $<>$ denotes the time average.

In the absence of light, from the balance of torques of ${\overrightarrow{\Gamma }}_{el}$ and ${\overrightarrow{\Gamma }}_{lf}$, the pretilt angle $\widehat{\theta }$ produced by the low-frequency electric field ${\overrightarrow{E}}_{rf}$ can be determined, which is the only function of X [48]

$$2K\frac{d^2\widehat{\theta }}{d{X}^2}+{\epsilon }_0{ϵ}_a^{rf}{E}_{rf}^{2}sin\left(2\widehat{\theta }\right)=0.$$

(32)

If ${ϵ}_a^{rf}>0$ (the NLC with positive dielectric anisotropy), $d^2\widehat{\theta }/d{X}^2<0$, which indicates that the angle at the center ($X=0$) of the cell ${\theta }_0(=\widehat{\theta }{|}_{X=0})$ is at the maximum and decreases as it proceeds closer to the boundaries. The opposite happens for the case that ${ϵ}_a^{rf}<0$ (the NLC with negative dielectric anisotropy), the angle at the center of the cell ${\theta }_0$ is at the minimum and increases as it approaches the boundaries. As a result, NLC molecules should be anchored at the boundaries in a manner where ${\theta |}_{X=L/2}{=\theta |}_{X=-L/2}=0$ for the positive dielectric anisotropy, as was observed in Refs. [3,10,50], and where ${\theta |}_{X=L/2}{=\theta |}_{X=-L/2}=\pi /2$, as shown in Figure 15, for negative dielectric anisotropy. In addition, the dependence of ${\theta }_0$ on low-frequency voltage $V_{rf}(\approx E_{rf}L)$ is different for the NLC with positive and the negative dielectric anisotropy, when $V_{rf}$ is higher than the Fréedericks threshold $V_{fr}=\pi \sqrt{K/{\epsilon }_0\left|{ϵ}_a^{rf}\right|}.$ For the NLC with positive dielectric anisotropy [51], ${\theta }_0\approx (\pi /2)\left[1-{\left(V_{fr}/V_{rf}\right)}^3\right]$, and for the NLC with negative dielectric anisotropy, we have the following

$${\theta }_0\approx (\pi /2){\left(V_{fr}/V_{rf}\right)}^3.$$

(33)

Although the molecule anchoring at the boundaries is different, the physical processes for both the positive and the negative are the same. Following the procedure to deal with the propagation of optical beams in an NLC with positive dielectric anisotropy [3,10,50], we obtained the equations for that in the NLC with negative dielectric anisotropy and found that the two cases can be uniformly expressed. In the presence of an externally applied (low-frequency) electric field $E_{rf}$, the propagation of the slowly varying envelope $\Phi$ of the optical field linearly polarized along the X-direction (the extraordinary light) and propagating along the Z-direction in the NLC-cell can be described by the following system:

$$2ik\frac{\partial \Phi }{\partial Z}+{\nabla }_{XY}^2\Phi +k_0^2{ϵ}_a^{op}\left[{sin}^2\theta -{sin}^2{\theta }_0\right]\Phi =0,$$

(34)

$$2K\left(\frac{{\partial }^2\theta }{\partial Z^2}+{\nabla }_{XY}^2\theta \right)+{ϵ}_0\left({ϵ}_a^{rf}{E}_{rf}^2+{ϵ}_a^{op}\frac{{|\Phi |}^2}{2}\right)sin\left(2\theta \right)=0,$$

(35)

where $k=k_0{n}_0$ with $k_0$ being the vacuum wavenumber and $n_0$ being the (linear) refractive index of the extraordinary light. Although it looks the same in form for both cases, Equation (35) is in fact different for either the positive or the negative cases because the sign of ${ϵ}_a^{rf}$ is different. The term ${\partial }_Z^2\theta$ in Equation (35) was proven to be negligible compared to ${\nabla }_{XY}^2\theta$ [50]. Furthermore, we can write $\theta =\widehat{\theta }+(\widehat{\theta }/{\theta }_0)\Psi$, where $\Psi$ represents optically induced angle perturbation. By simple substitution into Equation (35), we have the following:

$$\begin{array}{ccc}& & 2K\frac{\widehat{\theta }}{{\theta }_0}{\nabla }_{\perp }^2\Psi +\frac{4K}{{\theta }_0}\frac{\partial \widehat{\theta }}{\partial x}\frac{\partial \Psi }{\partial x}-{\epsilon }_0{ϵ}_a^{rf}{E}_{rf}^{2}sin\left(2\widehat{\theta }\right)\left[1-cos(2\frac{\widehat{\theta }}{{\theta }_0}\Psi )+\frac{\Psi }{{\theta }_0}-cot\left(2\widehat{\theta }\right)sin(2\frac{\widehat{\theta }}{{\theta }_0}\Psi )\right]\hfill \\ & & +{\epsilon }_0{ϵ}_a^{op}\frac{{|\Phi |}^2}{2}\left[sin\left(2\widehat{\theta }\right)cos(2\frac{\widehat{\theta }}{{\theta }_0}\Psi )+cos\left(2\widehat{\theta }\right)sin(2\frac{\widehat{\theta }}{{\theta }_0}\Psi )\right]=0,\hfill \end{array}$$

(36)

which is the exact result after direct substitution without any approximation. Supposing that the beam width is far smaller than cell thickness L and noting that $\widehat{\theta }\approx {\theta }_0$ and ${\partial }_x\widehat{\theta }\approx 0$ in the middle of the cell, we can simplify Equation (36) into the following form:

$$\begin{array}{c}\hfill {\nabla }_{XY}^2\Psi -\frac{{\epsilon }_0{ϵ}_a^{rf}{E}_{rf}^{2}sin\left(2{\theta }_0\right)[1-2{\theta }_{0}cot\left(2{\theta }_0\right)]}{2K{\theta }_0}\Psi +{\epsilon }_0{ϵ}_a^{op}\frac{sin\left(2{\theta }_0\right)}{4K}{|\Phi |}^2=0,\end{array}$$

(37)

where we also made sucn an assumption of ${ϵ}_a^{op}{|\Phi |}^2/2\ll {ϵ}_a^{rf}{E}_{rf}^2$ [3]. Multiplied by ${\theta }_0{ϵ}_a^{op}/E_{rf}^2{ϵ}_0{ϵ}_a^{rf}{n}_0[1-2{\theta }_{0}cot\left(2{\theta }_0\right)]$, the equation can be re-expressed as follows:

$$\mathrm{sgn}\left({ϵ}_a^{rf}\right)W_{mL}^2{\nabla }_{XY}^2{n}_{nl}-n_{nl}+n_2{|\Phi |}^2=0,$$

(38)

where

$$\begin{array}{ccc}\hfill n_{nl}& =& \frac{{ϵ}_a^{op}sin\left(2{\theta }_0\right)\Psi }{2n_0},\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill n_2& =& \frac{{\left({ϵ}_a^{op}\right)}^2{\theta }_{0}sin\left(2{\theta }_0\right)}{4n_0{ϵ}_a^{rf}{E}_{rf}^2[1-2{\theta }_{0}cot\left(2{\theta }_0\right)]},\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill W_{mL}& =& \frac{1}{E_{rf}}{\left(\frac{2{\theta }_{0}K}{{ϵ}_0\left|{ϵ}_a^{rf}\right|sin\left(2{\theta }_0\right)[1-2{\theta }_{0}cot\left(2{\theta }_0\right)]}\right)}^{1/2},\hfill \end{array}$$

(41)

are the NRI, the Kerr coefficient and the nonlinear characteristic length for the NLC. When $W_{mL}\to 0$, Equation (38) above will be reduced to $n_{nl}=n_2{|\Phi |}^2$ for the local limit [39]. Therefore, our definition of $n_2$ in Equation (40) can guarantee that its expression form is consistent when the degree of nonlocality transitions from the nonlocal case to the local case. On the other hand, the propagation equation for optical beams, Equation (34), is simplified as

$$i\frac{\partial \Phi }{\partial Z}+\frac{1}{2k}{\nabla }_{XY}^2\Phi +k_0{n}_{nl}\Phi =0.$$

(42)

When ${ϵ}_a^{rf}>0$, the Kerr coefficient $n_2$ given by Equation (40) is positive, and Equation (38) becomes $W_{mL}^2{\nabla }_{XY}^2{n}_{nl}-n_{nl}+n_2{|\Phi |}^2=0$. This equation is that of the molucular reorientation for NLC with positive dielectric anisotropy and has been discussed extensively [3,10,19,20,50,51,52]. For the NLC with negative dielectric anisotropy, however, ${ϵ}_a^{rf}<0$, the Kerr coefficient $n_2$ is negative, and Equation (38) becomes $W_{mL}^2{\nabla }_{XY}^2{n}_{nl}+n_{nl}-n_2{|\Phi |}^2=0$, which has been rarely investigated so far. In this case (${ϵ}_a^{rf}<0$), by the introduction of the dimensionless transform $u=\Phi /{\Phi }_0$, $\Delta n=n_{nl}/N_{L0}$, $x\left(y\right)=X\left(Y\right)/W_0$, $z=Z/k{W}_0^2$, and $w_m=W_{mL}/W_0$, where ${\Phi }_0=\sqrt{n_0/\left|n_2\right|}/k{W}_0$, $N_{L0}=1/k_0^2{W}_0^2{n}_0$, and $W_0$ is the beam width in Lab coordinate system, Equations (42) and (38) can be transformed as the dimensionless form $i{\partial }_{z}u+(1/2){\nabla }_{xy}^{2}u+\Delta nu=0$ and $w_m^2{\nabla }_{xy}^2\Delta n+\Delta n+{\left|u\right|}^2=0$, which does be Equations (1) and (2) with $s=-1$.

6.2. Optical Nonlinearities of NLC with Negative Dielectric Anisotropy

By substitution of Equation (33) and $E_{rf}\approx V_{rf}/L$, the nonlinear characteristic length $W_{mL}$ given by Equation (41) for the NLC with negative dielectric anisotropy is reduced to the following:

$$W_{mL}=\frac{L}{{\pi }^{4/3}}\sqrt{\frac{{\gamma }^{5/3}}{sin\gamma -\gamma cos\gamma }},$$

(43)

where $\gamma ={\pi }^4{K}^{3/2}/{\left[{\epsilon }_0(-{ϵ}_a^{rf})\right]}^{3/2}{V}_{rf}^3$. $W_{mL}$ is, obviously, the function of the cell thickness L, the bias voltage $V_{rf}$ and the material parameters of the NLC. By setting $d{W}_{mL}\left(V_{rf}\right)/d{V}_{rf}=0$, we obtain ${\gamma }_{min}=2.35$ at which $W_{mL}/L$ takes the minimum $0.29$, which does not depend on the material parameters of the NLC and is confirmed in Figure 16. Since model Equations (42) and (38) are obtained with the assumption that the beam width in the lab coordinate system, $W_0$, is far smaller than NLC thickness L, we can assume that $W_0=0.1L$ so that the weakest degree of nonlocality ${\sigma }_{min}[={\left(W_{mL}\right)}_{min}/W_0]$ is $2.9$, and the degree of nonlocality ${\sigma }_0(=W_{mL}/W_0)$ is always higher than ${\sigma }_c$ (for an example, ${\sigma }_c{|}_{z=1}=0.63$ for the (1 + 2)-dimensional case). Therefore, the NLC with negative dielectric anisotropy can only behave with self-focusing nonlinearity, and the transition from the self-focusing to the self-defocusing when ${\sigma }_0$ decreases across the critical point ${\sigma }_c$, which was predicted by the system of Equations (1) and (2) with $s=-1$, cannot be observed in the system of Equations (38) and (42) for the negative dielectric anisotropy, because the condition ${\sigma }_0<{\sigma }_c$ cannot be realized in such a NLC. We also numerically investigated the propagation of the optical beam in the NLC with negative dielectric anisotropy described by Equations (34) and (35) without simplification and found that the beam can always sample the self-focuing nonlinearity indeed [53]. Consequently, the nonlinear refractive index due to reorientation in the pure NLC is always self-focusing, despite the sign of the anisotropy in the low-frequency domain.

It is a well-known fact [54,55,56] that only dark solitons can exist in nonlocally nonlinear media with the positively defining attenuating response and the negative Kerr coefficient ($n_2<0$) described by Equations (1) and (3) for $s=-1$, no matter how much the degree of nonlocality is. However, as discussed above, the nonlocally nonlinear system of Equations (1) and (2) with $s=-1$, that is, Equations (42) and (38) for the NLC with negative dielectric anisotropy, can exhibit self-focusing nonlinearity in high ranges of the degree of nonlocality $({\sigma }_0>{\sigma }_c)$ and support the bright solitons. It is the unique feature of a nonlocally nonlinear system with oscillatory responses. By the relation of the input optical power between the Lab coordinate system $P_0^{\left(\mathrm{Lab}\right)}$ and the dimensionless system $P_0$

$$P_0^{\left(\mathrm{Lab}\right)}=\frac{k}{\omega {\mu }_0}{\int \int |\Phi |}^{2}dXdY=n_0\sqrt{\frac{{ϵ}_0}{{\mu }_0}}{\Phi }_0^2{W}_0^2\int \int {\left|u\right|}^{2}dxdy=\sqrt{\frac{{ϵ}_0}{{\mu }_0}}\frac{1}{(-n_2)k_0^2}{P}_0,$$

we can have the theoretical value of the critical power $P_{tc}$ via the variational result, i.e., Equation (23) for $s=-1$:

$$\begin{array}{c}\hfill P_{tc}=\sqrt{\frac{{ϵ}_0}{{\mu }_0}}\frac{1}{(-n_2)k_0^2}{P}_c{|}_{s=-1}=\sqrt{\frac{{\epsilon }_0}{{\mu }_0}}\frac{4n_0{ϵ}_a^{rf}{V}_{rf}^2\left[1-2{\theta }_{0}cot\left(2{\theta }_0\right)\right]}{{(2\pi /\lambda )}^2{L}^2{\left({ϵ}_a^{op}\right)}^2{\theta }_{0}sin\left(2{\theta }_0\right)}\frac{{\pi }^2}{2N_2\left(1,{\sigma }_0\right)},\end{array}$$

(44)

where $n_0=n_{\Vert }{n}_{\perp }/\sqrt{n_{\Vert }^2{cos}^2{\theta }_0+n_{\perp }^2{sin}^2{\theta }_0}$, and $N_2(1,{\sigma }_0)=\frac{\pi }{16{\sigma }_0^2}\left[\frac{1}{{\sigma }_0^2}exp(-\frac{1}{2{\sigma }_0^2})\mathrm{Ei}\left(\frac{1}{2{\sigma }_0^2}\right)-2\right]$ with Ei being the exponential integral function [28]. By using the parameters for the NLC KY19-008 [57], $K=1\times {10}^{-11}$ N, ${ϵ}_a^{op}=0.74106$, ${ϵ}_a^{rf}=-5.3$, $L=80\mathsf{\mu }$m, and beam parameters including $\lambda =532$ nm and $W_0=4\mathsf{\mu }$m, we obtain $P_{tc}\approx 6.37$ mW when voltage bias $V_{rf}\approx 3.4$ V is applied.

7. Conclusions

We reviewed recent research on the nonlocally nonlinear system with oscillatory responses. This nonlinear system exhibits various new features, such as the nonlocality-controllable transitions of nonlinearities, unique modulational instabilities and new kinds of solitons. The unique features of nonlocally nonlinear system come from oscillatory responses without positive definiteness, which are quite different from the positively defining attenuating ones discussed so far. In particular, in nonlocally nonlinear media with oscillatory responses, we theoretically predicted that bright solitons can exist even when the Kerr coefficient is negative, where only dark solitons were supposed to exist if the response functions are positively defining and attenuating. We find that oscillatory responses can describe the interaction between the optical beam and the nematic liquid crystal with negative dielectric anisotropy, although self-defocusing nonlinearity cannot be exhibited in such a nematic liquid crystal. These new and unexpected behaviors found in the nonlocally nonlinear system with oscillatory responses are of significance at a fundamental level, especially for the relevance between the nonlocality and the focusing-defocusing nonlinearities, which appears to be promising for tailoring optical properties in materials.