Flowing Liquid Crystal Torons Around Obstacles

1. Introduction

Liquid crystals (LCs) represent a fascinating class of soft materials characterized by their fluidity combined with anisotropic properties, making them highly responsive to external fields and boundary conditions. This unique combination has led to widespread applications, particularly in display technologies, and has also spurred interest in fundamental research on the stability of their phases and topological structures [1]. Among these structures, solitons in chiral nematic LCs have emerged as a significant area of study due to their rich topological properties and potential for novel technological applications [2,3]. Besides skyrmions, which are two-dimensional configurations, many other topological-protected structures were realized in experiments, including torons, hopfions, skyrmions and solitonic macromolecules. We will refer to torons as the three-dimensional elementary skyrmion terminating at two-point defects to satisfy the uniform surface boundary conditions and match the topologically nontrivial skyrmion tube with the uniform far-field background and reserve skyrmion to describe the two-dimensional midplane structure.

Solitons in liquid crystals are localized, non-singular distortions of the director field that exhibit robust particle-like behaviour and complex interactions. These structures are of particular interest due to their ability to be manipulated by external stimuli, such as electric fields, which can drive their motion and induce collective behaviour [4,5]. The study of soliton dynamics in LCs has provided insights into fundamental aspects of topology in soft matter and has highlighted their potential in advanced material applications, including reconfigurable optical devices and microfluidic systems [4,6,7,8,9,10,11,12,13,13].

While much of the research on LC torons has focused on their behaviour under uniform external fields, little is known under non-uniform conditions such as static colloids or other obstacles trapped in the LC matrix. Furthermore, only a handful of studies addressed the interaction between torons and mass flows [14,15,16,17] remaining an overlooked problem of relevance in various fields. Understanding the dynamical response of LC torons in non-uniform or heterogeneous media is relevant for practical applications, where flowing torons may encounter various types of obstacles, such as geometrically patterned substrates, dispersed colloidal particles or defects in the LC matrix. This response will significantly alter the torons’s initial trajectory, in addition to its shape and stability, potentially driving novel dynamical phenomena.

In the theoretical study of torons, three approaches are commonly used in the literature. The first approach involves the numerical solution of continuum equations resulting from the minimization of the elastic free energy, which is typically implemented using finite difference or finite volume methods [4,18,19,20]. This approach has focused either on the dynamic evolution of torons or on the stable configuration of static torons. The second approach is based on analytical studies of the continuum equations, under simplified conditions required to make the problems tractable [21,22]. This method often provides valuable insights into the toron behaviour in idealized situations. The third approach, a more recent development, uses particle-based modelling, where torons are treated as individual particles [23]. Langevin dynamics is employed to simulate their motion, and potential functions for the toron interactions are inferred from experimental data or from simulations using the first method. In this study, we use the first approach combining it with the lattice Boltzmann method to simulate material flow, providing a comprehensive and versatile framework to study the interaction between topology and external mass flows that may be perturbed by an arbitrary array of obstacles.

In what follows, we study the flow of liquid crystal torons, driven by an external mass flow, around solid obstacles. The toron’s stability and trajectory are found to depend strongly on the impact parameter of the obstacle. When the impact parameter is less than $h<5\mu$m, the toron becomes unstable and disintegrates, in flows where torons would remain stable otherwise. However, when the impact parameter, is larger than $h\ge 5\mu$m, the stability of the toron is preserved while the toron’s trajectory is deflected by an amount that decays exponentially with the impact parameter. Notably, the toron’s velocity is found to increase as it approaches the obstacles due to flow constriction. We have also examined the dynamics and stability of torons flowing around one and between two obstacles, providing additional insights into how multiple obstacles may influence the dynamics of flowing torons.

2. Method

2.1. Hydrodinamic Equations

The dynamics of the liquid crystal (LC) director field are described by the Ericksen-Leslie model [24,25,26]. This model couples two dynamical equations: one for the material flow and another for the director field. These equations are particularly adequate to describe the behaviour of LCs in the nematic or cholesteric phases.

For the velocity field, we use the Navier-Stokes equation along with the continuity equation:

$$\begin{array}{ccc}& & \rho {\partial }_t{u}_{\alpha }+\rho u_{\beta }{\partial }_{\beta }{u}_{\alpha }={\partial }_{\beta }\left[-p{\delta }_{\alpha \beta }+{\sigma }_{\alpha \beta }^v+{\sigma }_{\alpha \beta }^e\right]+\rho g\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & {\partial }_{\alpha }{u}_{\alpha }=0,\hfill \end{array}$$

(2)

where the viscous stress tensor is defined as:

$$\begin{array}{ccc}\hfill {\sigma }_{\alpha \beta }^v=& & {\alpha }_1{n}_{\alpha }{n}_{\beta }{n}_{\mu }{n}_{\rho }{D}_{\mu \rho }+{\alpha }_2{n}_{\beta }{N}_{\alpha }+{\alpha }_3{n}_{\alpha }{N}_{\beta }\hfill \\ & & +{\alpha }_4{D}_{\alpha \beta }+{\alpha }_5{n}_{\beta }{n}_{\mu }{D}_{\mu \alpha }+{\alpha }_6{n}_{\alpha }{n}_{\mu }{D}_{\mu \beta }.\hfill \end{array}$$

(3)

In these equations, $\rho$ denotes the fluid density, g is the external acceleration that drives the fluid flow, p is the pressure, $\mathit{u}$ is the fluid velocity, $\mathit{n}$ is the director field describing the direction of alignment of the LC molecules, and ${\alpha }_n$’s are the Leslie viscosities. The kinematic transport, which describes the influence of the macroscopic flow field on the microscopic structure, is given by:

$$\begin{array}{c}\hfill N_{\beta }={\partial }_t{n}_{\beta }+u_{\gamma }{\partial }_{\gamma }{n}_{\beta }-W_{\beta \gamma }{n}_{\gamma }\end{array}$$

(4)

while the shear rate and vorticity tensors are defined as:

$$\begin{array}{c}\hfill D_{\alpha \mu }={\displaystyle \frac{1}{2}}\left({\partial }_{\alpha }{u}_{\mu }+{\partial }_{\mu }{u}_{\alpha }\right),W_{\alpha \mu }={\displaystyle \frac{1}{2}}\left({\partial }_{\mu }{u}_{\alpha }-{\partial }_{\alpha }{u}_{\mu }\right).\end{array}$$

(5)

The elastic stress tensor is:

$$\begin{array}{c}\hfill {\sigma }_{\alpha \beta }^e=-{\partial }_{\alpha }{n}_{\gamma }{\displaystyle \frac{\delta E}{\delta \left({\partial }_{\beta }{n}_{\gamma }\right)}},\end{array}$$

(6)

where E represents the Frank-Oseen elastic free energy:

$$\begin{array}{ccc}\hfill E=& & \int dV({\displaystyle \frac{K_{11}}{2}}{(\nabla ·\mathit{n})}^2+{\displaystyle \frac{K_{22}}{2}}{\left(\mathit{n}·[\nabla \times \mathit{n}]+q_0\right)}^2\hfill \end{array}$$

(7)

$$\begin{array}{ccc}& & +{\displaystyle \frac{K_{33}}{2}}{[\mathit{n}\times [\nabla \times \mathit{n}]]}^2).\hfill \end{array}$$

(8)

with $K_{11}$, $K_{22}$, $K_{33}$ the Frank elastic constants, and $q_0=2\pi /P$, with P the cholesteric pitch.

The second set of equations describes the evolution of the director field:

$$\begin{array}{cc}\hfill & {\partial }_t{n}_{\mu }={\displaystyle \frac{1}{{\gamma }_1}}{h}_{\mu }-{\displaystyle \frac{{\gamma }_2}{{\gamma }_1}}{n}_{\alpha }{D}_{\alpha \mu }-u_{\gamma }{\partial }_{\gamma }{n}_{\mu }+W_{\mu \gamma }{n}_{\gamma },\hfill \end{array}$$

(9)

where ${\gamma }_1={\alpha }_3-{\alpha }_2$ is the rotational viscosity determining the relaxation rate of the director, and ${\gamma }_2={\alpha }_3+{\alpha }_2$ is the torsion coefficient characterizing the contribution to the viscous torque from the velocity field gradients. The ratio ${\gamma }_2/{\gamma }_1$ is known as the alignment parameter, with $|{\gamma }_2/{\gamma }_1|>1$ indicating flow-aligning and $|{\gamma }_2/{\gamma }_1|<1$ flow-tumbling systems. The molecular field is given by

$$\begin{array}{c}\hfill h_{\mu }=-{\displaystyle \frac{\delta E}{\delta n_{\mu }}}.\end{array}$$

(10)

The simulations employed a hybrid numerical method. The velocity field is resolved using the lattice Boltzmann method [27,28], summarized in the next subsection, with the elastic and viscous stress tensors (except the term proportional to ${\alpha }_4$) introduced as force terms. The director field equation, Eq.(9), was solved using a predictor-corrector finite-difference algorithm. On solid boundaries, including the obstacles, infinite homeotropic anchoring and no-slip conditions were applied using the bounce-back condition [27].

The simulations started with the liquid at rest, and the directors mostly aligned perpendicular to the plates except near the toron, whose configuration was obtained by minimizing its free energy starting from the Ansatz of Ref. [14]. The material parameters were chosen to be close to those of MBBA at 22^∘C [29], except for the absolute viscosity (or equivalently, ${\alpha }_4$), which was doubled to ensure reasonable simulation times. The code was parallelized in CUDA-C, and the simulations were executed on GPUs [30]. A typical performance of the simulation is ∼328 MLUPS (Mega Lattice Updates Per Seconds) in a Nvidia Tesla V100 using double precision calculations.

2.2. Lattice Boltzmann Method

We obtain the velocity field by solving the discretized Lattice Boltzmann Equation (LBE) using the single relaxation time (SRT) collision approximation and Guo’s forcing scheme [31] as follows:

$$\begin{array}{cc}\hfill & f_i(\mathit{x}+{\mathit{e}}_i\Delta t,t+\Delta t)-f_i(\mathit{x},t)=\hfill \\ & -{\displaystyle \frac{\Delta t}{\tau }}(f_i(\mathit{x},t)-f_i^{\mathrm{eq}}(\mathit{x},t))+\Delta t\left(1-{\displaystyle \frac{1}{2\tau }}\right)S_i^F(\mathit{x},t),\hfill \end{array}$$

(11)

where $f_i$ is the particle distribution function, $f_i^{\mathrm{eq}}$ is the equilibrium distribution function, which is determined by:

$$\begin{array}{cc}\hfill f_i^{\mathrm{eq}}\left(\mathit{x}\right)& =w_i\rho \left(\mathit{x}\right)\left[1+{\displaystyle \frac{{\mathit{e}}_i·\mathit{u}\left(\mathit{x}\right)}{c_s^2}}\right.\hfill \\ & \left.+{\displaystyle \frac{{({\mathit{e}}_i·\mathit{u}\left(\mathit{x}\right))}^2}{2c_s^4}}-{\displaystyle \frac{\mathit{u}\left(\mathit{x}\right)·\mathit{u}\left(\mathit{x}\right)}{2c_s^2}}\right],\hfill \end{array}$$

(12)

where $c_s$ is the lattice speed of sound, defined as $c_s={\displaystyle \frac{\Delta x}{\sqrt{3}\Delta t}}$. $\Delta x$, $\Delta t$, and the reference density ${\rho }_0$ is set to unity (lattice units). $w_i$ and ${\mathit{e}}_i$ are the lattice weights and the discrete velocity vector in the $i^{\mathrm{th}}$ direction, respectively. In this work, we use the $D_3{Q}_{19}$ lattice arrangement, hence $w_i$ and ${\mathit{e}}_i$ are given by:

$${\mathit{e}}_i=\left(\begin{array}{ccccccccccccccccccc}0& 1& -1& 0& 0& 0& 0& 1& -1& 1& -1& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 1& -1& 0& 0& 1& 1& -1& -1& 0& 0& 0& 0& 1& -1& 1& -1\\ 0& 0& 0& 0& 0& 1& -1& 0& 0& 0& 0& 1& 1& -1& -1& 1& 1& -1& -1\end{array}\right),$$

(13)

$$w_i=\left(\begin{array}{ccccccccccccccccccc}{\displaystyle \frac{1}{3}}& {\displaystyle \frac{1}{18}}& {\displaystyle \frac{1}{18}}& {\displaystyle \frac{1}{18}}& {\displaystyle \frac{1}{18}}& {\displaystyle \frac{1}{18}}& {\displaystyle \frac{1}{18}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}& {\displaystyle \frac{1}{36}}\end{array}\right).$$

(14)

The macroscopic fluid density is defined as:

$$\rho \left(\mathit{x}\right)=\sum _i{f}_i\left(\mathit{x}\right).$$

(15)

The relaxation time ($\tau$) is related to the kinematic viscosity of the fluid ($\nu$) through:

$$\nu =c_s^2\left(\tau +{\displaystyle \frac{1}{2}}\right)\Delta t.$$

(16)

The forcing term in the LBE equation $S_i^F$ is defined as:

$$\begin{array}{cc}\hfill & S_i^F\left(\mathit{x}\right)=w_i\left[{\displaystyle \frac{{\mathit{e}}_i·\mathit{F}\left(\mathit{x}\right)}{c_s^2}}+\right.\hfill \\ & \left.{\displaystyle \frac{({\mathit{e}}_i{\mathit{e}}_i-c_s^2\mathit{I}):(\mathit{u}\left(\mathit{x}\right)\mathit{F}\left(\mathit{x}\right)+\mathit{F}\left(\mathit{x}\right)\mathit{u}\left(\mathit{x}\right))}{2c_s^4}}\right],\hfill \end{array}$$

(17)

where $\mathit{F}$ represents the resultant force acting on the flow field (external force and stresses arising from the liquid crystal distortions), and $\mathit{u}$ is the macroscopic fluid velocity, which is defined as:

$$\mathit{u}\left(\mathit{x}\right)={\displaystyle \frac{1}{\rho \left(\mathit{x}\right)}}\sum _i{\mathit{e}}_i{f}_i\left(\mathit{x}\right)+{\displaystyle \frac{\Delta t}{2\rho \left(\mathit{x}\right)}}\mathit{F}\left(\mathit{x}\right).$$

(18)

The bounce-back boundary condition is used to model the no-slip boundary condition over the solid nodes [27]. The idea behind the bounce-back boundary condition is to reflect the distribution functions hitting a solid node with the same magnitude in the reverse direction. The exact location of the solid boundary is assumed to be between two lattice points [27].

3. Results

We report the results of numerical simulations of a flowing toron driven by an external force, around a cylindrical pillar near its initial trajectory. The liquid crystal was initialized at rest, with alignment along the vertical direction, except at the toron, where we used the relaxed configuration obtained from the ansatz [14]. The system dimensions were set to $L_x\times L_y\times L_z=96\times 56\times 16$, and the material parameters are listed in Table 1. The pillar, with a radius of $R=3.13,\mu$m (or 5 $\Delta x$), homeotropic anchoring and no-slip boundary was placed at $x=L_x/2$ and $y=L_y/2+h$. The fluid was driven by an external body force of magnitude $g=312.5$ m/s² (or $2\times {10}^{-9}$ l.u.) in the x direction.

As expected, the pillar distorts the velocity field and perturbs the toron initial trajectory, as shown in Figure 1. In the midplane, Figure 1(a), the toron shrinks slightly as it flows around the pillar but later recovers its original configuration. To confirm the toron’s stability post-pillar interaction, we relaxed the final configuration in Figure 1 without the external force, confirming that it remains stable (not shown). The toron cross-section in Figure 1(b) reveals that the toron becomes asymmetric in the direction of the flow, consistent with previous observations in Poiseuille-like flows without the pillar [16]. The velocity field is perturbed by the toron, while the pillar introduces large velocity gradients in the flow field (Figure 1(c)). In the cross-sectional view, Figure 1(d), the flow exhibits a Poiseuille-like profile, except near the toron and the pillar.

The deflection of the initial trajectory of the toron depends on the impact parameter h, as shown in Figure 2. Trajectories are defected significantly when the impact parameter is less than cholesteric pitch. However, below a threshold, $h=3.75,\mu$m, or $0.43$P, in our simulations, the toron becomes unstable and disintegrates. The toron in Poiseuille flows already becomes unstable without the obstacle at higher speeds [17] and here it occurs due to the local increase in the velocity field around the obstacle. When the toron collides directly with the pillar ($h=0$), it disintegrates before flowing around the obstacle. For small but non-zero values of h, the toron flows around the pillar but later becomes unstable, as depicted in Figure 3. The inset on the left of Figure 2 depicts the deflection $\Delta y$ of the toron initial trajectory as a function of the impact parameter, h. As h increases, $\Delta y$ decays exponentially, reaching zero when the obstacle is sufficiently far from the toron so that the toron does not feel the change in the director field caused by the obstacle. The exponential decay is usually observed in the distortions in the liquid crystal orientation [32]. The inset on the right shows the toron velocity as a function of position, where we note that a stable toron accelerates near the pillar, with its velocity increasing by about $26\%$ before returning to its initial value, due to flow constriction by the pillar.

The free energy integrated over the entire volume, calculated using Eq.(7), is plotted in Figure 4 as a function of time for various impact parameters, h. At small h, the free energy rises sharply, reaching a peak when the toron loses stability and disintegrates. At larger impact parameters, h, the toron remains stable, and the free energy stabilizes, with slight perturbations occurring as the toron flows around the pillar, particularly evident in the orange curve ($h=5\mu$m) around $t=1$ s.

Until now, we considered the flow of a toron around a single pillar. To explore the effects of multiple obstacles, we simulated the toron flow between two pillars, which act as a flow constriction (see Figure 5). In this configuration, the toron accelerates by 44% as it flows between the pillars, as shown in Figure 5(c). Similar to the single-pillar case, the toron becomes unstable at small values of h. In our simulations, the toron remains stable for $h>6.25,\mu$m, which is almost twice the impact parameter required for toron stability in the single-pillar case. This is likely due to the toron’s size reduction when passing around the pillars, an effect that is amplified by the presence of two obstacles. Additionally, in this configuration, the toron does not deviate from its initial trajectory.

In addition to simulating the interaction between the toron and a cylindrical pillar with homeotropic anchoring, we investigated other obstacle geometries (not shown). For example, we replaced the pillar with a solid sphere, mimicking a colloid, and applied different anchoring conditions (both homeotropic and toron-like). However, no qualitative differences were observed. Thus, we conclude that the primary factor determining the flow of torons around obstacles, as well as their stability, appears to be the effect of the large velocity gradients introduced by the obstacles rather than their detailed form.

4. Conclusion

In this study, we investigated the interaction between flowing liquid crystal torons and solid obstacles. Through numerical simulations using a combination of the lattice Boltzmann method and finite differences, we demonstrated how the dynamics of flowing torons respond to the presence of obstacles in the flow domain. Our findings reveal that the toron’s stability and trajectory depend strongly on the impact parameter of the obstacle. At short impact parameters (smaller than half cholesteric pitch), the toron becomes unstable and disintegrates, while at larger impact parameters, the toron’s trajectory is deflected around the obstacle. The deflection of the trajectory and the velocity increase as the toron approaches the obstacle depend strongly on the impact parameter. In addition, we explored the dynamics of torons flowing between two obstacles, providing insights into how multiple obstructions influence the toron´s stability and motion. This study contributes to our understanding of toron interactions in complex environments which could be used in applications in reconfigurable optical devices and microfluidic systems.