Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Simon
Čopar
*^{a} and
Žiga
Kos
*^{abc}
^{a}Faculty of Mathematics and Physics, University of Ljubljana, Ljubljana, Slovenia. E-mail: simon.copar@fmf.uni-lj.si; ziga.kos@fmf.uni-lj.si
^{b}Department of Condensed Matter Physics, Jožef Stefan Institute, Ljubljana, Slovenia
^{c}International Institute for Sustainability with Knotted Chiral Meta Matter (WPI-SKCM^{2}), Hiroshima University, Higashi-Hiroshima, Japan

Received
16th May 2024
, Accepted 5th August 2024

First published on 8th August 2024

From incompressible flows to electrostatics, harmonic functions can provide solutions to many two-dimensional problems and, similarly, the director field of a planar nematic can be determined using complex analysis. We derive a closed-form solution for a quasi-steady state director field induced by an arbitrarily large set of point defects and circular inclusions with or without fixed rotational degrees of freedom, and compute the forces and torques acting on each defect or inclusion. We show that a complete solution must include two types of singularities, generating a defect winding number and its spiral texture, which have a direct effect on defect equilibrium textures and their dynamics. The solution accounts for discrete degeneracy of topologically distinct free energy minima which can be obtained by defect braiding. The derived formalism can be readily applied to equilibrium and slowly evolving nematic textures for active or passive fluids with multiple defects present within the orientational order.

Complex analysis is an important tool to generate solutions of field theories in two dimensions. As all complex analytical functions solve the Laplace equation, models with energy proportional to (∇ϕ)^{2} can be solved by finding analytical functions that satisfy the boundary conditions. Additionally, conformal mapping can be used to transform the problem to a domain where solution can be readily obtained. A well known example is 2D hydrodynamics of ideal incompressible fluids, with Kutta–Joukowsky theory giving the lift generated by an airfoil by transforming a solution on a unit circle to a realistic airfoil domain.^{24} Complex analysis is also often used to solve for two-dimensional potentials in electrostatics^{25} and solitonic solutions in ferromagnets.^{26} A similar approach can be taken for fluids with an orientational order parameter isomorphic to e^{inϕ}, where n describes the polar (n = 1), nematic (n = 2), or higher-order symmetry, and ϕ the polar angle of the orientational order. For example, in the approximation of a single elastic constant, alignment of the nematic liquid crystals restricted to the 2D Euclidean plane can be described by a free energy , where K is the elastic constant and = (cosϕ,sinϕ) is the director field. The equilibrium texture is solved by harmonic functions ϕ and singular points of ϕ correspond to topological defects (disclinations).

The conformal description of topological defects is regularly used to generate analytical models of liquid crystalline textures. Solutions for non-trivial planar geometries are often achieved by conformal mapping,^{27–33} or by stereographic projection for spherical domains,^{34,35} for describing defects dynamics in active nematics in Q-tensor formalism,^{36} and active nematic textures in the director formalism.^{37} Complex field approach can also be used to solve for a shape of nematic domains with a free boundary under different anchoring conditions.^{28,38} Some of the well known solutions are the core energy of a single nematic defect,^{39,40} the director field and the interaction energy of a pair of defects,^{41} and drag force and mobility tensor for moving defects.^{40,42} These solutions focus on topological defects, where the director field rotates by a multiple of π along a loop circumnavigating the defect core. Another singular solution of the planar Laplace equation is a type of defect, where the director polar angle changes logarithmically with the distance from the defect core.^{40} Such spiral solutions with a radially varying director naturally occur in, for example, ‘magic’ spiral problems.^{39,43,44} However, a solution for an arbitrary number of defects with a given orientation and boundary condition has to include a combination of topological defects and logarithmic singularities. No such general framework has yet been formulated.

In this paper, we develop a complete formalism for planar nematic textures involving an arbitrarily large set of point defects with arbitrary winding numbers and logarithmic singularities, quantified with the newly introduced spiral charge. We explore the expression for the director azimuthal angle ϕ in the form of

(1) |

Φ(z) = ϕ(z) + iψ(z) | (2) |

w(z) = e^{iΦ(z)} = e^{−ψ(z)}e^{iϕ(z)}. | (3) |

The nematic elastic energy can be expressed in the regime of one-elastic constant K as

(4) |

(5) |

(6) |

w(z) = e^{iϕ0}(z/ε)^{k+iμ} = (r/ε)^{k+iμ}e^{iθ(k+iμ)+iϕ0}, | (7) |

We can extract the complex polar angle Φ(z) = ϕ(z) + iψ(z) by taking a logarithm of eqn (7):

(8) |

(9) |

The singularity at the origin has a diverging elastic energy, the nematic region is truncated at a finite radius ε, representing the size of the defect core. In the far field, we bound the system at a large radius R, and integrate over the contour in Fig. 1a. The elastic energy

(10) |

(11) |

(12) |

(13) |

(14) |

The defect pair free energy is given by integrating eqn (6) along the contour shown in Fig. 2a. Considering the confinement radius R that is much larger than defect–defect separation d, yields the free energy that completely decouples the winding numbers and the spiral charges, F = F_{k} + F_{μ}:

(15) |

(16) |

The free energy is the basis for finding equilibrium director configurations at fixed boundaries, or computing forces and torques acting on defects and inclusions. The complete symmetry between winding and spiral charges with no cross terms indicates that the only way for these contributions to couple is through boundary conditions, which we will explore later.

The free energy is a positive definite quadratic form in terms of the spiral charges μ_{1} and μ_{2}, and in the absence of boundary conditions, reaches a global minimum when all spiral charges are zero. Applying constraints changes that. Consider an example with a k_{1} = +1 defect enforced by a circular inclusion with a radius ε_{1} enforcing radial director, an accompanying k_{2} = −1 defect with core radius ε_{2} pinned at the distance d along the x axis, and homogeneous far-field in the direction ϕ(R) = ϕ_{∞}.

The boundary condition on the far-field director (eqn (13) under the limit R ≫ d) yields the following constraint:

(17) |

(18) |

(19) |

An important aspect of this problem is that the length scales are hierarchically well separated. The size of the particles or defect cores must be significantly smaller than the separation between them, ε ≪ d, so that the boundary conditions on the circuits of radius ε around the defects can be accurately met. If the particles are too close together, spatial variation of the director field ϕ(z) caused by surrounding particles, cause a deviation from the ideal boundary condition, even though the conditions are still matched in an average sense. Such configuration can still be sufficient to create initial conditions for further numerical simulations. A similar condition applies to the size of the entire system that should be much larger than interparticle separation, R ≫ d. For more tightly confined systems, the method of mirror charges can be used to improve upon our approximation in particular cases that allow it, or an exact solution can be obtained if we can find a conformal mapping from the nematic domain to one where boundary conditions can be met, such as used by Tarnavskyy et al.,^{29–31} instead of just using point defects represented by simple zeroes and poles.

An important feature of the energy minimisation described above is that we obtain not only a single solution, but a family of local energy minima, enumerated by n, i.e. the number of π turns the director makes to match the boundary condition. Each solution does have a different free energy, but the discrete nature of the solutions means they are metastable, and more than one can be realised in experiments, depending on the initial condition. It may also be possible to switch between them with an appropriate manipulation technique.

(20) |

The free energy expression (eqn (20)) can be rewritten in the matrix form,

(21) |

To apply the boundary conditions, we take the same steps we took for the two defects, evaluating the director angle on a path around each defect, and neglecting the displacement from the core radius. Eqn (18) rewrites into

(22) |

The difference f_{i} between the boundary orientation b_{i} and the k-contribution of the rest of the defects equals the angle that must be accumulated by the μ-part, and may be complemented by an arbitrary multiple of π due to rotational symmetries of the director,

(23) |

(24) |

(25) |

We want to find a free energy minimum corresponding to a constrained orientation of some of the defects. To achieve this, free energy has to be reparameterized in terms of f_{i} terms rather then μ_{i}. As some of the defects or the outer boundary might not be necessarily fixed, some of the components f_{i} may be unconstrained “slack variables” that need to be minimized over to find the ground state. We compute the block-wise inverse of the constraint matrix

(26) |

μ = M^{−1}(f_{0}1 − f). | (27) |

(28) |

For any constraint, there is an infinite discrete set of local minima, parameterized by indices n_{i}, as shown for selected values of n_{i} in Fig. 3. The meaning of these indices reveals that a system of defects with m constraints has m − 1 integer-valued topological indices (excluding one due to the fact that adding the same multiple of π to the entire director changes nothing). Varying the index n_{i} changes the number of relative half-turns the director makes between i-th defect and all other constrained defects. A system of a large number of constrained 2D defects, such as a lattice of micropillars, colloidal inclusions or other topological obstructions with fixed boundary conditions, will thus always be multistable, with a multidimensional lattice of multistable states. Although there will be a practical limit on heavily wound spirals that are hard to create and stabilize, we expect at least a limited number of metastable states. Although this solution was derived for small circular defect-inducing boundaries, the topology and multiplicity of solutions will hold for inclusions of arbitrary shape with topology of a disk, as long as no additional singularities appear on the surface itself.

The values of the indices n_{i} can be set in advance to obtain solutions for each set of indices. Alternatively, we can view the final system of equations as a quadratic minimization problem with m − 1 discrete variables, N − m continuous variables and one zero-mode. In this view, constraining the director on each boundary quantizes one dimension of the parameter space.

(29) |

(30) |

The first term gives rise to pairwise central forces obeying an inverse distance law:

(31) |

(32) |

We illustrate the importance of orientational constraints on the forces acting on the defect pair in Fig. 4. Without orientational contraints, all spiral charges are zero and the forces are central (Fig. 4a). The magnitude and the direction of forces changes if both defects are constrained (Fig. 4b) or if one of them is constrained (Fig. 4c and d). One possibility to prescribe the orientation of a defect profile is to determine the anchoring of nematic molecules at the surface of a colloidal inclusion with winding number k = +1. In Fig. 5, we show how switching between homeotropic and tangential anchoring at the surface of one inclusion alters the forces on each defect within the system.

In addition to forces, the director field also imparts torques on the particles. Torques reflect energy minimization by changing the boundary condition vector b in eqn (23). However, rotation of the boundary condition is nontrivially connected to rotation of particles. For instance, a particle with k = 1 is invariant to rotation, and for a particle with k = 0, rotation of the particle by some angle results in rotation of the boundary condition by the same angle. Considering the free energy expression in eqn (28), the torque on i-th particle reduces to

(33) |

With forces and torques computable for every configuration, one can write down effective dynamic equations

(34) |

Fig. 6 Annihilation and braiding of a ±1/2 defect pair with a fixed orientation. Defects have equal core sizes (ε_{1} = ε_{2} = ε). (a) Non-central forces lead to spiral annihilation trajectories. The initial separation is set to d_{12} = 20ε. (b) Defects with the same separation of d_{12} = 20ε are initiated without any spiral charge (μ_{1} = μ_{2} = 0). By braiding the defect pair with a fixed orientation in the counter-clockwise direction around each-other, the +1/2 defect obtains a spiral charge of −π/log(20) and the −1/2 defect obtains a spiral charge of π/log(20), in line with eqn (36). |

Spiral charge can be modified by not only changing the orientations of defects, but also by changing defect positions. Eqn (23) tells us that the orientation of a defect is determined from the contributions from all other defects through the boundary condition b_{i}:

(35) |

(36) |

We show that movement of defects with a fixed orientation is directly linked to changes of the spiral charge, similarly to how anyonic braiding was shown in k-atic liquid crystals by modulation of the boundary condition.^{50} Braiding of defects was observed also for spiral waves in living cells^{51} and in active nematics it was associated with topological entropy and chaos.^{52–54} Harmonic director fields are of great interest also in three dimensions.^{55} While our theory considers purely two-dimensional fields, a future challenge would be to address the defect textures also in quasi-two-dimensional layers, where the director is allowed to point out of plane.

Construction of director structures can be extended from simple rational expressions describing finite sets of defects, to any analytical functions with appropriate positioning of poles and zeros to account for infinite sets of defects, such as periodic defect arrays. The formalism remains the same, but the tractability of the elastic energy integrals depends on the function. For example, linear chains of dipoles can be mapped to trigonometric functions while doubly periodic lattices involve Jacobi or Weierstrass elliptic functions.

In this manuscript, we only considered point defect monopoles, which only includes complexified angle functions Φ with logarithmic singularities. Our formalism can be extended to include higher multipoles with singularities of the form z^{−} for -th multipole (see ref. 37), which do not require addition of spiral contributions. This would allow to account for particles with more elaborate director profiles, described through a multipolar expansion, at the expense of more convoluted expressions in the interaction matrix M and boundary constraints.

Our construction is not directly meant to solve the director field on domains that involve nontrivial boundary shapes or larger inclusions. For such problems, see other works that employ conformal mapping to transform the domain to one that offers a solution in terms of simple analytical functions – such solutions are limited to geometries that have an analytically tractable conformal mapping,^{27,30,31,56} or can otherwise be solved numerically.^{28} Our approach is a good approximation for most systems of small, preferably spherical particles. What is lost in exactness of boundary conditions, is offset by the fact that a closed-form solution can be obtained by a routine linear algebra algorithm for an arbitrary number of constrained boundaries. Here, we assume an infinitely strong anchoring, with director at the boundary prescribed exactly. For handling of finite anchoring at arbitrarily shaped boundaries, see Chandler and Spagnolie.^{33} Finally, the method's ability to handle larger numbers of defects suggests a possibility of developing a numerical scheme capable of handling finite-sized boundaries of arbitrary shapes, akin to methods used in hydrodynamics and electrostatics. We leave this as a future challenge.

(A1) |

(A2) |

The first term in eqn (A1) is directly set to zero once the zero total charge condition is applied. We can rewrite the second term of the free energy expression

(A3) |

(A4) |

(A5) |

(A6) |

(A7) |

For k = N − 1, we can construct a simplex with N vertices in a k-dimensional space, where the Euclidean distances between vertices correspond to . Note that this is possible because the condition d_{ij} ≥ 2ε for i ≠ j leads to non-negative values of and also to the triangle inequality . The k-dimensional volume of the simplex is always non-negative and is given by the Cayley–Menger determinant^{57}

(A8) |

An unbounded system does not have an outer boundary and thus also has no boundary condition given there. Instead, the sum of spiral charges must equal zero. To solve for this, we must make f_{0} an unknown quantity, specifying the far-field orientation, instead of the overall rotation ϕ_{0}, which can be eliminated as an unknown and can be computed from eqn (24) after the calculation. Using eqn (27) together with the zero total charge condition, which can be written as 1·μ = 0, we can write a modified version of eqn (25):

(B1) |

Repeating the block-wise inverse eqn (26), we obtain,

(B2) |

F = 2πkMk + 2πfM^{−1}f − 2πr(fM^{−1}1)^{2}. | (B3) |

- G. P. Alexander, B. G.-G. Chen, E. A. Matsumoto and R. D. Kamien, Rev. Mod. Phys., 2012, 84, 497 CrossRef CAS .
- C. Meng, J.-S. Wu and I. I. Smalyukh, Nat. Mater., 2023, 22, 64–72 CrossRef CAS PubMed .
- D. S. Kim, S. Čopar, U. Tkalec and D. K. Yoon, Sci. Adv., 2018, 4, eaau8064 CrossRef CAS PubMed .
- C. Peng, T. Turiv, Y. Guo, Q.-H. Wei and O. D. Lavrentovich, Science, 2016, 354, 882–885 CrossRef CAS PubMed .
- P. C. Mushenheim, R. R. Trivedi, H. H. Tuson, D. B. Weibel and N. L. Abbott, Soft Matter, 2014, 10, 88–95 RSC .
- S. Hernàndez-Navarro, P. Tierno, J. A. Farrera, J. Ignés-Mullol and F. Sagués, Angew. Chem., Int. Ed., 2014, 53, 10696–10700 CrossRef PubMed .
- M. Kim and F. Serra, Adv. Opt. Mater., 2022, 10, 2200916 CrossRef CAS .
- I. Nys, B. Berteloot and K. Neyts, J. Mol. Liq., 2023, 386, 122472 CrossRef CAS .
- J. Jiang, K. Ranabhat, X. Wang, H. Rich, R. Zhang and C. Peng, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2122226119 CrossRef CAS PubMed .
- M. Jiang, Y. Guo, R. L. Selinger, O. D. Lavrentovich and Q.-H. Wei, Liq. Cryst., 2023, 1–9 Search PubMed .
- S. Yi, H. Chen, X. Wang, M. Jiang, B. Li, Q.-H. Wei and R. Zhang, arXiv, 2023, preprint, arXiv:2312.14735 DOI:10.48550/arXiv.2312.14735.
- P. Pieranski and M. H. Godinho, Liquid Crystals: New Perspectives, 2021, pp. 193–309 Search PubMed .
- N. Sebastián, M. Lovšin, B. Berteloot, N. Osterman, A. Petelin, R. J. Mandle, S. Aya, M. Huang, I. Drevenšek-Olenik, K. Neyts and A. Mertelj, Nat. Commun., 2023, 14, 3029 CrossRef PubMed .
- J. Prost, F. Jülicher and J.-F. Joanny, Nat. Phys., 2015, 11, 111–117 Search PubMed .
- B. Halperin and D. R. Nelson, Phys. Rev. Lett., 1978, 41, 121 CrossRef CAS .
- M. Mur, Ž. Kos, M. Ravnik and I. Muševič, Nat. Commun., 2022, 13, 6855 CrossRef CAS PubMed .
- S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti and V. Vitelli, Nat. Rev. Phys., 2022, 4, 380–398 CrossRef .
- S. J. DeCamp, G. S. Redner, A. Baskaran, M. F. Hagan and Z. Dogic, Nat. Mater., 2015, 14, 1110–1115 CrossRef CAS PubMed .
- T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212–216 CrossRef CAS PubMed .
- F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic and A. R. Bausch, Science, 2014, 345, 1135–1139 CrossRef CAS PubMed .
- A. T. Brown, Soft Matter, 2020, 16, 4682–4691 RSC .
- Y.-H. Zhang, M. Deserno and Z.-C. Tu, et al. , Phys. Rev. E, 2020, 102, 012607 CrossRef CAS PubMed .
- L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. Cristina Marchetti, Philos. Trans. R. Soc., A, 2014, 372, 20130365 CrossRef PubMed .
- C. Pozrikidis, Introduction to theoretical and computational fluid dynamics, Oxford University Press, 2011 Search PubMed .
- R. P. Feynman, R. B. Leighton and M. Sands, The Feynman lectures on physics; New millennium ed., Basic Books, New York, 2010, vol. II: Mainly electromagnetism and matter.
- A. Belavin and A. Polyakov, JETP Lett., 1975, 22, 245 Search PubMed .
- A. J. Davidson and N. J. Mottram, Eur. J. Appl. Math., 2012, 23, 99 CrossRef .
- R. M. W. van Bijnen, R. H. J. Otten and P. van der Schoot, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 168 CrossRef PubMed .
- O. Tarnavskyy, M. Ledney and A. Lesiuk, Liq. Cryst., 2018, 45, 641 CrossRef CAS .
- O. S. Tarnavskyy, A. M. Savchenko and M. F. Ledney, Liq. Cryst., 2020, 47, 851 CrossRef CAS .
- O. S. Tarnavskyy and M. F. Ledney, Liq. Cryst., 2023, 50, 21 CrossRef CAS .
- H. Miyazako and T. Nara, R. Soc. Open Sci., 2022, 9, 211663 CrossRef PubMed .
- T. G. J. Chandler and S. E. Spagnolie, A nematic liquid crystal with an immersed body: equilibrium, stress, and paradox, 2023.
- V. Vitelli and D. Nelson, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 021711 CrossRef CAS PubMed .
- A. T. Brown, Soft Matter, 2020, 16, 4682 RSC .
- F. Vafa, M. J. Bowick, M. C. Marchetti and B. I. Shraiman, arXiv, 2020, preprint, arXiv:2007.02947 DOI:10.48550/arXiv.2007.02947.
- A. J. Houston and G. P. Alexander, New J. Phys., 2023, 25, 123006 CrossRef .
- J. Rudnick and R. Bruinsma, Phys. Rev. Lett., 1995, 74, 2491 CrossRef CAS PubMed .
- P. G. de Gennes and J. Prost, Physics of Liquid Crystals [PDF], Clarendon Press, 1993 Search PubMed .
- M. Kleman and O. Lavrentovich, Soft Matter Physics: An Introduction, Springer, 2003 Search PubMed .
- X. Tang and J. V. Selinger, Soft Matter, 2017, 13, 5481 RSC .
- X. Tang and J. V. Selinger, Soft Matter, 2019, 15, 587 RSC .
- D. R. M. Williams, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1994, 50, 1686 CrossRef CAS PubMed .
- A. H. Lewis, D. G. A. L. Aarts, P. D. Howell and A. Majumdar, Stud. Appl. Math., 2017, 138, 438 CrossRef .
- D. J. G. Pearce and K. Kruse, Soft Matter, 2021, 17, 7408 RSC .
- S. Shankar, S. Ramaswamy, M. C. Marchetti and M. J. Bowick, Phys. Rev. Lett., 2018, 121, 108002 CrossRef CAS PubMed .
- N. Kumar, R. Zhang, J. J. De Pablo and M. L. Gardel, Sci. Adv., 2018, 4, eaat7779 CrossRef CAS PubMed .
- K. Thijssen, M. R. Nejad and J. M. Yeomans, Phys. Rev. Lett., 2020, 125, 218004 CrossRef CAS PubMed .
- R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel and J. J. De Pablo, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, E124–E133 CAS .
- A. Mietke and J. Dunkel, Phys. Rev. X, 2022, 12, 011027 CAS .
- J. Liu, J. F. Totz, P. W. Miller, A. D. Hastewell, Y.-C. Chao, J. Dunkel and N. Fakhri, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2104191118 CrossRef CAS PubMed .
- A. J. Tan, E. Roberts, S. A. Smith, U. A. Olvera, J. Arteaga, S. Fortini, K. A. Mitchell and L. S. Hirst, Nat. Phys., 2019, 15, 1033–1039 Search PubMed .
- S. A. Smith and R. Gong, Front. Phys., 2022, 10, 880198 Search PubMed .
- K. A. Mitchell, M. M. H. Sabbir, K. Geumhan, S. A. Smith, B. Klein and D. A. Beller, Phys. Rev. E, 2024, 109, 014606 CrossRef CAS PubMed .
- J. Binysh and G. P. Alexander, J. Phys. A: Math. Theor., 2018, 51, 385202 CrossRef .
- F. Vafa, D. R. Nelson and A. Doostmohammadi, Phys. Rev. E, 2024, 109, 064606 CrossRef PubMed .
- D. M. Y. Sommerville, An Introduction to the Geometry of n Dimensions, Dover Publications, New York, 1958 Search PubMed .

This journal is © The Royal Society of Chemistry 2024 |