Analytical model for the motion and interaction of two-dimensional active nematic defects

Abstract

We develop an approximate, analytical model for the velocity of defects in active nematics by combining recent results for the velocity of topological defects in nematic liquid crystals with the flow field generated from individual defects in active nematics. Importantly, our model takes into account the long-range interactions between defects that result from the flows they produce as well as the orientational coupling between defects inherent in nematics. Our work complements previous studies of active nematic defect motion by introducing a linear approximation that allows us to treat defect interactions as two-body interactions and incorporates the hydrodynamic screening length as a tuning parameter. We show that the model can analytically predict bound states between two +1/2 winding number defects, effective attraction between two −1/2 defects, and the scaling of a critical unbinding length between ±1/2 defects with activity. The model also gives predictions for the trajectories of defects, such as the scattering of +1/2 defects by −1/2 defects at a critical impact parameter that depends on activity. In the presence of circular confinement, the model predicts a braiding motion for three +1/2 defects that was recently seen in experiments, as well as stable and ergodic trajectories for four or more defects.

---

1 Introduction

Topological defects are present in many physical systems, including solids, superfluids and superconductors, magnetic materials, complex fluids, soft and biological materials, and the cosmos.^1–6 Although the physical properties and scales of all of these systems vary widely, the progenitor of topological defects in each case is broken symmetry. Topological defects manifest as singularities in the broken symmetry order parameter, and, as such, take the form of quantized points, lines, surfaces, or hypersurfaces. The quantization of topological defects thus allows for an alternative picture of the structure and dynamics of these systems; namely, viewing topological defects as interacting particles. For some systems, such as vortices in type-II superconductors or Skyrmions in chiral magnets, the “particle model” approach has been successfully used to describe the behavior of the system.^7,8 In other systems, such as dislocations in solids and disclinations in liquid crystals, only recently have the required theoretical advances been made that permit defects to be treated as particles by relating them to gradients of the continuum order parameter.^9–16

Here, we focus on topological defects in two-dimensional active nematic liquid crystals. Active nematic liquid crystals, or simply “active nematics,” are fluids in which anisotropic constituents individually generate forces leading to collective macroscopic flows.^17,18 Active nematics are inherently out-of-equilibrium and exhibit chaotic flows on large scales.^19–22 They also exhibit spontaneous positive–negative defect pair production with self-advecting positive defects.^11,23,24

Much of the theoretical and computational work surrounding active nematics has focused on a continuum model for the time dependence of the nematic tensor order parameter Q and the fluid velocity field v.^17,18,25 There has been recent interest, however, in the dynamics of the topological defects themselves for their role in colloidal assembly, microfluidic cargo transport, morphogenesis, logic operations, and optimal mixing.^5,6,26–37 The behavior of defects is typically inferred from solutions to the continuum equations for the order parameter, but the nonlinearity of the time evolution equations limits the progress that can be made in obtaining analytical results, though recent progress has been made by using low-order expansions of the equations of motion of Q.^38,39

There have been considerable efforts to understand the flow fields produced by individual defects in active nematics, both with and without substrate friction.^{12,20,23,40–44} In many of these studies, defect motion is assumed to be dominated by the advection of the defect by the flows they produce. Additionally, because of the reliance on exact solutions for the fluid flows, understanding the motion of multiple interacting defects through them remains a challenge. Other recent studies of defect dynamics on curved surfaces have made the simplifying assumptions that positive defects are self-propelled, while negative defects are passive and only interact through Coulomb potentials and surface curvature coupling.^45,46 These models introduce noise (an experimental reality) and allow for independent orientations of defects, however, they do not include the coupling of defect orientations with relative positions, or the advection and shear flow interactions inherent in active nematics.^47,48 Recent work has also explored the effect of orientational noise on the unbinding of defects²⁴ and developed a hydrodynamic description of many interacting defects.^11,12,35 The hydrodynamic descriptions are suitable for systems with large numbers of active defects in the turbulent regime, but do not give descriptions of the trajectories of individual defects.

In this work, we combine recent analytical results for the flow field around a defect in an active nematic with an analytical framework to calculate the velocity of defects given the time dependence of the order parameter.^12,16,43 Using these elements, we posit an analytical model for the velocity of defects in an active nematic with many defects. A critical component to our model is a linear perturbation approximation for the Q-tensor order parameter near defect cores. Using this, we show that the defect velocity results from two-body interactions with other defects, and may be decomposed into three parts: the Coulomb interaction between defects, advection from fluid flows generated by other defects, and deflection by shear flows generated by other defects. We analytically probe the interaction between two defects and show that the model is able to predict phenomena observed in continuum simulations and experiments, and predicted by other analytical models such as bound states between positive defects and deflections of positive defects by negative defects.^12,49 We also show that the model predicts an orientation dependent, effective attraction between negative defects and that the critical unbinding length scales proportional to the active length, both of which have not been shown in previous analytical models to our knowledge. We study the predicted trajectories for multiple defects in a circularly confined system and identify several potential behaviors, such as stable configurations, periodic dynamics, ergodic trajectories, and defect braiding similar to that observed in recent experiments.³⁶ We also compare the analytical model to a continuum model for the case of defect coarsening in a confined system. Finally, we summarize our results and discuss potential future extensions of the model and challenges that remain. Our work complements the previous efforts to understand active nematic defect motion by explicitly including the hydrodynamic screening length as a tuning parameter, and allowing for the study of multiple interacting defects through a key approximation.

2 Velocity of active nematic defects

In this section, we derive an analytical expression for the velocity of defects in an active nematic in the presence of other defects. We present the primary equations in the main text while relegating some of the details to the Appendix.

2.1 Defect velocity in an arbitrary flow field

The local orientational order in a nematic phase may be characterized by a unit vector field called the director, n̂(r). Nematic phases have apolar symmetry, so −n̂ ≡ n̂. In two-dimensional nematics, topological defects are points of orientational singularity, and, hence, are points where n̂ cannot be defined. To identify disclinations, the winding number or “topological charge” is computed:where C is a closed loop in the nematic and ϕ is the angle of n with respect to an arbitrary reference axis. Since n̂ and −n̂ represent the same state, it is possible for defects of charge m = n/2, where n is an integer, to exist. The lowest energy nontrivial defects have charge m = ±1/2.⁵⁰ Here, we will only consider these half integer defects.

For systems with many defects, it is common to represent the nematic configuration using a rank-two tensor, Q. In two dimensions, Q may be parameterized as Q = S[n̂ ⊗ n̂ − (1/2)I] where S is a scalar representing the degree of local order. The tensor order parameter characterization has the benefit that topological defects are no longer singularities and are instead regularized so that S varies continuously and S → 0 at the defect cores. The apolar symmetry is also inherently present in Q due to its quadratic dependence on n̂.

Recently, an exact expression for the velocity of defects in terms of derivatives of Q was derived in both two and three dimensions.^12,16 The derivation, of which we give more details in the Appendix, is based on the conservation of topological charge for systems with a continuous, broken symmetry order parameter, first explored by Halperin⁵¹ and then later developed by Mazenko.^52–54 In two dimensions it reads

where ε is the two-dimensional antisymmetric matrix, ∂_t ≡ ∂/∂t, ∂_j ≡ ∂/∂x_j, and summation on repeated indices is assumed. Importantly, the velocity of the defect only depends on derivatives of Q at the location of the defect core, r = r₀. While eqn (2) is exact, Q may be complicated and non-linear, leading to difficulties in obtaining analytical results for defect velocities. Further, the time dependence of Q must also be known in order to use eqn (2).

To obtain an approximate analytical expression for the velocity of defects, we assume the Beris–Edwards model for the time evolution of Q in the presence of a fluid flow u:⁵⁵

where γ is a rotational viscosity andis a generalized tensor advection. Here E is the strain rate tensor, Ω is the vorticity tensor and λ is a dimensionless parameter that characterizes the tendency of nematogens to either tumble or align due to shear flow.⁵⁶ The free energy, F, is chosen so that its bulk density is analytic in Q (e.g. Landau-de Gennes⁵⁰) and has a one-constant elastic density f_E = L|∇Q|². For the remainder of the paper we work in dimensionless units by scaling lengths by the nematic correlation length and times by the nematic relaxation time τ = γ/(Cξ²), where k is a dimensionless parameter that determines the size of defects and relative strength of elastic forces and C is a characteristic energy density scale of the free energy.

Eqn (3) may be substituted into eqn (2) to yield an equation for the defect velocity in terms of spatial derivatives of Q only. To approximate spatial derivatives at the location of the defect core, we use a linear core approximation for a nematic defect in the presence of other defects.^13,16 In the absence of flow, u = 0, the velocity of the jth defect reduces to

where r_ij = r_j − r_i is the vector pointing from defect i to defect j. Eqn (5) is precisely the Coulomb interaction between defects in passive nematic liquid crystals.⁵⁰ In what follows, we will set k = 4 unless otherwise specified. In the presence of a flow field, u ≠ 0, eqn (2) gives two additional contributions to the defect velocity. The first is advection by the flow:

v^Advection_j = u(r = r_j). (6)

Eqn (6) shows that defects are simply advected by the flow velocity at the location of the defect core, which is an expected and intuitive result. The second contribution comes from the deflection of defects in shear flow and has a more complicated form. To simplify the notation, we introduce the unit vector

p̂_j = (cos 2φ_j, sin 2φ_j) (7)

where ϕ(x,y) is the angle of the director, n̂(x,y), with respect to the x-axis and the index j refers to the jth defect. Thus φ_j is the angle of the director field just outside the defect core of the jth defect in the +x direction and p̂_j gives information about the “orientation” of the jth defect. We also write the shear rate tensor aswhich we assume to be traceless since we will only consider incompressible flows. Then the last contribution to defect velocity iswhere E₁ and E₂ are computed at r = r_j. We note that eqn (10) depends on the charge of the defect and, in particular, the y-component of the velocity changes sign for different signed charges. For the remainder of the paper we will set λ = 1.

The defect velocity in the presence of flow is the sum of eqn (5), (6), and (10). In Fig. 1 we show an illustration of each contribution. In the next section we use recent analytical calculations of the flow field created by a defect in an active nematic to give an analytical expression for the advection and shear contributions to defect velocity.

Fig. 1 Illustration of defect motion in two-dimensional nematics. Coulomb forces result from elastic interactions while advection and shear interactions only occur when a flow (blue arrows) is present.

2.2 Flow from active nematic defects

In this section we review the analytical results from ref. 43 for the flow field around a defect in an active nematic and combine them with the results of the previous section to obtain an expression for the velocity of a defect in an active nematic. We consider solutions of the incompressible Stokes equation with active nematic stress:where η is the fluid viscosity, Γ is the friction coefficient with the substrate, p is the fluid pressure, and [small alpha, Greek, tilde] is the strength of the active force. Dividing eqn (11) by Γ yields the parameters α ≡ [small alpha, Greek, tilde]/Γ, which characterizes the magnitude of the active force and hence determines the magnitude of the flows, and , which is the hydrodynamic screening length that characterizes the spatial extent of the flows. We note that (as is done in ref. 43) we have neglected the effect of elastic stresses in eqn (11). To include these terms in our model, we would require an analytical solution to eqn (11) that includes them, which, to our knowledge, has not been found. Nevertheless, it has been shown that the active stress dominates elastic stresses in other solutions of the Stokes equation for nematic defects.^23,40,41 In a passive nematic, it would likely be much more important to include these terms, as the elastic backflow has been shown to alter the interaction of defects when active stresses are not present.⁵⁷

We will assume that the total flow at the core of a defect is the sum of all flow fields created by defects. We justify this assumption by first noting that the director field in a system of many defects is determined by the Euler–Lagrange equation ∇²ϕ = 0. This equation is linear, so the full solution is the sum of each individual defect solution. In particular, we will assume the director structure is given by the minimal energy solution, so that there are no torques on defects due to elastic forces.^9,58 Physically, this corresponds to assuming that the nematic relaxation time is much smaller than the active force time scale. The flow field created by a defect is sourced by the term α∇·Q in the rescaled eqn (11). Away from the core of other defects, this term depends only on the director, and, under our linear core assumptions, is given by the sum of all other defects. We thus approximate the total active force by the sum of active forces produced by each defect. This approximation is better when defects are well separated compared to their core size, and will begin to break down as defect cores overlap. Nevertheless, these approximations and assumptions allow us to make analytical progress in treating interactions between multiple active nematic defects, since they allow us to treat the velocities as arising from two-body interactions.

In active nematics, +1/2 defects are much more motile than −1/2 defects. This is because they generate a nonzero flow at their core⁴³

On the other hand, −1/2 defects generate zero flow at their core, u⁻(0) = 0, though they still move due to interactions with other defects and flows. The flow field away from the core of ±1/2 defects is given by⁴³where r and θ are polar coordinates and φ gives the orientation of the defect, given by eqn (8). The functions β_i are plotted in Fig. 2(a) and the functional forms are given in the Appendix. As shown in Fig. 2(a), at large arguments β_i(x) ∼ 1/x, while at smaller arguments the behavior is more complicated and even nonmonotonic in the case of β_2–4.

Fig. 2 (a) Radial parts β₁, β₂, β₃, and β₄ of the flow generated by active nematic defects plotted vs. r/2[small script l]_h. Inset: Log–log plot of β₁vs. π/4 + r/2[small script l]_h. (b) Radial parts χ₁, χ₂, and χ₃ of the shear rate generated by active nematic defects plotted vs. r/2[small script l]_h.

E ₁ and E₂ may be computed from gradients of the velocity. For ±1/2 defects these are given by

where
are plotted in Fig. 2(b). For large arguments, χ_i(x) ∼ 1/x², and hence the shear flow interaction will be subdominant when distances are much larger than the hydrodynamic screening length. This is the intuitive result for “dry” active nematics where [small script l]_h → 0. In this case, the self-propulsion of +1/2 defects dominates other terms in the velocity.

We note that the limit [small script l]_h → ∞ is not a sensible limit to take in the model. This is because, in this limit, the substrate friction goes to zero and solutions for the two-dimensional Stokes equation in an infinite domain diverge. This does not contradict experiments because substrate friction always plays a role, or in continuum simulations that set Γ = 0 since the finite size of domains always acts as an effective screening length on the system. A more sensible limit to examine in our model is how the flow solutions behave in the limit r ≪ [small script l]_h. While the form of β_i and χ_i is very complicated in this limit, we use the first two terms of β₁ (given in the Appendix) to estimate that β₁(x) ∼ (π²/4)/(π/4 + x) for small x. Indeed, for small r/2[small script l]_h we find this to be a satisfactory approximation, as we show in the inset of Fig. 2(a).

The velocity of an individual defect is found by treating all other defects as point sources for the flow at the location of the defect of interest. Due to flows throughout the nematic, the actual nematic configuration may be more complicated; however, we do not consider this here and assume that all flows in the system are due to topological defects. In the limit of small relaxation time compared to the timescale for active forces, this is a good approximation, since any non-topologically protected deformation will quickly be relaxed.

3 Interaction between two active nematic defects

Here we analyze the interaction between active nematic defects with the simplifying assumption that only two defects exist in an infinite domain. In each case we set our coordinates so that the defect of interest is at polar coordinates (r,θ) and the other defect is at (0,0). We will also assume “optimal” orientation between defects.⁹ This set of orientations minimizes the elastic free energy of the configuration. Further, given a global phase φ₀, the orientation of each defect may be determined by the defect positions:where i₊ indexes positive defects, i₋ indexes negative defects, and φ_j− may be found similarly by summing over negative defects i₋ ≠ j₋.

3.1 Two +1/2 defects

Given the approximations and assumptions of the previous section, the velocity of a +1/2 defect in the presence of another +1/2 defect iswhere r̂ is the unit vector pointing to the defect of interest from the other defect, â = [cos(θ − 2φ₀),sin(θ − 2φ₀)], and b̂ = [cos(θ − 4φ₀),sin(θ − 4φ₀)]. Components of eqn (22) proportional to β_i result from advection while components proportional to χ_i result from the interaction with shear flows. We note that the other defect will have exactly the opposite velocity of the defect under consideration.

Eqn (22) gives several insights into the motion of two +1/2 active defects. We first note that the self-propulsion velocity is reduced compared to an isolated defect, particularly at small distances relative to the screening length (note that β₁ → π as r → 0). The motion is complicated by terms proportional to vectors â and b̂ which depend on the relative orientation of the defects. However, we may analyze this coupling by decomposing the velocity into terms parallel and transverse to the direction between defects:

where we have suppressed the arguments of β₁,β₂,χ₁,χ₂ for brevity. Eqn (23) and (24) are independent of θ, the angle between defects, and the angular dependence is given by the fixed global phase φ₀. Thus, since v_⊥ is proportional to the change of angle between defects, the model predicts that defects will rotate around one another as long as v_⊥ ≠ 0. We note that the independence of θ in the defect dynamics is a consequence of the unique rotational symmetry of the far-field director for total topological charge +1. As we will show in later sections, other defect interactions do not share this property.

Due to the complicated forms of the flow velocity, it is not immediately obvious whether defects will tend to repel (as they do in passive nematics) or attract. Eqn (23) predicts that both behaviors are possible, with φ₀ = π/4 acting as a separatrix for the two behaviors. In Fig. 3(a) we plot v_‖ as a function of r/[small script l]_h for φ₀ = π/4 and [small script l]_h = 10 at a variety of activity strengths. For α/[small script l]_h² > 1.74 there is a distance r₀ such that v_‖ = 0, indicating that defects will tend to rotate around each other at this distance indefinitely. In Fig. 3(b) and (c) we plot the critical distances as a function of α/[small script l]_h² (b) and φ₀ (c). For φ₀ < π/4 and , r₀ diverges at zero activity (as expected), but decreases rapidly with increasing activity. For α/[small script l]_h² > 1, the critical distance remains roughly constant and small compared to [small script l]_h. As a function of φ₀, r₀ remains roughly constant until φ₀ ≈ 7π/32, and then increases rapidly as φ₀ approaches π/4.

Fig. 3 (a) Velocity component of a +1/2 defect in the presence of another +1/2 defect parallel to the direction between defects, v_‖⁺ [eqn (23)], versus scaled distance between defects r/[small script l]_h. Colors represent different values of α/[small script l]_h² with [small script l]_h = 10, θ = 0, φ₀ = π/4, and λ = 1. The dashed line shows v_‖⁺ = 0; r₀ is defined as the point at which the velocity crosses this line. (b) Scaled critical distance r₀/[small script l]_hversus α/[small script l]_h² for various values of φ₀. (c) r₀/[small script l]_hversus φ₀ for various values of α/[small script l]_h². (d) Example trajectories for two +1/2 defects with initial positions at x = −0.5[small script l]_h and x = 0.5[small script l]_h for α/[small script l]_h² = 0.1, φ₀ = 5π/24 (blue) and α/[small script l]_h² = 0.25, φ₀ = π/8 (orange). (e) Distance between two +1/2 defects versus time step t/Δt for the parameters in (d) as well as at α/[small script l]_h² = 0.1, φ₀ = 5π/16 (yellow).

The behavior of r₀ may be attributed primarily to the self-advection part of eqn (23), which is the dominant term in the equation. It does not depend on the distance between defects, and it contributes most strongly at high activities and smaller φ₀. However, we note that the shear flow interaction contributes most strongly to the behavior near φ₀ = π/4 where the advection terms do not contribute to eqn (23). If λ = 0, r₀ → ∞ for all values of α at φ₀ = π/4.

In Fig. 3(d) we also plot trajectories of defects obtained from numerically integrating eqn (22) using a forward Euler method with timestep Δt = 0.001 and [small script l]_h = 10. In both computations the defects are initially placed at (−0.5[small script l]_h,0) and (0.5[small script l]_h,0). For one trajectory, we set α/[small script l]_h² = 0.1 and φ₀ = 5π/24. In this case, r₀ is larger than the initial separation, and the defects initially move away from each other to the stable distance. For the other trajectory, we set α/[small script l]_h² = 0.25 and φ₀ = π/8, corresponding to a critical distance that is smaller than the initial defect separation. Here the defects move toward each other in order to reach the stable distance. In Fig. 3(e) we plot the distance between defects for the computed trajectories illustrated in Fig. 3(d), as well as for trajectories computed for α/[small script l]_h² = 0.1 and φ₀ = 5π/16 > π/4. In this case, while the defects still rotate around one another, the distance between them increases without bound.

The effective attraction between positive defects caused by hydrodynamic effects can lead to dynamical states in which +1/2 defects continuously rotate around one another, or even merge to form a stable higher order +1 defect, similar to behaviors predicted with a different analytical model in ref. 49. While our calculations above predict this attraction for the simplest case of a system with only two defects, this behavior should, in principle, occur even when other defects are present in the system if the two defects are far from the other defects. Indeed, bound states of two +1/2 defects may be observed in continuum simulations at high activities, as illustrated in Fig. 4(a) where we show a time snapshot of the scalar order parameter and director field in a system with periodic boundary conditions, α/[small script l]_h² = 2, and [small script l]_h = 10 (for details of the continuum simulations, see the Appendix). We have marked the bound states of +1/2 defects in green boxes in Fig. 4(a), and for one of these boxes the subpanels provide several time snapshots illustrating two bound +1/2 defects rotating around one another and eventually merging. We note that the bound states are short-lived in the continuum simulation due to interactions with other defects. However, it has been shown in numerical simulations that such bound states may be longer lived and ordered if subjected to confinement.^59–61 Partially stable bound +1/2 defect pairs have also been observed experimentally under circular confinement⁶² and in recently discovered acoustically actuated active nematics with no confinement, where the activity may be much larger than in biologically based active nematics.⁶³

Fig. 4 (a) Time snapshot of an active nematic continuum simulation with periodic boundary conditions, α/[small script l]_h² = 2, and [small script l]_h = 10. The color represents the strength of the scalar order parameter S and the white lines indicate the local director. The green boxes highlight bound states of two +1/2 defects rotating around one another. The sub panels show four time snapshots of one of the bound states at t = −4Δt, t = −2Δt, t = 0, and t = 2Δt, where t = 0 is the time in the main snapshot. The circles indicate the +1/2 defects that rotate around one another. (b) Three time snapshots of a confined continuum simulation with only two +1/2 defects with α/[small script l]_h² = 0.2, [small script l]_h = 10, φ₀ = π/8, and domain radius R = 30. (c) Scaled stable distance between defects r₀/[small script l]_h in a continuum simulation for various values of α/[small script l]_h² with [small script l]_h = 10 and φ₀ = π/8 (orange points). The dashed line is a reproduction of the analytically derived curve from eqn (23) (solid line), with a vertical shift of 1.2.

By simulating a continuum system under circular confinement, we can perform a rough comparison to the scaling of r₀ with α shown in Fig. 3(b). In Fig. 4(b), we plot several time snapshots of a continuum simulation with α/[small script l]_h² = 0.2, [small script l]_h = 10, φ₀ = π/8, and a confining radius of R = 30. In the continuum simulations, two defects rotate around one another at a fixed distance. As predicted by the analytical equations, this distance depends on the activity. In Fig. 4(c) we plot the critical distance r₀versus α/[small script l]_h² for several continuum simulations. We also plot the analytically derived curve from Fig. 3(b), r₀(φ₀ = π/8)/[small script l]_h + 1.2, where the vertical shift in the curve is made to highlight the similarities and differences in the shape of the relationship. Both the continuum simulations and analytic model show a relatively rapid drop in r₀ for small α; however, at larger activities r₀ falls off more rapidly with α in the continuum simulations than in the analytic model, which we attribute to nonlinear effects produced by the close proximity of the defect pair in the continuum simulations. When the defects become close enough together that their cores begin to overlap, our approximations break down. In particular, in the continuum model, defects are able to merge to form higher order defects, as is visible in several of the boxed regions in Fig. 4(a). We surmise that the quantitative difference (the vertical shift between the analytically predicted behavior and the continuum simulations) arises due to effects of the boundaries on the continuum simulations. We use u = 0 boundary conditions for the flow field in the confined continuum simulations, but eqn (22) corresponds to defects in an infinite, not bounded, domain. These boundary conditions reduce the flow magnitude and spatial extent, and, hence, increase the critical distance between defects. The vertical shift factor, here found to be ≈1.2, likely depends on both [small script l]_h and the size of the computational domain, though a comprehensive study of this effect is outside the scope of the present work. Additionally, because defects can spontaneously nucleate for higher activities in the continuum simulations, we cannot probe the high activity scaling result. In the lower activity regime, however, we obtain good qualitative agreement between our analytic model and the behavior of defects in continuum simulations and experiments.

3.2 Two −1/2 defects

The velocity of a −1/2 defect in the presence of another −1/2 defect iswhere ĉ = [cos(3θ − 2φ₀), −sin(3θ − 2φ₀)], d̂ = [cos(5θ − 2φ₀),sin(5θ − 2φ₀)], and f̂ = [cos(7θ − 4φ₀), −sin(7θ − 4φ₀)]. As with the two +1/2 defects, it is instructive to decompose the velocity into components parallel and transverse to the direction between the defects:

Eqn (27) shows that the transverse component of the velocity will be zero if 4θ − 2φ₀ = nπ. In this case, two −1/2 defects can only move toward or away from one another. From eqn (26), we find that if 4θ − 2φ₀ = (2n + 1)π, the defects repel each other. Geometrically, this occurs when two −1/2 defects “point” towards one another if viewed as having a triangular shape. On the other hand, when the defects point away from one another for 4θ − 2φ = 2nπ, the hydrodynamic interactions compete with the Coulomb force on the defects and there is a critical distance r₀ at which v_‖ = 0. In Fig. 5(a) we plot this critical distance as a function of α/[small script l]_h² for [small script l]_h = 10. For two −1/2 defects, r₀ rapidly decreases as α/[small script l]_h² increases until α/[small script l]_h ≈ 1 and the value of r₀ saturates. As we did with the +1/2 defects, in Fig. 5(a) we plot r₀/[small script l]_h for several confined continuum simulations with two −1/2 defects. The dashed line in Fig. 5(a) is the solid curve in the figure with a vertical shift, r₀/[small script l]_h + 2, showing that apart from the vertical shift, the finite system continuum simulations behave qualitatively similarly to the curve derived analytically for an infinite domain. As with the +1/2 defects in the previous section, we suspect the quantitative difference to result from the confinement of the continuum simulations, where the flow speeds are reduced even more than in the case of two +1/2 defects. We expect that the value of the vertical shift, found to be 2 for the parameters used here, should be dependent on [small script l]_h and the domain size of the confined system.

Fig. 5 (a) Scaled critical distance r₀/[small script l]_h at which the velocity of −1/2 defects v = 0 for 4θ − 2φ₀ = 0 versus α/[small script l]_h². The solid line shows the calculated r₀/[small script l]_h from eqn (26). Orange points show r₀/[small script l]_hversus α/[small script l]_h² for two −1/2 defects in confined continuum simulations with [small script l]_h = 10 and domain radius R = 30. The dashed curve is a reproduction of the solid curve with a vertical shift of 2. (b) Example trajectories of two −1/2 defects with initial x = −0.5[small script l]_h and x = 0.5[small script l]_h at φ₀ = 0.001.

Although a critical distance between −1/2 defects exists, the velocity equations do not support stable, bound defect pairs. This is because the dynamics are explicitly dependent on θ, and, in particular, v_⊥ changes as the defects rotate around one another. Further, expanding around the fixed points of eqn (27) shows that 4θ − 2φ₀ = 2nπ are unstable fixed points, while 4θ − 2φ₀ = (2n + 1)π are stable fixed points. The model thus predicts that active −1/2 defects will tend to reorient toward a stable relative orientation when interacting. We illustrate this instability in Fig. 5(b) by plotting the computed trajectories of two −1/2 defects with initial positions (−0.5[small script l]_h,0) and (0.5[small script l]_h,0), for φ₀ = 0.001, [small script l]_h = 10 and α/[small script l]_h² = 0.5. Initially the defects move towards each other until they rotate enough for the interaction to become repulsive. They then repel each other along a line such that asymptotically approaches 4θ − 2φ₀ = π. This instability provides an explanation for why bound states of −1/2 defect pairs are not observed in continuum simulations and experiments.

3.3 ±1/2 defects

Unlike the previous two cases we have examined, ±1/2 defect pairs interact non-reciprocally.¹² That is, v⁺⁻ ≠ −v⁻⁺ and, thus, the predicted trajectories will not be symmetric. The velocities of the defects are given bywhere â* = [cos(θ − 2φ₀), −sin(θ − 2φ₀)], ĉ* = [cos(3θ − 2φ₀), sin(3θ − 2φ₀)], ĝ = [cos(5θ − 4φ₀), sin(5θ − 4φ₀)], ĥ = [cos(3θ − 4φ₀), −sin(3θ − 4φ₀)], and we note that we have written the velocities in coordinates where the origin is located at the other defect.

We first analyze the interaction by examining the limiting case where the global phase is φ₀ = π/2 and the axes are oriented so that θ = 0 for the +1/2 defect and θ = π for the −1/2 defect. In this case the comet head of the +1/2 defect is pointing “away” from the −1/2 defect [see Fig. 6(a)]. The velocity difference, and, hence, the rate at which the distance between defects changes, is

Given a value for α and [small script l]_h, eqn (30) predicts that there will be a critical distance r_C above which the defects will move away from one another and below which the defects will inevitably annihilate. We argue that in a large active nematic system with many defects, the average distance between defects 〈r_d〉 ∼ r_C since this is the distance at which +1/2 defects may unbind from −1/2 defects.

Fig. 6 (a) Scaled critical unbinding distance r_C/[small script l]_hversus α/[small script l]_h² for [small script l]_h = 10 and φ₀ = π/2. The inset shows a schematic of the orientation of the defects. (b) r_C/[small script l]_hversus α/[small script l]_h² for various screening lengths [small script l]_h plotted on a log–log scale. The lower dashed line indicates x^−1/2 scaling while the upper dashed line indicates x⁻¹ scaling.

In Fig. 6(a) we plot r_C/[small script l]_hversus α/[small script l]_h² for [small script l]_h = 10. Here r_C → ∞ as α → 0 as expected, since passive ±1/2 defect pairs will always annihilate. As shown by the log–log plot in Fig. 6(b), we find that r_C ∼ α^−1/2, a scaling that has been observed for 〈r_d〉 in continuum simulations.^20,64 In some sense, this scaling is expected because the active length scale (the length scale that arises due to the competition between active and elastic forces¹⁸) [small script l]_a ∝ α^−1/2; however, a naive approximation would lead to the expectation that dr/dt ∼ v − 1/r, where v is the speed of the +1/2 defect, and hence that r_C ∼ α⁻¹. Only by including all of the hydrodynamic effects can the α^−1/2 scaling be predicted. Interestingly, we find that the screening length [small script l]_h plays an important role in this scaling. In Fig. 6(b) we also plot r_C/[small script l]_hversus α/[small script l]_h² for [small script l]_h ∈ {1,2,5}. For small screening lengths and small activity, r_C ∼ α⁻¹, confirming the naive approximation. In this regime, the self propulsion is the dominant interaction term, and only at higher activities is the α^−1/2 scaling restored. The crossover between regimes as a function of α occurs for higher relative activities at smaller screening lengths. Thus, we predict that active nematic experiments and simulations performed in the “dry” limit (smaller screening length) for relatively small activities will measure 〈r_d〉 ∼ α⁻¹.

To gain some analytical insight we use the approximation of β₁ in the limit of r_C ≪ [small script l]_h [see Fig. 2(a) and the Appendix]. Setting β₁(r_C/2[small script l]_h) = (π²/4)/(π/4 + r_C/2[small script l]_h) in eqn (30) and neglecting the terms which are small in this limit gives

where we have written the right hand side explicitly in terms of α/[small script l]_h² to match the scaling in Fig. 6. From this approximation, we see that if either [small script l]_h or α is large then we recover r_C ∼ (α/[small script l]_h²)^−1/2, while if both [small script l]_h and α are small r_C ∼ (α/[small script l]_h²)⁻¹. Thus, both limits are able to be captured by this approximation, and we emphasize that the full model can interpolate between these two regimes.

One might predict that if the comet head of the +1/2 defect is oriented toward the −1/2 defect, then annihilation will always occur. This is precisely true if 2θ − 2φ₀ = 2nπ; however, due to the orientational coupling and self propulsion of the +1/2 defect, if 2θ − 2φ₀ ≠ 2nπ, the +1/2 defect may deflect away from the −1/2 defect. In this more general case, the complicated non-reciprocity of the ±1/2 defect velocities does not lend itself to a straightforward analytical analysis. Instead, we analyze the interaction by numerically integrating eqn (28) and (29) and observing the resulting trajectories. In Fig. 7(a) we plot an example trajectory with φ₀ = π/12. In this case, although the ±1/2 defect pair distance initially decreases, eventually it begins to grow with time, as can be seen by the three snapshot points plotted along the trajectory. We note that although −1/2 defects are typically thought of as “passive” particles, their motion is not only influenced by Coulomb interactions and can be pushed or pulled by the flows of other defects, which can be seen in the trajectory in Fig. 7(a).

Fig. 7 (a) Example trajectories of ±1/2 defects with initial positions y⁺ = 0.5[small script l]_h and y⁻ = −0.5[small script l]_h, φ₀ = π/12, and [small script l]_h = 10. The points are time snapshots of the positions of the defects at t/Δt = 0, t/Δt = 5000 and t/Δt = 10 000. The arrows on the points show the orientation of the defect. (b) Critical global phase φ_C above which +1/2 defects deflect from −1/2 defects versus α/[small script l]_h² for [small script l]_h = 10 at various initial defect separations r₀/[small script l]_h. Inset: The same data with the vertical axis scaled by (r₀/[small script l]_h)². (c) φ_Cversus r₀/[small script l]_h for various α/[small script l]_h². Inset: The same data with the vertical axis scaled by α/[small script l]_h².

We find that, if defects are initially oriented along the x-axis, there is a critical global phase, φ_C, depending on initial separation and activity, above which +1/2 defects will deflect and not annihilate. To study this systematically, we compute trajectories of ±1/2 defect pairs with varying initial distances and orientations. We consider the two defects “annihilated” if they come within one nematic correlation length of each other. In our normalized units this occurs when r < 1. In Fig. 7(b) we plot the critical angle φ_Cversus α/[small script l]_h² for [small script l]_h = 10 at various initial distances, and in Fig. 7(c) we plot φ_Cversus r₀/[small script l]_h for various activities. The insets of these figures show that φ_C ∼ r₀⁻²α⁻¹, or φ_C ∼ ([small script l]_a/r₀)². We can also map this critical angle to a critical impact parameter b_C by rotating our axes by 2φ_C (so that φ₀ → 0). Then, b_C = r₀ sin 2φ_C so that b_C ≈ 2[small script l]_a²/r₀. We sketch the rotation leading to this mapping in Fig. 8.

Fig. 8 Schematic depiction of rotating the system axes to define the impact parameter b in terms of φ₀.

Our analysis of the model's predictions for the critical unbinding lengths and angles of ±1/2 annihilation do not assume a particular method of nucleation. One future extension of the model would be to implement specific nucleation methods. However, once the defects have nucleated, our analysis should apply to situations where the ±1/2 pair are well separated or screened from other defects in the system. Therefore, we expect our analysis to apply to defects nucleated from the usual bend instability,¹⁸ as well as defects nucleated via recently observed disorder–order transitions in phototactic active nematics.⁶⁵

4 Multiple defects in circular confinement

We next consider the computation of the trajectories of many defects in a confined circular geometry. Confining the space in which the particles are allowed to move better facilitates comparisons between particle-based simulations and experiments or continuum simulations. The model presented in Section 2 could be used to simulate many defects in an infinite system. However, because +1/2 defects are able to self propel, these defects will eventually move off to infinity.

To simulate a confined circular geometry, we use the method of images for defects in a circular domain. For each defect located at polar coordinate (r,θ), a defect of equal charge is placed at (R²/r,θ), where R is the radius of confinement. The corresponding lowest energy director angle is given by

where the polar angles of each defect θ_i are subtracted so that φ₀ gives the director angle offset at R. Eqn (32) describes a director field which will always have a fixed director at R. If the total topological charge of the system is zero, then ϕ(r = R) = φ₀, a constant. The total topological charge does not have to be zero, however. Eqn (32) will give the lowest energy director regardless of the total topological charge. We can then obtain the orientation of each defect using .

In confined systems, it is typical to assume no-slip boundary conditions for the fluid velocity. However, to our knowledge, no analytical solution to eqn (11) is available for u(R) = 0 boundary conditions. Thus, here we neglect the effect of the boundaries on the flow solution, and limit our considerations to the case of small screening length, namely [small script l]_h = 1. In this parameter regime, it was shown in ref. 43 that the flow solutions for a single defect are only weakly perturbed by no slip boundary conditions when the analytical solution for an infinite domain was compared to bounded numerical solutions. Finally, we assume that defects in the system do not feel the flow solutions from the image charges but experience only the Coulomb force from elastic forces. For the parameters that we explore, this assumption only results in small quantitative changes to the defect motion. In what follows, we integrate the equations of motion for each defect using a forward Euler method with time step Δt = 0.001 and initialize the system with randomly generated positions of the defects.

We first examine multiple +1/2 defects in confinement. For two +1/2 defects, the results of Section 3 apply and the defects will rotate around one another if φ₀ ≠ 0,π/2. If φ₀ = 0,π/2, the defects will find stable positions with v = 0. For three +1/2 defects, we find that the defects tend to braid one another in a manner similar to that seen in recent experiments.³⁶ In Fig. 9(a) and (b) we show example trajectories for three defects with system radii R = 3 and R = 5, respectively. For the smaller system size, the defects maintain tighter “figure eight” trajectories, as seen in the experiments, until eventually settling into the stable figure eight trajectory shown by the dashed line in Fig. 9(a). For the larger system size the defects do not follow perfect braiding trajectories, but instead traverse figure eight pathways bounded by outer and inner figure eight loops. Here we do not observe the same stabilization of the trajectories that appears in the more confined system. Movies of the integrated analytical model are available in the Supplemental Material (1,2) (ESI†). The ability of our model to predict the braiding behavior seen in experiments highlights the importance of orientational coupling in the active nematic. Such behavior could not be reproduced with a simpler active Brownian particle model where only Coulomb interactions are present.

Fig. 9 (a) Example trajectories given by the analytical model for three +1/2 defects in a confined system with radius R = 3 and α = 2. Different colors depict different defects. The dashed line in the lower figure depicts the late-time stable trajectory followed by all three particles. (b) Example trajectories for three +1/2 defects with R = 5 and α = 2. The lines outside the circle depict the director field at the boundary, which also apply to the system in (a). (c) Example trajectories for four +1/2 defects with R = 3 and α = 2. The black dots indicate the stable final positions of the defects. (d) Example trajectories for four +1/2 defects with R = 5 and α = 2. (e) Example trajectories for five +1/2 defects with R = 3 and α = 3. No stable points or trajectories were observed for the number of integration steps we considered. (f) Example trajectories for five +1/2 defects with R = 5 and α = 3. (g) Example trajectories for six +1/2 defects with R = 3 and α = 3. The dashed lines in the lower figure depict the late-time stable trajectories followed by three particles each. (h) Example trajectories for six +1/2 defects with R = 5 and α = 3.

We also examine trajectories for four, five, and six defects. We plot example trajectories for all of these cases in Fig. 9(c)–(h) and show movies of their motion in the Supplemental Material (3–8) (ESI†). For four +1/2 defects we find that there is a stable configuration where the defects sit in a rectangular formation and v = 0 for each defect. At small confinement the system rapidly relaxes to this configuration regardless of the random initial positions, as shown in Fig. 9(c). For larger confinement, the defects eventually find the stable configuration, but only after a long transient motion of two bound states of two +1/2 defects that rotate around one each other, as illustrated in Fig. 9(d). We find that stable configurations exist for four +1/2 defects for all values of φ₀. We note that, unlike the case for two +1/2 defects, the value of φ₀ has no influence over the existence of a stable state. This is because a total topological charge of m = 1 is a special case since it is the only arrangement for which the relative angle of the director with respect to the boundary is constant. For any other total topological charge, the director angle relative to the boundary changes continuously, and hence any rotation of the director rotates the overall configuration. To emphasize this, in Fig. 9(b), (d), (f) and (h) we indicate the director boundary condition. We also note that although the flow solutions we assume here do not respect no-slip boundary conditions, the symmetry of the configuration dictates that a qualitatively similar configuration should still exist for the case of no-slip conditions.

For five defects we no longer find stable defect configurations. In Fig. 9(e) and (f) we plot example trajectories for five +1/2 defects and do not observe any closed defect trajectories, meaning that there is not a stable attractor for either system size. For six defects, when R = 3 we observe two periodic trajectories that contain three defects each, as shown in Fig. 9(g). The periodic trajectories may change their orientation for different random initial conditions or values of φ₀. Upon increasing the system size to R = 5, however, we no longer observe stable periodic trajectories, as illustrated in Fig. 9(h). We note that while it is beyond our scope here to extensively study the effect of α on these trajectories, we find that the qualitative results do not change upon a cursory survey of varying activities. That is, only the quantitative details of the points and trajectories vary when small changes are made to the activity, but the existence of stable points or trajectories is unaffected.

Finally, we compare active nematic defect coarsening in the analytical model and a continuum simulation. For the analytical model, we initialize a system of radius R = 50 with 15 + 1/2 defects and 15 − 1/2 defects at random positions. We then numerically integrate their velocities and remove ±1/2 pairs that come within r < 1 of each other. For the continuum simulation, we initialize the system with a director field that has defects at the same positions as the analytical model. For both computations, we use α/[small script l]_h² = 2, [small script l]_h = 1, k = 1, and φ₀ = 0. Importantly, for these parameters, the continuum simulation does not spontaneously nucleate defects. In Fig. 10(a) we plot the number of defects N versus dimensionless simulation time t for both the computations with logarithmic scaling on the t-axis. At early times, when the system is more dense, the coarsening dynamics are different, with faster annihilation of defects in the continuum simulation. However, at longer time scales, the two curves have roughly logarithmic scaling with time. We also find that at late times, +1/2 defects move out toward the boundaries while −1/2 defects are left in the bulk. The +1/2 defects congregate at polar angles θ = φ₀ ± π/2 due to the orientational coupling between the defects and director field. In Fig. 10(b) we plot a late-time snapshot of an analytically integrated computation while in Fig. 10(c) we plot a late-time snapshot of a continuum configuration. In both cases the charge separation occurs while the defects continue to coarsen. While a full understanding of confined active nematic defect coarsening would require many more realizations and an extensive study of the effect of model parameters, this example highlights a potential future application of the analytical model, which requires far less computational time and resources to simulate than the continuum model.

Fig. 10 (a) Comparison of the number of defects N versus time t for the numerically integrated analytical model and a continuum simulation with R = 50, α = 2, [small script l]_h = 1, k = 1, and φ₀ = 0. The analytical model was initialized with random positions of 15 + 1/2 and 15 − 1/2 defects, while the continuum simulation was initialized with a director field that has the same random positions of defects used in the analytical model. Time is given in the dimensionless units outlined in Section 2. (b) Late-time configuration snapshot of the numerically integrated analytical model. (c) Late-time configuration snapshot of the continuum model +1/2 defects are circled in orange while −1/2 defects are circled in blue.

5 Discussion

As we have demonstrated, our analytical model qualitatively reproduces features of active nematic defect interactions observed in both continuum simulations and experiments. It is also able to reproduce predictions from previous analytical models of defect motion, such as the mutual rotation of two +1/2 defects, and the deflection of +1/2 defects by −1/2 defects.^12,49 By including a well-defined screening length, the model can interpolate between wet and dry regimes of active nematics. Additionally, because we employ a linear perturbation approximation, the model may be used to investigate the motion of multiple defects interacting with one another.

One of the primary advantages of the analytical model is the ease with which the effect of different parameters may be investigated. For example, in future studies, the model can be used to probe the effect of screening length and tumbling parameter. Another advantage is the modular nature of the theory. One may, for example, investigate the effect of elastic interactions by changing the coefficient for the Coulomb force term; or one could impose different flow profiles by adding the appropriate terms to the velocity. Since the velocity of defects is derived explicitly from topological charge conservation via the methods of Halperin and Mazenko,^51–54 we do not need to interpret the results in terms of auxiliary fields.^38,39 Having an explicit form for the defect velocities is particularly amenable for computational implementation.

Our analytical model achieves its greatest accuracy under the separation of length scales a < 〈r_d〉 < [small script l]_a, where a is the core size, 〈r_d〉 is the average separation of defects, and [small script l]_a is the active length scale. Our linear perturbation approximation, which allows us to treat the flows of individual defects separately, breaks down when the cores of defects overlap. Additionally, we have neglected the effect of topologically trivial structures such as bend walls, which also generate flows in the system,⁴⁸ on the motion of defects. When 〈r_d〉 > a, the effect of defect core overlap will be minimal, and when [small script l]_a > 〈r_d〉, structures such as bend walls will be quickly relaxed due to elastic forces.

We make the approximation that the director field is always given by the minimal energy configuration. It has been shown in experiments that this is not always the case for active nematics,^47,66 and induced torques, transverse velocities, and orientational noise may play an important role in the motion of active nematic defects.^9,11,24,58 This presents an opportunity for a future extension of the analytical model, since it has been shown that the defect velocity formalism used here can capture the transverse motion induced by these torques.¹⁶ However, to our knowledge, the flows induced by the relative twisting of defects in active nematics are not yet well understood. If approximate equations for this can be derived, then the orientation of topological defects in active nematics may be treated independently of their positions, while still taking into account the intrinsic coupling of the director field and defect orientation, and a true particle model could be formulated using the present theoretical framework. Accurately treating the orientation of defects independently would allow for many additional features to be included in the model such as defect nucleation or orientational noise. It would also be possible to include periodic boundary conditions, which cannot be treated in the present model due to instantaneous jumps in the orientation when defects cross the periodic boundary. In each case, it would be possible to maintain the orientational coupling that we have shown is important when describing active nematic defect motion.

An interesting extension of the methodology presented here would be to treat line defects in three-dimensional active nematics. Although it has been shown that the defect velocity formalism we employ here can be generalized to line defects in three dimensions,^16,67 there has been much less work in understanding defect dynamics in three dimensional active nematics both computationally,^68,69 and, especially, experimentally.⁷⁰ Recent theoretical work has described some aspects of defect loop motion,^71,72 and it will be interesting to compare and combine the results of the methodologies used here to obtain a clearer picture of defect motion in active nematics.

6 Summary

We have presented an analytical model for active nematic defect motion based on recently derived theoretical techniques and solutions to the Stokes equation for active nematic defects. Our work builds upon previous efforts to understand defect motion by applying a linear perturbation approximation to the Halperin–Mazenko formalism of defect dynamics. Using this, explicit velocities for an arbitrary number of defects may be predicted. The model incorporates both the flow-induced long range interactions between defects and also the crucial orientational coupling between the defects. We demonstrate that our model provides analytical predictions for the stable interaction radius adopted by bound states of pairs of +1/2 defects, the effective attraction that appears between pairs of −1/2 defects, and the way in which the critical unbinding length between a +1/2 and −1/2 defect pair scales with activity. The model can also capture the dependence of the scattering trajectory of a +1/2 defect from a −1/2 defect on activity. It can be applied to the braiding motion of multiple confined +1/2 defects as well as the coarsening behavior that occurs for a mixture of equal numbers of +1/2 and −1/2 defects under confinement. It will serve as a useful tool for future studies.

Data availability

All data supporting the findings of this study are presented in the paper and are available from the corresponding author upon request.

Conflicts of interest

There are no conflicts to declare.

Appendices

Defect velocity calculations

Here we present more detailed calculations leading to the expressions in eqn (5), (6), and (10). For completeness, we start by briefly reviewing the arguments leading to eqn (2), though more detailed derivations may be found in ref. 12 and 16 for nematic defects and ref. 52–54 for more general systems of defects.

The topological defect density in two dimensions may be written as

where m_i is the winding number of the ith defect. Using the fact that zeros of Q correspond to topological defects in the nematic phase, the Dirac delta function may be rewritten over Q(r) and ρ may be written as

ρ(r,t) = D(r,t)δ[Q(r,t)] (34)

D = ε_{k[small script l]}ε_μν∂_kQ_μα∂_{[small script l]}Q_να (35)

where D is related to the Jacobian transformation from real space to order parameter (Q) space. Taking a time derivative of D reveals

∂_tD = ∂_k(2ε_{k[small script l]}ε_μν∂_tQ_μα∂_{[small script l]}Q_να) ≡ ∇·J_D. (36)

Thus, D is a conserved field with current J_D. Taking a time derivative of both expressions for ρ, rearranging terms, and enforcing delta function relations yields eqn (2).

To get analytical approximations for the velocity, we must approximate Q near defect cores. We assume a linear core structure:¹³

with n̂₀ = (cos φ₀,sin φ₀), n̂₁ = n̂₀·ε, and ñ_i = n̂_i + [small variant phi, Greek, tilde](n̂_i·ε). [small variant phi, Greek, tilde] is the director angle due to all other defects:

Thus, we approximate the director distortion near the defect core of interest as a linear perturbation from all other defects.

The Coulomb velocity in eqn (5) may be derived by substituting ∂_tQ_μα = −δF/δQ_μα = (1/k)∇²Q_μα into eqn (2). Using eqn (37) to compute the spatial derivative and setting x = y = 0 gives eqn (5). Given a fluid flow u(r = 0), substituting ∂_tQ_μα = −(u·∇)Q_μα into eqn (2) gives the advection component of the velocity, eqn (6), after using eqn (37) to compute the gradient of Q. To get the shear component of the defect velocity, eqn (10), we first note that upon substituting ∂_tQ_μα = S_μα into eqn (2), most of the terms are zero when x = y = 0 since Q = 0 at the core of the defect. The only term that contributes to the velocity is ∂_tQ_μα = λE_μα, and substitution of this into eqn (2) yields eqn (10).

Radial flow solutions

Here we give the functional form of the distance dependent part of the flow and shear flow vectors for a single ±1/2 defect. The flow solutions are equivalent to those derived in ref. 43, though we have rewritten them to decompose them into radial and angular parts. The radial parts written here are infinite series which are tabulated to find the results in the main text.

The radial parts of the flow solutions are given by

The radial parts of the strain rate solutions are given by

To obtain the small x approximation of β₁ we compute the first two terms giving β₁(x) = π − 4x + ⋯ Comparing to a function of the form f(x) = B/(A + x), we may write (for x < 1) f(x) = B/A − (B/A²)x + ⋯ Thus, assuming β₁(x) scales as f(x) for small x, we identify B = π²/4 and A = π/4, or β₁(x) ≈ (π²/4)/(π/4 + x). The inset of Fig. 2(a) shows empirically that this is a reasonable approximation for small x.

Continuum simulation details

The continuum simulations shown in the main text solve eqn (3) and (11) for the Q tensor and flow velocity, u, which are rendered dimensionless via the same scaling described in the main text. For the free energy, we use the Landau-de Gennes bulk free energy with a one-constant elastic term:For all simulations we set A = −0.5, so at equilibrium the liquid crystal is in the nematic phase with . To numerically solve the dynamical equations, eqn (3) is discretized in time with an implicit backward Euler scheme and a Newton–Raphson sub-iteration to solve the resulting nonlinear equations. Both eqn (3) and (11) are discretized in space on either a square domain or a circular domain and we solve the weak form of partial differential equations using the MATLAB/C++ package FELICITY.⁷³ At each time step, given a configuration of Q we first solve eqn (11) for u. We then use this solution to compute the next time step for Q. For all simulations we use a time-step of Δt = 0.1.

For the simulation snapshots shown in Fig. 4(a), we use periodic boundary conditions on a square domain. Periodic boundary conditions necessitate that the total topological charge remains neutral. On the other hand, for the continuum simulation snapshots shown in Fig. 4(b) and data shown in Fig. 4(c) and 5(b), we use Dirichlet boundary conditions on a circular domain, with u = 0, S = S_N, and varying director depending on the total topological charge. To induce two +1/2 defects, the angle of the director field at the boundary is given by ϕ(R) = arctan(y/x) + φ₀, while for two −1/2 defects it is given by ϕ(R) = −arctan(y/x) + φ₀.

For the continuum simulation where we compare defect coarsening, shown in Fig. 10, we use Dirichlet boundary conditions with ϕ(R) = 0, requiring a total topological charge of zero. To find the number of defects as a function of time, we integrate the magnitude of D, given in eqn (35):^14,74

Acknowledgements

We gratefully acknowledge the support of the U.S. Department of Energy through the LANL/LDRD program for this work. This work was supported by the U.S. Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No. 89233218CNA000001).

Notes and references

1. P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics, Cambridge University Press, 1995.
2. T. W. B. Kibble, Aust. J. Phys., 1997, 50, 697.
3. L. M. Pismen, Vortices in Nonlinear Fields, Oxford University Press, 1999.
4. M. Kleman and J. Friedel, Rev. Mod. Phys., 2008, 80, 61–115.
5. 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.
6. Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun and K. Keren, Nat. Phys., 2021, 17, 251–259.
7. S.-Z. Lin, C. Reichhardt, C. D. Batista and A. Saxena, Phys. Rev. B:Condens. Matter Mater. Phys., 2013, 87, 214419.
8. C. Reichhardt and C. Reichhardt, Rep. Prog. Phys., 2016, 80, 026501.
9. X. Tang and J. V. Selinger, Soft Matter, 2017, 13, 5481.
10. A. Skaugen, L. Angheluta and J. Viñals, Phys. Rev. B, 2018, 97, 054113.
11. S. Shankar and M. C. Marchetti, Phys. Rev. X, 2019, 9, 041047.
12. L. Angheluta, Z. Chen, M. C. Marchetti and M. J. Bowick, New J. Phys., 2021, 23, 033009.
13. C. Long, X. Tang, R. L. Selinger and J. V. Selinger, Soft Matter, 2021, 17, 2265.
14. C. D. Schimming and J. Viñals, Soft Matter, 2022, 18, 2234–2244.
15. V. Skogvoll, L. Angheluta, A. Skaugen, M. Salvalaglio and J. Viñals, J. Mech. Phys. Solids, 2022, 166, 104932.
16. C. D. Schimming and J. Viñals, Proc. R. Soc. A, 2023, 479, 20230042.
17. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143.
18. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 3246.
19. S. J. DeCamp, G. S. Redner, A. Baskaran, M. F. Hagan and Z. Dogic, Nat. Mater., 2015, 14, 1110–1115.
20. L. Giomi, Phys. Rev. X, 2015, 5, 031003.
21. A. Doostmohammadi, T. N. Shendruk, K. Thijssen and J. M. Yeomans, Nat. Commun., 2016, 8, 15326.
22. L. M. Lemma, S. J. DeCamp, Z. You, L. Giomi and Z. Dogic, Soft Matter, 2019, 15, 3264–3272.
23. L. Giomi, M. J. Bowick, X. Ma and M. C. Marchetti, Phys. Rev. Lett., 2013, 110, 228101.
24. S. Shankar, S. Ramaswamy, M. C. Marchetti and M. J. Bowick, Phys. Rev. Lett., 2018, 121, 108002.
25. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031921.
26. I. Lazo, C. Peng, J. Xiang, S. V. Shiyanovskii and O. D. Lavrentovich, Nat. Commun., 2014, 5, 5033.
27. C. Conklin and J. Viñals, Soft Matter, 2017, 13, 725–739.
28. M. M. Genkin, A. Sokolov, O. D. Lavrentovich and I. S. Aranson, Phys. Rev. X, 2017, 7, 011029.
29. X. Li, J. C. Armas-Pérez, J. P. Hernández-Ortiz, C. G. Arges, X. Liu, J. A. Martínez-González, L. E. Ocola, C. Bishop, H. Xie, J. J. de Pablo and P. F. Nealey, ACS Nano, 2017, 11, 6492–6501.
30. 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.
31. K. Copenhagen, R. Alert, N. S. Wingreen and J. W. Shaevitz, Nat. Phys., 2021, 17, 211–215.
32. N. Figueroa-Morales, M. M. Genkin, A. Sokolov and I. S. Aranson, Commun. Phys., 2022, 5, 301.
33. R. Zhang, A. Mozaffari and J. J. de Pablo, Sci. Adv., 2022, 8, eabg9060.
34. M. Serra, L. Lemma, L. Giomi, Z. Dogic and L. Mahadevan, Nat. Phys., 2023, 19, 1355–1361.
35. S. Shankar, L. V. D. Scharrer, M. J. Bowick and M. C. Marchetti, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2400933121.
36. F. L. Memarian, D. Hammar, M. M. H. Sabbir, M. Elias, K. A. Mitchell and L. S. Hirst, Phys. Rev. Lett., 2024, 132, 228301.
37. K. A. Mitchell, M. M. H. Sabbir, K. Geumhan, S. A. Smith, B. Klein and D. A. Beller, Phys. Rev. E, 2024, 109, 014606.
38. D. Cortese, J. Eggers and T. B. Liverpool, Phys. Rev. E, 2018, 97, 022704.
39. Y.-H. Zhang, M. Deserno and Z.-C. Tu, Phys. Rev. E, 2020, 102, 012607.
40. L. M. Pismen, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2013, 88, 050502.
41. L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. C. Marchetti, Philos. Trans. R. Soc., A, 2014, 372, 20130365.
42. L. M. Pismen and F. Sagués, Eur. Phys. J. E:Soft Matter Biol. Phys., 2017, 40, 92.
43. J. Rønning, C. M. Marchetti, M. J. Bowick and L. Angheluta, Proc. R. Soc. A, 2022, 478, 20210879.
44. J. Rønning, M. C. Marchetti and L. Angheluta, R. Soc. Open Sci., 2023, 10, 221229.
45. 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.
46. P. W. Ellis, D. J. G. Pearce, Y.-W. Chang, G. Goldsztein, L. Giomi and A. Fernandez-Nieves, Nat. Phys., 2018, 14, 85–90.
47. D. J. G. Pearce, J. Nambisan, P. W. Ellis, A. Fernandez-Nieves and L. Giomi, Phys. Rev. Lett., 2021, 127, 197801.
48. L. C. Head, C. Doré, R. R. Keogh, L. Bonn, G. Negro, D. Marenduzzo, A. Doostmohommadi, K. Thijssen, T. L’opez-Le’on and T. N. Shendruk, Nat. Phys., 2024, 20, 492–500.
49. F. Vafa, Soft Matter, 2022, 18, 8087–8097.
50. P. G. de Gennes, The Physics of Liquid Crystals, Oxford University Press, 1975.
51. B. I. Halperin, Physics of Defects, North-Holland Pub. Co., 1981.
52. F. Liu and G. F. Mazenko, Phys. Rev. B:Condens. Matter Mater. Phys., 1992, 46, 5963.
53. G. F. Mazenko and R. A. Wickham, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1997, 57, 2539.
54. G. F. Mazenko, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, 1574.
55. A. N. Beris and B. J. Edwards, Thermodynamics of flowing systems, Oxford University Press, 1994.
56. F. M. Leslie, J. Mech. Appl. Math., 1966, 19, 357.
57. G. Tóth, C. Denniston and J. M. Yeomans, Phys. Rev. Lett., 2002, 88, 105504.
58. A. J. Vromans and L. Giomi, Soft Matter, 2016, 12, 6490.
59. M. M. Norton, A. Baskaran, A. Opathalage, B. Langeslay, S. Fraden, A. Baskaran and M. F. Hagan, Phys. Rev. E, 2018, 97, 012702.
60. C. D. Schimming, C. J. O. Reichhardt and C. Reichhardt, Phys. Rev. E, 2023, 108, L012602.
61. C. D. Schimming, C. J. O. Reichhardt and C. Reichhardt, Phys. Rev. Lett., 2024, 132, 018301.
62. A. Opathalage, M. M. Norton, M. P. N. Juniper, B. Langeslay, S. A. Aghvami, S. Fraden and Z. Dogic, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 4788–4797.
63. A. Sokolov, J. Katuri, J. J. de Pablo and A. Snezhko, Synthetic active liquid crystals powered by acoustic waves, 2024, https://arxiv.org/abs/2403.17268.
64. E. J. Hemingway, P. Mishra, M. C. Marchetti and S. M. Fielding, Soft Matter, 2016, 12, 7943–7952.
65. A. Repula, C. Gates, J. C. Cameron and I. I. Smalyukh, Commun. Mater., 2024, 5, 37.
66. D. J. G. Pearce and K. Kruse, Soft Matter, 2021, 17, 7408.
67. Y. Zushi, C. D. Schimming and K. A. Takeuchi, Phys. Rev. Res., 2024, 6, 023284.
68. S. Čopar, J. Aplinc, I. C. V. Kos, S. Z. Žumer and M. Ravnik, Phys. Rev. X, 2019, 9, 031051.
69. N. Kralj, M. Ravnik and I. C. V. Kos, Phys. Rev. Lett., 2023, 130, 128101.
70. G. Duclos, R. Adkins, D. Banerjee, M. S. E. Peterson, M. Varghese, I. Kolvin, A. Baskaran, R. A. Pelcovits, T. R. Powers, A. Baskaran, F. Toschi, M. F. Hagan, S. J. Streichan, V. Vitelli, D. A. Beller and Z. Dogic, Science, 2020, 367, 112–1124.
71. J. Binysh, Ž. Kos, S. C̆opar, M. Ravnik and G. P. Alexander, Phys. Rev. Lett., 2020, 124, 088001.
72. A. J. H. Houston and G. P. Alexander, Phys. Rev. E, 2022, 105, L062601.
73. S. W. Walker, SIAM J. Sci. Comput., 2018, 40, C234–C257.
74. M. L. Blow, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2014, 113, 248303.

---

Footnote

† Electronic supplementary information (ESI) available: In all Supplemental Movies the points are the positions of defects, the arrows are the defect orientations p̂, and the black circles indicate the domain boundaries. Below we list the specific parameters for each movie. Supplemental Movie 1: three +1/2 defects, R = 3, α = 2, [small script l]_h = 1. Supplemental Movie 2: three +1/2 defects, R = 5, α = 2, [small script l]_h = 1. Supplemental Movie 3: four +1/2 defects, R = 3, α = 2, [small script l]_h = 1. Supplemental Movie 4: four +1/2 defects, R = 5, α = 2, [small script l]_h = 1. Supplemental Movie 5: five +1/2 defects, R = 3, α = 3, [small script l]_h = 1. Supplemental Movie 6: five +1/2 defects, R = 5, α = 3, [small script l]_h = 1. Supplemental Movie 7: six +1/2 defects, R = 3, α = 3, [small script l]_h = 1. Supplemental Movie 8: six +1/2 defects, R = 5, α = 3, [small script l]_h = 1. See DOI: https://doi.org/10.1039/d4sm00956h

---