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

To find the equation of motion for the defect core, we first express the variational equation (1) in an integral form as

$$\begin{equation} \int \mathrm{d}^2\boldsymbol{x} \; \delta \boldsymbol{\phi} \cdot \partial_t \boldsymbol{\phi} = -\delta {\cal F} , \end{equation} \tag{ 2 }$$

where $\delta \boldsymbol{\phi}$ is a fixed boundary condition perturbation. By construction, equation (2) is valid for any infinitesimal $\delta \boldsymbol{\phi}$. Here, we consider a specific type of perturbation:

$$\begin{equation} {\delta \boldsymbol{\phi}\left(\boldsymbol{x},t\right) = f\left(\boldsymbol{x},t\right) \left(\delta\boldsymbol{q} \cdot \nabla\right) \boldsymbol{\phi}\left(\boldsymbol{x},t\right),} \end{equation} \tag{ 3 }$$

where $\delta \boldsymbol{q}$ is an infinitesimal vector and $f(\boldsymbol{x},t)$ is a smooth interpolating envelope function equal to one at $\boldsymbol{x} = \boldsymbol{q}(t)$, which decays quickly to zero for $|\boldsymbol{x} - \boldsymbol{q}(t)| \gt 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 $\delta \boldsymbol{\phi}$ at the system boundaries and at the other defect cores. Discarding this prefactor would in fact lead to additional ${\cal 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 right hand side (r.h.s.) of equation (2) into contributions from inside and outside ϒ, corresponding respectively to the core and the bulk. Namely,

$$\begin{equation} \int_{\mathrm{core}}\mathrm{d}^2 \boldsymbol{x} \, \delta \boldsymbol{\phi} \cdot \partial_t \boldsymbol{\phi} + \int_{\mathrm{bulk}}\mathrm{d}^2\boldsymbol{x}\, \delta \boldsymbol{\phi} \cdot \partial_t \boldsymbol{\phi} = - \delta {\cal F}_{\mathrm{core}} -\delta {\cal F}_{\mathrm{bulk}}. \end{equation} \tag{ 4 }$$

We now re-express the bulk free energy variation as

$$\begin{equation} \delta{\cal F}_{\mathrm{bulk}} = \int_{\mathrm{bulk}}\mathrm{d}^2\boldsymbol{x} \, \delta \boldsymbol{\phi} \cdot \frac{\delta{\cal F}_{\mathrm{bulk}}}{\delta \boldsymbol{\phi}} - \oint_{\Upsilon} \mathrm{d} S_i \, \delta\boldsymbol{\phi} \cdot \frac{\partial{F}_{\mathrm{bulk}}}{\partial \left(\partial_{i} \boldsymbol{\phi}\right)}, \end{equation} \tag{ 5 }$$

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}_\textrm{bulk}(\boldsymbol{\phi},\nabla\boldsymbol{\phi})$ does not depend on higher derivatives of φ ⁴ . As $\delta\boldsymbol{\phi}$ is zero outside the core region, this surface term only includes a contribution from the boundary ϒ, while the preceding minus sign implies that $\mathrm{d}\boldsymbol{S}$ points to the outside of the core.

To compute $\delta{\cal F}_\textrm{core}$, we note that in the core $\delta\boldsymbol{\phi} = (\delta\boldsymbol{q} \cdot \nabla)\boldsymbol{\phi}$ as by construction $f(|\boldsymbol{x}-\boldsymbol{q}| \unicode{x2A7D} r_0) = 1$. Therefore, $\delta\boldsymbol{\phi}$ in the core corresponds to an infinitesimal translation. Since, from (i), ${F}_\textrm{core}$ is translationally invariant, it follows that

$$\begin{equation*} \delta{F}_{\mathrm{core}} = \frac{\partial {F}_{\mathrm{core}}}{\partial \phi_b}\left(\delta\boldsymbol{q}\cdot\nabla\right)\phi_b +\frac{\partial {F}_{\mathrm{core}}}{\partial \left(\partial_{j}\phi_b\right)}\left(\delta\boldsymbol{q}\cdot\nabla\right)\partial_{j}\phi_b + {\ldots} = \left(\delta\boldsymbol{q}\cdot\nabla\right) {F}_{\mathrm{core}} , \end{equation*}$$

where the dots stand for dependencies of ${F}_\textrm{core}$ on higher order derivatives of φ , if any. Hence, $\delta{\cal F}_\textrm{core}$ is an exact differential, and therefore, we obtain

$$\begin{equation} \delta {\cal F}_{\mathrm{core}} = \oint_{\Upsilon} \mathrm{d}\boldsymbol{S} \cdot \delta \boldsymbol{q} \, {F}_{\mathrm{core}} = \oint_{\Upsilon} \mathrm{d}\boldsymbol{S} \cdot \delta \boldsymbol{q} \, {F}_{\mathrm{bulk}}, \end{equation} \tag{ 6 }$$

where the second equality results from the fact that ϒ belongs to the matching region. Putting together equations (4)–(6), we find that

$$\begin{equation} \int_{\mathrm{core}}\mathrm{d}^2\boldsymbol{x} \,\partial_t \boldsymbol{\phi} \cdot \delta \boldsymbol{\phi} = \delta q_j \oint_{\Upsilon}\mathrm{d} {S}_i \, {T}_{{\mathrm{bulk}},ij} , \end{equation} \tag{ 7 }$$

where we have defined the canonical stress tensor of the bulk theory [52, 53]:

$$\begin{equation*} {T}_{{\mathrm{bulk}},ij} \equiv \frac{\partial{F}_{\mathrm{bulk}}}{\partial \left(\partial_{i} \boldsymbol{\phi}\right)}\cdot \left(\partial_{j}\boldsymbol{\phi}\right) - \delta_{ij} {F}_{\mathrm{bulk}}. \end{equation*}$$

Hence, the r.h.s. of equation (7) corresponds to the net momentum flux through the matching region, and is solely determined by the large-scale bulk physics.

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 equation (7) holds for an arbitrary choice of curve ϒ in the matching region. Choosing without loss of generality ϒ to be a circle of radius r₀ around the singularity, we thus expect that the final result will be independent of r₀. Using the rigid core assumption (iii), it is clear that in the core $\partial_t \boldsymbol{\phi} = -(\boldsymbol{v} \cdot \nabla)\boldsymbol{\phi}$ ⁵ , such that the l.h.s. of equation (7) is given by

$$\begin{equation} \int_{\mathrm{core}}\mathrm{d}^2\boldsymbol{x} \,\partial_t \boldsymbol{\phi} \cdot \delta \boldsymbol{\phi} = - v_i \delta q_j R_{ik}\left(t\right)R_{jl}\left(t\right)\int_{{\cal D}_{r_0/a}}\mathrm{d}^2\boldsymbol{y} \, \partial_{k} \boldsymbol{\bar{\phi}} \left(\boldsymbol{y}\right) \cdot \partial_{l} \boldsymbol{\bar{\phi}} \left(\boldsymbol{y}\right) , \end{equation} \tag{ 8 }$$

where we have used the change of variable $\boldsymbol{y} = \boldsymbol{R}^{-1}(t)(\boldsymbol{x} - \boldsymbol{q}(t))/a$ and ${\cal D}_{r_0/a}$ is the disk of radius $r_0/a$ centered at 0. Equation (8) defines the effective friction tensor of the defect dynamics:

$$\begin{equation} \qquad \zeta_{ij}\left(\frac{r_0}{a},t\right) \equiv R_{ik}\left(t\right)R_{jl}\left(t\right)\int_{{\cal D}_{r_0/a}}\mathrm{d}^2\boldsymbol{y} \, \partial_{k} \boldsymbol{\bar{\phi}} \left(\boldsymbol{y}\right) \cdot \partial_{l} \boldsymbol{\bar{\phi}} \left(\boldsymbol{y}\right) \equiv R_{ik}\left(t\right)R_{jl}\left(t\right) \bar{\zeta}_{kl}\left(\frac{r_0}{a}\right). \end{equation} \tag{ 9 }$$

We note that, since $\boldsymbol{\bar{\phi}}$ is fully determined by the core structure, it is independent of the macroscopic defect variables such at its position, velocity and orientation. The tensor $\bar{\boldsymbol{\zeta}}$ 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₀ and parametrizing y with the polar coordinates $(r_0,\varphi)$, we find that

$$\begin{equation} \frac{\mathrm{d} \bar{\zeta}_{ij}}{\mathrm{d} r_0} = r_0 \int_0^{2\pi} \mathrm{d}\varphi \, \left. \partial_{i} \boldsymbol{\bar{\phi}} \left(r_0,\varphi\right) \cdot \partial_{j} \boldsymbol{\bar{\phi}} \left(r_0,\varphi\right)\right|_{\mathrm{ssd}}, \end{equation} \tag{ 10 }$$

where the 'ssd' subscript indicates that the integrand is calculated from the static single defect solution of the bulk theory, since the r.h.s. 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 ${\cal F}_\textrm{core}$, we show in appendix A that in most relevant situations ζ is often isotropic and thus independent of the defect orientation (see, however, section 2.3.3 for a counter example). In the following, we therefore keep the time dependency of ζ implicit.

Gathering the results accumulated so far, and noting that equation (7) must be valid for all $\delta\boldsymbol{q}$, we find that the defect equation of motion takes the compact form

$$\begin{equation} \boldsymbol{\zeta}\left(\frac{r_0}{a}\right) \boldsymbol{v} = - \oint_{{\cal C}_{r_0}}\boldsymbol{T}^\mathrm{T}_{{\mathrm{bulk}}} \, \mathrm{d} \boldsymbol{S}, \end{equation} \tag{ 11 }$$

where ${\cal C}_{r_0}$ stands for the circle of radius r₀ centered 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 equation (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 equation (11) depend on the matching variable r₀, while the actual equation of motion of the defect must be independent of it, since r₀ is arbitrary. In fact, we show below that eliminating r₀ 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.

The simplest choice for ${\cal F}$ is certainly the Ginzburg–Landau free energy

$$\begin{equation} {\cal F}^{\mathrm{GL}} = \int \mathrm{d}^2 \boldsymbol{x} \, \left[ \frac{1}{2}|\nabla \boldsymbol{\phi}|^2 + \chi\left(1-|\boldsymbol{\phi}|^2\right)^2 \right] . \end{equation} \tag{ 12 }$$

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 \simeq \chi^{-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

$$\begin{equation} {\cal F}_{\mathrm{bulk}}^{\mathrm{GL}} = \frac{1}{2}\int \mathrm{d}^2\boldsymbol{x}\, |\nabla \theta|^2 . \end{equation} \tag{ 13 }$$

The dynamics of θ is thus ruled by the diffusion equation, while the stress tensor associated with (13) is given by $T_{\textrm{bulk},ij}^\textrm{GL}(\theta) = (\partial_{i}\theta)(\partial_{j}\theta)-\frac{1}{2}\delta_{ij}|\nabla \theta|^2$. In particular, we showed previously [45] that the orientation field gradient induced by a defect of charge s following a trajectory $\boldsymbol{q}(t)$ with velocity $\boldsymbol{v}(t)$ is given by

$$\begin{equation} \nabla\theta\left(\boldsymbol{x},t\right) = -\frac{s}{2}\boldsymbol{\epsilon} \int_{-\infty}^{t} \frac{\mathrm{d} t^{\prime} }{\left(t-t^{\prime}\right)} \left[\nabla+\boldsymbol{v}\left(t^{\prime}\right)\right]\mathrm{e}^{-\frac{|\boldsymbol{x}-\boldsymbol{q}\left(t^{\prime}\right)|^2}{4\left(t-t^{\prime}\right)}}, \end{equation} \tag{ 14 }$$

where ε denotes the two-dimensional antisymmetric Levi–Civita tensor

$$\begin{equation*} \boldsymbol{\epsilon} = \left(\begin{array}{cc} 0 & -1 \\ 1 & 0 \end{array}\right). \end{equation*}$$

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 equation (11) in the limit of a → 0 and $r_0 \to 0$, keeping $r_0 \gg a$. To obtain the r.h.s., we note that the presence of a defect leads to a divergence of $|\nabla \theta| \sim a^{-1}$ at the core, such that for r₀ small the contour integral is dominated by the discontinuous part of $\nabla \theta$. Although the general expression (14) is nonlocal in time, the discontinuous part $\nabla\theta_\textrm{d}$ of $\nabla\theta$ in the vicinity of a defect is formally determined by its instantaneous position and velocity [45], namely

$$\begin{equation} \nabla\theta_{\mathrm{d}}\left(\boldsymbol{r},t\right) \underset{r\to 0}{\simeq} s\boldsymbol{\epsilon}\left[\frac{\hat{\boldsymbol{r}}}{r}+\frac{\boldsymbol{v}\left(t\right)}{2}\ln\left(\frac{r}{\lambda\left(t\right)}\right)-\frac{1}{2}\left(\boldsymbol{v}\left(t\right)\cdot\hat{\boldsymbol{r}}\right)\hat{\boldsymbol{r}}\right], \end{equation} \tag{ 15 }$$

where we have defined $\boldsymbol{r} \equiv \boldsymbol{x} - \boldsymbol{q}(t)$, $r\equiv |\boldsymbol{r}|$ and $\hat{\boldsymbol{r}} \equiv \boldsymbol{r}/r$. The additional length scale $\lambda(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 $\nabla\theta = \nabla\theta_\textrm{d}+\nabla\theta_\textrm{c}$, with $\nabla\theta_\textrm{c}$ accounting for the continuous part of the gradient, we calculate in appendix B the r.h.s. of equation (11), which leads for $r_0\to 0$ to

$$\begin{equation} \left[\zeta_{ij}\left(\frac{r_0}{a}\right) + \pi s^2 \ln\left(\frac{\mathrm{e}^{\frac{1}{2}}\lambda\left(t\right)}{r_0}\right) \delta_{ij} \right]v_j = - 2\pi s\epsilon_{ij} \partial_j \theta_{\mathrm{c}}\left(\boldsymbol{q}\left(t\right),t\right) . \end{equation} \tag{ 16 }$$

It is clear that only the terms on the l.h.s. of equation (16) depend on the matching variable r₀. Hence, we conclude that the friction $\zeta_{ij} = \zeta \delta_{ij}$ while the quantity $r_0 \exp[-\zeta/(\pi s^2)]$ must remain independent of r₀. Furthermore, from its definition (9) the defect friction is independent of any macroscopic scale, which yields

$$\begin{equation*} \zeta\left( \frac{r_0}{a} \right) = \pi s^2 \ln\left( \frac{r_0}{\bar{a}} \right), \end{equation*}$$

where $\bar{a} \propto 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 $\boldsymbol{v} = \textbf{0}$. The equation of motion for the defect therefore reads

$$\begin{equation} \ln\left(\frac{\mathrm{e}^{\frac{1}{2}}\lambda\left(t\right)}{\bar a}\right)\boldsymbol{v}\left(t\right) = -\frac{2}{s}\, \boldsymbol{\epsilon}\, \nabla \theta_{\mathrm{c}}\left(\boldsymbol{q}\left(t\right),t\right). \end{equation} \tag{ 17 }$$

Equation (17) is formally similar to equation (11). However, whereas the latter depends on the arbitrary matching variable r₀ through ζ and the integration contour on the r.h.s. Equation (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 $\bar{a}$. Its value is set by the core free energy ${\cal F}_\textrm{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 $\bar{a} \sqrt{\chi}\simeq 1.126$ [16]. The coefficient $\lambda(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 equation (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 $\lambda(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 equation (17) can be explicitly determined in particular configurations.

We first study the simple case of a defect moving uniformly with velocity $\boldsymbol{v} = \bar{\boldsymbol{v}}$. An experimental realization of this case would for example consist of a defect subject to an imposed spin wave $\theta_\textrm{ext}(\boldsymbol{x}) \equiv \boldsymbol{k} \cdot \boldsymbol{x}$, with k the corresponding wavevector, induced by an external field. Under these conditions, the total angular field is given by the sum of $\theta_\textrm{ext}(\boldsymbol{x})$ and the uniformly moving defect solution to the diffusion equation, namely [31, 35, 45]

$$\begin{equation} \nabla\theta\left(\boldsymbol{r}\right) = \frac{s}{2} \mathrm{e}^{-\frac{1}{2} \bar{\boldsymbol{v}}\cdot\boldsymbol{r}} \boldsymbol{\epsilon} \left[\bar{\boldsymbol{v}} K_0\left(\frac{\bar{v}r}{2}\right)-\bar{v}\hat{\boldsymbol{r}} K_1\left(\frac{\bar{v}r}{2}\right)\right]+\boldsymbol{k}, \end{equation} \tag{ 18 }$$

where K₀ and K₁ are modified Bessel functions of the second kind while, as previously, s denotes the charge of the defect and $\boldsymbol{r} = \boldsymbol{x}-\boldsymbol{q}(t)$. Expanding equation (18) for r → 0 and keeping only nonvanishing terms, we find that $\nabla\theta_\textrm{c} = \boldsymbol{k}$ while the discontinuous part $\nabla\theta_\textrm{d}$ is given by (15) with $\lambda = 4 \bar{v}^{-1} \mathrm{e}^{-\gamma_\textrm{E}}$; and where $\gamma_\textrm{E} \approx 0.578$ denotes the Euler–Mascheroni constant.

Replacing these expressions in equation (17), we thus recover the classical result [16]

$$\begin{equation} \ln\left(4\frac{\mathrm{e}^{\frac{1}{2}-\gamma_{\mathrm{E}}}}{\bar{v}\bar{a}}\right)\bar{\boldsymbol{v}} = -\frac{2}{s}\,\boldsymbol{\epsilon} \,\boldsymbol{k}, \end{equation} \tag{ 19 }$$

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.

We now consider a pair of defects with opposite charges $s_\pm = \pm s$. To proceed further, we note from (17) that the defect mobility $\mu \equiv \zeta^{-1} \sim \ln^{-1}(\lambda/a)$. Hence, so long as the macroscopic scale λ remains much larger than the core radius, defects move relatively slowly with a speed $v = {\cal O}(\mu)$. 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 $\lambda(t)$. Denoting, respectively, $\boldsymbol{q}_\pm(t)$ and $\boldsymbol{v}_\pm(t)$ as the positions and velocities of the two defects, it then follows from (17) that

$$\begin{equation} \qquad\qquad \ln\left(\frac{\mathrm{e}^{\frac{1}{2}}\lambda_\pm}{\bar{a}}\right)\boldsymbol{v}_\pm = -\frac{2}{s_\pm}\boldsymbol{\epsilon}\nabla\theta_{\mp} \left(\boldsymbol{q}_\pm\left(t\right),t\right) = \mp\frac{\hat{\boldsymbol{q}}}{q} + \ln\left( \frac{\lambda_\mp}{2q} \right)\boldsymbol{v}_\mp + \left(\boldsymbol{v}_\mp\cdot\hat{\boldsymbol{q}}\right)\hat{\boldsymbol{q}} , \end{equation} \tag{ 20 }$$

where $\theta_{\pm}$ denote the orientation field landscapes generated by the positively and negatively charged defects, while the second equality was obtained by replacing $\theta_\pm$ with (15). $\lambda_\pm$ in equation (20) denote the scales associated with the ± defects, while we have also defined $\boldsymbol{q} \equiv \frac{1}{2}(\boldsymbol{q}_+ - \boldsymbol{q}_-)$. It is a straightforward exercise to show that (20) implies that $\boldsymbol{v}_+ = -\boldsymbol{v}_- = v \hat{\boldsymbol{q}}$, such that the speed of both defects follows

$$\begin{equation} v = -\frac{\mu_{\mathrm{eff}}}{q}, \qquad \mu_{\mathrm{eff}} \equiv \left[\ln\left( \frac{\mathrm{e}^{\frac{3}{2}}\lambda_+\lambda_-}{2\bar{a}q } \right)\right]^{-1}. \end{equation} \tag{ 21 }$$

As q(t) is the only macroscopic scale of the problem, we propose the ansatz: $\lambda_\pm(t) \propto q(t)$. Hence, it follows that $\dot{\mu}_\textrm{eff} \sim v \mu_\textrm{eff}^2 \sim \mu_\textrm{eff}^3$ 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 [45], while the resulting expression for $\nabla\theta$ is again formally given by (15) with $\lambda_+(t) = \lambda_-(t) = 2\sqrt{2}q(t) \mu_\textrm{eff}^{-\frac{1}{2}}\mathrm{e}^{-\frac{1}{2}(1 + \gamma_\textrm{E})}$. Combining this result with the expression of µ_eff given in (21) yields the following self-consistent condition

$$\begin{equation} \mu_{\mathrm{eff}}\mathrm{e}^{\mu_{\mathrm{eff}}^{-1}} = \frac{4 q}{\bar{a}}\mathrm{e}^{\frac{1}{2} - \gamma_{\mathrm{E}}} \end{equation} \tag{ 22 }$$

whose solution can be expressed in terms of the velocity as $\mu_\textrm{eff}^{-1} = \ln \left[4\mathrm{e}^{\frac{1}{2}-\gamma_\textrm{E}}/(v\bar{a})\right]$. 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 $\mu_\textrm{eff} \, (\propto v)$. Evaluating the solution of equation (22) up to terms ${\cal O}\left( \ln[\ln(q/\bar{a})]/\ln(q/\bar{a}) \right)$, we thus obtain

$$\begin{equation} \ln\left[ \frac{4 q}{\bar{a}}\mathrm{e}^{\frac{1}{2}-\gamma_{\mathrm{E}}}\ln\left(\frac{4\mathrm{e}^{\frac{1}{2}-\gamma_{\mathrm{E}}}}{\bar{a}} q\right)\right]v = -\frac{1}{q}. \end{equation} \tag{ 23 }$$

Equation (23) fully determines the dynamics of the defect pair.

Extensive discussions on scale- [4, 37, 47] or velocity- [18, 25, 27, 31, 35–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 = \pm 1$ defect pair annihilation in a vectorial field described by the Ginzburg–Landau model (equation (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 behavior of v(q) at large q that depends on the initial defect separation. Past this transient, v(q) exhibits a behavior independent of the initial defect separation, which we find in clear departure from the $\propto q^{-1}$ scaling predicted by the constant mobility approximation (dashed lines in figure 2). Conversely, a parameter-free comparison with equation (23) reveals excellent agreement with the measured relation v(q) (blue curves in figure 2).

Zoom In Zoom Out Reset image size

In the low mobility approximation, equation (20) can be generalized to an arbitrary number of defects as

$$\begin{align} s_\alpha^2 \ln\left(\frac{\mathrm{e}^{\frac{1}{2}}\lambda_\alpha}{\bar{a}}\right)\boldsymbol{v}_\alpha = \sum_{\beta \ne \alpha} s_\alpha s_\beta\left[ \frac{\hat{\boldsymbol{q}}_{\alpha \beta}}{q_{\alpha \beta}} - \ln\left( \frac{\lambda_\beta}{2q_{\alpha \beta}} \right)\boldsymbol{v}_\beta - \left(\boldsymbol{v}_\beta\cdot\hat{\boldsymbol{q}}_{\alpha \beta}\right)\hat{\boldsymbol{q}}_{\alpha \beta} \right] , \end{align} \tag{ 24 }$$

where $\boldsymbol{q}_\alpha$, $\boldsymbol{v}_\alpha$ and s_α denote the position, velocity and charge of the αth defect, respectively, while $\boldsymbol{q}_{\alpha \beta} \equiv \frac{1}{2}(\boldsymbol{q}_\alpha - \boldsymbol{q}_\beta)$. Rearranging the terms, we recast (24) into a more compact form:

$$\begin{align} \boldsymbol{v}_\alpha = \sum_\beta \boldsymbol{{\cal M}}_{\alpha \beta} \left(- \nabla_{\boldsymbol{q}_\beta} U_{\mathrm{c}}\right), \end{align} \tag{ 25 }$$

where $U_\textrm{C} \equiv -\sum_{\alpha \beta} s_\alpha s_\beta \ln|q_{\alpha \beta}|$ is the two-dimensional Coulomb potential, while the mobility matrix $\boldsymbol{{\cal M}}$ is defined through $\sum_\gamma \boldsymbol{{\cal M}}_{\alpha \gamma} \boldsymbol{{\cal Z}}_{\gamma \beta} = \delta_{\alpha \beta}\boldsymbol{{\cal I}}$, where the friction matrix $\boldsymbol{{\cal Z}}$ is defined as

$$\begin{equation*} \boldsymbol{{\cal Z}}_{\alpha \beta} \equiv s_\alpha s_\beta \left\{\begin{array}{ll} \ln\left(\frac{\mathrm{e}^{\frac{1}{2}}\lambda_\alpha}{\bar{a}}\right) \boldsymbol{{\cal I}} & \;\;\;\alpha = \beta \\ \ln\left( \frac{\lambda_\beta}{2q_{\alpha \beta}} \right)\boldsymbol{{\cal I}} + \hat{\boldsymbol{q}}_{\alpha \beta}\hat{\boldsymbol{q}}_{\alpha \beta} &\;\; \;\alpha \ne \beta \end{array} \right. , \end{equation*}$$

and where $\boldsymbol{{\cal I}}$ denotes the two-dimensional identity matrix.

Equation (25) reveals that the many-body defect dynamics is coupled via the collective (non-diagonal) mobility $\boldsymbol{{\cal M}}$. When $\lambda_\alpha \ne \lambda_\beta$, $\boldsymbol{{\cal M}}$ is moreover not symmetric, such that it introduces effective non-reciprocal couplings between the defect velocities. Although the Coulomb force is conservative, the center 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 equation (25) are set by the orientation field landscape, which is itself driven out of equilibrium by the motion of the defects (see equation (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–58].

Unlike the two-defect case, it is generally not possible to determine analytically the parameters $\lambda_\alpha(t)$ by calculating the exact solution for the angular field $\theta(\boldsymbol{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(\lambda/\bar{a})$, while the off-diagonal coefficients are ${\cal O}(1)$. At the leading order in $\ln(\lambda/\bar{a})$, one can thus replace λ by any macroscopic scale of the problem. Indeed, replacing $\lambda \to \lambda + \lambda^{\prime}$ does not change the result to the leading order, namely

$$\begin{equation*} \ln\left( \frac{\lambda + \lambda^{\prime}}{\bar{a}}\right) = \ln\left( \frac{\lambda}{\bar{a}} \right) + {\cal O}\left(1\right). \end{equation*}$$

For two defects, figure 2(b) shows that the approximation $\lambda \propto 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.

Although the terms of equation (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 equation (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 $\theta_\textrm{ssd}(r,\varphi)$—where $(r,\varphi)$ stands for the polar coordinates in the defect frame—is scale-free by construction, we have $\theta_\textrm{ssd}(r,\varphi) = \theta_\textrm{ssd}(\varphi)$. Hence, equation (10) becomes

$$\begin{equation} \frac{\mathrm{d} \bar{\zeta}_{ij}}{\mathrm{d} r_0} = \frac{1}{r_0}\int_0^{2\pi} \mathrm{d}\varphi \, \left(\partial_{\varphi}\theta_{\mathrm{ssd}}\left(\varphi\right)\right)^2 \, \hat{\varphi}_{i}\hat{\varphi}_{j} \equiv \frac{\bar{\nu}_{ij}}{r_0}, \end{equation} \tag{ 26 }$$

where $\hat{\boldsymbol{\varphi}}$ is the unit vector tangent to the unit circle. Equation (26) therefore implies that the friction matrix takes the form

$$\begin{equation} \bar{\zeta}_{ij}\left(\frac{r_0}{a}\right) = \bar{\nu}_{ij}\ln\left(\frac{r_0}{a}\right)+C_{ij}, \end{equation} \tag{ 27 }$$

where both $\bar{\boldsymbol{\nu}}$ and C are dimensionless and independent of a. We are now in a good position to expand equation (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)$, equation (27) simplifies to $\bar{\zeta}_{ij} = -\bar{\nu}_{ij}\ln(a) + {\cal O}(1)$, such that the friction tensor does no more involve the matching variable r₀ nor the constant C .

Moreover, replacing ${\zeta}_{ij} \sim \ln(a)$ in equation (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 equation (11) can be written as ${\boldsymbol{T}}_\textrm{bulk} = \boldsymbol{T}_\textrm{sb} + {\cal 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 $\delta {\cal F}_\textrm{bulk}/\delta\boldsymbol{\phi}_\textrm{sb} = \mathbf{0}$ (from (1)), the corresponding stress tensor is divergence-free: $\nabla\cdot\boldsymbol{T}_\textrm{sb} = \mathbf{0}$. Consequently, the line integral on the r.h.s. of equation (11) is independent of the chosen path, and we get

$$\begin{equation} \ln\left(\frac{\lambda}{a}\right)\nu_{ij}v_j = -\oint \mathrm{d} S_j \, T_{{\mathrm{sb}},ji}, \end{equation} \tag{ 28 }$$

where $\nu_{ij} = R_{ik}(t)R_{jl}(t) \bar{\nu}_{kl}$ includes the possible effects of defect anisotropy (see equation (9)) and the integral on the r.h.s. can be evaluated on any loop enclosing the defect core. The scale $\lambda \gg 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₀. Contrary to (11), equation (28) can thus be directly used without the need to eliminate r₀ by enforcing the matching condition.

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

$$\begin{equation*} \boldsymbol{\mathcal{Q}}\left(\rho,\theta\right) = \frac{\rho}{\sqrt{2}}\left(\begin{array}{cc} \cos2\theta & \sin2\theta \\ \sin2\theta & -\cos2\theta \\ \end{array}\right), \end{equation*}$$

such that topological defects carry a half-integer charge: $s = \pm \frac{1}{2}$. The bulk dynamics is described by the Frank–Oseen free energy [59], which takes the form ⁶

$$\begin{equation} {\cal F}_{\mathrm{bulk}} = \frac{1}{2} \int \mathrm{d}^2\boldsymbol{x} \, \left[ |\nabla \theta|^2 + \sqrt{2}\alpha \nabla \theta\cdot \boldsymbol{\mathcal{Q}}\left(1,\theta\right)\cdot \nabla \theta \right], \end{equation} \tag{ 29 }$$

where the parameter $\alpha \in [-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 $\alpha \ne 0$, the relaxational dynamics of θ follows a nonlinear equation, such that deriving the counterpart of equation (14) is generally difficult; even perturbatively in α. Here, we show how some progress can be made from equation (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]

$$\begin{equation} \partial_\varphi\theta_{\mathrm{ssd}}\left(\varphi\right) = 1 + \frac{1}{p} \sqrt{\frac{1-\alpha p^2\cos\left[2\psi\left(\varphi\right)\right] }{1 - \alpha\cos\left[2 \psi\left(\varphi\right)\right] } }, \end{equation} \tag{ 30 }$$

where $\psi(\varphi) \equiv \theta_\textrm{ssd}(\varphi)-\varphi$ and p is a constant fixed by the circuitation condition $\int_0^{2\pi}\mathrm{d}\varphi\, \partial_\varphi\theta_\textrm{ssd} = 2\pi s$. Equation (30) can in principle be solved iteratively at any order in α. At ${\cal O}(1)$, the solution is simply that of the isotropic theory, namely, $\theta_\textrm{ssd}(\varphi) = s\varphi + \theta_\textrm{c} + {\cal O}(\alpha)$, where θ_c is an integration constant. At linear order in α, the solution becomes

$$\begin{equation} \partial_\varphi\theta_{\mathrm{ssd}}\left(\varphi\right) = s + \frac{\alpha}{2}\frac{s\left(s-2\right)}{s-1}\cos\left(2\left(s-1\right)\varphi + 2\theta_{\mathrm{c}}\right)+{\cal O}\left(\alpha^2\right). \end{equation} \tag{ 31 }$$

Inserting this expression in equation (26) with $s = \pm \frac{1}{2}$ and integrating over ϕ, we find that $\nu_{ij} = \bar{\nu}_{ij} = \frac{\pi}{4} \delta_{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 section 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 $\theta_\textrm{sb}(\boldsymbol{x}) = \theta_0(\boldsymbol{x}) + \alpha \theta_1(\boldsymbol{x}) + {\cal O}(\alpha^2)$, where $\theta_0(\boldsymbol{x})$ solves the linear problem (α = 0) and $\theta_1(\boldsymbol{x})$ is calculated perturbatively. However, the formal expression of θ₁ still involves intricate integrals. We thus use the fact that the integration contour on the r.h.s. of equation (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 $\boldsymbol{q} = \frac{1}{2}(\boldsymbol{q}_+ - \boldsymbol{q}_-)$

$$\begin{align} \ln\left(\frac{\lambda}{a}\right) \dot{\boldsymbol{q}} = -\frac{1}{q}\left[\hat{\boldsymbol{q}} - \frac{4 \alpha}{3}\sin\left(2\theta_q\right)\boldsymbol{\epsilon}\hat{\boldsymbol{q}} \right], \end{align} \tag{ 32 }$$

where θ_q denotes the angle between q and the direction of the nematic order parameter at infinity parametrized by θ_c. Comparing equations (21) and (32), 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 figure 3 that varying θ_q amounts to changing the relative orientation of the two defects. The tangential force in equation (32) then reflects the energetic cost of misaligned defects due to elastic anisotropy, and favors configurations with $\theta_q = \frac{\pi}{2}$ or 0.

Zoom In Zoom Out Reset image size

To confirm the predictions of equation (32), we performed numerical simulations of the dynamics described by

$$\begin{eqnarray} \partial_t \mathcal{Q}_{ij} & = -\frac{\delta {\cal F}_{\mathrm{an}} }{\delta \mathcal{Q}_{ij}} + \frac{\delta_{ij}}{2} {\mathrm{Tr}}\left[\frac{\delta {\cal F}_{\mathrm{an}} }{\delta \boldsymbol{\mathcal{Q}}} \right],\\ {\cal F}_{\mathrm{an}} & = \int \mathrm{d}^2\boldsymbol{x} \, \left[ \frac{1}{2}\left(\delta_{kl} + \sqrt{2}\alpha \mathcal{Q}_{kl}\right)\left(\partial_k \mathcal{Q}_{ij}\right)\left(\partial_l \mathcal{Q}_{ij}\right) + \chi^2\left( 1- {\mathrm{Tr}}\left[\boldsymbol{\mathcal{Q}}^2\right] \right)^2 \right], \nonumber \end{eqnarray} \tag{ 33 }$$

for the full nematic order parameter $\boldsymbol{\mathcal{Q}}$, where $\chi \sim a^{-1}$ sets the size of the defect core. ${\cal F}_\textrm{an}$ is the simplest free energy that admits (29) as a bulk theory, which can be verified by setting ρ = 1 in equation (33). Figure 4 shows trajectories of the vector q obtained from initially prepared $\pm\frac{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 equation (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 equation (32) and simulations of the full model (33) without any fitting parameter. In qualitative agreement with the predictions of equation (32), we find that for $\alpha \ne 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 figure 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.

Zoom In Zoom Out Reset image size

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 ${\cal O}(\alpha^2)$. At this order, the orientation field obtained from (30) is given by

$$\begin{eqnarray} \qquad\quad \partial_\varphi\theta_{\mathrm{ssd}}\left(\varphi\right) = & s + \frac{\alpha}{2}\frac{s\left(s-2\right)}{s-1}\cos\left(2\left(s-1\right)\varphi + 2\theta_{\mathrm{c}}\right) \nonumber \\ & + \frac{\alpha^2\left(5 s^3-20 s^2+24 s-8 \right)}{16\left(s-1\right)^3}\cos\left(4\left(s-1\right)\varphi + 4\theta_{\mathrm{c}}\right) + {\cal O}\left(\alpha^3\right). \end{eqnarray} \tag{ 34 }$$

We use equation (26) to calculate the defect friction. Splitting it into isotropic and anisotropic parts, we obtain

$$\begin{equation} {\nu}_{ij} = \frac{1}{2}\int_0^{2\pi} \mathrm{d}\varphi \, \left(\partial_{\varphi}\theta_{\mathrm{ssd}}\left(\varphi\right)\right)^2 \left[ \delta_{ij} + \sqrt{2}\mathcal{Q}_{ij}\left(\varphi\right) \right]. \end{equation} \tag{ 35 }$$

For the isotropic part, we obtain after some algebra

$$\begin{equation*} \frac{1}{2}\int_0^{2\pi} \mathrm{d}\varphi \, \left(\partial_{\varphi}\theta_{\mathrm{ssd}}\left(\varphi\right)\right)^2 = s^2\pi\left(1+\alpha^2\frac{\left(s-2\right)^2}{8\left(s-1\right)^2}\right) + {\cal O}\left(\alpha^3\right). \end{equation*}$$

In addition, the anisotropic part of the friction is found to be

$$\begin{equation*} \frac{1}{\sqrt{2}}\int_0^{2\pi} \mathrm{d}\varphi \, \left(\partial_{\varphi}\theta_{\mathrm{ssd}}\left(\varphi\right)\right)^2 \mathcal{Q}_{ij}\left(\varphi\right) = \frac{3\pi \alpha^2}{8\sqrt{2}}\mathcal{Q}_{ij}\left(2\theta_{\mathrm{c}}\right)\delta_{s,\frac{1}{2}} + {\cal O}\left(\alpha^3\right). \end{equation*}$$

Hence, only $s = +\frac{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 $+\frac{1}{2}$ defect. Namely, it is straightforward to show from (34) and for $s = \frac{1}{2}$ that rotating θ_c by an angle ϑ amounts to rotate the defect solution by an angle 2ϑ. Hence, we define the $+\frac{1}{2}$ defect polarization as $\hat{\boldsymbol{p}} \equiv (\cos2\theta_\textrm{c} , \sin2\theta_\textrm{c})$ (see the blue arrows in figure 3), such that we finally obtain for the two types of nematic defects

$$\begin{eqnarray} &&\boldsymbol{\nu}_{\frac{1}{2}} = \frac{\pi}{4}\left(1+\frac{21}{8}\alpha^2\right) \hat{\boldsymbol{p}}\hat{\boldsymbol{p}} + \frac{\pi}{4}\left(1 - \frac{3}{8}\alpha^2\right)\left( \boldsymbol{I} - \hat{\boldsymbol{p}}\hat{\boldsymbol{p}}\right), \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} &&\boldsymbol{\nu}_{-\frac{1}{2}} = \frac{\pi}{4}\left(1+\frac{25}{72}\alpha^2\right) \boldsymbol{I}. \end{eqnarray} \tag{ 37 }$$

Comparing the coefficients of equation (36), we thus find that the $+\frac{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 $-\frac{1}{2}$-charged defect upward in all directions.

We now discuss the situation where the order parameter φ is two-dimensional and evolves in a three-dimensional 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 $\boldsymbol{\gamma}(l,t)$ the vector defining the defect line, and choose an adequate parametrization such that $|\boldsymbol{\gamma}^{\prime}(l,t)| = 1$. In three dimensions, the curve $\boldsymbol{\gamma}(l,t)$ is locally defined by the Frenet–Serret frame

$$\begin{equation} \hat{\boldsymbol{t}} = \boldsymbol{\gamma}^{\prime}, \qquad \hat{\boldsymbol{n}} = \kappa^{-1} \hat{\boldsymbol{t}}^{\prime}, \qquad \hat{\boldsymbol{b}} = \hat{\boldsymbol{t}} \times \hat{\boldsymbol{n}}, \end{equation} \tag{ 47 }$$

such that $\hat{\boldsymbol{t}}$ is tangent to γ , $\hat{\boldsymbol{n}}$ and $\hat{\boldsymbol{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 $\kappa a \ll 1$. Given the geometry of the problem, the natural choice for the volume ${\cal V}$ is an infinitesimal cylinder with top and bottom surfaces S₁ and S₂ orthogonal to the curve γ respectively in l and $l + \mathrm{d} l$, while the whole lateral surface $S_\textrm{L}$ lies in the matching region (see a sketch in figure 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

$$\begin{eqnarray} \zeta_{ij} & = \int_{{\cal V}}\mathrm{d}^3\boldsymbol{x} \, \left(\partial_i\boldsymbol{\phi}\right) \cdot \left(\partial_j\boldsymbol{\phi}\right) \nonumber \\ & \simeq \frac{1}{2} \mathrm{d} l \left(\delta_{ij} - \hat{t}_i \hat{t}_j\right)\int_{{\cal D}_{r_0}} \mathrm{d}^2\boldsymbol{x}_{\mathrm{n}}\, \left(\partial_k\boldsymbol{\phi}\right) \cdot \left(\partial_k\boldsymbol{\phi}\right) \nonumber\\ & \equiv \mathrm{d} l \left(\delta_{ij} - \hat{t}_i \hat{t}_j\right)\zeta\left(\frac{r_0}{a}\right), \end{eqnarray} \tag{ 48 }$$

at leading order in $\mathrm{d} l$, where ${\cal D}_{r_0}$ denotes the disk of radius r₀ orthogonal to $\hat{\boldsymbol{t}}$. The second equality is obtained by neglecting the part of $\nabla\boldsymbol{\phi}$ tangential to $\hat{\boldsymbol{t}}$ and assuming that the tensor coming from the integral on the r.h.s. is isotropic on the plane normal to $\boldsymbol{\gamma}(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 section 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₁ and S₂, respectively orthogonal to $\boldsymbol{\gamma}(l)$ and $\boldsymbol{\gamma}(l+\mathrm{d}l)$, gives the surface tension contribution:

$$\begin{equation} -\left[\int_{S_1} + \int_{S_2}\right] d S_j T_{ji} = \sigma\left(\frac{r_0}{a}\right) \hat{t}_i^{\prime} \mathrm{d} l = \sigma\left(\frac{r_0}{a}\right)\kappa \hat{n}_i\mathrm{d} l, \end{equation} \tag{ 49 }$$

where $\sigma = \int_{{\cal D}_{r_0}}\mathrm{d}^2\boldsymbol{x}_\textrm{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₀ and the core size a. In fact, we will show later that they are both logarithmically divergent in r₀, similarly to the situation studied in section 2.2 for point-defects.

Zoom In Zoom Out Reset image size

We now calculate the contribution of the lateral surface $S_\textrm{L}$. Contrary to domain walls, the bulk stress tensor does not vanish here because of the presence of the long-range modulations of θ which mediate interactions between the defect lines. Denoting $(\hat{\boldsymbol{r}},\hat{\boldsymbol{\varphi}},\hat{\boldsymbol{t}})$ as the local cylindrical frame centred in γ , the surface element on $S_\textrm{L}$ can be expressed as $\mathrm{d}{S}_i = - r_0 \mathrm{d}\varphi \mathrm{d} l \epsilon_{ijk}\hat{t}_j \hat{\varphi}_k$ where _ijk is the 3D totally antisymmetric Levi–Civita tensor. Hence,

$$\begin{equation} - \int_{S_{\mathrm{L}}} \mathrm{d} {S}_j \, {T}_{{\mathrm{bulk}},ji} = r_0\mathrm{d} l \epsilon_{klj} \hat{t}_k \int_{0}^{2\pi} \mathrm{d} \varphi \, \hat{\varphi}_l {T}_{{\mathrm{bulk}},ji}, \end{equation} \tag{ 50 }$$

where we have used the fact that the integration is done in the matching region to substitute T by its bulk counterpart. Combining equations (48)–(50), and denoting $\dot{\boldsymbol{\gamma}}_\perp \equiv (\boldsymbol{I} - \hat{\boldsymbol{t}}\hat{\boldsymbol{t}})\dot{\boldsymbol{\gamma}}$, we finally obtain

$$\begin{equation} \zeta\left(\frac{r_0}{a}\right)\dot{\gamma}_{\perp,i} = \sigma\left(\frac{r_0}{a}\right)\kappa\hat{n}_i - r_0 \epsilon_{klj} \hat{t}_k \int_{0}^{2\pi} \mathrm{d} \varphi \, \hat{\varphi}_l {T}_{{\mathrm{bulk}},ji}. \end{equation} \tag{ 51 }$$

Similarly to point-defects, equation (51) explicitly depends on the matching variable r₀. As r₀ is arbitrary, it must simplify when a bulk free energy is specified.

We now consider a bulk free energy analogous to that used in section 2.2:

$$\begin{equation} {\cal F}_{\mathrm{bulk}} = \frac{1}{2}\int \mathrm{d}^3\boldsymbol{x} \, |\nabla\theta|^2. \end{equation} \tag{ 52 }$$

The corresponding dynamical evolution of the orientation field $\theta(\boldsymbol{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₀, we obtain

$$\begin{eqnarray} &&\frac{\mathrm{d} \zeta}{\mathrm{d} r_0} = \frac{r_0}{2} \int_0^{2\pi}\mathrm{d}\varphi\, |\nabla\theta|^2,\qquad \frac{\mathrm{d} \sigma}{\mathrm{d} r_0} = r_0 \int_0^{2\pi}\mathrm{d}\varphi \, F_{\mathrm{bulk}} = \frac{\mathrm{d} \zeta}{\mathrm{d} r_0}, \end{eqnarray} \tag{ 53 }$$

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:

$$\begin{equation} \nabla\theta_{\mathrm{ssd}}\left(\boldsymbol{x}\right) = s\,\frac{\hat{\boldsymbol{\varphi}}}{r}, \end{equation} \tag{ 54 }$$

where r is the minimal distance between x and γ . Plugging (54) in equation (53) then straightforwardly leads to $\zeta \simeq \sigma \simeq \pi s^2 \ln(r_0/a)$. We now take a → 0 and perform a low mobility expansion as was done in section 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:

$$\begin{equation} \nabla\theta_{\mathrm{sd}}\left(\boldsymbol{x}\right) = \frac{s}{2} \int \mathrm{d} \boldsymbol{\gamma} \times \frac{\boldsymbol{x} - \boldsymbol{\gamma}}{|\boldsymbol{x} - \boldsymbol{\gamma}|^3}. \end{equation} \tag{ 55 }$$

In particular, in the specific case of a straight disclination line aligned along $\hat{\boldsymbol{z}}$, the expression (55) simplifies into (54). Since we integrate the stress tensor over a circle of radius $r_0 \to 0$ enclosing the line around γ , we write the phase field for $\boldsymbol{x} \to \boldsymbol{\gamma}$ as $\nabla\theta = \nabla\theta_\textrm{d} + \nabla\theta_\textrm{c}$. The discontinuous part is obtained by expanding (55) in the vicinity of γ , giving $\nabla\theta_\textrm{d}(\boldsymbol{x}) = \nabla\theta_\textrm{ssd}(\boldsymbol{x}) - \kappa s \ln(\kappa r_0) \hat{\boldsymbol{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

$$\begin{equation} \ln\left(\frac{L}{a}\right)\dot{\boldsymbol{\gamma}}_\perp = \kappa\ln\left(\frac{1}{\kappa a}\right)\hat{\boldsymbol{n}} - \frac{2}{s}\, \hat{\boldsymbol{t}} \times \nabla\theta_{\mathrm{c}}\left(\boldsymbol{\gamma}\right) . \end{equation} \tag{ 56 }$$

While the matching scale r₀ in the expression of the surface tension was simplified due to the curvature contribution to $\nabla\theta_\textrm{d}$, as previously discussed for slow defects the r₀ 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 ${\cal O}(1)$ implies that the line velocity is also ${\cal O}(1)$, while the interaction term on the r.h.s. of (56) is ${\cal O}(\ln^{-1}(a))$. Hence, the dynamics is dominated by the tension and, after taking $L \approx \kappa^{-1}$, we recover an equation similar to that ruling the motion of domain walls:

$$\begin{equation} \dot{\boldsymbol{\gamma}}_\perp = \kappa\hat{\boldsymbol{n}}, \qquad \left(\kappa = {\cal O}\left(1\right)\right). \end{equation} \tag{ 57 }$$

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 $\kappa = {\cal O}(\ln^{-1}(a))$, the interaction term in (56) becomes relevant. In the quasi-static approximation and assuming open boundaries, the term $\nabla\theta_\textrm{c}(\boldsymbol{\gamma}_\alpha)$ appearing in the equation of motion of the defect α can thus be generally expressed as $\nabla\theta_\textrm{c}(\boldsymbol{\gamma}_\alpha) = \sum_{\beta \ne \alpha} \nabla\theta_{\textrm{sd},\beta}(\boldsymbol{\gamma}_\alpha)$, where $\theta_{\textrm{sd},\beta}$ corresponds to the orientation profile generated by the defect β.

In the specific case where all disclination lines are straight and aligned along $\hat{\boldsymbol{z}}$, it is straightforward to check from (55) that (56) reduces to the description of the two-dimensional Coulomb gas [30]. Conversely, for straight disclination lines with arbitrary orientations, equation (56) takes the form

$$\begin{equation} \ln\left( \frac{L}{a} \right)\dot{\boldsymbol{\gamma}}_{\perp,\alpha} = \sum_{\beta \ne \alpha}\frac{s_\beta}{s_\alpha }\;\frac{ \hat{\boldsymbol{r}}_{\alpha \beta} \left(\hat{\boldsymbol{t}}_\alpha \cdot \hat{\boldsymbol{t}}_\beta\right) - \hat{\boldsymbol{t}}_\beta \left(\hat{\boldsymbol{t}}_\alpha \cdot \hat{\boldsymbol{r}}_{\alpha \beta}\right)}{r_{\alpha \beta}}, \end{equation} \tag{ 58 }$$

where $\boldsymbol{r}_{\alpha \beta} \equiv \frac{1}{2}(\boldsymbol{\gamma}_\alpha - \tilde{\boldsymbol{\gamma}}_{\alpha \beta})$ and $\tilde{\boldsymbol{\gamma}}_{\alpha \beta} \equiv \textrm{argmin}_{\boldsymbol{\gamma}_\beta}(|\boldsymbol{\gamma}_\alpha - \boldsymbol{\gamma}_\beta|)$. 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_\alpha s_\beta < 0$ or $s_\alpha s_\beta \gt 0$, respectively. Hence, due to the first term the effective interaction between disclination lines is always attractive regardless of their respective charges.

Another type of structure of interest are loop singularities. To study their dynamics, we consider a circular ring α centered at position $\boldsymbol{X}_\alpha$ with radius R_α . We define the associated dipole moment as

$$\begin{equation} \boldsymbol{m}_\alpha \equiv \frac{s_\alpha}{4}\oint \, \boldsymbol{\gamma}_\alpha \times \mathrm{d}\boldsymbol{\gamma}_\alpha = \frac{\pi s_\alpha R_\alpha^2}{2}\, \hat{\boldsymbol{b}}_\alpha, \end{equation} \tag{ 59 }$$

where $\hat{\boldsymbol{b}}_\alpha$ is the binormal unit vector orthogonal to the loop plane. Taking the time derivative of equation (59) and using (56), we find that the loop radius and orientation obey

$$\begin{eqnarray} \dot{R}_\alpha & = -\frac{1}{R_\alpha}-\frac{2}{s_\alpha \ln\left( R_\alpha/a \right)} \left(\hat{\boldsymbol{b}}_\alpha\cdot \nabla\theta_{\mathrm{c}}\right), \end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray} \dot{\hat{\boldsymbol{b}}}_\alpha & = \frac{2}{s_\alpha R_\alpha \ln\left( R_\alpha/a\right)} \hat{\boldsymbol{b}}_\alpha\times\left(\hat{\boldsymbol{b}}_\alpha\times\nabla\theta_{\mathrm{c}}\right), \end{eqnarray} \tag{ 61 }$$

while the dynamics of $\boldsymbol{X}_\alpha$ can be obtained by simply averaging (56) over the ring contour. Note that we have chosen to replace L with R_α in (60) and (61), as here $1/R_\alpha$ sets the natural scale for the line velocity. We first consider the case where the loop is subject to a uniform orientation gradient $\nabla\theta_\textrm{c} = \boldsymbol{k}$, as was done in section 2.2.2. It is then straightforward to check that $\dot{\boldsymbol{X}}_\alpha = \textbf{0}$. Moreover, we conclude from equation (61) that the ring anti-aligns (aligns) with k when $s_\alpha \gt 0$ ($s_\alpha < 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_\textrm{c}/\ln(R_\textrm{c}/a) = |s_\alpha|/(2|\boldsymbol{k}|)$, this expanding force overcomes the curvature-induced contraction and the ring expands indefinitely.

When $\nabla\theta_\textrm{c}$ results from other defect loops located far away from $\boldsymbol{X}_\alpha$, it is given at the leading order in the far-field approximation by a superposition of dipoles, namely,

$$\begin{equation} \nabla\theta_{\mathrm{c}} \simeq \sum_{\beta \ne \alpha} \frac{3 \left(\hat{\boldsymbol{r}}_{\alpha \beta}\cdot \boldsymbol{m}_\beta\right)\hat{\boldsymbol{r}}_{\alpha \beta}-\boldsymbol{m}_\beta}{r_{\alpha \beta}^3}, \end{equation} \tag{ 62 }$$

where $\boldsymbol{r}_{\alpha \beta} = \boldsymbol{X}_\alpha - \boldsymbol{X}_\beta$. In this case too, the loop is globally immobile but only contracts and rotates according to

$$\begin{eqnarray} \dot{R}_\alpha & = -\frac{1}{R_\alpha} - \sum_{\beta\ne \alpha}\frac{s_\beta}{s_\alpha}\frac{\pi R_\beta^2}{\ln\left(R_\alpha/a\right)}\frac{3 \left(\hat{\boldsymbol{r}}_{\alpha \beta}\cdot \hat{\boldsymbol{b}}_\alpha\right)\left(\hat{\boldsymbol{r}}_{\alpha \beta}\cdot \hat{\boldsymbol{b}}_\beta\right)-\hat{\boldsymbol{b}}_\alpha\cdot\hat{\boldsymbol{b}}_\beta}{r_{\alpha \beta}^3} , \end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray} \dot{\hat{\boldsymbol{b}}}_\alpha & = \sum_{\beta\ne \alpha}\frac{s_\beta}{s_\alpha} \frac{\pi R_\beta^2}{R_\alpha \ln\left(R_\alpha/a\right)} \frac{\hat{\boldsymbol{b}}_\alpha\times\left[\hat{\boldsymbol{b}}_\alpha\times\left(3 \left(\hat{\boldsymbol{r}}_{\alpha \beta}\cdot \hat{\boldsymbol{b}}_\beta\right)\hat{\boldsymbol{r}}_{\alpha \beta}-\hat{\boldsymbol{b}}_\beta\right)\right]}{r_{\alpha \beta}^3}. \end{eqnarray} \tag{ 64 }$$

In contrast to straight disclination lines (equation (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_\alpha R_\beta^2 \simeq r_{\alpha \beta}^3\ln(R_\alpha/a)$, which clearly breaks the far field assumption $r_{\alpha \beta} \gg R_\alpha, R_\beta$. Hence, in most cases distant loops will only rotate and self-annihilate over a finite time.