A unified field theory of topological defects and non-linear local excitations

Abstract

Topological defects and smooth excitations determine the properties of systems showing collective order. We introduce a generic non-singular field theory that comprehensively describes defects and excitations in systems with O(n) broken rotational symmetry. Within this formalism, we explore fast events, such as defect nucleation/annihilation and dynamical phase transitions where the interplay between topological defects and non-linear excitations is particularly important. To highlight its versatility, we apply this formalism in the context of Bose-Einstein condensates, active nematics, and crystal lattices.

Introduction

Topological defects are hallmarks of systems exhibiting collective order. They are widely encountered from condensed matter, including biological systems, to elementary particles, and the very early Universe^{1,2,3,4,5,6,7,8}. The small-scale dynamics of interacting topological defects are crucial for the emergence of large-scale non-equilibrium phenomena, such as quantum turbulence in superfluids⁹, spontaneous flows in active matter¹⁰, or dislocation plasticity in crystals¹¹. In fact, classical discrete modeling approaches such as point vortex models¹² and discrete dislocation dynamics¹³ describe turbulence and plasticity in terms of the collective dynamics of topological defects as interacting charged points (in 2D) or line defects (in 3D). In most of these theories, the interactions of topological defects are modeled through the linear excitations that they induce in the far fields. The physics of events on short time- and length scales, such as core energies, nucleation conditions, defect interaction, etc., are often introduced by ad-hoc rules, such as cut-off parameters, Schmidt stress nucleation criteria, and defect line recombination rules. However, the dynamics of these events play a vital role in the transitions between different dynamical regimes. This is the case, for example, in stirred Bose-Einstein condensates where different superfluid flow regimes are observed depending on the size and speed of the moving obstacle^{14,15,16,17,18,19}, and where there is a subtle interplay between vortices and shock waves. Active nematic fluids are characterized by a dynamic transition to active turbulence at a sufficiently large activity where the spontaneous flows are sustained by the creation and annihilation of orientational defects^20,21. During plastic deformation of polycrystals, grains are progressively fragmented, a process governed by the nucleation and patterning of dislocations²². A number of macroscopic criteria exist for the nucleation of topological defects in crystals^23,24,25. Due to the highly non-linear nature of this process, however, it still remains poorly understood.

In this paper, we present a formalism to describe the evolution of ordered systems from the dynamics of their topological defects and their interactions with smooth but localized excitations. The versatility of the approach allows us to gain insight into defect annihilation, the onset of collective behavior, and perspectives on defect structures. In particular, we apply the method to systems of increasing topological and dynamical complexity. First, we study the motion of isolated vortices in Bose-Einstein condensates, which, in addition to confirming that the method correctly identifies topological defects and their velocities, sheds light on changes in quantum pressure arising from the interplay between phase slips and shock waves. For active nematics, we observe that the onset of active turbulence as a melting of periodic arches is signaled by the formation of bound dipoles of nematic defects at the core of dislocations in the nematic arches. Similarly, bound dipoles of phase slips are also associated with the nucleation of dislocations in a crystal lattice.

The proposed approach builds upon the classical method introduced by Halperin and Mazenko (hereafter called the HM-method)^26,27 to track and derive analytical results for topological defects. Therefore, in the section “Classical description of topological defects”, we begin with preliminary details of homotopy theory for topological defects and how the HM-method can be used for O(n)-symmetric theories to track their location and kinematics. In the section “Non-singular defect fields”, we then develop a non-singular field theory as a generalization of the HM-method which constitutes our primary reduced defect field. The method is then applied to the aforementioned physical systems in the sections “Defect annihilation: vortices in Bose-Einstein condensates”, “Onset of collective behavior: active nematics”, and “Defect structures: solid crystals”. For the sake of readability, a rigorous derivation of the theoretical framework for arbitrary dimensions and details of the numerical simulations are reported in the Supplementary Notes. Conclusions and perspectives for further study are outlined in the section “Discussion”.

Classical description of topological defects

Collective order is typically described by an order parameter field representative of symmetries and carrying information about topological defects and smooth, localized excitations. Although the order parameters are well-established for conventional systems, one often needs to define them for more exotic systems^28,29. In this paper, we focus on well-known order parameters for systems with broken O(n) rotational symmetries, where n is the intrinsic dimension of the order parameter.

Homotopy theory provides a valuable identification and classification of topological defects¹. The fundamental idea of homotopy theory is that the order parameter can be mapped onto a particular topological space R and the homotopy group of R classify topological defects. For example, in the XY-model of ferromagnetism, and more generally for any system with O(2) broken symmetry, the order parameter is mapped by a 2D unit vector u onto $R={{{{\mathcal{S}}}}}^{1}$, the unit circle. On ${{{{\mathcal{S}}}}}^{1}$, we may define classes of closed circuits (loops), where loops of the same class are homotopic, i.e., they can be continuously deformed into each other. These classes, together with an appropriate binary operation, define the homotopy group of ${{{{\mathcal{S}}}}}^{1}$. This group is isomorphic to ${\mathbb{Z}}$ under addition since the difference between two loops that are not homotopic is how many times they have looped around the circle ${{{{\mathcal{S}}}}}^{1}$. Therefore, in regions of space where u is continuous and well-defined, a closed circuit $\partial {{{\mathcal{M}}}}$ in real space corresponds to a closed circuit in ${{{{\mathcal{S}}}}}^{1}$, and the topological charge s_top contained in $\partial {{{\mathcal{M}}}}$ is given as an integer by the isomorphism between homotopy group of ${{{{\mathcal{S}}}}}^{1}$ and ${\mathbb{Z}}$. This topological charge is obtained from ${{{\boldsymbol{u}}}}=(\cos \theta ,\sin \theta )$, by the contour integral

$${s}_{{{{\rm{top}}}}}=\frac{1}{2\pi }{\oint }_{\partial {{{\mathcal{M}}}}}d\theta ,$$

(1)

which is invariant under any smooth deformations of $\partial {{{\mathcal{M}}}}$. This also implies that by shrinking $\partial {{{\mathcal{M}}}}$ down to a point and given that s_top is a constant, there must be regions inside $\partial {{{\mathcal{M}}}}$ where u is undefined. These are the topological defects that have a s_top charge. Therefore, topological defects for $R={{{{\mathcal{S}}}}}^{1}$ in 2D are points with their charge determined by corresponding loop integration. On the other hand, such topological defects (with $R={{{{\mathcal{S}}}}}^{1}$) in three dimensions are lines.

In field theories of symmetry-breaking transitions, the ground state of the order parameter minimizes a free energy constructed from symmetry considerations. For broken rotational symmetries, the order parameter is a vector field Ψ, which in the ordered (ground) state has a constant magnitude ∣Ψ∣ = Ψ₀, meaning that the ground state manifold is ${{{{\mathcal{S}}}}}^{n-1}$, where n is the number of components of Ψ. The link between the order parameter Ψ, and 2D unit vector (director) field ${{{\boldsymbol{u}}}}\in {{{{\mathcal{S}}}}}^{1}$ is given by u = Ψ/∣Ψ∣ and topological defects are located at positions where u is undefined, which corresponds to ∣Ψ∣ = 0 as shown in Fig. 1a, b.

Fig. 1: Different types of excitations in a 2D vector field theory.

A + 1 defect is shown in a the order parameter field Ψ and b in the unit vector field u = Ψ/∣Ψ∣. Excitations of c the ground state can be categorized into d linear excitations with variations in the orientation of Ψ, (e) local non-linear excitations for which also the magnitude ∣Ψ∣ varies and f topological defects.

A description of topological defects as zeros of order parameters in O(n) models and their kinematics was proposed originally by Halperin and Mazenko in the context of phase-ordering kinetics^26,27 and extended to systems driven out of equilibrium, such as in stirred Bose-Einstein condensation^18,30,31, active nematics^32,33, and deformed crystals^34,35,36. Sticking to O(2)-symmetry in two dimensions and using the definition of a topological charge given in Eq. (1), it is possible to express the topological defect density in terms of the zeros of the order parameter Ψ²⁶ tracked by Dirac-delta functions as

$${\rho }_{{{{\rm{top}}}}}({{{\boldsymbol{r}}}})\equiv \mathop{\sum}\limits_{\alpha }{q}_{\alpha }{\delta }^{(2)}({{{\boldsymbol{r}}}}-{{{{\boldsymbol{r}}}}}_{\alpha })=D({{{\boldsymbol{r}}}}){\delta }^{(2)}({{{\boldsymbol{\Psi }}}}),$$

(2)

where q_α and r_α are, respectively, the charge and position of the topological defect α, δ⁽²⁾(Ψ) = δ(Ψ₁)δ(Ψ₂), and D(r) is the (signed) Jacobian determinant of the map Ψ,

$$\begin{array}{rcl}D=\frac{\partial ({\Psi }_{1},{\Psi }_{2})}{\partial (x,y)}\,=\,{\partial }_{x}{\Psi }_{1}{\partial }_{y}{\Psi }_{2}-{\partial }_{x}{\Psi }_{2}{\partial }_{y}{\Psi }_{1}\\ \,=\,\frac{1}{2}{\epsilon }^{ij}{\tilde{\epsilon }}^{mn}({\partial }_{i}{\Psi }_{m})({\partial }_{j}{\Psi }_{n}),\end{array}$$

(3)

where ϵ^ij are the components of the Levi-Civita tensor in real space. The Levi-Civita tensor $\tilde{\epsilon }$ in order parameter space is written with a tilde to emphasize that it is contracted with the order parameter Ψ. In the Cartesian space, both ϵ and $\tilde{\epsilon }$ are simply the Levi-Civita (permutation) symbols. Note that Eq. (2) is the usual scaling property of the delta function taking Ψ as input, apart from the sign of D carrying information of the charge q_α of the topological defects. This result was shown in ref. ²⁶ by considering as explicit ansatz a negative point defect, but can, in general, be justified using differential forms. Nominally, the D field in Eq. (3) is evaluated at the location of the defect only, because of the δ-function in Eq. (2).

Results

Non-singular defect fields

The δ-function in the topological charge density of Eq. (2) locates the topological defects at singular points where u is undefined. In O(2) models, however, even though the ground state manifold is ${{{{\mathcal{S}}}}}^{1}$, the topological excitations have a finite core over which the magnitude of the order parameter goes smoothly to zero. This feature is also seen in physical systems, for instance, in liquid crystals, where optical retardance is an order parameter that goes to zero at the core. This has been used to quantify the size and structure of the defect cores in liquid crystals³⁷. Motivated by this, we seek to generalize Eq. (2) in a way that will avoid singularities in the resulting charge density.

Since the equilibrium value Ψ₀ of ∣Ψ∣ is constant, the order parameter effectively resides in ${{{{\mathcal{D}}}}}^{2}$, the unit disk. We propose in this paper that the simplest generalization of s_top is to consider the relative area of ${{{{\mathcal{D}}}}}^{2}$ swept by Ψ on the circuit $\partial {{{\mathcal{M}}}}$. During an infinitesimal displacement along $\partial {{{\mathcal{M}}}}$, Ψ sweeps the infinitesimal area given by half of the parallelogram spanned by Ψ and dΨ. This (signed) area is given by $\frac{1}{2}{\tilde{\epsilon }}^{mn}{\Psi }_{m}d{\Psi }_{n}$, see Fig. 2. The complete area of ${{{{\mathcal{D}}}}}^{2}$ is $\pi {\Psi }_{0}^{2}$, and we define the charge s as the area swept by Ψ relative to the area of ${{{{\mathcal{D}}}}}^{2}$,

$$s=\frac{1}{\pi {\Psi }_{0}^{2}}{\oint }_{\partial {{{\mathcal{M}}}}}\frac{1}{2}{\tilde{\epsilon }}^{mn}{\Psi }_{m}d{\Psi }_{n},$$

(4)

where $\partial {{{\mathcal{M}}}}$ is defined as in Eq. (1). Naming s a “charge” suggests that it satisfies a global conservation law, which we shall prove shortly. The connection between s and s_top is made by recognizing that for a path $\partial {{{\mathcal{M}}}}$ in the far field of a topological defect, where ∣Ψ∣ = Ψ₀, s = s_top. To see this, note that if ∣Ψ∣ = Ψ₀, then the infinitesimal area swept by Ψ is simply $\frac{1}{2}{\Psi }_{0}^{2}d\theta$, which inserted into Eq. (4) gives Eq. (1). Closer to the core, however, the magnitude ∣Ψ∣ decreases and s is no longer an integer, which is why the associated defect density will give information about the core extent. Using Green’s theorem, we get

$$s=\frac{1}{2\pi {\Psi }_{0}^{2}}{\oint }_{\partial {{{\mathcal{M}}}}}{\tilde{\epsilon }}^{mn}{\Psi }_{m}{\partial }_{k}{\Psi }_{n}d{l}^{k}={\int}_{{{{\mathcal{M}}}}}{d}^{2}r\rho ({{{\boldsymbol{r}}}}),$$

(5)

where ρ(r) is the charge density of s, given by

$$\rho ({{{\boldsymbol{r}}}})=\frac{D({{{\boldsymbol{r}}}})}{\pi {\Psi }_{0}^{2}}.$$

(6)

Whereas ρ_top describes topological defects as point singularities in the physical space, ρ describes topological defects with a finite core size.

Fig. 2: A continuous field Ψ(r) containing defects with integer charges +1, −1, and +2.

The net integer topological charge contained in the circuits is given by the winding number of the unit vector field u in ${{{{\mathcal{S}}}}}^{1}$. The (signed) relative area gives the value of s for the circuits spanned by the order parameter Ψ in ${{{{\mathcal{D}}}}}^{2}$.

The time derivative of Eq. (6) gives a continuity equation

$${\partial }_{t}\rho +\nabla \cdot {{{\boldsymbol{J}}}}=0,$$

(7)

with the current density determined by the evolution of the order parameter

$${J}^{i}=-\frac{1}{\pi {\Psi }_{0}^{2}}{\epsilon }^{ij}{\tilde{\epsilon }}^{mn}({\partial }_{t}{\Psi }_{m})({\partial }_{j}{\Psi }_{n}).$$

(8)

Thus, ρ is a globally conserved quantity, and the change in s contained in a circuit $\partial {{{\mathcal{M}}}}$ is given by

$${\partial }_{t}s={\partial }_{t}{\int}_{{{{\mathcal{M}}}}}{d}^{2}r\rho ({{{\boldsymbol{r}}}})={\int}_{\partial {{{\mathcal{M}}}}}{{{\boldsymbol{J}}}}\cdot d{{{\boldsymbol{n}}}},$$

(9)

where dn is an infinitesimal surface area normal to the circuit $\partial {{{\mathcal{M}}}}$. Far away from defects, ∣Ψ∣ = Ψ₀ and the time evolution of Ψ is carried by its phase θ(r, t) through ${{{\boldsymbol{\Psi }}}}={\Psi }_{0}(\cos \theta ,\sin \theta )$ which can be inserted in Eq. (8) to show that J = 0. This means that linear perturbations of the ground state, which affect the orientation of Ψ only, are not described by the charge density ρ. However, it describes a certain type of local non-linear perturbations, where the magnitude is affected; see Fig. 1c–f. We will exemplify this distinction in the applications. Due to the standard continuity form of Eq. (7), we can connect it to a velocity field v through the charge flux ρv. Equation (7) only determines the current ρv up to an unknown divergence free contribution K, i.e., ${{{\boldsymbol{v}}}}=\frac{1}{\rho }({{{\boldsymbol{J}}}}+{{{\boldsymbol{K}}}})$, where ∇ ⋅ K = 0. However, when ρ ≠ 0, there exists a unique velocity field v^(Ψ) such that the evolution of Ψ can be written in a generic advection form ∂_tΨ + (v^(Ψ) ⋅ ∇ )Ψ = 0, equivalently expressed as

$$\left(\begin{array}{c}{\partial }_{t}{\Psi }_{1}\\ {\partial }_{t}{\Psi }_{2}\end{array}\right)+\left(\begin{array}{cc}{\partial }_{1}{\Psi }_{1}&{\partial }_{2}{\Psi }_{1}\\ {\partial }_{1}{\Psi }_{2}&{\partial }_{2}{\Psi }_{2}\\ \end{array}\right)\left(\begin{array}{c}{v}_{1}^{({{{\boldsymbol{\Psi }}}})}\\ {v}_{2}^{({{{\boldsymbol{\Psi }}}})}\end{array}\right)=0.$$

(10)

This equation can be inverted to uniquely determine v^(Ψ) if $\det ({\partial }_{i}{\Psi }_{n})=D({{{\boldsymbol{r}}}})\,\ne \,0$. To find v^(Ψ) where this condition holds true, i.e., the regions of interest where also ρ(r) ≠ 0 from Eq. (6), it is then possible to invert Eq. (10). However, it is easier to insert ∂_tΨ = − (v^(Ψ) ⋅ ∇ )Ψ into the expression J/ρ and see that it is the solution of Eq. (10). Thus, to fix the gauge on v, we set K = 0 to get v = v^(Ψ) and find

$${v}^{i}=\frac{{J}^{i}}{\rho }=-2\frac{{\epsilon }^{ij}{\tilde{\epsilon }}^{mn}({\partial }_{t}{\Psi }_{m})({\partial }_{j}{\Psi }_{n})}{{\epsilon }^{ij}{\tilde{\epsilon }}^{mn}({\partial }_{i}{\Psi }_{m})({\partial }_{j}{\Psi }_{n})},$$

(11)

where it is implied that repeated indices are summed over independently in the numerator and denominator. It should be noted that the velocity v only describes the velocity of the defect density ρ and is not, in general, the same as the advection velocity of the order parameter. We have only shown that if ρ ≠ 0 in some region then it is possible to write the evolution of Ψ in this way. If the actual evolution of Ψ is given as the advection v_D of a density field (i.e., including the term Ψ ∇ ⋅ v_D), then v ≠ v_D, because the compressible part of the advection will not directly translate into the motion of topological defects. However, if a localized topological defect moves without changing its core structure, i.e., with a frozen core, Eq. (11) will give this velocity in the region of the core, which we will show in the section “Defect annihilation: vortices in Bose-Einstein condensates”. While the expression for the current of D and the velocity equation (11) have previously been used in the HM-method, several important distinctions can be highlighted. Firstly, the derivation of the ρ field from the redefined charge, Eq. (6), shows that the field carries topological information and does not only serve as auxiliary transformation determinants of δ-functions. Secondly, the velocity field has previously only been rigorously shown to apply to topological defects. In contrast, this derivation also describes the velocity of ρ for other non-linear excitations. Thirdly, the fixing of the gauge K has not been adequately addressed in previous works to the authors’ knowledge. While the derivation above was done for a n = 2 order parameter in d = 2 spatial dimensions for simplicity, topological defects exist whenever d ≥ n. Equation (4) can be generalized to arbitrary dimensions by replacing the integrand with the volume of the n-sphere spanned by Ψ = (Ψ₁, …, Ψ_n) and normalizing by the volume ${V}_{n}{\Psi }_{0}^{n}$ of the n-sphere. We show in the Supplementary Notes the formal derivation, and here we state the result that the charge density becomes

$$\begin{array}{l}n=1:\quad {\rho }_{i}=\frac{{\partial }_{i}\Psi }{2{\Psi }_{0}}\\ n\ge 2:\quad {\rho }_{{i}_{1}\ldots {i}_{d-n}}=\frac{{D}_{{i}_{1},\ldots {i}_{d-n}}}{{V}_{n}{\Psi }_{0}^{n}}\end{array}$$

(12)

with

$${D}_{{i}_{1}\ldots {i}_{d-n}}=\frac{1}{n!}{{\epsilon }^{{\mu }_{1}\ldots {\mu }_{n}}}_{{i}_{1}\ldots {i}_{d-n}}{\tilde{\epsilon }}^{{\nu }_{1}\ldots {\nu }_{n}}({\partial }_{{\mu }_{1}}{\Psi }_{{\nu }_{1}})\ldots ({\partial }_{{\mu }_{n}}{\Psi }_{{\nu }_{n}}).$$

(13)

Generalizing the derivation of the defect kinematics, we find general expressions for the reduced defect velocity field

$${v}^{{\mu }_{1}}=-n\frac{{\delta }_{{\nu }_{1}^{{\prime} }}^{\left[{\nu }_{1}\right.}{\delta }_{{\nu }_{2}^{{\prime} }}^{{\nu }_{2}}\ldots {\delta }_{{\nu }_{n}^{{\prime} }}^{\left.{\nu }_{n}\right]}({\partial }_{t}{\Psi }_{{\nu }_{1}})({\partial }^{{\mu }_{1}}{\Psi }^{{\nu }_{1}^{{\prime} }})\mathop{\prod }\nolimits_{l = 2}^{n}({\partial }_{{\mu }_{l}}{\Psi }_{{\nu }_{l}})({\partial }^{{\mu }_{l}}{\Psi }^{{\nu }_{l}^{{\prime} }})}{{\delta }_{{\nu }_{1}^{{\prime} }}^{\left[{\nu }_{1}\right.}{\delta }_{{\nu }_{2}^{{\prime} }}^{{\nu }_{2}}\ldots {\delta }_{{\nu }_{n}^{{\prime} }}^{\left.{\nu }_{n}\right]}\mathop{\prod }\nolimits_{l = 1}^{n}({\partial }_{{\mu }_{l}}{\Psi }_{{\nu }_{l}})({\partial }^{{\mu }_{l}}{\Psi }^{{\nu }_{l}^{{\prime} }})},$$

(14)

$$\,{{\mbox{Special case}}}\,n=d:\quad {v}^{{\mu }_{1}}=-n\frac{{\epsilon }^{{\mu }_{1}{\mu }_{2}\ldots {\mu }_{n}}{\tilde{\epsilon }}^{{\nu }_{1}\ldots {\nu }_{n}}({\partial }_{t}{\Psi }_{{\nu }_{1}})\mathop{\prod }\nolimits_{l = 2}^{n}({\partial }_{{\mu }_{l}}{\Psi }_{{\nu }_{l}})}{{\epsilon }^{{\mu }_{1}\ldots {\mu }_{n}}{\tilde{\epsilon }}^{{\nu }_{1}\ldots {\nu }_{n}}\mathop{\prod }\nolimits_{l = 1}^{n}{\partial }_{{\mu }_{l}}{\Psi }_{{\nu }_{l}}},$$

(15)

where [ν₁ν₂…ν_n] is the antisymmetrization over the indices ν₁ν₂…ν_n. Equation (15) is the special case of n = d, where the velocity can be written in a simpler way. Still, Eq. (14) looks complicated due to the arbitrary number of dimensions and so we have summarized the most important cases of n ≤ d ≤ 3 in Supplementary Figure 2 of the Supplementary Notes. Thus, Eqs. (12) and (14) are the primary general expressions of the reduced defect field. The equations generalize the description of topological defects in the HM-method to include both topological defects and non-linear excitations.

There are two important notes to be made on the generalization beyond the case d = n = 2. Firstly, for n ≥ 2, the charge density is a rank (d − n) tensor that represents the defect density per n-dimensional volume-oriented normal to the manifold, e.g., how the charge density on a 2D surface is expressed in terms of the normal vector to the surface. The case of n = 1 is special because densities on one-dimensional manifolds are usually expressed in terms of the density along the manifold, i.e., the charge density per length along the curve. Secondly, in the case of n < d, the gauge K cannot be uniquely determined by looking at the evolution of Ψ alone. Therefore, another condition is required to obtain Eq. (14). This condition implies that topological defects live effectively on a d − n dimensional submanifold and will move perpendicular to this structure, e.g., how the motion of a line defect is given by the velocity normal to its tangent vector. Due to the difference in definitions of the integrals to yield the topological content, this translates to the velocity being parallel to the charge density for n = 1 and perpendicular to it for n ≥ 2. This velocity will be normal to topological structures in the case of topological lines or walls. While the systems of study in this manuscript exhibit ground state manifolds with ${{{{\mathcal{S}}}}}^{1}$ symmetries (n = 2), the generalization can be directly applied to systems with n = 1, where the defect density represents domain walls in interfacial systems such as viscous fingering³⁸, or with n = 3, such as the 3D Heisenberg model of ferromagnetism, where the defect density will show emergent magnetic monopoles³⁹. For further discussions, see the Supplementary Notes.

With the method at hand, we study phenomena involving both topological charges and non-linear local excitations through the reduced defect field and the information it conveys, such as the velocity of topological defects. This is done by considering progressively such phenomena in three representative systems with broken O(2) symmetry and featuring increasing complexity in terms of order parameters and collective behaviors. Both system-specific information and general behaviors will be outlined. As a starting point, we consider a Bose-Einstein condensate where the order parameter is isomorphic to ${{{\boldsymbol{\Psi }}}}\in {{{{\mathcal{D}}}}}^{2}$ so that the method can be directly applied.

Defect annihilation: vortices in Bose-Einstein condensates

Within the Gross Pitaevskii theory of a superfluid Bose-Einstein condensate (BEC), the condensed bosons are described by a macroscopic wavefunction ψ, and its evolution can be described by damped Gross Pitaevskii equation^18,40

$${\mathbb{i}}\hslash {\partial }_{t}\psi =(1-{\mathbb{i}}\gamma )\left(-\frac{{\hslash }^{2}}{2m}{\nabla }^{2}+g| \psi {| }^{2}-\mu \right)\psi ,$$

(16)

where g is an effective scattering parameter between condensate atoms, γ > 0 is an effective thermal damping coefficient and μ is the chemical potential. The complex condensate wavefunction ψ is isomorphic to a real 2D vector order parameter Ψ = (Ψ₁, Ψ₂) through $\psi ={\Psi }_{1}+{\mathbb{i}}{\Psi }_{2}$, the norm of which is given by the absolute value ∣ψ∣. In the equilibrium ground state, the phase of ψ (and therefore the direction of Ψ) is constant, and the magnitude is given by $| \psi | ={\Psi }_{0}=\sqrt{\mu /g}$. Topological defects in the orientational (unit vector) field correspond to quantized vortices captured by the charge density field

$${\rho }^{(\psi )}({{{\boldsymbol{r}}}})=\frac{gD({{{\boldsymbol{r}}}})}{\pi \mu }.$$

(17)

In this context, the D field (calculated from Ψ) has the physical interpretation of the generalized superfluid vorticity³¹. Linear perturbations of the ground state are phonons, which are characterized by traveling waves in the phase of the order parameter ψ, and will not be signaled in the defect density field ρ. Non-linear local perturbations, e.g., brought on by external stirring potentials or obstacles, will lead to a decrease in the magnitude of the order parameter near the obstacle^14,16,17,41, leading to an increase in the quantum pressure, defined as

$$P=-\frac{{\hslash }^{2}}{2m}\frac{{\nabla }^{2}| \psi | }{| \psi | }.$$

(18)

Such excitations are detected by ρ^(ψ), and mediate the nucleation or annihilation of topological defects. To showcase this, we simulate a BEC as dictated by Eq. (16) with an initial condition featuring two vortices at (x, y) = (±5, 0). Numerical details are reported in the “Methods” section. Dimensionless units are defined so that ℏ = m = g = μ = 1 and the damping coefficient is set to γ = 0.1. Figure 3 illustrates the defect density from Eq. (17) during the fast event of annihilating a vortex with an anti-vortex due to a small thermal drag.

Fig. 3: Annihilation of a vortex dipole in a Bose-Einstein condensate.

Snapshots of a defect density, b condensate phase arg(ψ), and c quantum pressure at different times from bottom to top: at t = 5, t = 60 (before annihilation), t = 105 (after) and t = 110. a Defect velocity is included prior before annihilation. Notice in (b) the large phase gradients after the annihilation due to the induced shock-waves which can also be seen in the (c) quantum pressure profiles. The plots in column (c) have saturated colorbars because of the singular pressure at the defect core.

The velocity field from Eq. (14) is plotted close to vortices and shows two exciting features. At the beginning of the simulations (t = 5), the non-uniform velocity over the vortex core indicates the early core deformation induced by the initial conditions. After this relaxation, however, vortices retain stationary or rigid cores and consistently feature a uniform velocity. After the annihilation event, we can see traces of their diffusive cores in the excitations produced by the vortex annihilation, as seen by the quantum pressure in the system, which is shown in Fig. 3c. We will see in the following that similar traces appear as precursory patterns for defect nucleation. Moreover, after having dealt with a system with only one broken symmetry, we now consider systems that have multiple rotational or translational symmetries.

Onset of collective behavior: active nematics

In this section, we consider the case of an active nematic system. This system is peculiar as we can construct the defect density from different order parameters. By applying the proposed formalism we can investigate the transition among different regimes and the interplay among defects. Interestingly, we will show that defects in one broken symmetry are the nucleation sites of defects for a separate order.

Within the hydrodynamic approach⁴², the nematic orientational order of active matter in two dimensions is described by a rank-2 symmetric and traceless tensor Q determined by the nematic director ${{{\bf{n}}}}=(\cos (\theta ),\sin (\theta ))$

$$Q=S\left(\begin{array}{cc}{n}_{1}{n}_{1}-\frac{1}{2}&{n}_{1}{n}_{2}\\ {n}_{2}{n}_{1}&{n}_{2}{n}_{2}-\frac{1}{2}\end{array}\right)\equiv \left(\begin{array}{cc}{\Psi }_{1}&{\Psi }_{2}\\ {\Psi }_{2}&-{\Psi }_{1}\end{array}\right),$$

(19)

where S is an order parameter which is 0 in the disordered phase. Q is thus related to the ${{{{\mathcal{D}}}}}^{2}$ order parameter Ψ field by ${{{\boldsymbol{\Psi }}}}=\frac{S}{2}(\cos (2\theta ),\sin (2\theta ))$. The evolution of the Q-tensor follows dissipative dynamics coupled with an incompressible Stokes flow with substrate friction⁴³. Details on the evolution equation and its numerical method are reported in the “Methods” section. The system is here initialized in a homogeneous nematic phase with small perturbations in the angle of the director field. These perturbations are enhanced by the active stress creating a striped phase that is further destabilized and eventually melts due to the creation of topological defects leading into active turbulence. The ground state corresponds to a constant magnitude $| {{{\boldsymbol{\Psi }}}}| \equiv {\Psi }_{0}=\sqrt{B}/2$ dependent on the parameter B, which is defined in the “Methods” section. Within the framework introduced in the section “Non-singular defect fields”, this gives the following expression for the defect density

$${\rho }^{(Q)}=\frac{4D({{{\boldsymbol{r}}}})}{\pi B},$$

(20)

which supports orientational defects with half-integer charge s_top = ± 1/2. In Fig. 4a, we show the nematic orientation θ in the colorbar to emphasize the breaking of translational symmetry and the formation of a (transient) striped order. The striped order arises from modulations in the nematic orientation which, to first order, do not change the magnitude of the order parameter Ψ. Thus, these are linear perturbations not signaled by ρ^(Q).

Fig. 4: Onset of active turbulence in a nematic liquid crystal mediated by the nucleation of topological defects.

a The angle of the nematic director at t = 240, prior to nucleation of half-integer defects from the unstable periodic arches, and d at t = 260, after nucleation. b The defect density ${\rho }^{({\eta }_{{{{\boldsymbol{k}}}}})}$ at t = 240, corresponding to the broken translational symmetry, shows the charge signature of dislocations with large core structures. The dislocation core harbors a bound dipole (inset) shown in (c) the defect density ρ^(Q) associated to the nematic order at t = 240, which splits into fully formed $\pm \frac{1}{2}$ defects after nucleation as shown in (f) ρ^(Q) at t = 260. Panel (e) shows the speed v = 〈∣v∣〉 of the two localized blobs in the charged defect density ρ^(Q) around the nucleation site. After the nucleation event indicated by the dashed line, these correspond to the speed of the ±1/2 defects.

The inset of Fig. 4a shows a dislocation in the periodic arches in the nematic director. To describe these defects, we represent the parameter Ψ as a complex field ψ = ∣Ψ∣e^iθ and decompose it into a slowly-varying amplitude field of the periodic arch mode as

$$\psi ({{{\boldsymbol{r}}}})={\psi }_{0}({{{\boldsymbol{r}}}})+{\eta }_{{{{\boldsymbol{k}}}}}({{{\boldsymbol{r}}}}){e}^{i{{{\boldsymbol{k}}}}\cdot {{{\boldsymbol{r}}}}}+{\eta }_{-{{{\boldsymbol{k}}}}}({{{\boldsymbol{r}}}}){e}^{-i{{{\boldsymbol{k}}}}\cdot {{{\boldsymbol{r}}}}},$$

(21)

where ψ₀(r), η_k, η_−k, are slowly-varying complex fields on the length scale a₀ of the director field modulations. k is the wave vector of the modulations which is ${{{\boldsymbol{k}}}}=\frac{2\pi }{{a}_{0}}{{{{\boldsymbol{e}}}}}_{x}$ due to the initial condition. We can extract the complex amplitude of a k mode by a demodulation of ψ,

$${\eta }_{{{{\boldsymbol{k}}}}}=\langle \psi {e}^{-i{{{\boldsymbol{k}}}}\cdot {{{\boldsymbol{r}}}}}\rangle ,$$

(22)

through the convolution with a Gaussian kernel denoted by 〈 ⋅ 〉, which filters out the small-scale variations, Eq. (40). The modulation length scale a₀ and the equilibrium value η₀ of ∣η_k∣ are found numerically to be a₀ = 10.6 and η₀ = 0.20 for the given parameters. From the order parameter η_k, we can construct the defect density ${\rho }^{({\eta }_{{{{\boldsymbol{k}}}}})}$ as for the complex wavefunction in the BEC. This field locates the dislocations from the nematic arches as shown in panel (b) at t = 240, just prior to the nucleation of nematic defects.

By also showing the reduced defect field ρ^(Q) associated with the rotational symmetry (Fig. 4c), we clearly notice that each dislocation detected by ${\rho }^{({\eta }_{{{{\boldsymbol{k}}}}})}$ is a source for the nucleation of a dipole of half-integer defects. The precursory pattern of the two bound defects prior to nucleation is similar to the pattern retained by the dipole annihilation in the BEC. However, for active nematics, the bound state is associated with a dislocation in the periodic arches, while for BECs it is a source of quantum pressure. We observe numerically that the melting of the smectic-like arches is mediated by the dissociation of the dislocations into dipoles of ± 1/2 nematic defects. This occurs very fast and simultaneously at various locations, such that the system quickly transitions to active turbulence. Notice also that the core size of the dislocations in the periodic arches is bigger than the core size of the ± 1/2 nematic defects that form in the transition. To quantify such nucleation events, we compute the defect velocity Eq. (11) associated to the charged defect density ρ^(Q) which is localized in well-defined blobs of opposite signs around a dislocation as illustrated in Fig. 4b, c. By averaging the speed around these blobs, we can track the defect speed v = 〈∣v∣〉 as function of time and show that prior to dissociation, the defects are in a bound state while afterwards they move apart as ± 1/2 defects with different speeds as shown in Fig. 4e. Notice that the − 1/2 defect slows down while the + 1/2 acquires a net speed related to its self-propulsion.

To summarize this part, our analysis offers an alternative perspective on the onset of active turbulence using the presence of competing symmetries. The transition to a turbulent state from a periodic arch state seems to be mediated by the dissociation of one type of topological defect into a different kind associated with changes in the global symmetries. In the following section, we study a system where the order parameters with O(n)-symmetry are found by decomposing a more complicated topological space.

Defect structures: solid crystals

We focus here on the study of defects and collective order in crystals. The ground state manifold of the crystal can be factorized in fundamental ${{{{\mathcal{S}}}}}^{1}$ spaces, which has a straightforward physical interpretation related to the crystal’s Bravais reference lattice reflecting the broken translational symmetry. As discussed below, this implies that a dislocation, i.e., a topological defect in the crystal, can be represented by bound vortices in the amplitudes of the fundamental periodic modes. Indeed, by applying the formalism introduced in the section “Non-singular defect fields”, analogies with previously discussed systems emerge, as well as peculiar features that will be discussed in detail.

In the conserved Swift-Hohenberg modeling of crystal lattices, commonly named phase-field crystal (PFC)^44,45, the order parameter is a weakly distorted periodic scalar field ψ(r), and can be approximated as

$$\psi ({{{\boldsymbol{r}}}})=\bar{\psi }+\mathop{\sum }\limits_{n=1}^{N}{\eta }_{n}{e}^{i{{{{\boldsymbol{q}}}}}^{(n)}\cdot {{{\boldsymbol{r}}}}},$$

(23)

where $\bar{\psi }$ and ${\{{\eta }_{n}\}}_{n = 1}^{N}$ are slowly varying (on the lattice unit length scale) amplitude fields, and N is the number of reciprocal lattice vectors ${\{{{{\boldsymbol{{q}}}^{(n)}}}\}}_{n = 1}^{N}$ taken into consideration. Disordered or liquid phases are described by η_n(r) = 0. For a perfect lattice, $\bar{\psi }({{{\boldsymbol{r}}}})={\psi }_{0}$ and η_n(r) = η₀ are constant, and an affine displacement r → r − u amounts to a phase change ${\eta }_{n}={\eta }_{0}{e}^{-{{{{\boldsymbol{q}}}}}^{(n)}\cdot {{{\boldsymbol{u}}}}}$. The displacement field u supports dislocations, which are line topological defects. For a path $\partial {{{\mathcal{M}}}}$ in real space circling one dislocation, the charge is given by the vector difference between the end and starting point, namely the Burgers’ vector b,

$${\oint }_{\partial {{{\mathcal{M}}}}}d{{{\boldsymbol{u}}}}=-{{{\boldsymbol{b}}}},$$

(24)

(minus sign by convention). The corresponding dislocation density tensor α_ij is defined through the integral of some 2D surface ${{{\mathcal{M}}}}$ bounded by $\partial {{{\mathcal{M}}}}$

$${\int}_{{{{\mathcal{M}}}}}{\alpha }_{ij}{n}^{i}dS={b}_{j},$$

(25)

where n is the normal vector to the surface element dS. By multiplying Eq. (24) with a reciprocal lattice vector q⁽ⁿ⁾ of the structure, we get

$${\oint }_{\partial {{{\mathcal{M}}}}}d({{{{\boldsymbol{q}}}}}^{(n)}\cdot {{{\boldsymbol{u}}}})=-2\pi {s}_{n},$$

(26)

where s_n is an integer by definition of the reciprocal lattice vector. This shows that the phase of an amplitude θ_n ≡ (−q⁽ⁿ⁾ ⋅ u) is a topological order parameter that has integer winding numbers, i.e., ${\theta }_{n}\in {{{{\mathcal{S}}}}}^{1}$.

The amplitude η_n acts as an order parameter in ${{{{\mathcal{D}}}}}^{2}$, i.e., ${\Psi }_{1}^{(n)}=\Re ({\eta }_{n})$ and ${\Psi }_{2}^{(n)}=\Im ({\eta }_{n})$. A topological description of dislocations using the HM-framework has been provided in two and three dimensions in refs. ^34,36. Here, we adopt an alternative and convenient description using the charge density from Eq. (12), which is a vector field for 3D lattices

$${\rho }_{i}^{({\eta }_{n})}=\frac{{D}_{i}^{(n)}}{\pi {\Psi }_{0}^{2}},$$

(27)

where

$${{{{\boldsymbol{D}}}}}^{(n)}=\left(\begin{array}{c}({\partial }_{y}{\Psi }_{1}^{(n)})({\partial }_{z}{\Psi }_{2}^{(n)})-({\partial }_{y}{\Psi }_{2}^{(n)})({\partial }_{z}{\Psi }_{1}^{(n)})\\ ({\partial }_{x}{\Psi }_{2}^{(n)})({\partial }_{y}{\Psi }_{1}^{(n)})-({\partial }_{x}{\Psi }_{1}^{(n)})({\partial }_{z}{\Psi }_{2}^{(n)})\\ ({\partial }_{x}{\Psi }_{1}^{(n)})({\partial }_{y}{\Psi }_{2}^{(n)})-({\partial }_{x}{\Psi }_{2}^{(n)})({\partial }_{y}{\Psi }_{1}^{(n)})\end{array}\right).$$

(28)

By contracting Eq. (25) with q^j, we can relate the dislocation density tensor with the defect charge density in a given amplitude³⁶

$${\alpha }_{ij}=\frac{2d}{N{\eta }_{0}^{2}}\mathop{\sum }\limits_{n=1}^{N}{D}_{i}^{(n)}{q}_{j}^{(n)},$$

(29)

where d is the spatial dimension. The amplitudes η_n used to calculate D⁽ⁿ⁾ are extracted from the phase-field ψ as in Eq. (22), and only the modes corresponding to the shortest reciprocal lattice vectors are used to calculate α_ij.

Next, we focus on two examples to highlight insights obtained from using this approach. We consider the nucleation of dislocations in a square lattice from the point of view of its precursory pattern formations and quantify the dislocation core size. Then, we consider the classical inclusion problem of a rotated spherical crystal embedded in another crystal with the same lattice symmetry, to show how the surface of the inclusion changes its topology as a function of the lattice misorientation.

Dislocations in 2D square lattices

A minimal PFC free energy which can be minimized by a square lattice reads^46,47

$${F}_{\psi }^{{{{\rm{sq}}}}}=\int\,{d}^{2}r\left(\frac{1}{2}{({{{{\mathcal{L}}}}}_{1}{{{{\mathcal{L}}}}}_{2}\psi )}^{2}+\frac{r}{2}{\psi }^{2}+\frac{1}{4}{\psi }^{4}\right),$$

(30)

where ${{{{\mathcal{L}}}}}_{X}=X+{\nabla }^{2}$ and r is a parameter. We recall that PFC energy functionals describe order–disorder (solid–liquid) phase transitions. The minimizer field ψ of (30), for certain model parameters, has a perfect square lattice symmetry with an accurate two-mode amplitude expansion

$$\psi =\bar{\psi }+\mathop{\sum }\limits_{n=1}^{2}{\eta }_{n}{e}^{i{{{{\boldsymbol{q}}}}}^{(n)}\cdot {{{\boldsymbol{r}}}}}+\mathop{\sum }\limits_{n=3}^{4}{\eta }_{n}{e}^{i{{{{\boldsymbol{q}}}}}^{(n)}\cdot {{{\boldsymbol{r}}}}}+\,{{\mbox{c.c.}}}\,,$$

(31)

where {q⁽ⁿ⁾} = {(1, 0), (0, 1), (1, 1), (1, − 1)} are the reciprocal lattice vectors of the square lattice with lengths 1 and $\sqrt{2}$. This sets the characteristic length a₀ = 2π of the system, which is the width of the square unit cell. At equilibrium, the amplitude field η_n goes to the equilibrium values η_1,2 → A_sq, η_3,4 → B_sq. The characteristic unit of stress is given by the elastic shear modulus $\mu =16{B}_{{{{\rm{sq}}}}}^{2}$⁴⁷. The dislocation density tensor can be factorized as ${\alpha }_{ij}={t}_{i}{{{{\mathcal{B}}}}}_{j}$, where ${{{\mathcal{B}}}}$ is a 2D Burgers vector density and t the tangent vector to the dislocation line. In two dimensions, we define t to point out-of-plane so that the Burgers vector density is given by

$${{{\mathcal{B}}}}=({\alpha }_{31},{\alpha }_{32}),$$

(32)

where α_ij can be computed by using q^(1,2). We initiate a perfect square lattice of 101 × 101 unit cells and use the sHPFC model of ref. ⁴⁸ to apply a local stress in the central region which causes the nucleation of a dislocation dipole. The PFC deforms gradually, trying to account for the externally imposed stress, increasing from linear to non-linear strains until nucleation of a pure ± a₀e_x dislocation dipole. Once formed, these dislocations move under the action of the Peach-Koehler force⁴⁹, namely they separate at large speeds due to the external stress and slow down as they reach the far-field regions of the crystal. Simulation details are given in the “Methods” section. Figure 5 shows the region of applied stress during the nucleation event. Like for the nucleation of nematic defects, the nucleation is singled by a precursory localized pattern formation in the Burgers vector density, which corresponds to a bound dipole of phase slips. While variations only in the phase of the complex amplitudes are associated with linear elastic perturbations, non-linear elastic strains cause a decrease in the equilibrium value of the amplitudes⁵⁰ and so produce a signal in the reduced defect density given by the expression of the dislocation density. Thus, the excitations visible in the dislocation density ${{{\mathcal{B}}}}$ prior to nucleation are due to non-linear elastic strains. From the signal profile, Fig. 5c, we observe that these large non-linear elastic strains can be connected to a bound dislocation dipole.

Fig. 5: Nucleation of a dislocation dipole in a square PFC model.

a The PFC at t = 1600 prior to the nucleation of (b) a dislocation dipole at t = 1800. Panels (c) and (d) show the x-component ${{{{\mathcal{B}}}}}_{x}$ of the dislocation density ${{{\mathcal{B}}}}$ at t = 1600 and t = 1800, respectively. The magnitude of ${{{{\mathcal{B}}}}}_{y}$ is, in both cases, two orders of magnitude smaller and not shown. e The average speed v = 〈∣v∣〉 at the nucleation site of positive charge (${{{{\mathcal{B}}}}}_{x}\, > \,0$) and negative charge (${{{{\mathcal{B}}}}}_{x}\, < \,0$) where the dashed line indicates the time of nucleation (see text).

From the defect density corresponding to η₁ for q = (1, 0), we can also determine the average speed v = 〈∣v∣〉 of dislocations with positive and negative charge before and after nucleation. The defect speed as a function of time is shown in Fig. 5e. Like for the nucleation of defects in the active nematic, we observe a speed build up prior to nucleation succeeded by a relaxation to a constant speed. Unlike the ±1/2 defects in active nematics, however, both dislocations are equally mobile in this case.

The Burgers vector density, in addition to describing the process of nucleation itself, provides us with useful information about the defect core. To extract the core size directly from the Burger vector density without free-tuning parameters, we consider a coarse-grained version of the PFC model, namely its amplitude expansion (APFC)^51,52. This approach gives access to phases and lattice deformation directly rather than through the demodulation of Eq. (22). It builds on the definition of a free energy functional F_η derived from the PFC free energy ${F}_{\psi }^{{{{\rm{sq}}}}}$ under the approximation of slowly-varying amplitudes. We simulate a square lattice hosting dislocations in a static, periodic configuration, and focus on a single defect therein. The expression for F_η, the choice of q⁽ⁿ⁾, and details of the simulation setup are given in the “Methods” section. For the given lattice structure, the extension of its core depends on the parameters ${r}^{{\prime} }$ and s in the free energy F_η. The parameter ${r}^{{\prime} }$ corresponds to a phenomenological temperature controlling a first-order order–disorder phase transition at ${r}^{{\prime} }={r}_{0}$ with r₀ the critical point and ordered (disordered) phase for ${r}^{{\prime} }\, < \,{r}_{0}$ (${r}^{{\prime} }\, > \,{r}_{0}$), and s is a constant scaling the elastic moduli^53,54. $\Delta r={r}_{0}-{r}^{{\prime} }$ is referred to as the quenching depth. These parameters affect the competition among gradient terms and the bulk energy terms in F_η. Figure 6a, b illustrates two different core sizes for the same dislocation obtained with different values for ${r}^{{\prime} }$ and s. They show the reconstructed densities obtained by computing Eq. (23) with the numerical solution for the amplitudes (first column), the Burgers vector density component ${{{{\mathcal{B}}}}}_{x}$ (second column), a plot of ${{{{\mathcal{B}}}}}_{x}(x,0)$ and ${{{{\mathcal{B}}}}}_{x}(0,y)$ (third column, empty symbols) with Gaussian fits (solid lines). The data fitting is obtained via $G\exp (-{x}^{2}/2{\sigma }_{x}^{2}-{y}^{2}/2{\sigma }_{y}^{2})$ with G, σ_x and σ_y fitting parameters (dashed lines), well reproducing its shape and allowing for an estimation of the core size. The definition here introduced for the Burgers vector density fully characterizes the loss of coherency at the dislocation core. Importantly, it realizes a spreading of the topological charge at the core similar to non-singular continuum theories based either on regularization of singularities⁵⁵ or within strain-gradient elasticity theories^56,57.

Fig. 6: Dislocation core size near melting by APFC modeling.

a, b Reconstructed density (left), ${{{{\mathcal{B}}}}}_{x}$ (center), and ${{{{\mathcal{B}}}}}_{x}$ along x and y direction for a relatively small and large core size, respectively obtained with a Δr = 10⁻⁴, s = 3.16 and b Δr = 10⁻¹, s = 1, with r₀ = 7.455 × 10⁻² the critical point. Symbols show values from APFC simulations; dashed lines correspond to Gaussian fits. The latter are exploited to quantify the size of the core in terms of the variance along x and y, namely σ_x and σ_y. c Periodic modes ${\eta }_{n}{e}^{i{{{{\bf{q}}}}}_{{{{\bf{n}}}}}\cdot {{{\bf{r}}}}}+\,{{\mbox{c.c.}}}\,$ for the density in panel (b). d Core size in terms of σ_x as a function of the order–disorder correlation length w, for various values of s and ${r}^{{\prime} }$ (the latter shown by different colors and symbols). e Comparison of σ_x and σ_y as function of w for Δr = 0.00464 and s ∈ [10⁻¹, 3.16].

The amplitude expansion defined in Eq. (23), and thus the density field ψ, correspond to the sum of plane waves (Fourier modes) which are periodic stripe phases similar to the one shown in Fig. 4. The dislocation in the crystal then corresponds to the superposition of defects in such stripe phases. Interestingly, dislocations do not necessarily correspond to a defect for all the coupled stripe phases. Indeed, by applying Eq. (1) to the phase of the amplitudes one gets −∮q⁽ⁿ⁾ ⋅ u = 2πq⁽ⁿ⁾ ⋅ b. At least for perfect dislocations, those having a translation vector of the lattice as Burgers vector, we have that q⁽ⁿ⁾ ⋅ b = 0, for some n. Therefore, at the dislocation core, a different ordered phase forms as some amplitudes may have non-singular phases and, in turn, do not vanish. This differs from the case of dislocations forming in pure stripe phases, e.g., in Fig. 4, where the single complex amplitude vanishes, pointing to a disordered phase. In Fig. 6c, the fields ${\eta }_{n}{e}^{i{{{{\boldsymbol{q}}}}}^{(n)}\cdot {{{\boldsymbol{r}}}}}$ entering the sum in Eq. (23) are reported. Three out of four stripe phases (n = 1, 3, 4) vanish at the core, while one (n = 2) features a small variation of its amplitudes with no topological content.

The defect core can then be interpreted as a transition region between two different ordered phases, one of which is present at the dislocation core only. To explore the analogy with phase interfaces, we compare its extension with the width of a solid–liquid (order–disorder) interface, w, which measures the correlation length for these phases. We find some analogies and differences in the dependence on the parameters entering the free energy. Traveling-wave solutions exist for solid–liquid (order-disorder) interfaces with amplitudes having hyperbolic tangent profiles, $\eta \propto ({\eta }_{0}/2)\{1-\tanh [(x-Vt)/w]\}$ with $w\propto \sqrt{g}/(1+\sqrt{1-8{r}^{{\prime} }/9{r}_{0}})$, V the interface velocity along its normal and g a parameter in the free energy which multiplies gradient terms and scales the elastic constants^54,58,59 (see also “Methods”). For a given set of parameters, we determine the specific amplitude profile and w by fitting the result of numerical calculations with the hyperbolic tangent profile mentioned above for an interface with normal along the x-axis (〈10〉 crystallographic direction, further details are reported in the “Methods” section, Measuring the size of the dislocation core through σ_x and σ_y from a Gaussian fit as in Fig. 6a, b, we find that it scales linearly with w when varying g, while a different scaling is observed when varying ${r}^{{\prime} }$, c.f. Fig. 6d. Here g is an energy scale associated with amplitudes gradients, similar to theories based on Ginzburg-Landau energy functionals⁵⁹. ${r}^{{\prime} }$, instead, affects the equilibrium values of the amplitudes, which are qualitatively different for an interface, where they all vanish in the disordered phase, and a defect, where some amplitudes are non-zero owing to a non-singular phase (see Fig. 6c). Also, for ${r}^{{\prime} }\,\ne\, {r}_{0}$, interfaces move, which affects the width w⁶⁰. A more detailed analysis would require finding a solution for the amplitudes’ profile at defects, which goes beyond the goals of this investigation and will be addressed in future work.

The evaluation of the Burgers vector density also allows for the characterization of anisotropies in the behavior of phases at the core as illustrated in Fig. 6d. σ_y/σ_x ≈ 0.75 throughout the whole range of parameters investigated here as also illustrated in Fig. 6e. This may be ascribed to the asymmetry introduced by the specific orientation of the Burgers vector. We conclude that, for systems described by order parameters as in the phase-field crystal model, as well as in descriptions exploited in previous sections, the defect density may be exploited to characterize the loss of coherency at defects.

Order transition for 3D crystal inclusions

Like the melting of translational order in the nematic liquid crystal through the nucleation of defects in the nematic field, the global translational order in a single crystal is also destroyed under large deformations and rotations. To highlight this, we use a full 3D PFC model corresponding to a cubic lattice for which the PFC density in the one-mode approximation reads as

$$\psi ({{{\boldsymbol{r}}}})={\psi }_{0}+\mathop{\sum}\limits_{{{{\boldsymbol{q}}}}\in {{{{\mathcal{R}}}}}_{{{{\rm{bcc}}}}}^{(1)}}{\eta }_{0}{e}^{i{{{{\boldsymbol{q}}}}}^{(n)}\cdot {{{\boldsymbol{r}}}}},$$

(33)

where ${{{{\mathcal{R}}}}}_{{{{\rm{bcc}}}}}^{(1)}$ are the reciprocal lattice vectors of the bcc Bravais lattice with unit length⁴⁷. This sets the length of the bcc unit cell as ${a}_{0}=2\pi \sqrt{2}$. We consider spherical inclusions with radius 17a₀, rotated at an angle θ_rot about the [1, 1, 1]-axis. The initial condition is relaxed by dissipative dynamics with an appropriate symmetry-conserving free energy; see further details in the “Methods” section. We choose three representative angles θ_rot and calculate the Frobenius norm $| \alpha | =\sqrt{{\alpha }^{ij}{\alpha }_{ij}}$ of α at each angle. Since ∣α∣ > 0, we plot its isosurface at half its maximum value $| \alpha {| }_{M}=\mathop{\max }\nolimits_{{{{\boldsymbol{r}}}}}(| \alpha | ({{{\boldsymbol{r}}}}))$ in Fig. 7 for three representative misorientation angles θ_rot. For small lattice misorientations, ∣α∣ ≪ 1, indicating only slight non-linear elastic excitations (and no fully formed dislocations) at the interface between the inclusion and the matrix. As expected, these non-linear strains are largest in the plane perpendicular to the rotation axis, since the rotation deformation field scales with distance from the rotation axis. Notably, we observe a three-fold symmetry in the profile of ∣α∣, which can be ascribed to the underlying crystallographic orientation. For larger values of θ_rot, the non-linear distortions increase and localize into a network of dislocations. Notice that such a defect network is determined directly by the Burgers vector density rather than through arbitrary reconstructions^61,62. The description breaks down at large misorientations, as witnessed by the decrease in the magnitude of the defect density field since there is no longer a global translational order. Indeed, large misorientations lead to the nucleation of grain boundaries which are fully described by accounting for the bicrystallography of the two crystals meeting at the interface rather than the deformation with respect to a reference lattice⁶³. Such a regime shift echoes the onset of active turbulence in the nematic liquid, where the description in terms of the order parameter ${\rho }^{({{{{\boldsymbol{\eta }}}}}_{{{{\boldsymbol{q}}}}})}$ also breaks down.

Fig. 7: Rotated inclusions in the bcc PFC model.

The panels show, for three representative rotation angles θ_rot the isosurface of the Frobenius norm of the coarse dislocation density tensor $| \alpha | =\sqrt{{\alpha }^{ij}{\alpha }_{ij}}$ at 50% of its maximal value $| \alpha {| }_{M}=\mathop{\max }\nolimits_{{{{\boldsymbol{r}}}}}(| \alpha | ({{{\boldsymbol{r}}}}))$, which is given in the panels.

Discussion

In-depth understanding and tailoring of collective behaviors require a unified description of defects associated with symmetry breaking and the non-topological excitations of ground states. Here, we proposed a systematic way of deriving reduced defect fields from order parameters associated with O(n) broken symmetries which captures topological defects, localized non-linear excitations, and their dynamics. This enables the non-singular description of defects and their interaction, accounting for precursory and resulting patterns involving non-topological excitations. In this way, short-scale interactions between topological defects may be more accurately described, since features such as core overlap and high-energy excitations become more prominent at shorter length scales. This paves the way for a more thorough characterization of defect interactions, particularly in cases where the defects get close or are annihilated, as in the applications shown above. Moreover, the proposed framework can be used to study concurrent symmetry breakings and order transitions. Applications to systems of general interest, such as superfluids, active nematics, and solid crystals, are shown to showcase the considered framework, while we envisage applications in many other contexts.

We have shown that the method accurately tracks topological defects since these appear as localized blobs in the defect density field. The associated current density and velocity field determine the kinematics of the defects, and its utility has been shown to extend beyond tracking the velocity of topological defects. For example, in the case of the motion of vortices in a BEC, the velocity field accounts for both the overall velocity of the defect and local variations associated with the early-stage rearrangements of the defect core evolving towards its stationary shape. Thus, the uniformity of the velocity field over the core extent tests whether the frozen-core approximation¹ is valid. For active nematics and solid crystals, the velocity formula is shown to track the dynamics of defect dipoles, during, and after the nucleation of topological defects, pointing at interesting analogies and differences between processes in different physical systems. The rigorous derivation of these fields given in the Supplementary Notes for any dimensions makes the equations readily applicable to tracking topological defects and localized excitations in general.

We have found interesting features and insights about the evolution of these systems with broken symmetries. After the annihilation of the vortex dipole in the BEC, the remaining shock wave produces a signal in the defect density field that echoes the charge density pattern of the dipole, remnant also of other similar observations during mass-driven vortex collision⁶⁴. In active nematics, the large cores of the dislocation in the translational order harbor a bound dipole of orientational defects associated with the rotational order. This picture presents the idea of a hierarchy of topological defects, where the defects associated with one symmetry can spontaneously dissociate into stable defects for a different symmetry and melt the former ordered state. This is a non-equilibrium transition that echos the equilibrium Kosterlitz-Thouless transition for melting of 2D crystals via the hexatic phase⁶⁵.

In the case of a 3D crystal, a rotated inclusion was shown to be described as a network of topological defects (dislocations) up to a point before these dissociated into other types of defects (grain boundaries) and the global orientational order was destroyed. The best topological description of polycrystalline materials is an open challenge, even though candidates, such as interacting disconnections⁶³, exist. Applying this formalism to such topologies is a fascinating avenue of research. Employing the APFC framework, where the periodic nature of crystal densities is inherently coarse-grained, we have shown that dislocation cores in Swift-Hohenberg theories emerge as transition regions from crystalline to pointwise stripe-like phases. When approaching the solid–liquid coexistence limit, analogies between the dislocation core size and the extensions of order-disordered interfaces have been found.

Finally, while the whole framework is presented for systems with one broken rotational symmetry, it is a powerful tool that can be generalized to systems with multiple broken symmetries and reveal hidden hierarchies of topological defects associated with each symmetry, laying the foundation for unified theories in systems characterized by collective behaviors.

Methods

Bose-Einstein condesates

The damped Gross Pitaevskii equation, Eq. (16), is solved by using a Fourier pseudo-spectral integration scheme which is described in detail in ref. ³⁰. We use a periodic grid of size [− 32, 32] × [− 32, 32] with spatial discretization Δx = Δy = 0.25. To initialize the dipole we use the ansatz $\psi ={\Pi }_{\alpha = 1}^{2}\chi (| r-{r}_{\alpha }| ){e}^{i{q}_{\alpha }{\theta }_{\alpha }}$, where r_α is the position of the vortex labeled α, ${\theta }_{\alpha }=\arctan [(y-{y}_{\alpha })/(x-{x}_{\alpha })]$, and

$$\chi (r)=\left\{\begin{array}{ll}r,\quad &r \,< \,1\\ 1,\quad &r \,> \,1\end{array}\right..$$

(34)

This order parameter is then evolved in imaginary time, t → iτ, with γ = 0 to lower the energy and find a better estimate for the core structure of the vortices we use as the initial condition.

Active nematic liquid crystals

The evolution of the Q-tensor follows dissipative dynamics coupled with an incompressible Stokes flow⁴³

$${\partial }_{t}{Q}_{ij}+{{{\boldsymbol{v}}}}\cdot \nabla {Q}_{ij}-{Q}_{ik}{\Omega }_{kj}+{\Omega }_{ik}{Q}_{kj}=\lambda {{{{\mathcal{W}}}}}_{ij}+{\gamma }^{-1}{H}_{ij},$$

(35)

$$(\Gamma -\eta {\nabla }^{2}){v}_{i}={\partial }_{j}(\alpha {Q}_{ji})-\nabla p,\quad \nabla \cdot {{{\bf{v}}}}=0,$$

(36)

where v is the flow velocity that advects the nematic structure, p is the fluid pressure, Γ is the friction with a substrate, η is the viscosity and αQ is the active stress. The vorticity tensor 2Ω_ij = (∂_iv_j − ∂_jv_i) rotates the nematic structure, λ is the flow alignment parameter which aligns the nematic orientation in the direction of shear

$${{{{\mathcal{W}}}}}_{ij}={E}_{ij}+({E}_{ik}{Q}_{kj}+{Q}_{ik}{E}_{kj})-{Q}_{lk}{E}_{kl}({\delta }_{ij}+{Q}_{ij}),$$

with the trace less strain rate 2E_ij = (∂_iv_j + ∂_jv_i − δ_ij∂_kv_k). The molecular field

$${H}_{ij}=K{\nabla }^{2}{Q}_{ij}+A(B-2{Q}_{kk}^{2}){Q}_{ij}.$$

(37)

controls the relaxation to equilibrium with γ as the rotational diffusivity. We have here assumed a single Frank elastic constant K, treating splay and bend distortions similarly. The second term in the molecular field is a relaxation to a homogeneous nematic state. The parameter A is the quench depth and B sets the value of the order parameter ${S}_{0}=\sqrt{B}$ in the homogeneous state. We discretize the above equations on a [−64, 64] × [−64, 64] grid with spatial discretization Δx = Δy = 0.5, and solve the system using pseudo-spectral methods. The parameters are set to K = Γ = γ = 1, A = λ = η = 0.5, B = 2 and α = − 1.4. The initial state is $S=\sqrt{2}$ with the angle of the director θ being uniformly distributed in the interval (−0.05, 0.05). We solve the equations for the flow field, Eq. (36), in Fourier space and evolve the equation for the Q tensor, Eq. (35), using the same scheme as for the BEC.

2D square lattice PFC

To simulate the PFC dynamics, we use the sHPFC model proposed in ref. ⁴⁸, namely

$${\partial }_{t}\psi =\Gamma {\nabla }^{2}\frac{\delta {F}_{\psi }^{{{{\rm{sq}}}}}}{\delta \psi }-{{{\boldsymbol{v}}}}\cdot \nabla \psi ,$$

(38)

coupled to a momentum equation for ∂_tv

$${\rho }_{0}{\partial }_{t}{{{\boldsymbol{v}}}}=\langle {\tilde{\mu }}_{{{{\rm{c}}}}}\nabla \psi -\nabla \tilde{f}\rangle +{\Gamma }_{S}{\nabla }^{2}{{{\boldsymbol{v}}}}+{{{{\boldsymbol{f}}}}}^{({{{\rm{ext}}}})}.$$

(39)

〈 ⋅ 〉 is a convolution with a Gaussian kernel given by

$$\langle \tilde{X}\rangle =\int\,d{{{{\boldsymbol{r}}}}}^{{\prime} }\frac{\tilde{X({{{{\boldsymbol{r}}}}}^{{\prime} })}}{2\pi {w}^{2}}\exp \left(-\frac{{({{{\boldsymbol{r}}}}-{{{{\boldsymbol{r}}}}}^{{\prime} })}^{2}}{2{w}^{2}}\right),$$

(40)

which filters out variations on length scales smaller than w. The quench depth in Eq. (30) is set to r = − 0.3 and the average density to $\bar{\psi }=-0.3$. Parameters are set to Γ = 1, ρ₀ = Γ_S = 2⁻⁶, and an initial velocity field v = 0. We solve the system of coupled equations with a Fourier pseudo-spectral method. The spatial grid of the simulation is set to Δx = Δy = a₀/7. Further details can be found in ref. ⁴⁸.

In the simulation reported in Fig. 5, the perfect lattice is indented by an applied external force density given by a Gaussian profile ${{{{\boldsymbol{f}}}}}^{({{{\rm{ext}}}})}={f}_{0}\frac{(y-{y}_{0})}{{a}_{0}}\exp (-\frac{{({{{\boldsymbol{r}}}}-{{{{\boldsymbol{r}}}}}_{0})}^{2}}{2{w}^{2}}){{{{\boldsymbol{e}}}}}_{x}$. Above a critical strength f₀ = 3.5μ/a₀ and width w = a₀, this force causes the nucleation of a dislocation dipole.

2D square lattice APFC

The evolution of the amplitudes as delivered by the APFC model can then be directly expressed as

$$\frac{\partial {\eta }_{n}}{\partial t}=-{\left\vert {{{{\bf{q}}}}}^{(n)}\right\vert }^{2}\frac{\delta {F}_{\eta }}{\delta {\eta }_{n}^{* }},$$

(41)

with F_η the free energy depending on {η_n} that can be derived by substituting (31) in ${F}_{\psi }^{{{{\rm{sq}}}}}$ and integrating over the unit cell⁵⁹. By assuming constant $\bar{\psi }$ it reads

$${F}_{\eta }=\int\,{d}^{2}r\left(g\mathop{\sum }\limits_{n=1}^{N}| {{{{\mathcal{G}}}}}_{n}{\eta }_{n}{| }^{2}+W(\{{\eta }_{n}\})+C(\bar{\psi })\right),$$

(42)

with ${{{{\mathcal{G}}}}}_{n}=({\nabla }^{2}+2i{{{{\bf{q}}}}}^{(n)}\cdot \nabla )$, g a coefficient that controls elastic constants, $W(\{{\eta }_{n}\})={r}^{{\prime} }\Phi /2+(3/4){\Phi }^{2}-(3/4)\mathop{\sum }\nolimits_{n = 1}^{N}| {\eta }_{n}{| }^{4}+{f}^{{{{\rm{s}}}}}(\{{\eta }_{n}\})$, ${r}^{{\prime} }=r+3{\bar{\psi }}^{2}$, $\Phi =\mathop{\sum }\nolimits_{n = 1}^{N}| {\eta }_{n}{| }^{2}$, and f^s({η_n}) a symmetry-dependent polynomial in the amplitudes. For the square symmetry as encoded in Eq. (30) and the choice q⁽¹⁾ = (1, 0), q⁽²⁾ = (0, 1), q⁽³⁾ = (1, 1), q⁽⁴⁾ = (−1, 1) and {q⁽ⁿ⁾} = {−q⁽ⁿ⁻⁴⁾} for n = 5, …, 8, we have ${f}^{{{{\rm{s}}}}}(\{{\eta }_{n}\})=2\bar{\psi }({\eta }_{1}{\eta }_{2}{\eta }_{3}^{* }+{\eta }_{1}{\eta }_{2}^{* }{\eta }_{4})+3({\eta }_{1}^{2}{\eta }_{3}^{* }{\eta }_{4}+{\eta }_{2}^{2}{\eta }_{3}^{* }{\eta }_{4}^{* })+\,{{\mbox{c.c.}}}\,$, with $\{{\eta }_{n}^{* }\}=\{{\eta }_{n-4}\}$ for n = 5, …, 8 as ψ is a real function. Therefore, one may consider just η_n with n = 1, …, 4 as variables. $C(\bar{\psi })$ is a constant depending on $\bar{\psi }$⁵⁹, set here to $\bar{\psi }=-0.3$ as set in the corresponding PFC modeling of the 2D square lattice. ${r}^{{\prime} }$ corresponds to a phenomenological temperature. With r₀ the solid–liquid critical point, the solid crystalline phase is favored for ${r}^{{\prime} }\, < \,{r}_{0}$.

We simulate a stationary system hosting dislocations with the APFC model exploiting the (FEM) numerical approach with adaptive grid refinement outlined in refs. ^66,67. The semi-implicit integration scheme adopted for numerical simulations can be found therein. We consider dislocations with spacing L = 50a₀ arranged in a periodic, 2D matrix with alternating Burgers vectors $\pm {a}_{0}\hat{x}$. The system is initialized by setting the displacement field of dislocation known from classical continuum mechanics⁴⁹ in the phase of amplitudes, −q⁽ⁿ⁾ ⋅ u, and let relaxed according to the amplitudes evolution law (41). We can consider a system 2L × 2L by exploiting periodic boundary conditions.

In the section “Defect structures: solid crystals”, we characterize the extension of the core of dislocations through the field D⁽ⁿ⁾ as entering the definition of the dislocation density tensor α, Eq. (29). We compare the size of the defects extracted with the aid of Gaussian fits (see Fig. 6a, b) with the extension of a solid–liquid interface, w, computed numerically as the average of interface width for single amplitudes. This is obtained by initializing the solid phase with a straight interface having normal along the x-axis and letting the system evolve by Eq. (41) until reaching a steady state. Then, a fit of each amplitude with a function ${\phi }_{i}={\bar{A}}_{i}[1-\tanh (x-{\bar{x}}_{i})/{\bar{w}}_{i}]$, representing a traveling wave solution for a solid–liquid interface^54,58,60, is performed with ${\bar{A}}_{i}$, ${\bar{x}}_{i}$ and ${\bar{w}}_{i}$ parameters and the solid–liquid interface thickness extracted as $w=\mathop{\sum }\nolimits_{i = 1}^{4}{\bar{w}}_{i}/4$.

3D bcc lattice PFC

Numerical simulations reported in the section “Defect structures: solid crystals” are obtained by solving the classical PFC equation encoding dissipative dynamics,

$${\partial }_{t}\psi ={\nabla }^{2}\frac{\delta {F}_{\psi }^{{{{\rm{bcc}}}}}}{\delta \psi },$$

(43)

where ${F}_{\psi }^{{{{\rm{bcc}}}}}$ is a free energy functional that produces a stable bcc lattice, given by

$${F}_{\psi }^{{{{\rm{bcc}}}}}=\int\,{d}^{3}r\frac{1}{2}{({{{{\mathcal{L}}}}}_{1}\psi )}^{2}+\frac{r}{2}{\psi }^{2}+\frac{1}{4}{\psi }^{4}.$$

(44)

As parameters, we use r = − 0.3 and ψ₀ = −0.325 with spatial discretization Δx = Δy = Δz = a₀/7 and exploiting a Fourier pseudo-spectral integration scheme. We consider a 51 × 51 × 51 cubic crystal as matrix in which we embed a spherical inclusion with radius 17a₀ rotated at an angle θ_rot about the [1, 1, 1]-axis. This initial condition is obtained just by a rotation of grid points inside the inclusion. This leaves a sharp (and unphysical) interface which is regularized by letting this initial condition relax as dictated by Eq. (43) for 300 time steps with Δt = 0.1.