1. Introduction

Nematic liquid crystals (NLC) exhibit huge optical nonlinearities, owing to large refractive index anisotropy, coupled with the optically-induced collective molecular reorientation. They behave in a fluid-like fashion, but display a long-range order that is characteristic of crystals [1, 2]. Thanks to the optically nonlinear, saturable, nonlocal and nonresonant response, NLC have been the subject of considerable study in recent years, from both theoretical [3, 4] and experimental points of view [5–10].

Nematic liquid crystals possess properties of both liquids and solids. They contain rod-like molecules which display orientational alignment without positional order. Several types of long-range order are observed in thermotropic NLC, the simplest one being when the position of molecules is arbitrary but their orientation is the same. A useful property of NLC is the ability to change optical properties under the action of an external electric field, which enables the macroscopic reorientation of the angle of the director (a unit pseudovector pointing along the predominant direction of molecules). In other words, the light incident on NLC modifies the electric permittivity tensor, leading to the reorientational nonlinearity.

The average alignment of the molecules is associated with the director, and in our work it is prescribed to be parallel to the top and bottom bounding surfaces. As it is well known [1], at a critical value of the strength of applied electric field, a static distortion of the nematic occurs, and this phenomenon is referred to as the Freedericksz transition.

In an earlier publication [11] we investigated the propagation of laser beams in NLC, both in time and in 3 spatial dimensions, using an appropriately developed theoretical model and a numerical procedure based on the split-step fast Fourier transform technique. We extend the analysis here to the counterpropagating (CP) beams. We find that the spatial solitons [12] exist in a narrow region of beam intensities, similar to the case of copropagating beams [11], but at lower values of the control parameters. Below this region the beams diffract, above the region the beams display periodic and even chaotic behavior. Novelty here is that for higher intensities a transversal motion of the beams (the undulation, or the beam shift) is observed. We also consider the propagation of broader CP Gaussian beams, which offers rich opportunities for observing a complex pattern-forming dynamical behavior. The propagation and interactions of more complex beam structures, such as optical vortices, are also studied.

2. The model

As mentioned, a useful property of NLC is to change optical properties under the influence of an external electric field, producing a reorientation of the director tilt angle θ. To describe the evolution of slowly-varying CP beam envelopes A and B, linearly polarized along the x axis and propagating along the z axis, we utilize the paraxial wave equations [4, 5, 8, 9]:

where A and B are the forward and the backward propagating beam envelopes, 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, 8]:

with θ₀ (V) being the orientation distribution due to the applied voltage V far from the input interface. θ_in is the director orientation at the boundaries z=0 and z=L, where L is the propagation distance and z̅ is the relaxation distance. The temporal evolution of the angle of reorientation is given by the diffusion equation [2, 4]:

where γ is the viscous coefficient and K is Frank’s elastic constant. Here θ is the overall tilt angle, owing to both the light and the voltage influence. Using the rescaling z=zk${{\mathrm{x}}_0}^2$, x=xx₀, y=yx₀, and t=tτ, the equations are transformed into a dimensionless form:

where τ is the relaxation time and x₀ is the transverse scaling length. Equations (5), (6) and (7) form the basis of the model. By solving these equations we will be describing the beam propagation in both space and time. We develop a novel numerical procedure, utilizing our prior experience in treating the propagation of CP beams in nonlinear media [13]. The temporal equation for the angle of reorientation (Eq. (7)) and the spatial propagation equations for the beam envelopes (Eqs. (5) and (6)) are solved together in a system of nested loops, including an iterated convergence loop. Starting from a given distribution of θ (t; x, y, z), the incident fields A and B are propagated (from their entrance sides) along z, to obtain a new distribution of the fields. Then the distribution of θ is integrated for a time step at each spatial point. The fields now do not correspond to the new distribution of θ, and have to be propagated again along z. This two-step procedure is iterated until stable self-consistent distributions of the reorientation angle and the beam envelopes are obtained, which means that θ, A and B are advanced for a time step dt. The convergence in the self-consistency loop signifies that a solution at the time t+dt is found, and a new temporal step can start.

3. Results and discussion

Numerical studies of Eqs. (5)–(7), for CP beams propagating in electrically biased plane-oriented NLC, are performed in different conditions and for different beam configurations. The geometries are similar to the ones used in experiments [5–10]. The propagation of narrow and broad Gaussian beams, as well as vortices, is investigated. Even though the planar NLC are often made thin, we consider the propagation in the bulk. All the pictures and the movies in the transverse (x, y) plane are presented at the exit face of the crystal (z=L), and those in the (y, z) plane are at the x=0 plane (in the middle of the crystal). All the pictures and the movies are presented for one beam (forward) only. In all the figures the times are indicated, except where the steady state is displayed.

In all the 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, the birefringence ε_a=0.5. All of these data are consistent with the values reported in the experimental investigations [5–8].

First we consider the behavior of CP Gaussian beams in NLC, as the beam intensity is increased. The effect of the input intensity variation on the CP Gaussian beams propagation is presented in Fig. 1. For smaller intensities (Fig. 1(a)) self-focusing is too weak to keep the beam tightly focused, so it is not passing through unchanged, as a soliton. By increasing the beam intensity (Fig. 1(b)), at a certain value a stable solitonic propagation is achieved. For still higher intensities we observe transversal motion of the beam, in the form of transverse jumps [13]. One (Fig. 1(c)), and two consecutive jumps (Fig. 1(d)) are noted in Fig. 1. This beam displacement is due to a force exerted on the beam, caused by the refractive index change in the medium, and originating from the presence of the other beam. Such transverse motion of beams corresponds to the small-amplitude undulations [5–10], observed in experiments. For further increase in the intensity we see unstable dynamical behavior, Fig. 1(e), in the form of convective instabilities.

 

Fig. 1. The forward beam propagation, shown in the (y, z) plane (the first and third rows) and in the (x, y) plane (the second row), for different input intensities. a), f) and k) Diffracting beam, I=6×10⁺⁹ V²/m². b), g) and l) Soliton, I=7×10⁺⁹ V²/m². c), h) and m) One beam shift, I=8×10⁺⁹ V²/m². d) (632 KB), i) (334 KB) and n) (1.381 MB) Two beam shifts, I=1×10⁺¹⁰ V²/m². e) (131 KB), j) (229 KB) and o) (760 KB) Convective instability, I=9×10⁺¹⁰ V²/m². The first and second rows depict intensity distributions, the third row θ(y,z). For all the simulations FWHM=4 µm, L=0.5 mm and ε_a=0.5.

Download Full Size | PDF

In Fig. 2 we present examples of intensity distributions from Fig. 1, visualized by iso-surfaces in 3D, for different input beam intensities, for which stable solitonic propagation, undulation, filamentation and irregular dynamics occur. Several isosurfaces of the same beam are superimposed onto the same picture, represented by different colors. Figure 2(c) clearly displays the convective nature of the modulational instability occurring in the system.

 

Fig. 2. Beam propagation, shown in 3D, for different input intensities: a) I=7×10⁺⁹ V²/m² at t=15 τ, using two intensity levels (0.1 and 0.2 of I); b) I=1×10⁺¹⁰ V²/m² at t=3.5 τ and t=4.8 τ, the same intensity levels.; c) I=9×10⁺¹⁰ V²/m² at t=1 τ, using eight intensity levels (in relative units: 0.0001, 0.001, 0.01,0.1, 0.5, 1.0, 3.0 and 5.0). Parameters: FWHM=4 µm, L=0.5 mm and ε_a=0.5.

Download Full Size | PDF

In Fig. 3 we compare the behavior of copropagating (CO) and CP beams. In the CO case we launch two overlapping beams from the same side of the cell, and look for the formation of vector solitons. The transversal motion of the beams found in CP geometry has no counterpart in the CO geometry. For narrow CP beams (4 µm), the intensity needed for soliton existence (to pass through the medium without diffraction) is about three times lower than in the CO case [11]. This is not difficult to understand: each of the CP beams produces a larger change of the refractive index and contributes more to beam focusing in its own half of NLC (where it enters). The change in the index needed for soliton formation is caused by total intensity, hence less intensity of each beam is required. The same conclusion holds for broader beams, with the input FWHM=20 µm (Fig. 3). As expected, the same (or similar) behavior occurs for CO and CP cases, but at different input beam intensities. For the CP beams, the soliton existence region is wider than for the CO beams. For higher intensities we again observe different kinds of instabilities: soliton breathing (the second column), stable filamentation (the third and fourth columns), and long-lasting convective instabilities (the fifth column).

 

Fig. 3. Comparison between the copropagating (the first row (708 KB, 591 KB, 296 KB)) and the counterpropagating (the second row (521 KB, 406 KB, 184 KB)) beams in NLC, for different input intensities, indicated in the figures. For all the simulations input FWHM=20 µm, L=0.5 mm and ε_a=0.5.

Download Full Size | PDF

 

Fig. 4. Intensity distributions at the output face of the crystal for two values of intensities: I=1×10⁺⁹ V²/m² (the first row (71 KB, 1.428 MB, 1.303 MB, 1.749 MB)) and I=7×10⁺⁹ V²/m² (the second row (116 KB, 108 KB, 873 KB, 1.123 MB, 1.290 MB)), for different input FWHM of beams (indicated in each figure). In all the simulations L=0.5 mm and ε_a=0.5.

Download Full Size | PDF

Further, we consider the behavior of broader CP Gaussian beams in NLC (Fig. 4). We utilize broader Gaussian beams to display modulational instabilities and pattern formation of CP beams. We vary the input FWHM of Gaussian beams, for two different values of the input beam intensity. Interestingly, by increasing FWHM we see more regular patterns developing after an irregular transient dynamical phase, unlike the case of CP broad Gaussians in photorefractive crystals [13]. There one sees ordered patterns first, which become increasingly irregular. This behavior is the consequence of the long-range nonlocal interactions that are characteristic of NLC, in contrast to the local interactions, characteristic of the isotropic photorefractive medium. Nonlocal interactions tend to suppress modulational instabilities of broader beams and produce smoother patterns. Still, the transverse symmetry of developing initially axially symmetric beams need not remain axially symmetric. Pattern formation often involves symmetry-breaking mechanisms, and the resulting patterns could be of a different symmetry. Another feature of broad CP beams in NLC is relatively long time needed to achieve the steady state.

Besides varying beam parameters, we also vary the birefringence ε_a, which characterizes the strength of the nonlinearity. In Figure 5 we show the propagation for ε_a=0.8, analogous to Fig. 1. One can see that the instabilities develop similarly to the case ε_a=0.5, but at different, and smaller, input intensities. In particular, the values of intensity where the stable soliton propagation is observed are different for the two cases. For higher ε_a this value is smaller, as expected (for ε_a=0.8 the threshold intensity for the stable soliton existence is 2.7×10⁺⁹ V²/m², whereas for ε_a=0.5 the soliton appears at I=7×10⁺⁹ V²/m²).

 

Fig. 5. The intensity distributions for ε_a=0.8 and for different input intensities: a) I=2×10⁺⁹ V²/m², b) I=2.7×10⁺⁹ V²/m², c) I=5×10⁺⁹ V²/m² (757 KB), d) I=7×10⁺⁹ V²/m² (518 KB). In all simulations FWHM=4 µm and L=0.5 mm.

Download Full Size | PDF

Finally, we consider the interaction of CP vortices with opposite topological charges (±1) in NLC. Vortex beams carry angular momentum, proportional to the charge, and as a rule provide for interesting transverse dynamics. The input beam intensities are varied for two values of the input FWHM (Fig. 6). We find stable vortex propagation for lower intensities (first and forth columns), as well as standing filaments in the form of dipoles and quadrupoles for higher intensities (second and fifth columns). By further increasing the input beam intensities, an irregular long-lasting behavior is observed (third and sixth columns).

 

Fig. 6. CP vortices in the (y, z) plane (the first row) and in the (x, y) plane (the second row), for different input FWHM: 8 µm (first, second and third columns (125 KB, 82 KB)) and 23 µm (fourth, fifth and sixth columns (94 KB, 169 KB)), and for different intensities, noted in each figure. The first and fourth columns depict stable vortex propagation. The second and fifth columns display the filamentation of vortices. The third and sixth columns display dynamical instabilities. In all the simulations L=0.5 mm and ε_a=0.5.

Download Full Size | PDF

4. Conclusions

We report, for the first time to our knowledge, on the behavior of CP self-trapped beam structures in NLC. We demonstrate the existence of CP vector solitons in a narrow region of beam intensities. For narrow CP Gaussian beams we discover the transversal motion of the beams, corresponding to the undulation of the beams. At higher input intensities, we see filamentation and irregular dynamics of CP beams. The consideration of broader Gaussian beams offers opportunities for observing complex long-lasting dynamical behavior. By increasing FWHM of input beams regular patterns, after an irregular transient dynamical behavior, are displayed. For CP vortices with the same topological charge, stable vortex propagation, standing and rotating filaments, and an irregular behavior are observed, depending of the beam parameters.