## Abstract

Dynamical and steady-state behavior of beams propagating in nematic liquid crystals (NLCs) is analyzed. A well-known model for the beam propagation and the director reorientation angle in a NLC cell is treated numerically in space and time. The formation of steady-state soliton breathers in a threshold region of beam intensities is displayed. Below the region the beams diffract, above the region spatiotemporal instabilities develop, as the input intensity and the material parameters are varied. Curiously, the only kind of solitons we could demonstrate in our numerical studies was the breathers. Despite repeated efforts, we could not find the solitons with a steady profile propagating in the NLC model at hand.

© 2009 Optical Society of America

## 1. Introduction

By now, nematic liquid crystals (NLCs) are common materials found in many consumer electronic devices. They behave in a fluid-like fashion, but display long-range order that is characteristic of crystals [1]. They contain rod-like molecules which exhibit orientational alignment without positional order. NLCs exhibit strong optical nonlinearities, owing to large refractive index anisotropy. Another useful property of NLCs is the ability to change optical characteristics under the action of an external electric field, which enables the macroscopic reorientation of the director tilt angle *θ*. The light incident on a NLC modifies the electric permittivity tensor, leading to the reorientational nonlinearity. For these reasons the propagation of self-focused beams [2] in NLCs have been subjected to increased scrutiny in recent years, in both experimental [3, 4] and theoretical [5–8] aspects of research.

We investigate the propagation of single self-trapped laser beams in a bulk NLC cell. Starting from the standard equations describing nonlocal and nonlinear (NN) interactions in NLCs, we develop a numerical model that exhibits rich dynamical behavior in three spatial dimensions and time. We treat the full system of partial differential equations (PDEs) in both space and time, using an appropriately modified split-step FFT procedure. We vary different beam parameters: FWHM of the input beam and the optical and static permittivity anisotropies of liquid-crystal molecules. We present breathing soliton transverse profiles at the cell exit (x,y) plane, the optical field intensity I(0,y,z) and the optically induced molecular reorientation *θ*̂ (0,y,z) along the propagation axis, in the middle of the crystal. We demonstrate the formation of soliton breathers in a threshold region of beam intensities. Despite repeated efforts and careful numerics, we could not observe steady spatial solitons, propagating with an unchanging transverse profile. The only kind of solitons we could identify in the model at hand was the breathers. We observe self-focusing and modulational instabilities (MI). The effects of self-focusing and breathing we describe resemble the physics of self-written waveguides, first noted in photosensitive glass [9–11] and later developed for polymers [12–14]. Even though rich dynamical behavior is observed, we confine our attention mostly to spatial effects, *i.e.* the effects observed in 3 spatial dimensions, after temporal steady state is reached in the system. To ascertain the reliability of our numerical procedure and the reality of effects we report upon, in the end a comparison is presented with the steady-state numerical method, based on a successive over-relaxation (SOR) algorithm.

The paper is organized as follows. Section 2 contains the description of the NN model adopted and of the numerical method employed. Section 3 presents results of our numerical simulations, grouped in a number of subsections. Section 4 brings a summary of the paper.

## 2. The model

The distortion of molecular orientation in NLC can be described by the reorientation angle *θ* of the director in the transverse plane. In the presence of an external low frequency electric field the spatial evolution of a slowly-varying beam envelope *A*, linearly polarized along the *x* axis and propagating along the *z* axis, is well captured by the dimensionless paraxial wave equation [15]:

where Δ_{x,y} is the transverse Laplacian. The temporal evolution of the angle of reorientation *θ* is described by the dimensionless diffusion equation [4, 16, 17]:

with the boundary conditions:

Here *γ* is the viscous coefficient and *K* is Frank’s elastic constant, *τ* is the director relaxation time, and

Also, *k*=*k _{0}n_{0}* is the wave vector in the medium,

*x*the transverse scaling length, and Δ

_{0}*ε*=

^{OPT}*n*is the optical permittivity anisotropy of the liquid-crystal molecules.

_{e}^{2}-n_{0}^{2}*E*=

^{DC}*V*/

*D*is the applied field strength (

*V*is the applied bias voltage,

*D*is the cell thickness) and Δε

^{DC}is the static permittivity anisotropy of the liquid-crystal molecules.

*θ*is the overall tilt angle (the total orientation of the molecules with respect to the

*z*axis), owing to both light and voltage:

*θ*=

*θ*

_{0}+

*θ*̂, where the angle

*θ*accounts for the molecular orientation induced by the static electric field only, while the quantity

_{0}*θ*̂ corresponds to the optically induced molecular reorientation. Hard boundary conditions on the molecular orientation at the transverse borders of the liquid-crystal cell are introduced by relation (3).

Equations (1) and (2) form the basis of our model. It is a closed NN model that describes the propagation of self-trapped beams in NLCs well. By solving these equations we will be describing the beam propagation in both space and time. We develop a novel numerical procedure, based on FFT, utilizing our prior experience in treating the beam propagation in NLCs [5]. The novelty is that we treat concurrently the system of coupled PDEs in both space and time. Both the spatial and the temporal integrations are based on the split-step FFT method, because both Eqs. (1) and (2) are of the parabolic type. The only difference is that the spatial equation is a paraxial wave equation, whereas the temporal equation is a diffusion equation. Our procedure is uniquely suited for observing the dynamical states in slowly-varying photosensitive systems, where the fast optical fields are slaved to the slow change in the photorefractive nonlinearity.

In all simulations the following data are kept fixed: the diffraction length *Ld*=*k*x_{0}
^{2}=76.6 µm, the propagation distance (*i.e.* the cell length) *L*=20 *Ld*=1.5 mm, the transverse scaling length x_{0}=2 µm, the laser wavelength λ=514 nm, the elastic constant *K*=12 10^{-12} N, the viscous coefficient γ=0.08 kg/ms, the ordinary refractive index *n _{0}*=1.53, the cell thickness

*D*=75 µm, and the bias voltage

*V*=1V. All of these data are consistent with the values reported in experimental investigations [4, 16]. The initial beam widths FWHM are varied between 2 µm and 5 µm. The values for the optical and static permittivity anisotropies of the liquid-crystal molecules are varied between 0.3 and 0.8, and between 10 and 14.5, respectively. The intensity is varied between

*I*=1.0×10

^{9}V

^{2}/m

^{2}and

*I*=5×10

^{12}V

^{2}/m

^{2}.

The main difference between the model used in this paper and the model used previously [5, 18] is in the form and the treatment of the temporal equation for the angle of reorientation. Now the temporal evolution of the angle of reorientation, owing to both light and voltage, is given by a more complete Eq. (2), in which *θ*=*θ*
_{0}+*θ*̂. Previously *θ _{0}* was given by a simple approximate formula or held constant,

*θ*=

_{0}*π/4*. Here, the initial distribution

*θ*is determined in the beginning, independent of the solution of the full system of Eqs. (1) and (2), using boundary conditions and a SOR algorithm for solving PDEs. Such a procedure is preferable when one deals with a steady-state boundary-value problem, as is the case of

_{0}*θ*

_{0}. Only after

*θ*

_{0}is determined the main integration of Eqs. (1) and (2) begins. This provides for a more realistic physical modeling of the system.

Our model, composed of a system of two PDEs, represents a typical example of a NN system. Nonlocality, provided by the temporal diffusion Eq. (2), exerts great impact on the beam propagation in NLCs [19–23]. Spatiotemporal nonlocality means that the response of the medium at a particular point and in a given moment, is not determined exclusively by the wave intensity at that point and in that moment, but also depends on the wave intensity in its vicinity and at earlier moments. Nonlocality also exerts great influence on the self-focusing and MI of propagating beams [20–23].

In the end, to stress again, we restrict our attention mostly to the spatial effects, even though we treat the system in 3+1 spatial and temporal dimensions. A number of movies are still provided, depicting the approach to steady state. Although it is common that people refer to the spatial effects along the propagation axis as the “dynamical effects”, we carefully distinguish between the dynamical and the spatial effects.

## 3. Results of numerical simulation

#### 3.1 Propagation of Gaussian beam in three spatial dimensions and time

### 3.1.1 Influence of input beam intensity

We consider first the behavior of propagating beams in NLCs with Gaussian input intensities. We increase the input beam intensity until spatial solitons – stable beams propagating without change in the transverse profile – are located. The effect of the input intensity variation on the beam propagation is presented in the movies in Fig. 1.

For smaller intensities self-focusing is too weak to keep the beam tightly focused, so that it cannot get through localized. The beam just diffracts and diminishes in intensity (not shown). By increasing the beam intensity (Figs. 1(a) and 1(b)) we achieve a threshold for stable solitonic propagation, however with breathing; the beam preserves its new shape, but its characteristic width and maximum intensity breathe periodically as it propagates. For still higher intensities (Fig. 1(c)) we observe irregular behavior of the beam, so that it is not passing through unchanged, as a soliton. Slight width-modulations, present at the onset of soliton formation, become more pronounced as the intensities are increased. In that range of intensities we see periodic spot size variation during propagation – the breathing of a spatial soliton. The width modulation pitch increases with the intensity increase. Eventually, the soliton is completely broken by the joint action of the longitudinal MI and the transversal instability, developing at the intensities above 5×10^{11} V^{2}/m^{2} [24].

In Fig. 1 we display how the steady-state beam intensities *I*(0,y,z) and the angle reorientation distributions *θ*̂ (0,y,z) are reached, as functions of the propagation distance, for different input beam intensities. The beam intensity *I* is proportional to |*A*|^{2}. The region of intensities in the figure corresponds to the threshold region in which the breathing solitons (BSs) are found. Figure 1 represents typical behavior during the propagation of input Gaussian beams in NLCs. A comparison between the beam intensity and the reorientation angle distributions shows that they change in unison, as they should; Eq. (2) is solved independently at every *z* position, using the current propagating value of *A* from Eq. (1). The reorientation angle distribution follows closely the intensity distribution, but owing to broader nonlocality for the given parameters, it is always wider. The profiles exhibit similar behavior to the experimental profiles reported in [4, 17, 25, 26] (cf. Figs. 6 in [4], 3 in [17], 5 in [25] and 5 in [26]).

In addition to Gaussians, we tried other initial beam profiles: hyper-Gaussians, hyperbolic secants, etc. Invariably, we obtained BSs whose profile is different from the input beams. We also launched the obtained BS profiles as initial beams, to again find breathing beams. The soliton breathing (maximum intensity breathing periodically as the soliton propagates) was the only way of soliton propagation that we could identify in this NN system. We did not observe steady spatial solitons. This, of course, doesn’t mean that such solitons do not exist in the system - it only means that we could not find them.

### 3.1.2 Influence of the input beam width FWHM

When the same sequence of increasing input intensity simulations is repeated for different input widths, for fixed parameters Δε^{OPT}=0.4 and Δε^{DC}=14.5, similar behavior is found. In Fig. 2(a) we show the cases of BS propagation, for FWHM=2 µm, 2.5 µm, 3 µm, 3.5 µm, 4 µm, 4.5 µm and 5 µm. For each input beam width, appropriate input beam intensity can be found to establish the existence of a breather. We also note that the changes in FWHM cause the changes in the period of oscillation. In Fig. 2(b) we depict the intensities of BSs as functions of the input widths (FWHM ranging as in Fig. 2(a)). For smaller FWHM higher input intensities are needed for the emergence of breathers. The curve fitted through the points is *y*=*1*/*x ^{4}*. Such a functional dependence is predicted by Snyder and Mitchell [27] for accessible (highly nonlocal) solitons in nonlocal Kerr media. It is clear from Figs. 1–3 that for the range of parameters considered in this paper we are close to the strongly nonlocal regime. This also helps explain the lack of steady solitons in our numerical simulations, as well as in the corresponding experiments [21, 25, 26]. It is typical of highly nonlocal solitons to display breathing behavior.

For the intensity of soliton breathing, (input FWHM=3.5 mm, *I*=8.6×^{10}10 V^{2}/m^{2}, Fig. 2(a)), we present in Fig. 3 the intensity and the optically induced molecular reorientation in the middle of the crystal, *I*(0,0,*z*) and *θ*̂ (*0*,*0*,*z*) and the corresponding FWHM of its transverse profiles, as functions of the propagation distance. In Fig. 3 we also depict the transverse profiles of BS and the induced tilt angle at the exit (*x*, *y*) plane. We show how the BS changes its width and intensity, and how the molecular reorientation changes during *z* propagation. We invariably observe soliton breathing - the beam preserves its new shape but its characteristic width and maximum intensity breathe as it propagates.

### 3.1.3 Influence of the optical and the static permittivity anisotropy

If we consider the change in Δε^{OPT} and Δε^{DC} for the input Gaussian beam propagation, we observe again stable beams in a narrow threshold region of control parameters. Below the threshold region the beams diffract, above it spatiotemporal instabilities are observed, as the input intensity and the optical and static permittivity anisotropies of the liquid crystal molecules are increased. As before, we only observe spatial soliton breathers.

If the same sequence of increasing input intensity simulations is repeated for different birefringences (0.3–0.8), and for the fixed input beam width FWHM=4µm, a similar behavioral pattern is observed (Fig. 4). We change Δε^{OPT} and repeat the sequence of increasing input intensities for two values of Δε^{DC}, Δε^{DC}=10 and Δε^{DC}=14.5.

In Fig. 4 we show the intensities of BSs as functions of Δε^{OPT} and the molecular reorientation induced by the electric field only *θ _{0}*(

*x*,

*y*), as functions of the transverse variables, for the two mentioned values of Δε

^{DC}. Inset in Fig. 4(A) depicts soliton breathing propagation in

*z*, for Δε

^{OPT}=0.4 and Δε

^{OPT}=0.3 (Δε

^{DC}=14.5). We can see that if we increase birefringence (for fixed ε

^{DC}) the threshold input intensities for soliton propagation are decreasing; for smaller birefringence the appropriate input beam intensity for the emergence of BSs appears at a higher intensity. For smaller Δε

^{DC}=10 we see similarly an increase in the intensity needed for the existence of stable BSs. We conclude that the changes in FWHM lead to the changes in the intensity and in the period of oscillation along

*z*direction.

### 3. 1. 4. Influence of the change in boundary condition θ(x=-D/2)=θ(x=D/2)=30^{0}

If we repeat the same sequence of increasing input intensity simulations for different input widths FWHM, for fixed parameters Δε^{OPT} and *L*, but change boundary conditions, we find similar behavior as before. For each input beam width the appropriate input beam intensity can be found for the existence of BSs. In Fig. 5 we show the cases of BS propagation, for different FWHM (2.5µm–4µm).

#### 3.2. Propagation of input Gaussian beams using a steady-state procedure

As an independent check on our time-dependent solution procedure, we considered the steady state propagation in NLCs, and found similar behavior - the existence of BSs as the only stable propagation mode of localized beams in NLC model at hand. The propagation in the steady state is achieved by putting *∂θ*/*∂t*=0 in Eq. (2) from the beginning and by solving Eqs. (1) and (2) as a time-independent system of PDEs. The temporal solution procedure is now eliminated, being substituted by a SOR algorithm. However, the results are very similar to the ones obtained in the limit of long time evolution in the original time-dependent procedure.

### 3. 2. 1. Influence of input beam intensity

For the steady state case we again examine the beam propagation, similar to the previous analysis. We consider the behavior of propagating single input Gaussian beams in NLCs. We increase beam intensity until spatial solitons are located. The effect of the input Gaussian intensity variation on the single beam propagation is presented in Fig. 6. Again, for smaller intensities, self-focusing is too weak to keep the beam tightly focused - it can not pass through localized. One should note the behavior of the beams with the smallest and the next to the smallest input intensities of *I*=1.0×10^{+9} V^{2}/m^{2} and *I*=1.0×10^{+10} V^{2}/m^{2}. While the beam of the smallest intensity (inset) smoothly diminishes in intensity with propagation, the beam next to it is initially greatly reduced in intensity, but then it recovers in an irregular fashion.

By increasing the beam intensity we achieve regular breathing and stable solitonic propagation. For FWHM=4 µm we see stable soliton propagation for an input intensity of *I*=3.5×10^{+10} V^{2}/m^{2}.

The comparison between the maximum beam intensity and the reorientation angle values shows that they change in unison, the amplitude of the angle of reorientation oscillation being small. This again comes out because of the nonlocal nature of nonlinearity. In Fig. 1 the comparison between the optical intensity distribution and the reorientation angle distribution shows the same, that the amplitude of the reorientation angle oscillation is small.

### 3. 2. 2. Comparison between the steady-state propagation procedure and the propagation in time.

We also compare directly results between the steady state propagation procedure and the propagation in time procedure, after steady state is reached. We find that the agreement is very good for shorter propagation distances (~10L_{D}). For longer propagations small quantitative discrepancies appear, however the qualitative behavior in both numerical procedures remains the same. When one compares the spatial distributions obtained by the two procedures, the differences are hardly discernible. It is reassuring that two very different numerical procedures provide very similar results over a broad range of parameters. Still, we believe that the propagation in time procedure provides more reliable and physically more transparent results.

A comparison between the time-dependent and the steady-state procedure is displayed in Fig. 7. Black curve in Fig. 7 represents the steady-state propagation result, for a set of parameters. Three profiles at the noted maximum, middle, and minimum points of the steady-state curve are used as the input intensity profiles in the time-dependent procedure. The resulting curves after steady state is reached are represented in red, green, and blue. One can note essentially the same qualitative breathing behavior, only displaced in phase. The discrepancy noted at larger distances probably stems from the fact that the initial profile in the steady-state procedure was Gaussian, while the initial profiles in the time-dependent procedure were established BS profiles.

## 4. Summary

In this paper we studied a model for the laser light propagation in a cell containing a liquid crystal in the nematic phase. We investigated the behavior of beams in time and in three spatial dimensions, using an appropriately developed theoretical model and a numerical procedure based on the fast Fourier transform. Starting from the standard equations describing nonlocal and nonlinear interactions in NLCs, we developed a numerical model that exhibits rich dynamical behavior in three spatial dimensions and time. However, our attention here was mostly focused on the behavior in steady state.

We presented beam propagation in NLCs for different values of the beam and material parameters. We demonstrated the formation of stable solitons in a narrow threshold region of beam intensities for fixed parameters and display soliton breathing.

We examined the influence of input intensities, input FWHM, Δε^{OPT}, Δε^{DC}, and boundary conditions. For different combination of these parameters, we always obtained consistently similar results, which present typical behavior during the propagation of input Gaussian beams in NLCs. The typical stable soliton propagation mode is the soliton breathing. Also, we considered the behavior of propagating single input Gaussian beams using a steady state procedure, and found similar results. We could not locate steady (not breathing) spatial solitons in any of our numerical simulations, despite repeated efforts. We conclude that the soliton breathing is the predominant stable way of soliton propagation in NLCs.

## Acknowledgments

This study is supported by the Qatar National Research Foundation project NPRP 25-6-7-2 and by the Ministry of Science and Technological Development of the Republic of Serbia, under the project OI 141031.

## References and links

**1. **I. C. Khoo, *Liquid Crystals: Physical Properties and Nonlinear Optical Phenomena* (Wiley, New York, 1995).

**2. **Y. S. Kivshar and G. P. Agrawal, *Optical solitons* (Academic Press, San Diego, 2003).

**3. **G. Assanto, M. Peccianti, K. Brzdakiewicz, A. De Luca, and C. Umeton, “Nonlinear wave propagation and spatial solitons in nematic liquid crystals,” J. Nonlinear Opt. Phys. Mater. **12**, 123–134 (2003).
[CrossRef]

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

**5. **A.I. Strinić, D.V. Timotijević, D. Arsenović, M.S. Petrović, and M.R. Belic, “Spatiotemporal optical instabilities in nematic solitons,” Opt. Express **13**, 493–504 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=OPEX-13-2-493.
[CrossRef] [PubMed]

**6. **P. D. Rasmussen, O. Bang, and W. Krolikowski, “Theory of nonlocal soliton interaction in nematic liquid crystals,” Phys. Rev. E **72**, 066611 1–7 (2005).
[CrossRef]

**7. **G. D’Alessandro and A. A. Wheeler, “Bistability of liquid crystal microcavities,” Phys. Rev. A **67**, 023816 1–12 (2003).

**8. **J. F. Henninot, M. Debailleul, and M. Warenghem, “Tunable non-locality of thermal non-linearity in dye doped nematic liquid crystal,” Mol. Cryst. Liq. Cryst. **375**, 631–640 (2002).
[CrossRef]

**9. **T. M. Monro, C. Martijn de Sterke, and L. Poladian, “Investigation of waveguide growth in photosensitive germanosilicate glass,” J. Opt. Soc. Am. B. **13**, 2824–2832 (1996).
[CrossRef]

**10. **S. LaRochelle, V. Mizrahi, G. I. Stegeman, and J. E. Sipe, “Growth dynamics of photosensitive gratings in optical fibers,” Appl. Phys. Lett. **57**, 747–749 (1990).
[CrossRef]

**11. **T. M. Monro, C. Martijn de Sterke, and L. Poladian, “Catching light in its own trap,” J. Mod. Opt. **48**, 191–238 (2001).

**12. **A. S. Kewitsch and A. Yariv, “Nonlinear optical properties of photoresists for projection lithography,” Appl. Phys. Lett. **68**, 455–457 (1996).
[CrossRef]

**13. **A. A. Sukhorukov, S. Shoji, Y. S. Kivshar, and S. Kawata, “Self-written waveguides in photosensitive materials,” J. Nonlinear Opt. Phys. Mater. **11**, 391–407 (2002).
[CrossRef]

**14. **K. D. Dorkenoo, F. Gillot, O. Crégut, Y. Sonnefraud, A. Fort, and H. Leblond, “Control of the refractive index in photopolymerizable materials for (2+1)D solitary wave guide formation,” Phys. Rev. Lett. **93**, 143905 1–4 (2004).
[CrossRef]

**15. **M. Peccianti, A. De Rossi, G. Assanto, A. De Luca, C. Umeton, and I. Khoo, “Electrically assisted self-confinement and waveguiding in planar nematic liquid crystal cells,” Appl. Phys. Lett. **77**, 7–9 (2000).
[CrossRef]

**16. **J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, “Simulations and experiments on self-focusing conditions in nematic liquid-crystal planar cells,” Opt. Express **12**, 1011–1018 (2004), http://www.opticsinfobase.org/oe/abstract.cfm?URI=OPEX-12-6-1011.
[CrossRef] [PubMed]

**17. **J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, “Time-dependence of soliton formation in planar cells of nematic liquid crystals,” IEEE J. Quantum Electron. **41**, 735–740 (2005).
[CrossRef]

**18. **A. Strinić, D. Jović, M. Petrović, D. Timotijević, N. Aleksić, and M. Belić, “Counterpropagating beams in nematic liquid crystals,” Opt. Express **14**, 12310–12315 (2006), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-25-12310.
[CrossRef] [PubMed]

**19. **W. Krolikowski, O. Bang, N. Nikolov, D. Neshev, J. Wyller, J. Rasmussen, and D. Edmundson, “Modulational instability, solitons and beam propagation in spatially nonlocal nonlinear media,” J. Opt. B **6**, S288–S294 (2004).
[CrossRef]

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

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

**22. **M. Peccianti, C. Conti, and G. Assanto, “Optical modulational instability in a nonlocal medium,” Phys. Rev. E **68**, 025602(R) 1–4 (2003).
[CrossRef]

**23. **G. Assanto, M. Peccianti, and C. Conti, “One dimensional transverse modulational instability in nonlocal media witha reorientational nonlinearity,” IEEE J. Sel. Top. Quantum Electron. **10**, 862–869 (2004).
[CrossRef]

**24. **A. I. Strinic and M. R. Belic, “Beam propagation in nematic liquid crystals,” Acta Phys. Pol. A , **Vol. 112(5)**, 877–883 (2007).

**25. **M. Peccianti, C. Conti, G. Assanto, A. De Luca, and C. Umeton, “Nonlocal optical propagation in nonlinear nematic liquid crystals,” J. Nonlinear Opt. Phys. Mater. **12**, 525–538 (2003).
[CrossRef]

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

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