1. Introduction

Nematic liquid crystals (NLC) are ubiquitous materials, found in many consumer electronic devices, probably including the screen on which this very paper is being read.

NLC are known to exhibit enormous optical nonlinearities, owing to large refractive index anisotropy, coupled with the optically-induced collective molecular reorientation. These qualities make them suitable for investigation using lasers of relatively low power and low-cost detection equipment. For that reason they have been the subject of considerable study in recent years, from both experimental and theoretical points of view.

Research into nonlinear optical properties of NLC has been mostly limited to transverse effects in thin films, i.e., to one-transverse-dimensional (1D) geometries. Interest into instabilities has been rekindled by the work of Braun et al. [1] on filamentation and undulation. The first (1+1)D solitons were observed by Karpierz et al. [2] in a homeo-tropically- aligned NLC thin film waveguide. Lately emphasis has been shifted to 2D transverse geometries, that is, to the effects in bulk NLC. Warenghem et al. [3] and Derrien et al. [4] utilized the thermal nonlinearity of dye-doped NLC to demonstrate (2+1)D self-confinement in capillaries. Among these, of special interest are the localized (but nonlocal!) self-trapped optical wave structures, including most prominently the spatial solitons. Recently, Peccianti et al. [5–7] and Assanto et al. [8,9] have demonstrated stable (2+1)D solitons in bulk undoped NLC, taking advantage of the reorientational molecular response in the presence of a voltage-induced pre-tilt. Such response is saturating and nonlocal in space, as well as nonlocal in time, hence allowing for incoherent excitations [6].

We investigate the behavior of optical propagation in time and in 3 spatial dimensions, using an appropriately developed theoretical model and a numerical procedure based on the fast Fourier transform (FFT). We find that solitons exist in a narrow region of beam intensities. Below this region the beams diffract, above the solitons filament and breathe. At higher intensities and larger birefringence a modulational instability (MI) of self-trapped beams is observed, in the form of long-lasting convective instability propagating along the beam. We also consider the propagation and interactions of more complicated beam structures, such as dipoles, arrays of beams, and optical vortices. We obtain good agreement with experimental results [7–9].

A nematic liquid crystal possesses some properties of both liquids and solids. It contains rod-like molecules which exhibit orientational alignment without positional order. They behave in a fluid-like fashion, but display a long-range order that is characteristic of crystals [10]. Several types of long-range order are observed in thermotropic LC, the simplest being the one in which the position of molecules is arbitrary but their orientation is nearly the same. The unique property of NLC is their ability to change optical properties under the action of an external electric field, which enables the macroscopic reorientation of the director θ of NLC. In other words, the light incident on a LC modifies the electric permittivity tensor, leading to the reorientational nonlinearity.

The average alignment of the molecules is associated with the director (a unit pseudovector pointing along the preferred direction), and in our work it is prescribed to be parallel to the top and bottom bounding surfaces. This form of alignment is called homogeneous, whereas when the average molecular alignment is normal to the boundary, it is called homeotropic. As is well documented (see, e.g., de Gennes & Prost, 1993) [11], at a critical value of the strength of magnetic or electric field, a static distortion of the nematic occurs and this phenomenon is referred to as the Freedericksz transition.

2. The model

We assume NLC to be a birefringent and positive uniaxial medium, allowing for self-focusing in the bulk. We also assume a scalar representation of the optical field. These assumptions make possible simple theoretical modeling and numerical handling that captures the essential features of experimental investigations. The distortion of molecular orientation in NLC is interpreted by a reorientation angle θ of the director, measured from the propagation direction. In the presence of an externally applied (low frequency) voltage the evolution of a slowly-varying beam envelope A, linearly polarized along x axis and propagating along z axis, is well described by the paraxial wave equation [5]:

where k=k₀n₀ is the wave vector in the medium and ε_a =$n_e^{\mathit{2}}$ -$n_{\mathit{0}}^{\mathit{2}}$ is the birefringence of the medium. The rest distribution angle θ_rest in the presence of a low-frequency electric field is modeled by [5]:

with θ₀ (V) being the orientation distribution due to the applied voltage far from the input interface, θ_in is the director orientation at the boundary z=0, and z̄ is the relaxation distance. Temporal evolution of the angle of reorientation is given by the diffusion equation [12, 13]:

where γ is the viscous coefficient and K is Frank’s elastic constant. Here θ is the overall tilt angle, due to both light and voltage.

Using the rescaling z=zk${{\mathrm{x}}_0}^2$, x=xx₀, y=yx₀ and t=tτ, we transform the equations into a dimensionless form:

where τ is the relaxation time and x₀ is the transverse scaling length. Equations (4) and (5) form the basis of the model.

By solving these equations we will describe beam propagation in both space and time. We develop a novel numerical procedure, based on the split-step FFT method, utilizing our prior experience in treating beam propagation in photorefractive media [14, 15]. The novelty is that we treat concurrently a system of coupled partial differential equations, in both space and time. Upon discretization, the diffraction and diffusion are treated in the inverse transverse space, whereas the nonlinearities are treated in the direct space. The temporal equation for the angle of reorientation and the spatial propagation equation for the beam envelope are solved together in a system of nested loops, including an iterated convergence loop. The procedure is as follows.

The object is to determine the spatial distributions of θ and A at each moment of time that will satisfy Eqs. (4) and (5). Starting from a given distribution of θ, the incident field A at z=0 is integrated along z, to obtain a distribution of A. Then the distribution of θ is integrated for a time step at each z point. Note that this entails the solution of a partial differential equation. The field A now does not correspond to the new distribution of θ, and it has to be propagated again along z. This two-step procedure is iterated until stable self-consistent distributions of θ and A are obtained. Then the distributions are updated and the temporal loop advanced for a time step. The convergence in both temporal and self-consistency loops signifies that a steady-state solution is found. However, this need not always be the case. Time-dependent, dynamical states can also be observed, when the temporal loop refuses to converge. Our procedure is uniquely suited for observing the dynamical states in slowly-varying physical systems, where the fast optical fields are slaved to the slow change in the nonlinearity.

In all simulations the following data are kept constant: the diffraction length Ld=k${{\mathrm{x}}_0}^2$=79 µm, the propagation distance L=6.3 Ld=0.5 mm, the transverse scaling length x₀=2 µm, the laser wavelength λ=514 nm, the relaxation distance z̄=40 µm, the elastic constant K=0.7 10^-11 N, the viscous coefficient γ=0.08 kg/ms, the ordinary refractive index n₀ =1.53, the director orientation at the boundary θ_in =π/2, the orientation distribution θ₀ =π/4, and the initial beam widths FWHM=4 µm. Two values for the birefringence are used, ε_a=0.5 and ε_a=0.8, and the intensity is varied between I=0.5×10⁺¹⁰ V²/m² and I=5×10⁺¹³ V²/m². All of these data are consistent with the values reported in experimental investigations [5–9]. We also adopt the geometry of NLC cell and beam propagation, as well as the boundary conditions from these references.

3. Numerical results

Numerical studies of partial differential equations describing beam propagation in electrically-biased planarly-oriented NLC are performed in different conditions and for a variety of beam configurations. Single Gaussian, dipole, vortex and beam-array propagation are investigated. Stable self-confined solutions (optical solitons) are demonstrated, as well as their modulational instabilities. These phenomena are also observed in experiments [5–9, 16–18]. We have found that solitons can exist in a narrow range of beam powers. Below this threshold range the self-focusing of beams to various degrees is seen, and above the filamentation, breathing, and modulational instabilities of solitons are observed.

All the pictures and movies in the (x,y) plane are presented at the exit face of the crystal (z=L), and all the pictures and movies in the (y,z) plane are for the x=0 plane (in the middle of the crystal).

3.1 Solitons

Liquid crystal cells are often made very thin. However, NLC can be treated as bulk medium when the transverse dimensions of the cell are much larger than the size of the input beam. Such a configuration is convenient for studying 2D solitons. The basic physical mechanism for the existence of 2D solitons is self-focusing, and for LC it can be understood as follows. As the guided field reorients LC molecules, the refractive index of the medium increases most where the field is most intense, which in turn modifies the guided field itself. The perturbation of the reorientation angle is connected to the change in the index of refraction. Owing to large anisotropy of LC, the guided mode changes its profile significantly, resulting in a lensing effect and self-focusing. Because of the long-range order in LC, the reorientation effect depends not only on the local value of the electric field but also on the entire field profile across the waveguide [10].

Nonlocality has a great impact on the beam propagation in NLC [7, 16, 17]. 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 this moment, but also depends on the wave intensity in its vicinity and at earlier moments. Nonlocality also exerts great influence on the modulational instability (MI). MI represents an exponential growth of a weak spatial perturbation of the wave as it propagates [16]. The gain causes the amplification of sidebands, which leads to the breakup of wave and the appearance of transverse localized structures through the process of filamentation. Saturable nonlocal nonlinearity tends to suppress the exponential growth of instability, decreasing the growth rate and the width of the instability band. But, it can not eliminate the instability completely. The first experimental observation and study of MI in NLC is given in [18].

 

Fig. 1. Soliton propagation, shown in the (y,z) plane, for different input intensities: a) I=0.5×10⁺¹⁰ V²/m², b) movie for I=1.8×10⁺¹⁰ V²/m² (134 KB) c) I=3×10⁺¹⁰ V²/m², d) I=4×10⁺¹⁰ V²/m², e) I=5×10⁺¹⁰ V²/m², f) I=1×10⁺¹¹ V²/m², g) I=5×10⁺¹¹ V²/m², h) I=5×10⁺¹² V²/m², i) I=5×10⁺¹³ V²/m². For all simulations FWHM=4 µm, L=0.5 mm and ε_a=0.5. Stationary situation is depicted, at t=5 τ.

Download Full Size | PDF

We first consider the behavior of propagating single Gaussian beams in NLC. We increase beam intensity until spatial solitons – stable beams propagating without change in the profile -are located. The effect of input intensities variation on the single Gaussian beam propagation is presented in Fig. 1. For smaller intensities (Fig. 1(a)) self-focusing is too weak to keep the beam tightly focused, so that it can not get through localized. By increasing beam intensity (Fig. 1(b)) we achieve stable propagation i.e. a soliton. For FWHM=4µm we have stable soliton propagation for input intensities in the range of I=1.8×10⁺¹⁰ V²/m² to I=3×10⁺¹⁰ V²/m², depicted in Fig. 1(c). It should be mentioned that, with Gaussian input, a slight beam width modulation is always present, as the soliton forms. These width modulations become more pronounced as the intensities are increased. In that range of intensities we see a periodical spot size variation along propagation – the breathing of soliton. The pitch of width modulation increases with intensity increase (Fig. 1(d–f)). Eventually, the soliton is broken up (Fig. 1(h), 1(i)) by joint action of the longitudinal modulation instability, just described, and the transversal one, mixing in at the intensities above I=5×10⁺¹¹ V²/m², (Fig. 1(g)). Transverse modulational instabilities also cause filamentation, examples of which will be given later (Fig. 6(a)).

If the same sequence of increasing input intensity simulations is repeated for larger birefringence, ε_a=0.8, a similar behavioral pattern is observed, but this time for different intensities, e.g., stable solitons emerge at a smaller intensity I=1.0×10⁺¹⁰ V²/m² (Fig. 2). A comparison between the beam and angle reorientation distributions in the same figure reveals that they change in unison, the angle reorientation distribution being wider. This comes because of the nonlocal nature of nonlinearity. In Fig. 3 a comparison between the optical intensity distribution and the reorientation angle distribution shows this clearly (for the case from Fig. 2(c) and 2(i)). The profiles are made in the (x,y) plane, at the same moment. They exhibit similar behavior to the one reported in [8].

 

Fig. 2. Soliton propagation in the (y,z) plane, for different input intensities and ε_a=0.8. Upper row I(y,z); lower row θ(y,z) at t=10 τ. a) and g) I=0.5×10⁺¹⁰ V²/m², b) and h) I=1×10⁺¹⁰ V²/m², c) and i) I=1.8×10⁺¹⁰ V²/m², d) and j) I=3×10⁺¹⁰ V²/m², e) (253 KB) and k) (509 KB), movies for I=4×10⁺¹⁰ V²/m², f) and l) I=5×10⁺¹² V²/m², all in V²/m². For all simulations FWHM=4 µm and L=0.5 mm.

Download Full Size | PDF

 

Fig. 3. Soliton at the exit (x,y) plane, for t=5τ. a) I(x,y); b) θ(x,y). The input intensity I=1.8×10⁺¹⁰ V²/m², FWHM=4 µm, L=0.5 mm, and ε_a=0.8. Here, x and y axes represent actual points as used in simulation; physical lengths spanned are the same for x and y (NLC film thickness) and amount to 116 µm.

Download Full Size | PDF

The next figure (Fig. 4) presents the temporal development of the angle of reorientation (reorientation angle versus time) at a higher intensity. We note a gradual broadening, characteristic of diffusion-type dynamics. At the beginning one can notice transient transverse patterns evolving, which have no physical significance [19]. What is important here is the final steady-state distribution, which also gives an indication of the refractive index perturbation.

 

Fig. 4. Snapshots of θ(x,y) for soliton propagation, shown in the (x,y) plane, at different moments. The input intensity I=5×10⁺¹¹ V²/m², FWHM=4 µm, L=0.5 mm and ε_a=0.8. Here, x and y are presented in the same way as in Fig. 3.

Download Full Size | PDF

 

Fig. 5. Comparison of beam propagation, shown in the (y,z) plane, for different input intensities and different ε_a. Upper row ε_a=0.5, lower row ε_a=0.8. a) and e) I=0.5×10⁺¹⁰ V²/m², b) and f) I=1.8×10⁺¹⁰ V²/m², c) and g) I=4×10⁺¹⁰ V²/m², d) and h) I=5×10⁺¹² V²/m². For all simulations FWHM=4 µm, L=0.5 mm, t=5 τ.

Download Full Size | PDF

In Fig. 5 a comparison of the propagation with two different ε_a is depicted, ε_a=0.5 (Fig. 5(a–d)) and ε_a=0.8 (Fig. 5(e–h)). One can see that there exist similar phases in the propagation for different ε_a, which happen at different intensities (i.e., the soliton instabilities develop similarly). The values of intensity where the stable soliton propagation is observed are different for the two cases. For higher ε_a this value is smaller and similar phases in filamentation develop at lower intensities.

In Fig. 6(a) and Fig. 6(b) we show examples of beam distribution, visualized by isosurfaces in 3D, for an input beam intensity for which modulation instability and strong filamentation occur. An interesting feature is a long-lasting convective instability seen at the head of the beam, as it bores an optical path through the crystal. Such an instability is more pronounced at higher intensities and larger birefringence. In our numerical simulations they last for long time (a few τ, tens of s), and we believe that they can be observed experimentally. One can follow the propagation of solitons in 3D (for the case I=5×10⁺¹¹ V²/m²) with a movie presented in the same figure. These two figures represent 3D view at different moments (for t=1.4 τ, and t=2.95 τ). Several isosurfaces of the same soliton are merged onto the same picture, represented by different colors, the smaller intensities are lighter, the higher are darker.

 

Fig. 6. (a) Beam filamentation, shown in 3D, for input intensity I=5×10⁺¹¹ V²/m² at t=1.4 τ. FWHM=4 µm, L=0.5 mm and ε_a=0.5. (A) Movie of the view in the (y,z) plane (236 KB), (B) 3D view using four intensity levels (in relative units: 0.01, 1, 10, 30). (C) 3D view for the highest relative intensity (30).

Download Full Size | PDF

 

Fig. 6. (b) Continuation of a), the situation at t=2.95 τ, when the convective instability is about to end. FWHM=4 µm, L=0.5 mm and ε_a=0.5. (A) (y,z) plane, (B) 3D view, with three levels of intensity (in relative units: 1, 10, 100).

Download Full Size | PDF

3.2 Dipoles

Nonlocal nature of nonlinearity in NLC allows for interesting interaction possibilities between two or more beams, or between the components of more complicated multi-component beams. Here we consider only the interaction between two simple Gaussian in-phase beams, and between the (out-of-phase) components of a dipole beam. Below we also address the interaction of arrays of solitons that are in and out of phase.

Let us consider first two symmetric in-phase solitons, as in Fig. 7 [8]. The interaction takes place where the θ perturbations overlap. In a local medium the index overlap coincides with the field-intensity overlap. In NLC, which represents an example of strongly nonlocal medium, the perturbation diffuses out of the excitation regions (Fig. 2), providing for interaction between solitons over wider distances. In such a case [8, 20, 21], the interaction tends to be attractive and weakly phase dependent. Repulsion can only be observed in the beam-crossing region and for very small angles between them.

 

Fig. 7. Dual beam propagation, shown in the (y,z) plane, for different distances between beams, which are in phase. Distances between the components are: (a) 8 µm, (b) 12 µm, (c) 16 µm, (d) 40 µm. Other data: input intensity I=1.8×10⁺¹⁰ V²/m², FWHM=4 µm, L=0.5 mm and ε_a=0.5.

Download Full Size | PDF

Figure 7 displays the interaction of two in-phase solitons. We notice attraction and collapse for small distance between the components, which is also experimentally confirmed [8]. Once established, the solitons attract each other in the beginning, after that they fuse and collapse into a single filamented beam. Fusion can be noticed in the reduction of the period and the amplitude of beam interlacing, owing to the partial inelasticity of the interaction [8]. For smaller input separation between components, two solitons collapse into one beam faster (Figs. 7(a), (b), and (c)). After collapsing, the two beams continue to propagate together for a long distance, behaving like a single beam. For initial distance between component beams equal or larger than the relaxation distance (40 µm here), there is no attraction and fusion between the components (Fig. 7(d)).

Figure 8 presents the comparison between the propagation of dual beams with components that are in phase (upper row), and out of phase (lower row). In the case of components out of phase, we also notice attraction, but the fusion and collapse of beams is prevented. Owing to the singularity (I=0) between the beam components, the crossing is avoided, and no interlacing occurs.

 

Fig. 8. Movies of the dual beam propagation for different phase relations between beams. (a) (144 KB), (b) (42 KB) and (c) (2.486 MB) Components in phase, (d) (382 KB), (e) (89 KB) and (f) (753 KB) phase difference of π between the components. (a) and (d) Intensity in the (y,z) plane. (b) and (e) Intensity in the exit (x,y) plane. (c) and (f) Phase distributions in the (x,y) plane. Other data as in Fig. 7(a).

Download Full Size | PDF

3.3 Arrays

The behavior of dipoles suggests an idea how to construct more stable arrays of solitons, similar to the photorefractive nonlinear optics [15, 22, 23]. We investigate the propagation of 5-by-5 arrays, which in one case are all in-phase, and in the other case are out-of-phase, with the alternating phase distributed in a chess-board pattern (Fig. 9). For the in-phase case we observe strong attraction and straight fusing of the beams, as expected. However, in the out-of-phase case, after an initial attraction, the beams refuse to fuse, although strong finite-size effects and large lattice deformation are observed. After a longer propagation, the beams nonetheless collapse to a complicated transverse pattern. Considering the stability of propagating arrays, the out-of-phase case offers better results, but not for much. The reason for small improvement in the propagation characteristics of arrays, we believe, is complicated intricacies of the beam interaction in NLC. Because the interaction in NLC is of long-range order, out-of-phase arrays also collapse and the lattice is destroyed, similar to the case with no input phase difference between the lattice components, but over a longer propagation distance.

 

Fig. 9. Movies of the propagation of 5×5 arrays, L=1mm, initial distance between the array components 20 µm. (a) (1.506 MB) and (b) (1.237 MB) Array in phase; (c) (1.493 MB) and (d) (859 KB) Array out-of-phase, with the chess-board arrangement of input phase. (a) and (c) Intensity in the (y,z) plane; (b) and (d) Intensity in the (x,y) plane. Other data are the same as in Fig 7.

Download Full Size | PDF

3.4 Vortices

Beams carrying angular momentum are convenient for investigating the propagation and interaction characteristics of different arrangements. We confine our attention to the propagation of the basic Gauss-Laguerre mode, carrying unit topological charge.

In Fig. 10 an example of vortex propagation is depicted. For an input intensity I=5×10⁺¹¹ V²/m² and FWHM=8 µm a stable self-focused structure is obtained. For smaller intensity, the vortex could not get through the medium, i.e., it was diffracted. For the same input intensity and smaller vortex width, FWHM=4 µm, a breathing, almost stable structure is seen. For smaller intensity and the same FWHM=4 µm, vortex did not propagate through. For the same input intensity and two times broader vortex (FWHM=8 µm), we finally got a stable structure.

 

Fig. 10. Movies of the stable vortex propagation, (a) (1.205 MB) Intensity in the (y,z) plane. (b) (1.532 MB) Intensity, and (c) (1.522 MB) phase in the exit (x,y) plane. I=5×10⁺¹¹ V²/m², FWHM=8 µm, L=0.5 mm and ε_a=0.5.

Download Full Size | PDF

4. Conclusions

In conclusion, we investigated spatiotemporal optical propagation in nematic liquid crystals. A time-dependent model for the beam propagation and the director reorientation in NLC is numerically treated for various beam configurations. Starting from the standard equations describing nonlocal and nonlinear interactions in NLC, we developed a numerical model that exhibits rich dynamical behavior in three spatial dimensions and time. Solitons can propagate through the medium only in a narrow region of beam intensities. For lower intensities the beams diffract, for higher intensities the solitons filament. At higher intensities and higher birefringence dynamical effects are seen, in the form of soliton breathing and convective MI instabilities. For multiple beams – dipoles and arrays – generally attractive interaction is observed, resulting in collapse and fusion of beams. We find our numerical results to be in good qualitative agreement with experimental data.