Dynamical theory of topological defects II: universal aspects of defect motion

We study the dynamics of topological defects in continuum theories governed by a free energy minimization principle, building on our recently developed framework (Romano et al 2023 J. Stat. Mech. 083211). We show how the equation of motion of point defects, domain walls, disclination lines and any other singularity can be understood with one unifying mathematical framework. For disclination lines, this also allows us to study the interplay between the internal line tension and the interaction with other lines. This interplay is non-trivial, allowing defect loops to expand, instead of contracting, due to external interaction. We also use this framework to obtain an analytical description of two long-lasting problems in point defect motion, namely the scale dependence of the defect mobility and the role of elastic anisotropy in the motion of defects in liquid crystals. For the former, we show that the effective defect mobility is strongly problem-dependent, but it can be computed with high accuracy for a pair of annihilating defects. For the latter, we show that at the first order in perturbation theory, anisotropy causes a non-radial force, making the trajectory of annihilating defects deviate from a straight line. At higher orders, it also induces a correction in the mobility, which becomes non-isotropic for the +1/2 defect. We argue that, due to its generality, our method can help to shed light on the motion of singularities in many different systems, including driven and active non-equilibrium theories.

In strongly ordered systems, defects are usually described over macroscopic scales as quasi-particles in interaction with the surrounding order parameter phase field [16].The derivation of this reduced description from the dynamics of the full order parameter has been the subject of many studies [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35], while the ubiquity and dynamical characteristics of defects in active matter have led to a recent revival of interest in this problem [36][37][38][39][40][41][42][43].One difficulty in carrying out such coarse-graining procedure is that defects are intrinsically microscopic structures, such that their description a priori requires the knowledge of the order parameter dynamics over microscopic scales.Hence, many studies have relied on matched asymptotics by solving the field theory in the vicinity and far away from the defect, while the continuity of the full solution is imposed in an intermediate matching region (for a pedagogical introduction, see Ref. [16]).On the other hand, when the defect core size is made truly microscopic, there is no guarantee that it is faithfully described by the phenomenological theories expressed in terms of smoothly varying fields.
The effective equations of motion governing the dynamics of defects generally take a similar functional form [16], suggesting a certain degree of universality insofar as the details of the microscopic core structure only set the value of certain coefficients in the reduced description.Existing results are, however, mostly restricted to idealized cases (except, e.g., Ref. [44]), and thus omit important features present in real systems.For example, while most approaches consider the limit of slow defects, we have shown that significant memory effects emerge due to the dependencies of the order parameter landscape on the past position and velocity of the defect [45].Another feature of liquid crystals generally ignored is the effect of elastic anisotropy, which introduces higher order nonlinearities to the order parameter field theory.Anisotropy in the elastic response of the medium, on the other hand, is responsible for qualitative changes in the dynamics of defects [46][47][48][49] that cannot be accounted for in the single Frank constant approximation.
Here, we propose a new approach allowing to derive the dynamical equations of motion for defects from any dissipative field theory that satisfies a minimization principle.Importantly, this approach is formally valid at all orders in the defect velocity and for any free energy functional describing the dynamics.In Sec.2.1, we demonstrate that, under a set of rather weak assumptions, the defect equation of motion takes a universal form as the details of the core structure set the value of a unique length scale in the expression of the mobility, while the effective force moving the defect is fully determined by the large scale physics.
To illustrate the power of the approach, we apply it to multiple scenarios involving either point or line defects.We show in particular that the main features of the physics of topological defects can be captured by means of a low mobility expansion, which leads to substantial simplifications in the equations of motion.Whereas most of existing results correspond to the leading order contribution to this expansion, we show in Sec.2.2 how in simple cases improved approximations can be obtained by considering higher order corrections.Applying the method to a theory describing nematic liquid crystals in two dimensions, we moreover quantify in Sec.2.3 how elastic anisotropy spontaneously rotates pairs of annihilating defects and affects their mobilities.Section 3 is finally devoted to defect lines.A derivation of the Allen-Cahn equation [50] for domain walls is firstly given in Sec.3.1, while the dynamics of disclination lines and loops emerging in three-dimensional phases with broken U (1) symmetry are discussed in Sec.3.2.

Derivation of the defect dynamics from free energy variations
Throughout this section, we study a two-dimensional system described by an order parameter ϕ(x, t) whose dynamics minimizes a free energy F = d 2 x F (ϕ, ∇ϕ) as described by Equation ( 1) corresponds to a deterministic version of model A in the Halperin-Hohenberg classification [51], and thus serves as a general form to describe any relaxational dynamics without conservation law.For example, in dynamics with broken polar or nematic orientational order, ϕ corresponds to a vectorial or rank-2 tensor field: where ρ(x, t) and θ(x, t) respectively set the magnitude and orientation of order.On the other hand, phases with broken Z 2 symmetry will be described by a scalar order parameter.
In what follows, we work in a parameter regime for which the system is strongly ordered far away from the defects.In practice, this implies that the free energy F reduces to a known functional F bulk that depends only on the slow modes of the dynamics.The free energy F bulk then describes the dynamics of ϕ everywhere except in the vicinity of defects, where it is captured by a a priori unknown free energy F core .For instance, as will be detailed further in Sec.2.2, for systems with orientational order the strongly ordered limit corresponds to the case where the norm ρ of ϕ is a fast mode, while F bulk depends only on the orientation θ and its derivatives.
In the remaining of this section, we focus on the case of point defects such as those occurring in two-dimensional phases with orientational order, while the case with discrete symmetry will be addressed in Sec. 3. We consider a configuration with an arbitrary number of defects, and derive the equation of motion of a specific defect whose position and velocity are denoted respectively as q(t) and v(t).For convenience, we will use a generalized scalar product simply defined as the sum of squared components: Below, we present a detailed derivation of the equation of motion for defects which applies to an arbitrary free energy F satisfying the following assumptions: (i) Translational invariance: the free energy density F core varies in space only through the field ϕ, so that it is translationally invariant.
(ii) Microscopic core size: the dynamics of the core is associated with a length scale a → 0, which plays the role of the core size.As |ϕ| varies from zero at the centre of the core to one over a distance a, gradients of ϕ in the core are of order a −1 .
(iii) The rigid core assumption: at the leading order in a, the shape of the core is independent of its position or velocity.Up to O(a) terms, ϕ(x, t) can thus be expressed in the reference frame of the core in terms of a fixed function φ(y) with y = R −1 (t)[x − q(t)]/a and where the rotation matrix R(t) parametrizes the direction of the core.The relationship between ϕ and φ depends on the nature of the order.For vectorial and nematic (rank-2 tensor) orders, we respectively have (iv) Existence of a matching regime: lastly, we assume the existence of an intermediate region at distance ∼ r 0 from the core where F bulk and F core coincide and both describe the dynamics.This region shall be well-separated from the other scales of the problem, such that a ≪ r 0 ≪ L where L stands for a macroscopic scale, e.g. the system size or the typical inter-defect distance (see a sketch in Fig. 1).
Assumption (i) is natural so long as the system is not externally driven by a spatially dependent field ‡, while (ii) ensures a proper scale separation between the defect core size and the macroscopic dynamics.(iv) reflects the continuity condition between the core and bulk physics, and plays a central role in the matched asymptotic methods [16].Assumption (iii), in fact, follows from (i) and (ii).Indeed, for vanishing a the field ϕ at the core can be generally expanded in powers of a/L, with L a problem-dependent macroscopic scale.The leading order of this expansion is by construction independent of the relative positions and velocities of other defects, as they contribute to terms at least O(a/L).Moreover, from (i) the order parameter at the core is, up to rotations and translations, uniquely determined by F core , which eventually leads to (iii).‡ (i) should still be verified in the presence of weak or smoothly varying fields.Formally, we require that the external field varies on length scales ≫ a and timescales ≫ a 2 /D, where D is the effective diffusivity associated with the dynamics of ϕ.
Figure 1: Schematics of the three scales a, r 0 , and L involved in the derivation.The microscopic scale a sets the defect core size (red region).The macroscopic scale L is given by the system size, or the typical distance between defects.In the intermediate matching region parametrized by r 0 (yellow), it is assumed that both the core and bulk theories hold.
2.1.1.The force applied on a defect core To find the equation of motion for the defect core, we first express the variational equation ( 1) in an integral form as where δϕ is a fixed boundary condition perturbation.By construction, Eq. ( 2) is valid for any infinitesimal δϕ.Here, we consider a specific type of perturbation: where δq is an infinitesimal vector and f (x, t) is a smooth interpolating envelope function equal to one at x = q(t), which decays quickly to zero for |x − q(t)| > r 0 , such that it is zero at the boundaries of the system and at the positions of other defects.The introduction of the envelope function f is done to isolate the core for which we wish to derive the equation of motion.It is not strictly necessary, but simplifies the derivation by allowing us to discard the effect of the variation δϕ at the system boundaries and at the other defect cores.Discarding this prefactor, would in fact lead to additional O(a/L) contributions to the final equation, which are subdominant in the limit of well separated scales.Defining a closed curve Υ within the matching region, we split the l.h.s. and r.h.s. of Eq. (2) into contributions from inside and outside Υ, corresponding respectively to the core and the bulk.Namely, We now re-express the bulk free energy variation as where the second term on the r.h.s. is a surface contribution retained after integrating by parts, and where we assumed that the free energy density F bulk (ϕ, ∇ϕ) does not depend on higher derivatives of ϕ §.As δϕ is zero outside the core region, this surface term only includes a contribution from the boundary Υ, while the preceding minus sign implies that dS points to the outside of the core.
To compute δF core , we note that in the core δϕ = (δq • ∇)ϕ as by construction f (|x−q| ≤ r 0 ) = 1.Therefore, δϕ in the core corresponds to an infinitesimal translation.Since, from (i), F core is translationally invariant, it follows that where the dots stand for dependencies of F core on higher order derivatives of ϕ, if any.Hence, δF core is an exact differential, and therefore, we obtain where the second equality results from the fact that Υ belongs to the matching region.Putting together Eqs.(4,5,6), we find that where we have defined the canonical stress tensor of the bulk theory [52,53]: Hence, the r.h.s. of Eq. ( 7) corresponds to the net momentum flux through the matching region, and is solely determined by the large-scale bulk physics.

The defect friction tensor
As we show now, the l.h.s. of (7) weakly depends on the specific form of the core free energy.We first note that Eq. ( 7) holds for an arbitrary choice of curve Υ in the matching region.Choosing without loss of generality Υ to be a circle of radius r 0 around the singularity, we thus expect that the final result will be independent of r 0 .Using the rigid core assumption (iii), it is clear that in the core ∂ t ϕ = −(v • ∇)ϕ∥, such that the l.h.s. of Eq. ( 7) is given by where we have used the change of variable y = R −1 (t)(x − q(t))/a and D r 0 /a is the disk of radius r 0 /a centred at 0. Equation (8) defines the effective friction tensor of the defect dynamics: § The generalization to cases where F bulk depends on higher derivatives of ϕ is straightforward and does not affect the final equation (7).∥ The time derivative of ϕ also includes a contribution from the rotation matrix R(t).However, this contribution is subdominant in the limit r 0 /a ≫ 1.
We note that, since φ is fully determined by the core structure, it is independent of the macroscopic defect variables such at its position, velocity and orientation.The tensor ζ is thus a fixed function which specifies the structure of the defect friction, and can be evaluated at leading order in a from the static single defect solution.Namely, differentiating (9) w.r.t.r 0 and parametrizing y with the polar coordinates (r 0 , φ), we find that where the 'ssd' subscript indicates that the integrand is calculated from the static single defect solution of the bulk theory, since the right hand side of this equation is evaluated in the matching region.As already noted in a number of works [4,27], the functional shape of the friction tensor is only determined by the bulk theory, while the core theory only enters through an integration constant when (10) is integrated on both sides.This integration constant plays the role of a phenomenological parameter that captures the microscopic features of the core.The time dependency of ζ results from the fact that an anisotropic core structure may lead the defect to experience different friction strengths in different directions.Although ζ is determined by the a priori unknown core free energy F core , we show in Appendix A that in most relevant situations ζ is often isotropic and thus independent of the defect orientation (see, however, Sec.2.3.3 for a counter example).In the following, we therefore keep the time dependency of ζ implicit.

The general equation of motion for defects
Gathering the results accumulated so far, and noting that Eq. ( 7) must be valid for all δq, we find that the defect equation of motion takes the compact form where C r 0 stands for the circle of radius r 0 centred at the defect core.Equation (11) highlights that the defect equation of motion, up to a constant factor in the mobility, is universal as it does not depend on the core physics so long as the latter satisfies translational and rotational invariance.This equation moreover bears a transparent physical interpretation, since it simply states that the momentum flux through the boundary of the core is, up to frictional effects, entirely dissipated into the motion of the defect, which is a natural consequence of (iii).The r.h.s. of Eq. ( 11) corresponds to the Ericksen force defined by Eshelby [19], which moreover takes the same formal expression as the Peach-Koehler force acting on dislocations [19,54].Lastly, it is important to note that both sides of Eq. ( 11) depend on the matching variable r 0 , while the actual equation of motion of the defect must be independent of it, since r 0 is arbitrary.In fact, we show below that eliminating r 0 self-consistently allows in some cases to fix the functional form of the mobility.We now illustrate how (11) can be used to explicitly derive the dynamics of defects.We start by showing how standard results are recovered for systems described by the archetypal Ginzburg-Landau free energy (12).We then address nonlinear problems such as when defects dynamics evolve in a medium featuring elastic anisotropy.
2.2.Defects in the classical Ginzburg-Landau framework 2.2.1.Elimination of the matching scale The simplest choice for F is certainly the Ginzburg-Landau free energy Note that we work in time units such that the coefficient in front of the elastic term in (12) has been set to one.The scale of the defect core is then given by a ≃ χ −1/2 .Over scales much larger than a, the dynamics of the order parameter is then fully captured by that of its orientation θ, which is associated with the bulk free energy The dynamics of θ is thus ruled by the diffusion equation, while the stress tensor associated with ( 13) is given by In particular, we showed previously [45] that the orientation field gradient induced by a defect of charge s following a trajectory q(t) with velocity v(t) is given by where ϵ denotes the two-dimensional antisymmetric Levi-Civita tensor Using the linearity property of the diffusion equation, the total orientation field landscape induced multiple defects corresponds to the linear superposition of single defect solutions (14).
Since the matching region is assumed to be well separated from the other scales of the theory, we now explicitly evaluate both sides of Eq. ( 11) in the limit of a → 0 and r 0 → 0, keeping r 0 ≫ a.To obtain the r.h.s., we note that the presence of a defect leads to a divergence of |∇θ| ∼ a −1 at the core, such that for r 0 small the contour integral is dominated by the discontinuous part of ∇θ.Although the general expression ( 14) is nonlocal in time, the discontinuous part ∇θ d of ∇θ in the vicinity of a defect is formally determined by its instantaneous position and velocity [45], namely where we have defined r ≡ x − q(t), r ≡ |r| and r ≡ r/r.The additional length scale λ(t) leads to a continuous contribution to (15), but was included for dimensional consistency.This quantity formally depends on the whole knowledge of the past history of the defect trajectory, and it is generally not possible to evaluate it directly from (14).
For now, we thus retain it as a phenomenological constant, and we will show in the following sections how it can be determined or approximated.Denoting ∇θ = ∇θ d + ∇θ c , with ∇θ c accounting for the continuous part of the gradient, we calculate in Appendix B the r.h.s. of Eq. ( 11), which leads for r 0 → 0 to It is clear that only the terms on the l.h.s. of Eq. ( 16) depend on the matching variable r 0 .Hence, we conclude that the friction ζ ij = ζδ ij while the quantity r 0 exp[−ζ/(πs 2 )] must remain independent of r 0 .Furthermore, from its definition (9) the defect friction is independent of any macroscopic scale, which yields where ā ∝ a is a phenomenological constant that plays the role of the core effective radius, and depends on the details of the core physics.It is straightforward to show that this expression of ζ satisfies (10), while the same result could have been obtained directly by calculating the integral in (10) using the solution (15) with v = 0.The equation of motion for the defect therefore reads ln e Equation ( 17) is formally similar to Eq. ( 11).However, whereas the latter depends on the arbitrary matching variable r 0 through ζ and the integration contour on the r.h.s., Eq. ( 17) determines the defect motion independently of the choice of matching region.The only quantity that depends on the core physics in (17) is the parameter ā.Its value is set by the core free energy F core , which sets the profile of the full order parameter ϕ at the core.For the Ginzburg-Landau free energy (12), it has been shown that ā√ χ ≃ 1.126 [16].The coefficient λ(t) in the expression of the effective friction, on the other hand, is a macroscopic scale fixed by the history of the defect trajectory [45].The r.h.s. of Eq. (17) shows that defects are essentially moved by gradients of the orientation field [19].This gradient is in general generated by other defects, or by specific anchoring conditions at the system boundaries.Its expression and that of λ(t) can in principle be obtained by solving for the dynamics of the orientation field θ with the appropriate boundary conditions at the defect cores.Below, we show that the terms of Eq. ( 17) can be explicitly determined in particular configurations.

An isolated defect moving at constant velocity
We first study the simple case of a defect moving uniformly with velocity v = v.An experimental realization of this case would for example consist of a defect subject to an imposed spin wave θ ext (x) ≡ k • x, with k the corresponding wavevector, induced by an external field.Under these conditions, the total angular field is given by the sum of θ ext (x) and the uniformly moving defect solution to the diffusion equation, namely [31,35,45] where K 0 and K 1 are modified Bessel functions of the second kind while, as previously, s denotes the charge of the defect and r = x − q(t).Expanding Eq. ( 18) for r → 0 and keeping only nonvanishing terms, we find that ∇θ c = k while the discontinuous part ∇θ d is given by ( 15) with λ = 4v −1 e −γ E ; and where γ E ≈ 0.578 denotes the Euler-Mascheroni constant.
Replacing these expressions in Eq. ( 17), we thus recover the classical result [16] ln 4 e whose solution determines the velocity of the defect in response to an imposed spin wave.Consistently with the general result (17) and as noted by a number of authors [4, 18, 25, 27, 31, 35-37, 47, 55], the defect is subject to nonlinear friction.As we show below, this feature is rather generic, while the exact expression of the friction depends on the problem of interest.

Self-consistent solution for a pair of slowly moving defects
We now consider a pair of defects with opposite charges s ± = ±s.To proceed further, we note from ( 17) that the defect mobility µ ≡ ζ −1 ∼ ln −1 (λ/a).Hence, so long as the macroscopic scale λ remains much larger than the core radius, defects move relatively slowly with a speed In what follows, we thus treat µ as a small parameter, which allows us to derive an approximate expression for λ.
At first order in µ, the angular field created by a moving defect takes the universal form given by (15), regardless of its trajectory [45].The only trajectory-dependent quantity is then the parameter λ(t).Denoting, respectively, q ± (t) and v ± (t) as the positions and velocities of the two defects, it then follows from (17) that ln e where θ ± denote the orientation field landscapes generated by the positively and negatively charged defects, while the second equality was obtained by replacing θ ± with (15).λ ± in Eq. ( 20) denote the scales associated with the ± defects, while we have also defined q ≡ 1 2 (q + − q − ).It is a straightforward exercise to show that (20) implies that v + = −v − = v q, such that the speed of both defects follows As q(t) is the only macroscopic scale of the problem, we propose the ansatz: λ ± (t) ∝ q(t).
Hence, it follows that μeff ∼ vµ 2 eff ∼ µ 3 eff such that, at first order in the mobility, the angular field generated by a defect moving according to (21) can be calculated from (14) treating µ eff as a constant parameter.The details of this calculation are presented in Ref. [45], while the resulting expression for ∇θ is again formally given by (15) with Figure 2: Comparison of trajectories between simulation (data points) and the selfconsistent approximation (23) (solid blue line, SC) initial defect separation q(0) = 32 and q(0) = 64 (ā ≈ 0.36).Panel (a): Defect velocity as function of its position over the full simulation range.Panel (b): Zoom of (a) closer to the annihilation point.The yellow line is obtained by replacing q with a fixed scale (FS) q = 32 in the nested logarithm of ( 23).This approximation differs from the numerical data close to annihilation, but improves at larger distances.The dashed grey curves indicate the inverse distance scaling predicted by the Coulomb interaction.
. Combining this result with the expression of µ eff given in (21) yields the following self-consistent condition whose solution can be expressed in terms of the velocity as µ −1 eff = ln 4e . As there is no other macroscopic scale but the defects' degrees of freedom, we recover an expression for the friction similar to that of (19).The condition (22), however, only holds perturbatively in µ eff (∝ v).Evaluating the solution of Eq. ( 22) up to terms O (ln[ln(q/ā)]/ ln(q/ā)), we thus obtain Equation ( 23) fully determines the dynamics of the defect pair.
Extensive discussions on scale- [4,37,47] or velocity- [18,25,27,31,[35][36][37]47, 55] dependent defect mobility can be found in the literature.The derivation outlined above shows that, in fact, the scale λ entering the expression of the mobility is primarily determined by the relevant scales of the background orientation field.In the general case, moreover, λ may depend on the past configurations of the system, leading to nontrivial memory effects [45].
We show in Figure 2 comparisons between the relation v(q) obtained from numerical simulations of a s = ±1 defect pair annihilation in a vectorial field described by the Ginzburg-Landau model (Eq.( 12)), and theoretical predictions with different approximation schemes for the defect mobilities (details on numerical simulations can be found in Appendix D).Since the simulations are initialized with the static defect profiles, the order parameter field first relaxes to the dynamical solution of (1).This effect induces a transient behaviour of v(q) at large q that depends on the initial defect separation.Past this transient, v(q) exhibits a behaviour independent of the initial defect separation, which we find in clear departure from the the ∝ q −1 scaling predicted by the constant mobility approximation (dashed lines in Fig. 2).Conversely, a parameter-free comparison with Eq. ( 23) reveals excellent agreement with the measured relation v(q) (blue curves in Fig. 2).

Many defects and the collective mobility
In the low mobility approximation, Eq. ( 20) can be generalized to an arbitrary number of defects as where q α , v α and s α denote the position, velocity and charge of the αth defect, respectively, while q αβ ≡ 1 2 (q α − q β ).Rearranging the terms, we recast (24) into a more compact form: where U C ≡ − α<β s α s β ln |q αβ | is the two-dimensional Coulomb potential, while the mobility matrix M is defined through γ M αγ Z γβ = δ αβ I, where the friction matrix Z is defined as and where I denotes the two-dimensional identity matrix.Equation ( 25) reveals that the many-body defect dynamics is coupled via the collective (non-diagonal) mobility M. When λ α ̸ = λ β , M is moreover not symmetric, such that it introduces effective non-reciprocal couplings between the defect velocities.Although the Coulomb force is conservative, the centre of mass of the system is thus generally not immobile.A similar effect was reported in the context of active nematics [40].Here, we observe that this effect also arises in the absence of active drive, so that it is generic to the relaxation dynamics of systems with a large number of defects.To rationalize this result, we note that both the interaction and the mobility in Eq. ( 25) are set by the orientation field landscape, which is itself driven out of equilibrium by the motion of the defects (see Eq. ( 14)).During relaxation to equilibrium, the effective interactions between defects are therefore mediated by a nonequilibrium medium, which is a well-identified mechanism for the generation of non-reciprocal couplings [56][57][58].
Unlike the two-defect case, it is generally not possible to determine analytically the parameters λ α (t) by calculating the exact solution for the angular field θ(r, t).We however note that, since the λ α are given by the macroscopic scales of the system (for example the mean inter-defect separation), the mobility in ( 25) is dominated by its diagonal coefficients which are of order ln(λ/ā), while the off-diagonal coefficients are O(1).At the leading order in ln(λ/ā), one can thus replace λ by any macroscopic scale of the problem.Indeed, replacing λ → λ + λ ′ does not change the result to the leading order, namely For two defects, Figure 2(b) shows that the approximation λ ∝ q matches well with the self-consistent solution (23) for large enough defect separation (compare the blue and yellow curves).As we show in the following section, this approximation is moreover of great use when dealing with more complex problems.

Defects in nonlinear systems
2.3.1.The large ln(a) expansion Although the terms of Eq. ( 11) can be explicitly calculated in simple scenarios, calculations quickly become intractable if the bulk free energy contains higher order nonlinearities.Here, we thus show that Eq. ( 11) admits a powerful expansion in the inverse of ln(a) related to the low mobility expansion discussed above.We first note that, as the static single defect solution for the orientation field θ ssd (r, φ)-where (r, φ) stands for the polar coordinates in the defect frame-is scale-free by construction, we have θ ssd (r, φ) = θ ssd (φ).Hence, Eq. ( 10) becomes where φ is the unit vector tangent to the unit circle.Equation ( 26) therefore implies that the friction matrix takes the form ζij where both ν and C are dimensionless and independent of a.We are now in a good position to expand Eq. ( 11) using ln(a) as a large parameter.Such expansion-although it does not involve a ratio of two scales-shall be seen as a formal way to take the limit of slowly moving defects.We note that at the leading order in ln(a), Eq. ( 27) simplifies to ζij = −ν ij ln(a) + O(1), such that the friction tensor does no more involve the matching variable r 0 nor the constant C.Moreover, replacing ζ ij ∼ ln(a) in Eq. ( 11), one observes that the defect velocity is at least of order ln −1 (a).Hence, at the leading order the stress tensor on the r.h.s. of Eq. ( 11) can be written as T bulk = T sb + O(ln −1 (a)), where the subscript 'sb' indicates that T is calculated from the static solution of the bulk theory.We also note that since the static solution satisfies δF bulk /δϕ sb = 0 (from (1)), the corresponding stress tensor is divergence-free: ∇ • T sb = 0. Consequently, the line integral on the r.h.s. of Eq. ( 11) is independent of the chosen path, and we get where ν ij = R ik (t)R jl (t)ν kl includes the possible effects of defect anisotropy (see Eq. ( 9)) and the integral on the r.h.s.can be evaluated on any loop enclosing the defect core.
The scale λ ≫ a on the l.h.s. was included for dimensional consistency, and gives a subdominant contribution in ln(a).At the leading order, λ can thus formally be replaced by any relevant macroscopic scale of the theory.We finally note that the ln(a) expansion provides a formulation that is independent of the matching variable r 0 .Contrary to (11), Eq. ( 28) can thus be directly used without the need to eliminate r 0 by enforcing the matching condition.

Elastic anisotropy and misalignment-induced forces
To illustrate how ( 28) can be employed to address more complex scenarios, we now consider the dynamics of a pair of topological defects in a nematic liquid crystal featuring elastic anisotropy.The order parameter for this case is therefore the nematic tensor such that topological defects carry a half-integer charge: s = ± 1 2 .The bulk dynamics is described by the Frank-Oseen free energy [59], which takes the form ¶ where the parameter α ∈ [−1; 1] is nonzero whenever splay and bend deformations incur different elastic costs.Although elastic anisotropy commonly occurs in liquid crystals, the theoretical understanding of defect dynamics in this context remains limited [46,48].This feature results from the fact that, for α ̸ = 0, the relaxational dynamics of θ follows a nonlinear equation, such that deriving the counterpart of Eq. ( 14) is generally difficult; even perturbatively in α.Here, we show how some progress can be made from Eq. ( 28) for slow defects and in the presence of weak anisotropy.
To calculate the defect mobility, we use (26) and consider the static single defect solution associated with the free energy (29).This solution can be expressed in an implicit form as [60] where ψ(φ) ≡ θ ssd (φ) − φ and p is a constant fixed by the circuitation condition 2π 0 dφ ∂ φ θ ssd = 2πs.Equation ( 30) can in principle be solved iteratively at any ¶ Equation ( 29) corresponds to the free energy density 2 expressed in terms of the director field n, where 1  2 order in α.At O(1), the solution is simply that of the isotropic theory, namely, , where θ c is an integration constant.At linear order in α, the solution becomes Inserting this expression in Eq. ( 26) with s = ± 1 2 and integrating over φ, we find that ν ij = νij = π 4 δ ij .Therefore, at linear order in α oppositely charged defects have equal isotropic mobilities, while corrections are expected at higher order in perturbation [46,47] (see also Sec. 2.3.3).
The evaluation of the force term in ( 28) is less straightforward and calculation details are presented in Appendix C. In summary, we express the static orientation field created by the two defects as θ sb (x) = θ 0 (x) + αθ 1 (x) + O(α 2 ), where θ 0 (x) solves the linear problem (α = 0) and θ 1 (x) is calculated perturbatively.However, the formal expression of θ 1 still involves intricate integrals.We thus use the fact that the integration contour on the r.h.s. of Eq. ( 28) is arbitrary, and integrate over the bisector of the segment formed by the two defects.This choice conveniently allows to express all integrals independently of any model parameters, such that they reduce to numerical constants.After numerically computing these coefficients, we finally obtain for the vector q = 1 2 (q where θ q denotes the angle between q and the direction of the nematic order parameter at infinity parametrized by θ c .Comparing Eqs. ( 32) and ( 21), we see that elastic anisotropy generates, at the linear order in α, a tangential force that essentially aligns q along the background nematic order orientation for α < 0, or orthogonal to it for α > 0. To give a physical interpretation of this force, we show in Fig. 3 that varying θ q amounts to changing the relative orientation of the two defects.The tangential force in Figure 4: Trajectories of the vector q obtained from numerical integration of Eq. ( 32) (full lines) and direct simulations of the annihilation of isolated ± 1 2 defect pairs in the full theory (33) (symbols).Panel (a) shows trajectories at fixed θ q = 0.7 for various values of α, while panel (b) corresponds to α = 0.2 and different values of θ q (rad).In (a,b) the dashed lines correspond to the expected trajectories at α = 0.
Eq. ( 32) then reflects the energetic cost of misaligned defects due to elastic anisotropy, and favors configurations with θ q = π 2 or 0. To confirm the predictions of Eq. (32), we performed numerical simulations of the dynamics described by for the full nematic order parameter Q, where χ ∼ a −1 sets the size of the defect core.F an is the simplest free energy that admits (29) as a bulk theory, which can be verified by setting ρ = 1 in Eq. (33). Figure .4 shows trajectories of the vector q obtained from initially prepared ± 1 2 defect pairs at various values of α and initial orientation θ q (details on numerical simulations are given in Appendix D).Note that the mobility in Eq. ( 32) only affects the speed of the defects, and not their trajectories.Hence, studying the trajectories we can compare the results obtained from the integration of Eq. ( 32) and simulations of the full model (33) without any fitting parameter.In qualitative agreement with the predictions of Eq. ( 32), we find that for α ̸ = 0 they annihilate following curved trajectories, as a result of the tangential component of the force between misaligned defects.We observe that the magnitude of the curvature increases with α.Quantitatively, the curves shown in Fig. 4 reveal a good agreement between theory and simulations, whereas deviations appear at larger α and close to the annihilation point, which we interpret as the breakdown of our perturbative approach and the slow defect approximation, respectively.

Anisotropic mobility of + 1
2 defects We showed previously that at the linear order in α the defect mobility is not affected by elastic anisotropy.Data obtained at stronger anisotropy from numerical simulations [46] or in experiments [47], however, suggest that higher order corrections will affect the mobilities of the defects.In this section, we therefore calculate the leading order contribution of elastic anisotropy to the defect mobility, which arise at O(α 2 ).At this order, the orientation field obtained from (30) is given by We use Eq. ( 26) to calculate the defect friction.Splitting it into isotropic and anisotropic parts, we obtain For the isotropic part, we obtain after some algebra In addition, the anisotropic part of the friction is found to be Hence, only s = + 1 2 defects have an anisotropic friction tensor.As we argue in Appendix A, we expect this result to hold at all orders in α.
We now note that θ c is defined in the previous section as the background orientation of the nematic order.Conversely, by rotational invariance of the problem θ c may define the direction of the comet-shaped + 1 2 defect.Namely, it is straightforward to show from (34) and for s = 1  2 that rotating θ c by an angle ϑ amounts to rotate the defect solution by an angle 2ϑ.Hence, we define the + 1  2 defect polarization as p ≡ (cos 2θ c , sin 2θ c ) (see the blue arrows in Fig. 3), such that we finally obtain for the two types of nematic defects Comparing the coefficients of Eq. ( 36), we thus find that the + 1 2 defect experiences a larger friction in the direction longitudinal to its comet-like shape, whereas elastic anisotropy reduces the strength of the friction in the transverse direction.Conversely, elastic anisotropy renormalizes the friction of the − 1 2 -charged defect upwards in all directions.

Domain walls and disclination lines
We now illustrate how the variational approach outlined in the previous section can be extended to describe defects with different geometries.We start by discussing the simple case of domain walls in the presence of Z 2 broken symmetry, and then extend our derivation to disclination lines.

Domain walls: the Allen-Cahn equation
Domain walls form in systems described by a scalar order parameter ϕ(x, t) that accounts for spontaneously broken Z 2 symmetry.Considering units for which ϕ equilibrates to the values φ = ±1, they correspond to thin interfaces separating homogeneous regions where ϕ takes opposite signs (Fig. 5).A minimal continuous description for this class of system is given by the following free energy where V (a, ϕ) is an arbitrary bulk potential that depends on a microscopic scale a setting the thickness of the domain walls.In particular, a common example is the Landau potential The derivation below relies on the same assumptions (i-iv) that were extensively used in Sec.2.1.Importantly, we will work in the limit a → 0 ensuring a separation of scales between the domain wall thickness and any macroscopic length.We consider a domain wall described by the curve γ(l, t) on the plane, where l sets the curve parametrization and t is the time.In what follows we will denote as γ ′ and γ derivatives of γ with respect to l and t, respectively.The derivation of the equation of motion for an infinitesimal segment [γ(l, t); γ(l + dl, t)] roughly follows that outlined in Sec.2.1.However, as the object we now describe is one-dimensional, we need to choose an appropriate volume V to delimit the matching region.Given a point γ(l, t), we define V as a tetragon enclosing the domain wall (as sketched in Fig. 5).The edges S 1 and S 2 cross the interface orthogonally respectively at γ(l, t) and γ(l + dl, t), and extend on both sides up to a distance r 0 ≫ a.In addition, the other two edges S 3 and S 4 close the curve and are both fully in the matching region.
Following similar steps as in Sec.2.1, and noting that inside of V, ∂ t ϕ = − γ(l, t)•∇ϕ, we obtain γj (l, t) where T is the stress tensor associated with F and the dS vectors point to the outside of V.As the integrals over S 3,4 are evaluated in the bulk of the system where the field is homogeneous and takes values ±1 (Fig. 5), the stress tensor identically vanishes on both of these integration regions, which do not contribute to the force.To calculate the contributions from the integrals over S 1,2 , we assume that the curve γ is sufficiently smooth such that its local curvature satisfies κ(l, t)a ≪ 1. Denoting t(l, t) and n(l, t) the local tangential and normal vectors to γ, it is clear that In what follows, we thus only retain the dominant contribution ∂ nϕ ≡ n • ∇ϕ to ∇ϕ, and neglect ∂ tϕ ≡ t • ∇ϕ.Moreover, as the edges S 1,2 are orthogonal to γ in l, we have dS 1 = −dx n t(sl, t) and dS 2 = dx n t(l + dl, t), where hereafter x n denotes the running coordinate normal to the interface.Then, where F is the free energy density associated with (38), while the first contribution on the r.h.s. was eliminated noting that the stress tensor of a scalar theory is always symmetric.After expressing dS 2,j T ji in a similar way, we obtain where the ±∞ integration boundaries follow from the fact that r 0 ≫ a and F = 0 in the bulk phases.The constant σ represents the free energy per unit length of the interface, and thus corresponds to the surface tension.We now calculate the friction tensor in (39).At the leading order in dl, where we have again used the fact that ∂ tϕ is subdominant and that ∇ϕ vanishes in the bulk.Defining the friction as ζ = ∞ −∞ dx n (∂ nϕ) 2 , we finally recover the Allen-Cahn equation [4] where κ = | t′ | is the local curvature of the interface.As expected, we thus find that the transverse velocity of the interface is set by the local curvature.The surface tension and friction are moreover fully determined by the interface profile, which can be obtained from (38) by solving ∂ 2 n nϕ − ∂ ϕ V (a, ϕ) = 0. Multiplying both sides by ∂ nϕ and integrating across the interface, it is straightforward to show that 1 which is the well known relationship between surface tension and friction for domain walls [4].Defining γ⊥ (l, t) ≡ n • γ, and choosing the arc-length parametrization for γ, we end up with the compact formula The generalization of Eq. ( 45) to three dimensions, where domain walls take the from of two-dimensional surfaces, is straightforward.Parametrizing the domain wall with the coordinates l, any point Λ(l, t) on the manifold then evolves as where n is the vector normal to the interface (assuming that the surface is closed, n points to the outside) and H is its mean curvature.Equation ( 46) is known as the Allen-Cahn equation [50], and states that the domain wall velocity is only determined by its local mean curvature since the friction and the surface tension are proportional to each other.

The general equation of motion for disclination lines
We now discuss the situation where the order parameter ϕ is two-dimensional and evolves in a threedimensional space.In this scenario, defects take the form of one-dimensional manifolds around which the circulation of the order parameter orientation is topologically constrained, and that are usually referred to as vortex-or disclination lines [61].Similar structures are for example found in superfluid Helium [62] or in the displacement field of sheared glasses [63].Similarly to point defects, the motion of disclination lines has been mostly studied in the context of the classical Ginzburg-Landau theory via matched asymptotics [28] or using a variational approach based on the Rayleigh dissipation function [22].Here, we first show how to obtain the equation of motion for disclination lines arising in a system whose dynamics minimizes an arbitrary free energy, while the application to the Ginzburg-Landau framework is presented in a second part.
As before, we denote by γ(l, t) the vector defining the defect line, and choose an adequate parametrization such that |γ ′ (l, t)| = 1.In three dimensions, the curve γ(l, t) is locally defined by the Frenet-Serret frame such that t is tangent to γ, n and b are respectively the normal and binormal unit vectors, and κ denotes the local curvature of the curve.As before, we assume a separation of scales between the core and bulk physics, such that we work in the limit where the core radius a satisfies κa ≪ 1.Given the geometry of the problem, the natural choice for the volume V is an infinitesimal cylinder with top and bottom surfaces S 1 and S 2 orthogonal to the curve γ respectively in l and l + dl, while the whole lateral surface S L lies in the matching region (see a sketch in Fig. 6).The application of the free energy variation method in this case combines the results obtained in the previous sections for point-defects and domain walls.We first evaluate the defect friction under the assumption (iii) of a rigid core.When the local radius of curvature of γ is large as compared to the size of the core, we obtain at leading order in dl, where D r 0 denotes the disk of radius r 0 orthogonal to t.The second equality is obtained by neglecting the part of ∇ϕ tangential to t and assuming that the tensor coming from the integral on the r.h.s. is isotropic on the plane normal to γ(l).This assumption is motivated by observing that, in the limit of large loop curvature, the integral in the second line of (48) has properties analogous to the friction integral studied for topological defects in Sec. 2, which we have shown to lead to isotropic friction tensor in most situations.Analogously to domain walls, the integration of the stress tensor over the surfaces S 1 and S 2 , respectively orthogonal to γ(l) and γ(l + dl), gives the surface tension contribution: where σ = Dr 0 d 2 x n F .At this stage, we can already note that, contrary to the case of domain walls, both the friction and surface tension depend on the matching variable r 0 and the core size a.In fact, we will show later that they are both logarithmically divergent in r 0 , similarly to the situation studied in Sec.2.2 for point-defects.
We now calculate the contribution of the lateral surface S L .Contrary to domain walls, the bulk stress tensor does not vanish here because of the presence of the longrange modulations of θ which mediate interactions between the defect lines.Denoting (r, φ, t) as the local cylindrical frame centred in γ, the surface element on S L can be expressed as dS i = −r 0 dφdlϵ ijk tj φk where ϵ ijk is the 3D totally antisymmetric Levi-Civita tensor.Hence, where we have used the fact that the integration is done in the matching region to substitute T by its bulk counterpart.Combining Eqs.(48,49,50), and denoting γ⊥ ≡ (I − tt ) γ, we finally obtain Similarly to point-defects, Eq. ( 51) explicitly depends on the matching variable r 0 .As r 0 is arbitrary, it must simplify when a bulk free energy is specified.

Isotropic systems: the large ln(a) expansion
We now consider a bulk free energy analogous to that used in Sec.2.2: The corresponding dynamical evolution of the orientation field θ(x, t) is given by the three-dimensional diffusion equation.As this equation is linear, it can formally be solved using the approach presented in [45], albeit leading to a more complicated solution than in two-dimensions.Here, we instead note that taking the derivative of the friction coefficient and the surface tension with respect to r 0 , we obtain where the integrals are evaluated at the leading order in a from the static straight line defect solution.This solution simply corresponds to an extension of the point-defect solution into the third dimension: where r is the minimal distance between x and γ.Plugging (54) in Eq. ( 53) then straightforwardly leads to ζ ≃ σ ≃ πs 2 ln(r 0 /a).We now take a → 0 and perform a low mobility expansion as was done in Sec.2.3.1.Therefore, we calculate the force term in (51) via the static line defect solution which for general curves is given by the Biot-Savart law: In particular, in the specific case of a straight disclination line aligned along ẑ, the expression (55) simplifies into (54).Since we integrate the stress tensor over a circle of radius r 0 → 0 enclosing the line around γ, we write the phase field for x → γ as ∇θ = ∇θ d + ∇θ c .The discontinuous part is obtained by expanding (55) in the vicinity of γ, giving ∇θ d (x) = ∇θ ssd (x) − κs ln(κr 0 ) b/2, while the continuous part θ c accounts for other singularities or externally imposed boundary conditions.After evaluating the integral of the stress tensor explicitly, we obtain ln L a γ⊥ = κ ln While the matching scale r 0 in the expression of the surface tension was simplified due to the curvature contribution to ∇θ d , as previously discussed for slow defects the r 0 dependency of the friction can be replaced by any arbitrary macroscopic scale L. Equation ( 56) can be interpreted as the generalization of ( 17) (in the large ln(a) limit) to the case where the singularities extend along a third dimension.Importantly, it includes an additional term on the r.h.s.leading to the collapse of disclination lines with finite curvature, while the Peach-Koehler-type force describes how distinct lines interact.
We note that since the friction coefficient and the surface tension share the same scaling with ln(a), the importance of the interaction term to the dynamics of disclination lines depends on their curvature.Indeed, assuming κ to be O(1) implies that the line velocity is also O(1), while the interaction term on the r.h.s. of ( 56) is O(ln −1 (a)).Hence, the dynamics is dominated by the tension and, after taking L ≈ κ −1 , we recover an equation similar to that ruling the motion of domain walls: Expressions similar to (57) have been derived by several authors who focused on the dynamics of isolated line defects [28,61].
On the other hand, in the presence of several line defects with low curvature (namely, if κ = O(ln −1 (a)), the interaction term in (56) becomes relevant.In the quasi-static approximation and assuming open boundaries, the term ∇θ c (γ α ) appearing in the equation of motion of the defect α can thus be generally expressed as ∇θ c (γ α ) = β̸ =α ∇θ sd,β (γ α ), where θ sd,β corresponds to the orientation profile generated by the defect β.
In the specific case where all disclination lines are straight and aligned along ẑ, it is straightforward to check from (55) that (56) reduces to the description of the twodimensional Coulomb gas [30].Conversely, for straight disclination lines with arbitrary orientations, Eq. ( 56) takes the form ln L a γ⊥,α = where r αβ ≡ 1 2 (γ α − γαβ ) and γαβ ≡ argmin γ β (|γ α − γ β |).The second contribution to the r.h.s. of (58) aligns or anti-aligns the direction of the lines α and β (parametrized by their tangent vector) whenever s α s β < 0 or s α s β > 0, respectively.Hence, due to the first term the effective interaction between disclination lines is always attractive regardless of their respective charges.

Disclination loops
Another type of structure of interest are loop singularities.To study their dynamics, we consider a circular ring α centred at position X α with radius R α .We define the associated dipole moment as where bα is the binormal unit vector orthogonal to the loop plane.Taking the time derivative of Eq. ( 59) and using (56), we find that the loop radius and orientation obey while the dynamics of X α can be obtained by simply averaging (56) over the ring contour.Note that we have chosen to replace L with R α in (60,61), as here 1/R α sets the natural scale for the line velocity.We first consider the case where the loop is subject to a uniform orientation gradient ∇θ c = k, as was done in Sec.2.2.2.It is then straightforward to check that Ẋα = 0.Moreover, we conclude from Eq. ( 61) that the ring anti-aligns (aligns) with k when s α > 0 (s α < 0).In either case, the presence of the externally imposed gradient leads to a positive contribution to the growth of the ring radius in (60).Hence, beyond a threshold radius satisfying R c / ln(R c /a) = |s α |/(2|k|), this expanding force overcomes the curvature-induced contraction and the ring expands indefinitely.
When ∇θ c results from other defect loops located far away from X α , it is given at the leading order in the far-field approximation by a superposition of dipoles, namely, where r αβ = X α − X β .In this case too, the loop is globally immobile but only contracts and rotates according to In contrast to straight disclination lines (Eq. ( 58)), the interaction force between loops scales as the inverse of the cube of their separation.In particular, the alignment dynamics of disclination loops resembles that of magnetic dipoles, such that loops with equal and opposite charge will anti-align and align, respectively.Equation ( 63) further shows that the contraction of the ring is also affected by the presence of nearby loops.However, given two loops α and β, the requirement for the interaction contribution to compete with the curvature term is R α R 2 β ≃ r 3 αβ ln(R α /a), which clearly breaks the far field assumption r αβ ≫ R α , R β .Hence, in most cases distant loops will only rotate and self-annihilate over a finite time.

Discussion
The reduced particle-field description of the large-scale dynamics of the field theory (1) requires a considerably reduced number of degrees of freedom.Our previous work [45] was devoted to determining the perturbation of the order parameter introduced by a collection of moving defects.Here, we have focused on the other facet of the problem, namely, determining how defects are set into motion by the order parameter landscape.As we have demonstrated, for a large class of systems the defect dynamics takes a universal form (Eqs. ( 11) and ( 51)) since it is largely determined by the large-scale features of the theory, and should thus be qualitatively insensitive to the microscopic details of the system of interest.One important requirement for (11,51) to hold is that the underlying dynamics is described by a translationally invariant free energy functional.We note, however, that coupling the order parameter to a smoothly varying external field should not introduce major difficulties to the derivation presented in Sec.2.1.As a number of nonequilibrium (active) field theories take the form of an equilibrium-like order parameter dynamics coupled to an external flow [64,65], density [66], or chemical [67,68] field, some progress could be achieved in these cases via the approach outlined in Sec.2.1.When the order parameter evolution does not satisfy this structure, however, one cannot a priori exclude that the functional form of the defect equations of motion depends on the details of the microscopic scale physics [69,70].
In Sec.2.2, we have shown how the nonlinearity of the defect friction leads to a generic dependency on the length scale λ(t) (Eq.( 17)).Hence, a quantitative characterization of the many-body defect dynamics seems out of reach, since it would in principle require to evaluate this scale, which depends on the full history of the defect motion [45].As this task is typically unfeasible-except for relatively simple configurations-we have focused on slowly moving defects and showed how memory effects disappear in this case.Within this approximation, λ(t) can be replaced by any other relevant macroscopic scale of the problem, such that the equation of motion for the defect will be fully determined.
Taking the slow defect approximation moreover allowed us to address the case of defects evolving in a medium exhibiting elastic anisotropy in Sec.2.3.Working perturbatively in the anisotropy parameter α, we were able to quantify its linear order contribution to the bending of defect trajectories due to the energetic cost of mismatching.Interestingly, we also discovered that the + 1 2 nematic defect mobility becomes anisotropic, as it depends on the relative orientation of the comet shape of the defect relatively to the background orientation field.We moreover argue in Appendix A that this property may extend to −1 defects in polar systems.As backflows also substantially affect the dynamics of defects [46,71], observing the effect of elastic anisotropy experimentally may only be achieved when hydrodynamic effects are negligible, such as in Langmuir monolayers [47].
The low mobility expansion also made the derivation of the closed form of the equation of motion for disclination lines in Sec.3.2 more straightforward.The Peach-Koehler force on the r.h.s. of Eq. ( 56) takes a similar form as that proposed for disclination lines in three-dimensional nematics [72], suggesting an interesting connection between the two settings.
Our results moreover indicate that the proportionality of the friction coefficient and the surface tension observed for domain walls extends to diclination lines.Contrary to point defects in two dimensions, we have shown that line defects freely rotating in three dimensions will always annihilate regardless of their charge.Furthermore, the study of closed loops done in Sec.3.2.3reveals that they align their moment similarly to magnetic dipoles, while in the absence of an externally imposed phase gradient they always self-annihilate due to surface tension.
Further extensions of the formalism presented in this work can include taking into account the roles of backflow [71], fluctuations [73], or curved geometry [74].Another interesting extension of the method we propose would be to study the motion of dislocations [54] by using a field theoretical treatment of elasticity [75], or more sophisticated descriptions [76,77].While some of these developments may be technically challenging, they should not entail any additional conceptual difficulty.

Appendix A. Symmetry of the mobility matrix from the defect symmetries
We claim in the main text that the mobility of a defect is generally isotropic.We present here an argument supporting this claim, based on the symmetries of the defects.
An isolated s-charged defect is n-fold symmetric, with n = |s − 1| and n = 2|s − 1| for polar and nematic order parameters, respectively.When n ̸ = 1, the symmetry is nontrivial so that the order parameter field resulting from the presence of the defect must be invariant under 2π/n-rotations.In practice, this property implies that the defect mobility matrix µ satisfies: where R n corresponds to the rotation matrix with angle 2π/n.For n > 2 and noting that µ is a symmetric matrix, this relation is satisfied only if µ = µI, with I the identity.On the other hand, for n = 1 or 2 (A.1) is always satisfied.We thus conclude that the mobilities of defects with charge s = 1 in polar and s = − 1 2 in nematic systems must be isotropic.Inversely, the mobilities of + 1 2 -charged nematic and −1-charged polar defects may, in principle, be anisotropic.Focusing on nematic systems and defining p as the polarization of the + 1  2 defect (see Sec. 2.3.3),we express its mobility as where, as in the main text, r 0 and a denote the matching and core length scales.It is also worth noting that the polarization p is defined respectively to the orientation of the background nematic field (see Eq. ( 36)), and is therefore not an independent degree of freedom of the dynamics.As discussed in Sec.2.3.3,µ ∥ ̸ = µ ⊥ generally occurs in two-dimensional nematics with elastic anisotropy.
In polar systems, the orientation of the 2-fold symmetric −1 defect can be defined with an appropriate director, which will lead to a mobility matrix of the form of (A.2).In principle, a similar line of thought can be applied to other 1 and 2-fold symmetric defects.However, these are usually unstable due to their charge (like the s = 2 polar defect), and can thus be observed only with specific choices of boundary conditions [78].

Appendix B. Integration of the stress tensor for the isotropic case
In this appendix, we derive the dominant term on the r.h.s. of Eq. ( 11) for the linear Ginzburg-Landau theory in the limit r 0 → 0. As sketched in the main text, we split the near-field orientation gradient into a discontinuous and continuous parts: ∇θ = ∇θ d + ∇θ c .Substituting this expression into T bulk , and keeping only the nonvanishing contribution in the limit r 0 → 0, we have Using the expression of the discontinuous part given in (15), the mixed term part becomes where in the second equality we have used the identity ϵ il ϵ mj = δ im δ jl − δ ij δ lm , while in the next one was obtained by using the continuity of ∇θ c to bring it outside of the integral for vanishing r 0 .For the remaining contribution, using (15) we find after straightforward calculations Equation.( 16) is finally obtained after rearranging the terms of this equation.

Appendix C. Calculation of the pairwise force between defects in an anisotropic medium
In this appendix, we calculate the integral on the r.h.s. of Eq. ( 28) for a pair of defects evolving in a system described by the bulk free energy (29) perturbatively in the anisotropy parameter α.Without loss of generality, we consider two oppositely charged defects at positions q ± = ±q x.The static equation of motion for the orientation field θ deriving from (29) reads where − sin 2θ cos 2θ cos 2θ sin 2θ .
We assume α to be small, so that we write the solution of Eq. (C.1) at linear order as θ(x) ≃ θ 0 (x) + αθ 1 (x), where θ 0 solves the isotropic (α = 0) problem.Namely, θ 0 (x, q) = θ d (x, q) + θ c , (C.2) where θ c denotes the continuous part of the solution which has to be constant in the static case, while θ d (x, q) = 1 2 arg(x − q x) − 1 2 arg(x + q x).In turn, the first order perturbation θ 1 is solution of Equation (C.3) can in principle be solved by inverting the Laplacian.
Here, we however adopt an alternative approach.We first note that, since the bulk theory is conformal, the solution of (C.3) can generally be written as θ 1 (x, q, θ c ) = ϑ 1 (y, θ c ), where y = x/q and ϑ 1 (y, θ c ) is the solution of (C.3) for two defects at positions ± x.From the following property of the Q and Q tensors, (C.6)These two equations can be formally solved using the Green's function of the Laplacian.As will appear clear below, the solutions are importantly independent of q and θ c .
We are now able to evaluate the integral of the stress tensor of Eq. (28).Since this integral is independent of the integration contour, we choose it to be the straight line along y passing by x = 0.The r.h.s. of Eq. ( 28) can then be written as We moreover also expand T sb as T sb = T 0 + αT 1 so that T 0 is the stress tensor of the unperturbed theory and T 1 the corresponding linear order perturbation.The integral for T 0 is straightforward, and gives the usual Coulomb interaction.To calculate the which corresponds to a pair of defects with charges s = ±n −1 at positions ±q 0 on the complex plane.To reduce the importance of finite size effects, in the anisotropic case the initial distance between the defects was taken to be 2|q 0 | = 64, i.e. much less than the linear dimensions of the whole system.The results presented in Sec.2.2.3 (Fig. 2) were obtained from simulations of (D.2) with n = 1, which thus reduces to the classical Ginzburg-Landau theory.In turn, simulations corresponding to Sec. 2.2.3 (Fig. 4) were performed with n = 2 and varying the anisotropy parameter α.As in both cases (D.3) is not a solution of Eq. (D.2) on the torus, we initially let the dynamics evolve over a simulation time corresponding to roughly 10% of the time need by the defect pair to annihilate before starting the data acquisition.Defect tracking was performed by computing the circuitation of the phase of f around four neighbouring boxes of the numerical grid.

Figure 3 :
Figure 3: Configurations of a pair of nematic defects corresponding to θ q = 0 (a), θ q = π 4 (b), and θ q = π 2 (c).The grey lines indicate the local orientation of the nematic field.The polarization p of the + 1 2 defect is shown with the blue arrows, while red triangles are used to indicate the position of the − 1 2 defect.

Figure 5 :
Figure 5: A domain wall (grey-shaded region) separating two homogeneous regions with φ = ±1.The choice of boundary to calculate the variation is shown in red.The surfaces S 1,2 define the infinitesimal segment, while S 3,4 show the matching region.

Figure 6 :
Figure 6: Choice of boundary for defect lines.