Finding exact analytical soliton profile solutions is only possible for certain types of non-linear media. In most cases one must resort to numerical techniques to find the soliton profile. In this work we present numerical calculations of spatial soliton profiles in nematic liquid crystals. The nonlinearity is governed by the optical-field-induced liquid crystal director reorientation, which is described by a system of coupled nonlinear partial differential equations. The soliton profile is found using an iterative scheme whereby the induced waveguide and mode profiles are calculated alternatively until convergence is achieved. In this way it is also possible to find higher order solitons. The results in this work can be used to accurately design all-optical interconnections with soliton beams.
©2010 Optical Society of America
Spatial optical solitons of the bright (dark) type occur in a self-(de)focusing medium when the nonlinear self-focusing exactly balances the natural diffraction of the beam . Another point of view is that the optical beam creates a self-induced waveguide. In order to be a soliton, the beam profile must be equal to the mode of its induced waveguide [2, 3]. For certain types of nonlinearity, in which the optical field and the refractive index are straightforwardly related, it is possible to find analytical solutions. For example, an analytical solution is presented in  for a threshold type nonlinearity. When the nonlinearity is more complicated, one must resort to numerical techniques to find the soliton solution. Calculations of soliton mode profiles in nonlocal nonlinear media have been presented in a number of publications, such as [5–7].
Liquid crystals (LCs) prove to be an ideal testbed for nonlinear optical phenomena. Molecules in a nematic liquid crystal have no positional order but do posses an orientational degree of order. This leads to an anisotropy in the macroscopic properties of the material. In the continuum theory, the material properties (optical and electrical) are determined by the average orientation of the molecules, denoted as the director n̄. The material exhibits a number of optical nonlinear effects  of which the optical induced director orientation has been investigated extensively over the last decade. Due to the dielectric anisotropy at low frequencies Δεs = ε ∥ - ε ⊥ the biasing electric field causes a torque which tends to align the director along the electric field lines (for positive dielectric anisotropy). In a similar way an optical electric field will also cause a torque which is proportional to Δε = n 2 ∥ - n 2 ⊥, with Δn = n ∥ - n ⊥ being the birefringence. Due to the reorientation of the director, a light beam will normally experience an increase in refractive index with increasing intensity. This results in a self-focusing mechanism for optical beams which can be used to generate spatial optical solitons in a number of configurations, either with [9, 10] or without [11, 12] bias voltage. Often these spatial solitons in nematic liquid crystals are referred to as nematicons. Both experiments [10, 13] and theoretical calculations [14, 15] have demonstrated that the nonlinearity is highly nonlocal and the nonlocality depends on the thickness of the liquid crystal layer. Snyder and Mitchell have demonstrated  that spatial solitons in highly nonlocal media with a parabolic nonlinear response have a Gaussian profile. These solitons are highly stable which means that a deviation from the exact soliton profile or optical power will not lead to loss of confinement, but merely results in breathing of the beam. This is why they are called accessible solitons in . Liquid crystal orientational nonlinearity is highly nonlocal, but does not exhibit a perfect parabolic response. Therefore it is sufficient to launch a Gaussian beam into the cell to generate a soliton experimentally. However, it is only possible to observe breathing solitons as demonstrated by Conti in  and recently by Striníc et al. in , since this Gaussian profile does not match the soliton profile exactly. Breathing means that the width of the beam varies in a periodic way along the propagation direction. In literature, many articles can be found in which soliton profiles are calculated based on (semi-)analytical models [19–22]. The calculation of the nonlinear response of the liquid crystal is governed by a number of partial differential equations and previously mentioned articles all start from a simplified model, which in our opinion cannot describe accurately the whole reorientation dynamics. Moreover also the effect of the longitudinal anisotropy component is neglected, which leads to walk-off of the beam . In  the effect of walk-off is taken into account in the accurate description of the reorientational nonlinearity and soliton profiles.
In this work we will use an iterative numerical technique to find the exact soliton profiles in nematic liquid crystals. Section 2 describes the liquid crystal model and mode calculation model. Section 3 then describes some calculation examples of the liquid crystal behavior, while Section 4 shows some mode calculation examples. Section 5 then deals with the calculation of the exact soliton profiles for both zero and higher order soliton modes.
2. Numerical simulation methods
2.1. Liquid crystal simulation
One part of the simulation program consists of a finite element program to model the behavior of the liquid crystal. The model is able to simulate variations in order parameter and instead of working with a vector that describes the orientation of the director with fixed order parameter, the model works with the Q tensor. This tensor contains information on both the orientation and the order parameter at a certain position in the liquid crystal. The simulation model is based on the minimization of the Landau-de-Gennes free energy functional , defined as
Here, fD is the distortion free-energy density which is of the form
With the definition of the Q tensor that we used, the L coefficients are related to the elastic constants as , , , and . The bulk free-energy density is expanded in a power series near the transition temperature and contains the thermotropic constants A, B and C:
The electrostatic energy density takes the form
The optical electric energy density is of the same form, but it is important to take into consideration the complex nature of the optical electric field as the simulations are in the steady state regime (frequency domain).
This term accounts for the influence of the optical fields, which is responsible for the optical nonlinearity. The Q tensor model  allows for the simulation of order variations of the liquid crystal. For a uniaxial material the relation between the Q tensor and the director n̄ may be defined as Q = S(3n̄ ⊗ n̄ - I)/2, with S the local order parameter. Approximately, the dielectric and optical tensor can be described by the following equation: , in which S 0 is the equilibrium order parameter. Our description of the LC behavior is based on a very general model, because it incorporates the whole reorienation and order parameter calculation and is applicable to a wide range of configurations due to the versatility of the finite element implementation. One source of errors is the fact that the bulk energy density is truncated to include at most fourth order terms. Otherwise the accuracy of the modeling is only limited to the number of elements in the mesh.
Previous publications on solitons in LCs have assumed the order parameter constant [9, 14, 18–22]. Such an assumption is valid for uniaxial arrangements of molecules only, which may be adequately represented by a director field. However it is known that the order parameter of the liquid crystal can change drastically when rapid variations in orientation occur, when strong electric fields are present  or under the influence of surfaces [25, 28]. Therefore it is interesting to investigate to which extent a strong optical field (as in the case of spatial solitons) influences the order parameter.
2.2. Finite element anisotropic mode solver
In order to calculate the optical modes of the induced waveguide, a full-vectorial finite element mode solver has been used which can handle the full anisotropy of the dielectric tensor . The mode solver is based on the solution of the variational form of the curl-curl equations of the electric field
It has been initially developed for the analysis of optical waveguides, either in liquid crystal or with a liquid crystal cladding. The finite element mesh is different for the LC calculation and the mode calculation to allow for better tuning of the density of the mesh and to minimize the calculation time. The number of elements in the mode calculation mesh is more than 2000, which leads to an error on the effective index smaller than 10-5 as was shown in .
3. Liquid crystal behavior
3.1. Configuration with bias voltage
Figure 1 describes the first configuration that is investigated in this work. The LC simulation window is 150 μm × 50 μm and the optical beam is launched along the z axis in the middle of the simulation window. The beam is a circular Gaussian beam with a waist of 3 μm. At the top and bottom of the simulation window the liquid crystal orientation is fixed at a pretilt angle of 2° (θ = 2°, φ = 90°). A voltage of 1 V is applied across the layer so the dominant electric field component lies along x. The parameters of the liquid crystal E7 are used in accord with our previous publications .
Figure 2 shows the x and y component of the director, nx and ny respectively. It is clear that the component along y is two orders of magnitude smaller than the x component. This means that the elements εxy and εyz in the optical tensor ε̿ are almost zero and the element εyy exhibits only very small variations. Therefore the induced waveguide will not support any modes which are y polarized. The mode calculations that follow demonstrate that the x polarized modes have only a small y component.
The variation of the order parameter is not shown here, but calculations have revealed that the variations in order parameter are small for the optical power densities used in these simulations. The Q tensor method that is used for the simulations can only predict order variations accurately for temperatures close to the supercooling temperature. For a temperature of 2°C below the transition temperature a maximum change in order parameter of 3 ∙ 10-6 was found for an optical power of 5 mW. Increasing the optical power up to 50 mW leads to a maximum change of 5 ∙ 10-5. The change in the permittivity tensor can roughly be described by which leads us to a maximum variation of 3 ∙ 10-5 for 50 mW optical power. These variations are 3 orders of magnitude smaller than the variations due to the optical director orientation and can thus be neglected in the further calculations.
3.2. Twist configuration
The second configuration that is investigated in this work is the twisted nematic configuration, similar to the configuration used in . The simulation window is now 25 μm by 75 μm and the director is parallel to the glass plates at top and bottom and twists over 180° from bottom to top (θ= 0°, φ(x = 0μm) = −90°, φ(x = 25μm) = 90°). A Gaussian beam is injected in the middle of the simulation window with a polarization along the y direction and with an optical power of 56 mW. Figure 3 shows the result of the calculation in terms of the y component of the director, together with the resulting refractive index profile. The optical beam creates only very minor distortions in the orientation compared to the zero optical power situation. In fact, the refractive index at the beam center does not increase due to reorientation, it is only at the sides that the refractive index increases. In other words, the width of the higher index region increases. In this case the effect of reorientation on the refractive index is much smaller than in the biased configuration, which also means that variations of order parameter may play a role, at least for temperatures close to the nematic to isotropic transition temperature (which is important for liquid crystals with a low transition temperature such as 5CB).
4. Mode calculations
4.1. Biased configuration
From experimental data  and numerical calculations based on beam propagation methods  it is known that for an applied voltage of 1.0 V an optical power of roughly 4 mW is required for spatial soliton-like behavior. We have calculated the modes of a waveguide induced by a 3.5 mW x polarized Gaussian beam. Figure 4 shows the electric field components of the fundamental mode solution. The mode is mainly x polarized. The contour plots of the first 8 modes (i.e. the ones with the largest effective index) are shown in Fig. 5. All the modes are mainly x polarized due to the anisotropy of the induced waveguide. Moreover, due to the εxz terms the optical field exhibits a variation in phase along the x axis. Due to this phase variation, the beam will enter the cell without walk-off angle. Without this phase variation a beam incident in the cell exhibits a transverse shift. Therefore, it propagates at an angle through the cell until it reaches a certain height and is bent back into the bulk by the gradient in refractive index profile. This results in a sinusoidal-like undulation throughout the cell [23, 31]. By launching the beam under a certain angle – or equivalently by introducing a phase variation along the x direction – undulations can be avoided.
Another interesting observation is that the induced waveguide is highly multimode. Only 8 modes are shown, but our simulations reveal that the number of guided modes is larger than 100. The fact that the induced waveguide is highly multimode was shown by a number of experimental results in the past. In 2001 it was shown experimentally by Peccianti et al.  that incoherent light can be guided in a self-induced waveguide. Confinement of incoherent light is a sign of a highly multimodal waveguide. Additionally, in  it has been demonstrated experimentally that a first order soliton can be generated in a liquid crystal. On the other hand, from the comparison with the harmonic oscillator, it is obvious that the waveguide is multimode . The investigation of higher order solitons in highly nonlocal media has been presented in .
4.2. Twist configuration
Mode calculations for the twisted configuration reveal that, in contrast to the biased configuration, only one mode mode is guided (in two dimensions) for powers up to 100 mW. Next to that, the refractive index profile of Fig. 3 also reveals that the nonlinearity is now only slightly nonlocal.
These observations show that properties of nematicons are very different depending on the configuration used, spanning from highly nonlocal highly multimode self-induced waveguides, to single mode waveguides with small nonlocality.
5. Soliton calculations
5.1. Fundamental mode
In order to find the zero order soliton beam profile, the following iterative procedure is followed. First, a Gaussian beam with a waist of 3 μm is used to calculate the director distributions for different optical input powers. Next, the optical modes are calculated for the different director orientations that have been obtained. From these mode calculations only the zero order mode is of interest. The mode profiles for different optical powers are then compared to the initial input profile by calculating the normalized covariance along x and y (denoted as Cx and Cy). The mode profile for the optical power that results in the largest covariance value is then used as an input for the second iteration and the procedure is repeated in further iterations.
5.1.1. First iteration
Figure 6 compares the mode profiles for different optical powers with the input profile. It is clear that the shape of the mode profile along x and y is very similar to the input Gaussian profile, which is due to the high nonlocality of the nonlinearity. Although similar, it is not exactly Gaussian because the nonlocality is not perfectly parabolic . The width of the profile decreases for increasing optical input power. The reduction is slightly more pronounced along y than x.
The correlation coefficients in Fig. 7 are actually shown as 1/(1 − Cx) and 1/(1 − Cy) as this visualizes better the correlation between the profiles. It is clear that the correlation reaches a maximum value for a certain optical power. However, the ideal optical power is different along x and y direction. The mode is thus not perfectly circular, but slightly elliptical. The ellipticity arises from the fact that the configuration is different along the x and the y direction. Along the x direction the LC is limited by the boundaries, while the configuration is much larger along the y direction (modeled by periodic boundary conditions). In thinner cells, this ellipticity will even be more pronounced as observed in simulation results (not shown in this manuscript).
5.1.2. Second iteration
For the input field of the second iteration, we chose to take the mode profile for 4.26 mW because this resulted in the best correlation along the y direction in the first iteration (Fig. 7). With this input profile, the director profile was calculated for different optical powers, together with the resulting mode profiles. The resulting correlation with the input profile is shown in Fig. 8. Remarkably the maximum correlation now appears at roughly the same power level for both x and y directions. Furthermore, the maximum correlation is strongly increased compared to the first iteration. These results show that two iteration steps are sufficient to find a self-consistent soliton solution. Further iterations increase the correlation but ultimately the correlation is limited by the number of datapoints.
The calculation gives a field profile that is invariant along the propagation direction for a particular value of the optical power. This is however not the only zero order soliton solution that can be found because there is a family of zero order solitons for different incident beamwidths, basically following the relation P ~ 1/w 4 0 which is valid for accessible solitons . In this equation P is the critical soliton power and w 0 is the beam waist. Our algorithm is suitable for finding these different soliton solutions and this can be achieved by using Gaussian input profiles with different waist.
5.2. First order soliton
The algorithm must be adapted slightly to find higher order soliton solutions. The initial input field is still the same Gaussian profile, but in the further calculation steps, the profile of the nth mode is used as input.
In order to demonstrate that higher order mode solutions can be found with our algorithm, we will search for the soliton solution with two lobes along the y direction, similar to the second mode profile shown in Fig. 5. The mode profiles along the y direction for the first iteration are shown in Fig. 9, together with the correlation between these profiles and the input profile. Again the maximum correlation occurs for a different optical power along the x and y direction.
For the second iteration step, the mode profiles along the x and y direction are shown in Fig. 10, while the correlation with the input profile is shown in Fig. 11. Similar to the result for the fundamental mode, the correlation factor improves compared to the first iteration and the maximum correlation for x and y occurs at the same optical power. With this in mind we can state that we have found a first order soliton solution. The required optical power is 7 mW, which is higher than the required power for the fundamental soliton solution.
We have presented calculations of liquid crystal behavior in the presence of strong optical fields with a numerical model that incorporates the full description of orientation and order parameter. With an iterative scheme, we have numerically determined the spatial soliton profiles in a classic nematic liquid crystal soltion geometry for the fundamental and first order mode. The fundamental solution is not purely Gaussian, because a Gaussian input beam always leads to breathing solitons. Our technique is suitable to find all higher order soliton solutions. Moreover, for each soliton solution, a series of solutions can be found with different width. On the other hand, our numerical model is also able to find soliton solutions in other configurations, such as the twisted nematic geometry. The results in this work can be used to design all-optical interconnections with soliton beams in a more accurate way.
Jeroen Beeckman is postdoctoral fellow of the Research Foundation - Flanders (FWO) and Pieter Vanbrabant is PhD Fellow of the same institution. The project is a result of collaboration within the framework of the Photonics@be program of the Belgian Science Policy.
References and links
1. Y. Kivshar and G. Agrawal, Optical Solitons - From Fibers to Photonic Crystals (Academic Press, San Diego, 2003).
2. A. Snyder, D. Mitchell, and Y. Kivshar, “Unification of linear and nonlinear-wave optics,” Mod. Phys. Lett. B 9, 1479–1506 (1995). [CrossRef]
3. M. Mitchell, M. Segev, T. Coskun, and D. Christodoulides, “Theory of self-trapped spatially incoherent light beams,” Phys. Rev. Lett. 79, 4990–4993 (1997). [CrossRef]
6. C. Rotschild, O. Cohen, O. Manela, and M. Segev, “Solitons in nonlinear media with an infinite range of non-locality: First observation of coherent elliptic solitons and of vortex-ring solitons,” Phys. Rev. Lett. 95, 213904 (2005). [CrossRef] [PubMed]
8. I. Khoo, Liquid Crystals: Physical Properties and Nonlinear Optical Phenomena (Wiley-Interscience, New York, 1994).
9. 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]
10. J. Henninot, J. Blach, and M. Warenghem, “Experimental study of the nonlocality of spatial optical solitons excited in nematic liquid crystal,” J. Opt. A: Pure Appl. Opt. 9, 20–25 (2007). [CrossRef]
11. K. Jaworowicz, K. A. Brzdakiewicz, M. A. Karpierz, and M. Sierakowski, “Spatial solitons in twisted nematic layer,” Mol. Cryst. Liq. Cryst. 453, 301–307 (2006). [CrossRef]
12. M. Peccianti, C. Conti, G. Assanto, A. De Luca, and C. Umeton, “Routing of Anisotropic Spatial Solitons and Modulational Instability in Liquid Crystals,” Nature 432, 733–737 (2004). [CrossRef] [PubMed]
13. 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]
14. 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). [CrossRef] [PubMed]
16. A. Snyder and D. Mitchell, “Accessible Solitons,” Science 276, 1538–1541 (1997). [CrossRef]
20. A. Minzoni, N. Smyth, and A. Worthy, “Modulation solutions for nematicon propagation in nonlocal liquid crystals,” J. Opt. Soc. Am. B 24, 1549–1556 (2007). [CrossRef]
21. H. Ren, S. Ouyang, Q. Guo, W. Hu, and C. Longgui, “A perturbed (1+2)-dimensional soliton solution in nematic liquid crystals,” J. Opt. A: Pure Appl. Opt. 10, 025102 (2008). [CrossRef]
22. H. Zhang, D. Xu, and L. Li, “An approximate solution for describing a fundamental soliton in nonlocal nonlinear media,” J. Opt. A: Pure Appl. Opt. 11, 125203 (2009). [CrossRef]
24. C. Conti, M. Peccianti, and G. Assanto, “Spatial solitons and modulational instability in the precense of large birefringence: The case of highly nonlocal liquid crystals,” Phys. Rev. E 72, 066614 (2005). [CrossRef]
25. R. James, E. Willman, F. A. Fernández, and S. E. Day, “Finite-Element Modeling of Liquid-Crystal Hydrodynamics With a Variable Degree of Order,” IEEE T. Electron Dev. 53, 1575–1582 (2006). [CrossRef]
26. P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, International Series of Monographs on Physics (Oxford University Press, Oxford, 1995).
27. R. Barberi, F. Ciuchi, G. Durand, M. Iovane, D. Sikharulidze, A. Sonnet, and E. Virga, “Electric Field Induced Order Reconstruction in a Nematic Cell,” Eur. Phys. J. E Soft Matter 13, 61–71 (2004). [CrossRef] [PubMed]
29. J. Beeckman, R. James, F. Fernandez, W. De Cort, P. Vanbrabant, and K. Neyts, “Calculation of Fully Anisotropic Liquid Crystal Waveguide Modes,” J. Lightw. Technol. 27, 3812–3819 (2009). [CrossRef]
30. U. Laudyn, M. Kwasny, and M. Karpierz, “Nematicons in chiral nematic liquid crystals,” Appl. Phys. Lett. 94, 091110 (2009). [CrossRef]
31. J. Beeckman, K. Neyts, X. Hutsebaut, C. Cambournac, and M. Haelterman, “Simulation of 2-D Lateral Light Propagation in Nematic-liquid-crystal Cells with Tilted Molecules and Nonlinear Reorientational Effect,” Opt. Quantum Electron. 37, 95–106 (2005). [CrossRef]
32. M. Peccianti and G. Assanto, “Incoherent Spatial Solitary Waves in Nematic Liquid Crystals,” Opt. Lett. 15, 1791–1793 (2001). [CrossRef]
33. X. Hutsebaut, C. Cambournac, M. Haelterman, A. Adamski, and K. Neyts, “Single-component higher-order mode solitons in liquid crystals,” Opt. Commun. 233, 211217 (2004). [CrossRef]