Class: Nematic Liquid Crystal
A liquid crystal is a state of matter between a solid and a liquid. It is characterized by the orientation of the molecules, which is ordered, but not the position. The molecules are rodlike and the orientation is described by a unit vector, the nematic director. The nematic liquid crystal is the simplest form of liquid crystal, and is characterized by the nematic director being the only order parameter. The nematic liquid crystal is used to describe the behavior of many biological systems, such as the cytoskeleton, and is also used in many technological applications, such as in liquid crystal displays.
In this class, we simulate an active nematic liquid crystal using [framework].
Variables and parameters
The primary variables are the symmetric traceless tensor \(Q\) and the velocity field \(\\mathbf{u}\)
The NematicLiquidCrystal class takes the same keyword as the BaseSystem class in addition to
Keyword  Definition  Default value 

alpha 
The activity parameter. Negative is extensile  \(1\) 
K 
Frank elastic constant  \(1\) 
A 
Parameter in front of \(\text{Tr}(Q^2)\) in the freeenergy  \(1\) 
B 
Parameter in front of \(\text{Tr}(Q^2)^2\) in the freeenergy  \(1\) 
C 
Parameter in front of \(\text{Tr}(Q^3)\) in the freeenergy. 3D only  \(0\) 
gamma 
Rotational diffusion  \(1\) 
Gamma 
Friction  \(0\) 
eta 
Viscosity  \(1\) 
These parameters are discussed in more detail in the model section.
Note on the tensor parameters
The \(Q\) tensor is given by the nematic director as
in \(d\) dimensions. To take advantage of its symmetric nature we have saved \(Q\) as a vector field, which in two and three dimensions takes the forms
respectivly. We can translate from tensor indexes to the right value stored in the vector by using the function
This returns the element \(Q_{ij}\) of a symmetric and traceless tensor field. In addition to this we also have the function
so that we can optimally store the antisymetric tensors as well. In two dimensions these only have one independent component, which is stored as a scalar field, while in three dimensions it is stored as
In order to calculate the director field \(\mathbf n\) and the amount of order \(S\) in two dimensions we use that we can map the orderparameter to the complex field \(\psi = Q_{xx} + iQ_{xy} =Se^{2i\theta}/2\), where \(\theta\) is the angle of the director field. In three dimensions we use that \(S\) is given by the largest eigenvalue as \(S = 3\lambda/2\) with the director being the coresponding eigenvector ^{1}. This is taken care of in the function
Wich returns S, n
, where n
is the director field.
Note that a nematic liquid crystall can be biaxial and given as
where \(P\) is given by the difference between the smallest eigenvalues and \(\mathbf m\) and \(\mathbf l\) is the corresponding eigenvectors.
Model
We model the active nematic using a set of coupled differential equations, namely the EdvardBeris equation coupled to the Stokes equation ^{2} ^{3} ^{4} ^{5}
The EdwardBeris equation (evolve_nematic
)
Here \(2\Omega_{ij} = \partial_i u_j  \partial_j u_i\) is the vorticity tensor, \(P\) is the pressure, \(\gamma\) is the rotational friction coefficient, \(\sigma^p\) is the passive stress, \(\Gamma\) is friction with a substrate, \(\eta\) is viscosity and the active stress is given by \(\sigma^a = \alpha Q\). The vorticity tensor is calculated by the function
Note that the velocity has to be updated before this function is called. The calculation of the pressure and velocity is described furhter down. Since the active stress is simply proportional to \(Q\) we have not included any function to calculate it, but calculate the force directly with the function
The molecular field \(H\) is given as
The last term is here to make it trace less. For the free energy we use
where it is assumed that there is a single Frank elastic constant \(K\). Note that the last term only exists in three dimensions since the trace of \(Q^3\) vanishes in two dimensions. The molecular field is then given as
For the passive stress we use \(\(\sigma_{ij} = Q_{ik}H_{kj}H_{ik}Q_{kj}  K (\partial_i Q_{kl})(\partial_jQ_{kl}).\)\) Note that we use the convention that the dot product between a vector and a tensor contracts the last component when calculating the divergence of this. The first two terms is the asymmetric stress, while the second term is the Ericksen stress. Terms due to flow allingment and anisotropic viscosity are not included. The molecular field and the passive stress are calculated by the functions
The linear and nonlinear part of the evolution equation for \(Q\), eq. is given as
The evolution of this is handled by the function
Disipative dynamics
Note that if we set the velocity field to zero the dynamics become \(\(\partial_t Q= \frac{K}{\gamma} \nabla^2 Q_{ij} +\frac{A}{\gamma}(B  2Q^2_{kk})Q_{ij}.\)\) This is used to relax the initial system before starting the simulation. The linear and nonlinear part of this equation are
An evolver for this dissipative dynamics is included as
The velocity field
For a given orderparameter \(Q\) the velocity field is calculated in Fourier space using eq (??). We start by finding an expression for the pressure by taking the divergence of this eq. (??) and then using the incompressibility condition giving \(\(\nabla^2 P = \nabla \cdot \mathbf F,\)\) where \(\mathbf F = \nabla \cdot \sigma^a(Q) +\nabla \cdot \sigma^p(Q)\) is the active and passive forces. This is solved in Fourier space as
The above equation can be inverted in order to find all the modes of the pressure except the zero mode, i.e the pressure is determined up to a constant. We set this constant to zero. Once the pressure is found we obtain the velocity from
(\Gamma + \eta k^2){{\mathbf u}_{\mathfrak f}} =  i\mathbf k {{P}_{\mathfrak f}} + {{F}_{\mathfrak f}}.
Note that when \(\Gamma = 0\) we need to set the zero mode of the velocity by hand. This is set to zero. The pressure and velocity are calculated/updated by the two functions
Note that calc_pressure_f
only returns the Fourier transform of the
pressure. The function conf_u
updates both the velocity field self.u
and its Fourier transform self.u_f
.
Minimum of the free energy
When starting a simulation it is often interesting to start from a configuration that is the minimum of the free energy pluss some perturbations or with a vortex dipole/fillament. From the free energy we see that the minimum energy is given by a homogeneous nematic, and it is inedependent of the direction the nematogens are pointing. Assuming that the unitvector \(\mathbf n\) is homogeneous we can rewrite the free energy in terms of the parameter \(S\).
In two dimmensions
the free energy is only given by powers of \(\text{Tr}(Q^2)\) which is \(S^2/2\) in terms of \(S\). The freeenergy is therfore for a homogeneous two dimentional nematic given as
The minimum of this is given as \(S =\sqrt{B}\) when \(B >0\) and \(S = 0\) if \(B<0\).
In three dimensions
In three dimensions we have that \(\text{Tr}(Q^2) = 2 S^2/3\) and \(\text{Tr}(Q^3)= 2 S^3/9\). Using this we find that there are a local minima at
when \(B > 3 C^2/(16A^2)\).
Topological defects and active turbulence
Because of the headtail symmetry if the nematic director the topological defects in the nematic phase can have half integer winding number. We can see this by maping the \(Q\) tensor to a complex field. This is done by writing the nematic director as \(\\mathbf{n} = \cos{\theta} \hat x + \sin{\theta} \hat y\), with \(\hat x /\hat y\) being the unit vectors in \(x /y\) direction, and mapping the \(Q\) tensor, see eq. (??), to the complex field
Using the same arguments as for the BEC we find that the allowed winding numbers \(\(k = \int_C \nabla \theta \cdot d\\mathbf{l} = 2\pi q\)\) with \(q\) being a halfinteger. The defects of lowest absolute charge is the \(\pm 1/2\) defects, which are depicted below.
Liquid crystal disclination dipole: The nematic director (headless vectors) around a defect dipole. The \(+1/2\) defect is marked with red, while the \(1/2\) defect is in blue.
For tracing the defect nodes one can use the function
If dt_Q
is given this finds the defects velocity and if polarization
is given the polarization of the \(+1/2\) defects are
found.
This polarization is given by
where \(\mathbf{r}_+\) is the defects position. The field \(\mathbf{e}_+\) can be found by the function
Note that this function does not include the normalization to avoid division by zero.
Initial States
A homogeneous nematic with random noise can be implemented with the function
This gives a state where the nematogens are aligned with the \(x\)axis. The noise is in the angle \(\theta\). If the activity is high enough this state will after a while start to spontaneously generate topological defects. In addition there is possible to insert defect dipoles using the function
which works similarly as the one implemented for the BEC. This function can be used either to initialize a homogeneous state with a dipole, or it can be used to insert a dipole into an already existing nematic.
Spatially varying activity
The activity \(\alpha\) is can be spatially varying. This can be used to make active channels as shown in the following figure.
Nematic liquid crystal active channel: Illustration of an active channel. width= 20
and \(d = 2\).
This simple channel with the activity \(\alpha_0\) inside and \(\alpha = 0\) outside is included as the function
Which sets the activity to \(\(\alpha = \alpha_0 \left[1 1/2 \left(\tanh([xw/2]/d) \tanh([x+w/2]/d) \right)\right].\)\) \(w\) is here the width of the channel, \(\alpha_0\) is the activity in the channel and \(d\) is the width of the interface between the channel and the bulk. More complicated structures can be created if one wish.
Three dimensions
In three dimensions, the \(Q\)tensor can be characterized as
where \((\mathbf n, \mathbf m, \mathbf l )\) is an orthornmal triad. It is five parameters: \(S\), the two defining angles of \(\mathbf n\), \(P\) and the angle of \(\mathbf m\) in the plane orthogonal to \(\mathbf n\). \(S\) can always be determined from the highest eigenvalue \(\lambda_{\textrm{max}}\) of \(Q\) by ^{1}
Topological defects
Topological defects in nematic liquid crystals are called disclinations and are characterized the orientation of the rodlike particles having rotated after following a path around the dislocation. From ^{6},
We can write this as
to reduce the number of sums preformed.
In two dimensions, where \(\mathbf n = (\cos \theta,\sin \theta)\), we have
where \(\psi_1,\psi_2\) are the components of an \(\mathcal D^2\) order parameter. We get only one component of \(D_{\gamma i}\), which is \(D_{33}\) which is
= \epsilon_{\mu \nu} \epsilon_{kl} \partial_k Q_{\mu 1} \partial_l Q_{\nu 1}
+ \epsilon_{\mu \nu} \epsilon_{kl} \partial_k Q_{\mu 2} \partial_l Q_{\nu 2}
We have \(Q_{\mu1} = \frac{1}{2} \psi_\mu\) and \(Q_{\mu 2} = \frac{1}{2} \epsilon_{\mu q} \psi_q\), so
D_{33} = \frac{1}{4} \epsilon_{\mu \nu} \epsilon_{kl} (\partial_k \psi_\mu )(\partial_l \psi_\nu)
+ \frac{1}{4} \epsilon_{\mu \nu} \epsilon_{kl} (\partial_k \epsilon_{\mu q} \psi_q) (\partial_l \epsilon_{\nu r} \psi_r)
And using that
we get
This is the same determinant as we would get using the coarse grain density of ^{7}, only with \(\psi_0\), so, the disclination density should be
In three dimenstions, the story is more complicated, because we have a tensor \(\rho_{\gamma i}\). This tensor contains two pieces of information, namely which direction the disclination is pointing, and around which axis \(\boldsymbol \Omega\), near the disclination, the rods are rotating. In that way, it is similar to a dislocation density in a crystal structure, only that it allows for the orientation the "Burgers vector" to be any direction. It can be written like this ^{6}
where the unit vectors are \(\boldsymbol t\) is tangent vector and \(\boldsymbol\Omega\) is the rotational vector. From this, we see that
so \(\omega\) is the quantity we should integrate to find the nodes of the defects. The unitvectors \(\boldsymbol t\) and \(\boldsymbol \Omega\) is found as the eigenvectors of the matrices \(\rho^T \rho\) and \(\rho \rho^T\) respectivly. Since the eigenvectors are determined up to a sign one have to make sure that \(\boldsymbol t\) is continous along the defect and impose the condition \(\text{sign}(\boldsymbol \Omega \cdot \boldsymbol t) = \text{sign}(\text{Tr}(\rho))\).
From Ref.^{6}, we have
replacing the delta function, which we may generalize to
Inserting topological defects
How do we insert topological defects of a given character?
We can generate an initial state of \(Q_{ij}\) by writing
and then simply impose an orientation field corresponding to an angle field on the \(\mathbf n\) fields.

Schimming, C. D. (2022). Theoretical and Computational Methods for Mesoscopic Textures in Nematic Liquid Crystals with Anisotropic Elasticity. PhD Thesis. The University of Minnesota. https://hdl.handle.net/11299/241713 ↩↩

Marchetti, M. C., Joanny, JF., Ramaswamy, S., Liverpool, T. B., Prost, J., and Rao, M. and Simha, R. A. (2013). Hydrodynamics of soft active matter. Reviews of Modern Physics. 85, 3, 1143. https://doi.org/10.1103/RevModPhys.85.1143 ↩

Genkin, M. M., Sokolov, A., Lavrentovich, O. D. and Aranson, I. S. (2017). Topological defects in a living nematic ensnare swimming bacteria. Physical Review X. 7, 1,011029. https://doi.org/10.1103/PhysRevX.7.011029 ↩

Nejad, M. R., Doostmohammadi, A. and Yeomans, J. M. (2021). Memory effects, arches and polar defect ordering at the crossover from wet to dry active nematics. Soft Matter. 17, 9, 25002511. https://doi.org/10.1039/D0SM01794A ↩

Angheluta, L., Chen, Z., Marchetti, M. C. and Bowick, Mark J. (2021). The role of fluid flow in the dynamics of active nematic defects. New Journal of Physics. 23, 3, 033009. https://doi.org/10.1088/13672630/abe8a8 ↩

Schimming, C. D. and Viñals, J. (2023). Kinematics and dynamics of disclination lines in threedimensional nematics. Proceedings of the Royal Society A. 479, 2273, 20230042. https://doi.org/10.1098/rspa.2023.0042 ↩↩↩

Skogvoll, V., Rønning, J., Salvalaglio, M., Angheluta, L. (2023). A unified field theory of topological defects and nonlinear local excitations. npj Comput Mater, 9, 122. https://doi.org/10.1038/s41524023010776 ↩