Introduction

In nematic liquid crystals, the four distortion modes -splay, bend, twist, and saddle splay -can each contribute differently to the elastic distortion energy, 1,2 a phenomenon hereafter referred to as ''anisotropic elasticity''. Even though the origin of this anisotropic elasticity can be traced to the relative alignment of elongated nematogens, and it is well documented, there still remain many open questions related to the effects of anisotropic elasticity on the equilibrium and nonequilibrium properties of defected nematics. A better understanding of the role of anisotropy on the motion and interaction of disclinations is fundamental to modeling biologically inspired and synthetic active matter systems.

In common thermotropic liquid crystals comprising small rod like molecules, the contrast between splay, twist, and bend elastic constants is small, and the so called one constant (''isotropic elasticity'') approximation has been successful in a wide variety of applications. More recently, however, attention has shifted to systems comprised of more complex nematogens which exhibit large elastic anisotropy. Chief among them, we mention lyotropic chromonic liquid crystals [3][4][5][6][7] and nematic micellar systems. 8,9 Novel behavior has been uncovered which is a direct result of elastic anisotropy, such as spontaneously broken chiral symmetry due to confinement, [8][9][10][11][12] or the existence and motion of topological solitons. [13][14][15] Complex anisotropic effects have also been observed recently in studies of disclination line reconnection in three dimensions. 16 In contrast with two dimensions, disclination lines in three dimensions only have a topological charge of 1/2, and can annihilate despite having the same charge sign. An apparent asymmetry in the motion of wedge disclination segments (of effective charge AE1/2) seems to be eliminated through twist in anisotropic media, thus restoring the implied topological symmetry.

The topology of defected configurations in two and three dimensional nematic phases is well understood, including the case of biaxial ground states. 17 In two dimensions, the orientation y(x) (see Fig. 1) of the nematic director n ˆis a harmonic function of position x in the one constant (isotropic) approximation. Well known singular solutions are associated with disclination point sources. 2,18 Configurations comprising many disclinations can be described by linear superposition, and results have been given for a number of cases of interest, including, for example, binding-unbinding transitions in active matter, 19 or defect interactions in complex twisted configurations obtained by conformal mapping techniques. 20 In contrast, little is known about nematic director n ˆor tensor order parameter Q configurations corresponding to defected configurations in elastically anisotropic media, both in two and three dimensions. A key result in two dimensions was obtained by Dzyaloshinskii. 21,22 When the splay K 1 and bend K 3 elastic constants are different, he found an analytic-albeit only implicit-solution for the equilibrium nematic orientation y corresponding to an isolated disclination. The solution is independent of distance from the core, but depends on the azimuthal angle around the disclination. More generally, the Euler-Lagrange equations that follow from the Frank free energy are nonlinear and challenging to solve analytically.

While it is possible to study both equilibrium and transient configurations of nematics containing disclinations in the director representation, with the Frank free energy governing elastic distortion, and Leslie-Ericksen hydrodynamics, it is often the case that a Q tensor order parameter representation and the Landau-de Gennes theory are used instead. Virtually all studies of nematic active and biological matter use this representation as it eliminates the need for defect core regularization (especially in three dimensions), and hence it permits a more convenient computational treatment of disclinations and their motion. Unfortunately, this choice has the effect in practice of restricting these studies to the one constant approximation. Elasticity in the tensor order parameter representation is incorporated in a phenomenological series expansion in powers of order parameter gradients, eqn (14) below. For small distortions, Frank elastic constants can be related to the coefficients of the expansion as shown in eqn (15). In order to capture splay-bend anisotropy, one must resort to at least cubic terms in gradients of the order parameter. At this order, however, the Landau-de Gennes energy is known to become unbounded for any choice of parameters. 23,24 In principle, the requirement of a bounded free energy could be accomplished by consideration in the expansion defining F el of terms at least of fourth order in Q. 25 However, it is also possible to have a bounded free energy, only third order in Q, by constraining the eigenvalues of Q to lie within their physically admissible range. 23 The resulting singular potential method sidesteps the need to choose between fourteen possible fourth order invariants 26 (in addition to choosing among six possible third order invariants).

Building into the theory the constraint that the eigenvalues of Q must remain within the physically admissible range can be accomplished by an appropriately defined singular potential. 23,[27][28][29][30] The drawback of this theory is that the determination of the energy needs to be done entirely numerically at a significant computational cost relative to simple evaluations of the Landau-de Gennes energy. Two complementary issues are investigated below in relation to elastically anistropic nematic phases, both in the tensor order parameter representation. First, we build on the singular potential method analysis of ref. 29 to quantitatively describe both bialixiality and anisotropy of disclination cores. We use the method to compute the optical retardance, G = S À P, near a disclination core, where S and P are the uniaxial and biaxial order parameters respectively. Exactly at the disclination core, S = P, in agreement with experiments 31 and earlier calculations. 29 We then consider a Fourier decomposition of the optical retardance Gðr; jÞ ¼ P n G n ðrÞ cosðnjÞ and show that as the core is approached G 0 B r, as happens in elastically isotropic systems. We also show that G 1 for a +1/2 disclination and G 3 for a À1/2 disclination are nonzero in the region of r B 1. However, they vanish as r 2 as the core is approached. Hence, the uniaxial and anisotropic far field leads to an anisotropic and biaxial region as the core is approached. At even smaller distances, the configuration becomes both uniaxial and isotropic, as judged from the azimuthal Fourier transform of G.

Second, we focus on the interaction of a pair of disclinations of opposite sign (a disclination dipole), and examine the nature of their screening at distances much larger than their separation. For isotropic elasticity, the orientation angle far from the disclination pair behaves as y = q 1 + q 2 À d(q 1 À q 2 )sin j/(2r) where q 1,2 = AE1/2 are the charges of the disclinations separated by distance d, r is the radial distance from the pair, and j is the azimuthal angle measured relative to the separation distance vector. For two disclinations of opposite charge, the distortion is screened and decays algebraically as 1/r, modulated by sinj in angular dependence. In the anisotropic case, the far field dependence contains an additional term of the form AEd sin(3j)/r which has the same decay with distance, but a different angular dependence. As a consequence, disclination interactions in elastically anisotropic nematics are qualitatively different than their isotropic counterparts, and the implications of these findings on current phenomenology involving multiple defect interactions and motion need to be reexamined.

Nematic director and Q tensor representations

In the director representation, local order in the nematic phase is described by a director field, the unit vector n(x). This field corresponds to the local average orientation direction of the constituent molecules, with configurations being invariant under the transformation n -Àn. The Frank free energy considers distortions away from a uniform ground state, and contains all scalar combinations of gradients of n to second order that respect n -Àn 2 ,

with K 1 , K 2 , K 3 , K 24 the elastic constants that correspond to splay, twist, bend, and saddle splay distortion modes respectively. In two dimensions, the twist and saddle-splay terms are manifestly zero. We introduce an anisotropy parameter e = (K 3 À K 1 )/(K 3 + K 1 ), dimensionless lengths % x = x/x where x is a characteristic length scale defined in eqn (18) in relation to the Q-tensor representation, and a dimensionless free energy

Dropping the overlines for simplicity one finds,

The minimizer of eqn (2) for a single point disclination in an infinite medium and for arbitrary e has been given by Dzyaloshinskii, though only implicitly as an integral equation. 21,22 The nematic director n = (cos y, sin y) is determined by the orientation field y, which is found to be independent of the distance r from the point defect, and depends only the azimuth j, i.e. y(j) (see Fig. 1). The Euler-Lagrange equation for the minimizer of the Frank free energy (2) is

d 2 y dj 2 cos 2ðy À jÞ þ 2 dy dj À dy dj 2 ! sin 2ðy À jÞ " # :

(3)

In the isotropic limit of e = 0, the director orientation is multivalued y iso (j) = qj, where q = AE1/2 is the disclination charge. 31 A perturbative solution in e can be found by expanding,

where the first order correction is nonlinear in j 31 y c ¼ qð2 À qÞ 4ð1 À qÞ 2 sinð2ð1 À qÞjÞ:

This expression also follows directly from Dzyaloshinskii's solution -see Appendix D for details.

In order to capture both the magnitude of local order and biaxiality, a tensor order parameter representation is commonly introduced. It is a coarse-grained, statistical measure of nematic alignment. In three dimensions it is defined as

Here r(p) is the probability density function of molecular orientation p defined on S 2 , the unit sphere, and ds is the surface measure on the sphere. We have denoted by I the rank three identity tensor. Because of nematic symmetry, one has r(p) = r(Àp). By definition, Q is traceless and symmetric. Its three eigenvectors n, m, l form an orthonormal basis, so that Q may be written as,

S and P can be written in terms of the three eigenvalues, l 1 Z

The eigenvectors corresponding to l 1 and l 2 are n and m respectively. The scalar order parameter S describes the degree to which molecules are aligned along the director n, while P describes biaxiality, or the difference in alignment along the two remaining axes. A Landau-de Gennes free energy expansion is introduced in terms of scalar contractions of Q (the ''bulk'' terms), supplemented by terms in gradients of Q (the ''elastic'' terms). For small distortion and fixed S, the elastic terms in the Landau-de Gennes free energy may be mapped onto the Frank elastic free energy exactly. In order to include bend-splay anisotropy, one must expand the elastic energy at least to third order in gradients of Q. It is well known, however, that at this order the free energy is unbounded below. 23,24 A possible remedy involves consideration of gradient terms of fourth order in Q. 25 It is also possible to maintain a third order theory, and avoid choosing among fourteen possible fourth order terms allowed by symmetry, by introducing the Ball-Majumdar singular bulk potential method. 23,30 A bulk free energy F b ½Q ¼ E½Q À TDS½Q is defined where E is the bulk energy, T is the temperature, and DS is the entropy relative to the isotropic phase. The energy is chosen to be of the Maier-Saupe form

where k is a positive constant that characterizes alignment strength. The entropy may be written in terms of the molecular probability distribution function,

where n is the number density of nematogens, k B is Boltzmann's constant, and the probability density function of molecular orientation r is allowed to be a function of position for an inhomogeneous configuration. In order to find an explicit expression of DS in terms of Q, r is determined so that it maximizes DS subject to the constraint (6). The solution is,

with partition function Z given by:

where K is a tensor of Lagrange multipliers arising from the constraint (6). By substituting eqn (9) into eqn (6) we may relate the multipliers K to Q as a mean field consistency condition,

Substituting eqn (9) into eqn (8) and using eqn (11) to simplify, the entropy may be written in terms of Q as,

where: is a double index contraction.

For the elastic free energy in our present study, we include only one term of third order in Q to allow for bend-splay anisotropy,

where . .

. is a triple index contraction from inner indices to outer indices, and L i are the elastic constants. Written in index notation this equation reads,

We recall that the mapping to the Frank free energy coefficients in the case of a uniaxial and constant S nematic phase is given by: 32

The total free energy in the singular potential method is the sum

Rotational relaxation dynamics of the nematogens is considered through

with g a rotational diffusion constant. We introduce dimensionless variables,

where the length and time scales are,

Dropping the overlines for simplicity, the dimensionless equation of motion for Q is,

with the transpose of a rank-3 tensor being defined as (rQ) T klj = q j Q kl . Hereafter, all distances and times will be dimensionless.

The partition function defined on the unit sphere ( 10) must be evaluated numerically, as well as the self consistency condition (11) to find K = K(Q). Stationary solutions of eqn (19) are found by using the Newton-Rhapson relaxation method for the case of configurations with one isolated disclination. For the case of a disclination pair, however, the Newton-Rhapson method is not computationally efficient due it to its slow convergence for large systems. Instead we discretize eqn (19) in time by using a Crank-Nicolson method. We then use the same Newton-Rhapson method to solve for each subsequent time step, and iterate in time until q t Q is sufficiently small. We have implemented this singular potential method in a new finite element formulation, based on the framework deal.ii, that allows for efficient paralellization. Large three dimensional configurations can be efficiently studied at high resolution (in the scale of x). The Appendices provide additional numerical details.

Boundary conditions in a finite domain need to be discussed separately. Given the variational derivative of the energy

; we impose Neumann boundary conditions by requiring that the normal component at the outer boundary NÁqf/q(rQ) = 0, where N is the outward pointing normal. This reduces to the familiar Neumann boundary condition on Q in the isotropic limit, but more generally, it is the natural boundary condition to use for a fully anisotropic system.

A single disclination in the Q tensor representation

We present first the results of a high resolution numerical study of Q for a single disclination in an elastically anisotropic medium (L 3 a 0). We show that the singular potential method can quantitatively describe the biaxial core region around the disclination, and that the stationary configuration reduces to the Dzyaloshinskii solution away from the core where the nematic configuration becomes uniaxial. The thin film approximation for

so that the tensor is described by three independent components, not just two as in a strictly two dimensional case, and hence biaxiality can be accommodated. The xz, yz, zx, and zy components of the right-hand side of eqn (19) are manifestly zero because K and Q can be simultaneously diagonalized, 30 and q z Q = 0. Hence, any configuration initialized in the thin film approximation will remain as such without further constraint on the equation of motion. Additionally, the thin film approximation restricts all eigenvectors to lie in the x-y plane or along the z-axis. For a configuration with directors initialized in the x-y plane, the only way for the director to escape into the z-direction is for Q zz to become equal to the larger of the other two eigenvalues, creating the so-called ''pancake'' configuration. This does not happen in our configurations, though a clarifying visualization for how this manifests in the x-y plane for disclinations can be found in ref. 33 The biaxial core region has been extensively studied in the one constant approximation, 34,35 and in a more general case that included all possible terms in gradients up to second order in Q. 36 Strong biaxiality develops in the core region of the disclination. For a Landau-de Gennes bulk energy, a purely uniaxial configuration is shown not to be stable; although uniaxial far from the core, the three eigenvalues of Q become distinct as the core region is approached, and two of them eventually cross at the disclination line. 36 The core structure of Q has also been recently characterized experimentally in lyotropic chromonics, 31 enabled by a large size of their core (tens of microns). A biaxial region has been confirmed in the optical retardance, albeit with a strong angular dependence due to elastic anisotropy. This angular dependence of the retardance has been shown to be in agreement with results of the singular potential method. 29 A stationary solution of eqn (19) in the thin film approximation has been obtained in a two dimensional circular

; with an isolated AE1/2 disclination near its center maintained by appropriate Dirichlet boundary conditions on the outer boundary. We choose dimensionless values of the parameters k = 8.0, L 2 = 4.58, L 3 = 4.5. k has been chosen so that the system is below the supercooling limit as in the experiments of ref. 31 and simulations of ref. 29, which corresponds to an equilibrium value of S to be S 0 = 0.6751. L 3 is chosen to be as large as possible while maintaining numerical stability, while L 2 is chosen to maintain e = 0.4 through eqn (15), consistent with ref. 29 and 31. The most notable effect of taking a different e value would be to change the director profile far from the disclination core, as can be seen from eqn (5). The effect of taking L 3 larger while keeping a fixed e value is to increase the higher Fourier mode amplitude.

The computational domain is discretized with quadrilateral elements, initially with 12 cells. It is then globally refined 5 times, and further refined at distances R ¼ 8; 4; 2; 1; The director n and scalar order parameters S and P are determined by calculating the eigenvalues and corresponding eigenvectors of the Q tensor at each point in the computational domain. This is done with the eigh method from the Numpy numerical package, which calculates the eigensystem of a symmetric matrix. 37 We find that the stationary disclination cores are located at (x disc , y disc ) = (0, 0) and (0.868, 0) for the À1/2 and +1/2 disclinations respectively. The quantity Gðr 0 ; j 0 Þ ¼ ðS À PÞ is computed as it is proportional to the optical retardance in the experiments. 31 Primed variables are polar coordinates referred not to the center of the computational domain, but to the actual disclination center (x disc , y disc ) defined as the location where S = P. To probe the effect of anisotropy, an angular Fourier transform is introduced,

The Fourier coefficients are calculated with the rfft real Fourier transform method from the Numpy numerical package. The cosine coefficients in eqn (20) are 2/N times the real part of the discrete transform modes, where N is the number of grid points at each r 0 . 37 Fig. 2a and c show the director angle y vs. the azimuth j 0 plotted at several fixed distances from the disclination centers. At large distances, the director angle approaches the Dzyaloshinskii perturbative solution of eqn (3) calculated relative to the domain center (as is appropriate for the boundary conditions), but plotted as a function of j 0 at several values of r 0 . Explicitly, if y DZ (j) is the solution to eqn (3), the solid line in Fig. 2a is given by y DZ atan2 r 0 sin j 0 þ y disc ; r 0 cos j

Þ for r 0 = 10. For small values of r 0 the director angle approaches a straight line in the diagram, the isotropic solution y ¼ 1 2 j 0 .

As r 0 increases, however, the angle tends towards the Dzyaloshinskii uniaxial solution. In order to further probe the biaxial core region, Fig. 2b and d show the two dominant angular Fourier modes G n (r 0 ). The figures also show a fit to a power law with distance. The zeroth Fourier modes goes to zero linearly, while the higher Fourier modes appear to decrease quadratically as the disclination center is approached. The determination of this dependence has been made possible by the high spatial resolution of our numerical method. Neither prior research nor the experimental work could make this determination.

The singular potential method with L 3 a 0 predicts a compact biaxial core, with amplitudes of the angular Fourier components of G vanishing faster with distance to the defect center than the zeroth order component. Therefore the director angle approaches qj 0 as is the case for an elastically isotropic medium. Furthermore, the dominant dependence of the eigenvalues is also linear as the core center is approached, in agreement with earlier isotropic results. Both results suggest that the isotropic and linear core approximation is a reasonable approximation even in anisotropic media.

A disclination dipole

The complicating factor that remains, and to which we turn next, is that in two or multi defect configurations, the tensor field is not a superposition of configurations corresponding to isolated single defects. Therefore it remains to be seen whether interaction leads to a more complicated core structure in multi disclination systems.

The Euler-Lagrange equations corresponding to the Frank energy (2) in Cartesian coordinates read,

Consider now a pair of disclinations a distance d from each other, which are mutually aligned or anti aligned. We seek a perturbative solution for the director field to first order in e. 38,39 The solution in the isotropic limit of e = 0 can be written as

where q 1 , q 2 are the corresponding disclination charges, and we have introduced polar coordinates (r i , j i ) centered at each defect position (x i , y i ) (see Fig. 3 for a diagram of the relevant coordinates). The constant term rotates the director everywhere by p/2, a transformation under which eqn (2) is invariant. For q 1 and q 2 half integers of opposite sign, this solution and the corresponding one without the constant term are so-called ''isomorphs'', characterized by whether the line connecting the two defects is parallel or perpendicular to the far-field director. For example, with q 1 = +1/2 and q 2 = À1/2, eqn ( 22) is the perpendicular isomorph. By expanding yðx; yÞ ¼ y iso ðx; yÞ þ ey c ðx; yÞ þ O e 2 À Á ; and substituting into Eq. ( 21) we find a Poisson equation for the first Fig. 2 (a) and (c) Director angle y as a function of the azymuth j 0 at various distances from the core for +1/2 and À1/2 disclinations respectively, computed from the equilibrium Q tensor. The solid line is y DZ atan2 r 0 sin j 0 þ y disc ; r 0 cos j 3 Diagram showing a disclination pair in polar coordinates. Here (r i , j i ) are polar coordinates centered on the disclination with charge q i , and (r, j) are polar coordinates centered on the midpoint between the two disclinations.

order correction y c :

We point out that the other isomorph merely changes the righthand side -and therefore the solution -by a sign. In what follows, we find an approximate solution to eqn (23) in various regions which can then be compared against numerical results.

For concreteness, we choose q 1 = +1/2 and q 2 = À1/2. Near one of the disclinations, (x 1 , y 1 ), one may rewrite j 2 and r 2 in terms of j 1 and r 1 . In this region, r 1 /d { 1 so that we Taylor expand the right-hand side to find,

A particular solution y p,1 c can be found as given by

By comparing it with eqn (5), we note that the term independent of r 1 corresponds to the correction for an isolated disclination in an anisotropically elastic medium, while the term due to pairwise disclination interaction is new and goes linearly in r 1 close to q 1 . A similar calculation for the region close to q 2 yields a particular solution,

Again we obtain a term independent of r 2 which is identical to eqn (5), and an interaction term which is linear in r 2 .

Finally, in the far-field, one may rewrite the equation to first order in polar coordinates whose origin is midway between the two defects (r, j). Expanding the inhomogeneous term in d/r { 1 yields,

A particular solution to second order is given by,

The dependence on 3j and proportional to d/r at long distances is unexpected. Consider the isotropic solution eqn (22), and express it in terms of the midpoint polar coordinates,

If q 1 + q 2 = 0 the constant terms identically vanishes (charges mutually screen), and the dipolar term has the expected dependence in d/r sin j from a multipolar expansion. However, anisotropic elasticity changes charge screening, and it introduces a new term that, while also decaying as d/r at long distances, has a different angular dependence.

A general solution which matches the particular solutions in the inner and far field regions would also require the general solution to Laplace's equation. Far from the disclination pair, one would have,

The inner solutions include the components n = 1, n = 2, n = 3, and (although much smaller in magnitude as we will argue below) n = 4 components. Hence, we would expect those Fourier modes to be present in the far field in order to match at the near-field far-field boundary, giving an approximate far-field solution of:

We will not pursue this analytic expansion further. Rather we will argue that this dependence is consistent with our numerical solutions for weak elastic anisotropy shown below.

5 Numerical solutions for a disclination pair

Director representation

Eqn ( 23) is a Poisson equation in which the source term is singular at the location of the two disclinations. We have modified a preexisting deal.II library program to solve it. 40,41 The actual linear system is solved with the conjugate gradient method with Trilinos ML algebraic multigrid as a preconditioner. 42 As was the case with the Q tensor, we take as outer boundary condition a zero normal component of the configurational force, where here the configurational force is qf n /q(ry) with f n the Frank elastic energy density. Because the solution is found perturbatively, the boundary conditions must be specified order by order (see Appendix C, eqn (64) for details). We solve on a circular domain radius R = 5500 and defect spacing d = 60. These dimensions have been chosen to correspond with the Q-tensor configuration solution shown later. We also solve eqn (23) inside a modified circular domain that excludes the singular points in its right hand side. We cut out two small discs around each disclination, and impose Dirichlet boundary conditions on the circumference of each discs. For simplicity, we prescribe y c = 0 on these internal boundaries which corresponds to y = y iso from eqn (22). We choose the cutout radius r cutout = 10 because, as evidenced in Fig. 2c and a, an isolated disclination in the Q-tensor formulation becomes uniaxial with approximately constant-S at approximately r = 10. The choice of domain is motivated by the comparison carried out below with a full numerical solution in the Q representation with the same value of the anisotropy parameter e. In the Q-tensor formulation, the configuration with two disclinations is not stationary, and hence allowing an unconstrained configuration relax leads to disclination annihilation. This would prevent us from determining the constrained equilibrium configuration corresponding to two immobile disclinations.

Fig. 4a shows a colormap of y c , both in the far field and near field limits. Near the disclination cores one may clearly see the n = 1 and n = 3 mode contributions from eqn ( 25) and ( 26) around the +1/2 and À1/2 disclinations respectively. The far field appears to have six fold symmetry, consistent with a contribution from n = 3. In order to quantify the contribution from the various Fourier components to y c , we decompose the far field numerical solution into angular Fourier modes,

and fit each mode A n (r) by a polynomial in 1/r, with a degree consistent with eqn (30). For example, A 3 is allowed to have degree 1 and 3 in 1/r, while A 2 is only allowed to have degree 2. Fig. 4b shows the angular Fourier coefficients and the corresponding fits. Both the n = 1 and n = 3 Fourier modes are consistent with the prediction, while the n = 2 and n = 4 modes deviate somewhat from the expected quadratic and quartic behavior. The linear dependence of the n = 3 mode matches the prediction from eqn (30) in both magnitude and sign.

The effect of adding cutouts to the integration domain around disclination cores is to suppress the near field n = 1 and n = 3 mode contributions, as can be seen in Fig. 4d. This reduction translates in the far field into a small reduction in the magnitude of the n = 3 mode, and a noticeable reduction in the amplitude of the n = 1 mode.

In agreement with the perturbative calculation of Section 4, these numerical results show a different angular dependence of the director angle that arises from disclination interactions in an anisotropic medium. The n = 3 Fourier mode decays at the same rate with distance as the n = 1 mode arising from the isotropic solution, although it is a factor of e/2 in magnitude smaller. Depending on the value of the anisotropy parameter, this term could introduce a significant deviation relative to the isotropic interaction terms, and must therefore be considered in, for example, disclination ensemble dynamics in elastically anisotropic media. Note also that the sign of the n = 3 far field term changes under the transformation to a different disclination pair isomorph. Hence, it is possible that the effective contribution from elastic anisotropy could be smaller in an ensemble of defects containing a distribution of isomorphs. This journal is © The Royal Society of Chemistry 2024 5.2 Q Tensor representation

With our choice of elastic terms, eqn (13), elastic anisotropy is determined by the coefficients L 2 and L 3 while the Frank elastic anisotropy is solely determined by e. Given eqn (14), we focus on L 2 = 0 and find that L 3 = 0.3065 for e = 0.1, a regime in which eqn (23) should hold. We note that the results are essentially identical for any other L 2 value, supposing that L 3 is chosen to maintain e = 0.1. This is because the L 2 term in eqn ( 14) may be decomposed into gradients of the scalar order parameters and director. Since the disclinations are cut out, the scalar order parameter remains constant and uniform. The contribution from L 2 to the director is to introduce twist anisotropy which, in two dimensions, is manifestly zero. We consider a disc of radius R = 5500, defect spacing d = 60, and defect cutout radius r cutout = 10. The Maier-Saupe constant k = 8.0, which corresponds to an equilibrium value of S 0 = 0.6751.

Because of the large size of the computational domain, a direct solution of the minimization problem (eqn (19) with q t Q = 0) is difficult. We instead iterate eqn (19) in time until a stationary configuration is reached. As initial condition we choose,

where R is a rotation matrix about the z ˆaxis by angle y c , which is the numerical solution to eqn (23) with disclination cutouts fixed at zero. We define

and n ¼ cos y iso sin y iso 0 ½ T . Fig. 5 shows y c as calculated from the Q tensor representation compared to y c from eqn (23) within the cutout domain. y c is well-defined in this case because the director remains in the x-y-plane, as has been verified.

Isolated disclination motion far from a dipole

To give a suggestion for a potential experimental avenue which may be explored to verify the far-field dipole director profile, we derive the equation of motion of an isolated disclination under the influence of a dipole using the Halperin-Mazenko formalism developed in ref. 43. The calculation is done in 2D, though the results are similar to a previous calculation done in 3D. 33 For this, we assume a disclination director profile of qj, and a scalar order parameter which decreases linearly to zero at the disclination core. Further, we assume that the director field of the isolated disclination superposes with the ambient director field created by the dipole, and neglect distortions to the dipole profile that would arise from interactions with the isolated disclination.

For a given ambient director angle field y produced by the dipole, the velocity of a test disclination is determined by the disclination density current which is derived in Appendix E. For a +1/2 disclination, the defect velocity is

with r > = q y x ˆÀ q x y ˆ, while for a À1/2 disclination, it is

Because y is small in the far field, the contribution from the L 3 term in eqn (33) gives a nearly uniform contribution to the velocity in the Àx ˆdirection for both the isotropic and anisotropic parts of the dipole director profile. By contrast, the first term in eqn (33) gives qualitatively different behavior from these two parts. To see this, we calculate the following explicitly:

This field is plotted in Fig. 6 with n = 1 for the isotropic contribution and n = 3 for the anisotropic contribution. For an isotropic dipole profile, one would expect a disclination in the upper half plane to move mostly in the azimuthal direction, while the anisotropic dipole profile would tend to cause the disclination's path to fluctuate in the radial direction. We speculate that this fluctuation is measurable, and should vary linearly with e. A material for which e is tunable, such as the biopolymer suspension in ref. 44, could give a quantitative measure of the magnitude of this fluctuation in e.

Conclusions

We have presented an analysis of the radial and angular dependencies of the orientation order parameter around both an isolated disclination and a disclination dipole in an elastically anisotropic nematic. In the former case, a singular potential theory in the Q tensor order parameter representation of the nematic shows that the order parameter approaches isotropy near the core: the eigenvalues of the Q tensor become axisymmetric, in agreement with the elastically isotropic case. We provide a scaling law which shows that the zeroth order angular Fourier of the retardance goes to zero linearly with the radial distance r 0 , while the next order Fourier mode decreases quadratically.

For the case of a disclination dipole, we have presented analytical perturbative solutions in the director representation in the limit of weak anisotropy (small elastic constant e). Solutions are given for the nematic orientation angle both near one of the disclinations in the dipole, and in the far field. Particularly noteworthy is the far field dependence in which the n = 1 angular Fourier mode of the isotropic limit is supplemented by an n = 3 mode as a leading order term due to anisotropy. The predictions agree very well with numerical calculations in both the director field and Q-tensor representations of the nematic. We speculate that the difference in the dipole director profile due to anisotropy can be experimentally observed through the motion of a test disclination under the influence of the dipole.

Conflicts of interest

There are no conflicts of interest to declare.

Appendices

A Numerical method for an isolated disclination in the director representation

The numerical solution of eqn (3), the one dimensional profile of the director y, as a function of polar angle j is computed by using the finite element framework deal.II. 40,41 The equation is solved by iteration with a Newton-Rhapson method on the domain j A [0,2p]. The endpoints are fixed at 0 and 2pq to maintain azimuth continuity. The equation residual is defined as,

A Gateaux derivative is introduced,

with

Also define:

so that we may write the residual as:

An iteration in Newton-Rhapson method then reads: dR y ðnÞ dy ðnÞ ¼ ÀR y ðnÞ y ðnþ1Þ ¼ y ðnÞ þ ady ðnÞ (41) with damping parameter a r 1. To solve with the finite element method, we take the inner product with a test function Z and integrate by parts:

The test functions are zero on the boundaries so that the surface integrals vanish. Approximating dy ¼ P j dy j Z j with test functions Z j given by piecewise polynomial Lagrange elements, and enforcing eqn (42) for each test function Z i gives a linear system in dy j . We iterate until the L 2 norm of the residual is less than some desired threshold. For the simulations run in this paper, the domain is broken into 2 10 evenly-spaced segments, we use first degree Lagrange elements, and the residual L 2 norm tolerance is set to 10 À10 . We use the UMFPACK direct sparse matrix solver since, in one dimension at this size, performance is not an issue.

B Numerical method in the Q-tensor representation

In order to solve eqn (19) numerically we also use the deal.II finite element framework. 40,41 This library has the benefit of implementing adaptive mesh refinement, as well as being massively paralellization via MPI, allowing for very large scale computations. To solve all linear systems in our implementation, we use the Trilinos linear algebra library via deal.II. 45 The code developed is available in the GitHub repository. 46 To integrate eqn (16), consider that the variation of the free energy is given explicitly by:

where f is the free energy density. Here we take nÁqf/q(rQ) = 0 as a boundary condition which corresponds to zero normal configuration force at the boundary. Additionally, to ensure that q t F r 0 always, we must take:

One may understand this as taking the time evolution in the direction of the variation dQ where the variation is chosen to make dF negative definite. To simplify the exposition, take T Q = Àqf/qQ and T rQ = qf/q(rQ). Finally, T = T Q + rÁT rQ . These are given explicitly by:

We note that the divergence is contracted over the k index.

To discretize eqn (19) in time, we use a Crank-Nicolson method:

where Q 0 and Q are the Q-configurations at the previous and current timesteps respectively, dt is the timestep, and T and T 0 are evaluated at Q and Q 0 respectively. Because T is nonlinear, we define a residual:

To solve for the configuration when R = 0, we use a Newton-Rhapson method. The Gateaux derivative then reads:

Explicitly, this yields:

where dL ij is given by:

The Taylor series expansion of L about Q involves the directional derivative in the direction of dQ. Since Q and dQ are restricted to the submanifold of traceless, symmetric tensors, this directional derivative can be accomplished by differentiating L with respect to the degrees of freedom of Q and dotting into the degrees of freedom of dQ. This set of degrees of freedom is arbitrary, but we note that the space of traceless, symmetric tensors is five-

This journal is © The Royal Society of Chemistry 2024 Soft Matter, 2024, 20, 2900-2914 | 2913

Changing coordinates to correspond with Fig. 3 the first term of the disclination current may be represented as:

with n = 1 for the isotropic contribution and n = 3 for the anisotropic contribution.