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

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.

The small-scale dynamics of interacting topological defects are crucial for the emergence of large-scale nonequilibrium phenomena, such as quantum turbulence in superfluids [9], spontaneous flows in active matter [10], or dislocation plasticity in crystals [11]. In fact, classical discrete modeling approaches such as point vortex models [12] and discrete dislocation dynamics [13] 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 [22]. 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 Section I A, 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 Sec. II A, we then develop a non-singular field theory as a generalization of the HMmethod which constitutes our primary reduced defect field. The method is then applied to the aforementioned physical systems in Secs. II B-II D. 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 [28]. Conclusions and perspectives for further study are outlined in the final section III.

A. 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 wellestablished for conventional systems, one often needs to define them for more exotic systems [29,30]. 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 [1]. 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 = S 1 , the unit circle. On 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 S 1 . This group is isomorphic to Z under addition since the difference between two loops that are not homotopic is how many times they have looped around the circle S 1 . Therefore, in regions of space where u is continuous and well-defined, a closed circuit ∂M in real space corresponds to a closed circuit in S 1 , and the topological charge s top contained in ∂M is given as an integer by the isomorphism between homotopy group of S 1 and Z. This topological charge is obtained from u = (cos θ, sin θ), by the contour integral which is invariant under any smooth deformations of ∂M. This also implies that by shrinking ∂M down to a point and given that s top is a constant, there must be regions inside ∂M where u is undefined. These are the topological defects that have a s top charge. Therefore, topological defects for R = S 1 in 2D are points with their charge determined by corresponding loop integration. On the other hand, such topological defects (with R = 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 |Ψ| = Ψ 0 , meaning that the ground state manifold is S n−1 , where n is the number of components of Ψ. The link between the order parameter Ψ, and 2D unit vector (director) field u ∈ 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. 1 (a-b). 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,31,32], active nematics [33,34], and deformed crystals [35][36][37]. Sticking to O(2)-symmetry in two dimensions and using the definition of a topological charge given in Equation (1), it is possible to express the topological defect density in terms of the zeros of the order parameter Ψ [26] tracked by Dirac-delta functions as (2) where q α and r α are, respectively, the charge and position of the topological defect α, δ (2) (Ψ) = δ(Ψ 1 )δ(Ψ 2 ), and D(r) is the (signed) Jacobian determinant of the map Ψ, where ij are the components of the Levi-Civita tensor in real space. The Levi-Civita tensor˜ 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˜ are simply the Levi-Civita (permutation) symbols. Note that Equation (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. [26] by considering as explicit ansatz a negative point defect, but can, in general, be justified using differential forms. Nominally, the D field in Equation (3) is evaluated at the location of the defect only, because of the δ-function in Equation (2).

II. RESULTS
A. Non-singular defect fields The δ-function in the topological charge density of Equation (2) locates the topological defects at singular points where u is undefined. In O(2) models, however, even though the ground state manifold is 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 [38]. Motivated by this, we seek to generalize Equation (2) in a way that will avoid singularities in the resulting charge density.
Since the equilibrium value Ψ 0 of |Ψ| is constant, the order parameter effectively resides in D 2 , the unit disk. We propose in this paper that the simplest generalization of s top is to consider the relative area of D 2 swept by Ψ on the circuit ∂M. During an infinitesimal displacement along ∂M, Ψ sweeps the infinitesimal area given by half of the parallelogram spanned by Ψ and dΨ. This (signed) area is given by 1 2˜ mn Ψ m dΨ n , see Fig. 2. The complete area of D 2 is πΨ 2 0 , and we define the charge s as the area swept by Ψ relative to the area of D 2 , where ∂M is defined as in Equation (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 ∂M in the far field of a topological defect, where |Ψ| = Ψ 0 , s = s top . To see this, note that if |Ψ| = Ψ 0 , then the infinitesimal area swept by Ψ is simply 1 2 Ψ 2 0 dθ, which inserted into Equation (4) gives Equation (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 where ρ(r) is the charge density of s, given by Whereas ρ top describes topological defects as point singularities in the physical space, ρ describes topological defects with a finite core size.
The time derivative of Equation (6) gives a continuity equation with the current density determined by the evolution of the order parameter Thus, ρ is a globally conserved quantity, and the change in s contained in a circuit ∂M is given by where dn is an infinitesimal surface area normal to the circuit ∂M. Far away from defects, |Ψ| = Ψ 0 and the time evolution of Ψ is carried by its phase θ(r, t) through Ψ = Ψ 0 (cos θ, sin θ) which can be inserted in Equation (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. 1(cf). We will exemplify this distinction in the applications. Due to the standard continuity form of Equation (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., v = 1 ρ (J + 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 This equation can be inverted to uniquely determine v (Ψ) if det(∂ i Ψ n ) = D(r) = 0. To find v (Ψ) where this condition holds true, i.e. the regions of interest where also ρ(r) = 0 from Equation (6), it is then possible to invert Equation (10). However, it is easier to insert ∂ t Ψ = −(v (Ψ) · ∇)Ψ into the expression J /ρ and see that it is the solution of Equation (10). Thus, to fix the gauge on v, we set K = 0 to get v = v (Ψ) and find 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, Equation (11) will give this velocity in the region of the core, which we will show in Sec. II B. 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, Equation (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 nsphere spanned by Ψ = (Ψ 1 , ..., Ψ n ) and normalizing by the volume V n Ψ n 0 of the n-sphere. We show in the Supplementary Notes [28] the formal derivation, and here we state the result that the charge density becomes (13) Generalizing the derivation of the defect kinematics, we find general expressions for the reduced defect velocity field where [ν 1 ν 2 ...ν n ] is the antisymmetrization over the indices ν 1 ν 2 ...ν n . Equation (15) (12) and (14) are the primary general expressions of the re-duced 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 volumeoriented 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 Equation (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 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 [39], or with n = 3, such as the 3D Heisenberg model of ferromagnetism, where the defect density will show emergent magnetic monopoles [40]. For further discussions, see the Supplementary Notes [28].
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 Ψ ∈ D 2 so that the method can be directly applied.

B. 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 wave function ψ, and its evolution can be described by damped Gross Pitaevskii equation [18,41] 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 Ψ = (Ψ 1 , Ψ 2 ) through ψ = Ψ 1 +iΨ 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 |ψ| = Ψ 0 = µ/g. Topological defects in the orientational (unit vector) field correspond to quantized vortices captured by the charge density field In this context, the D field (calculated from Ψ) has the physical interpretation of the generalized superfluid vorticity [32]. 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,42], leading to an increase in the quantum pressure, defined as Such excitations are detected by ρ (ψ) , and mediate the nucleation or annihilation of topological defects. To showcase this, we simulate a BEC as dictated by Equation (16) with an initial condition featuring two vortices at (x, y) = (±5, 0) Numerical details are reported in Sec. IV. 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 Equation (17) during the fast event of annihilating a vortex with an anti-vortex due to a small thermal drag.
The velocity field from Equation (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. 3(c). 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. 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.

C. 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 [43], 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 n = (cos(θ), sin(θ)) where S is an order parameter which is 0 in the disordered phase. Q is thus related to the D 2 order parameter Ψ field by Ψ = S 2 (cos(2θ), sin(2θ)). The evolution of the Q-tensor follows dissipative dynamics coupled with an incompressible Stokes flow with substrate friction [44]. Details on the evolution equation and its numerical method are reported in Sec IV. 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 |Ψ| ≡ Ψ 0 = √ B/2 dependent on the parameter B, which is defined in Sec. IV. Within the framework introduced in Sec. II A, this gives the following expression for the defect density which supports orientational defects with half-integer charge s top = ±1/2. In Fig. 4(a), 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) . The inset of Fig. 4(a) 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 where ψ 0 (r), η k , η −k , are slowly-varying complex fields on the length scale a 0 of the director field modulations. k is the wave vector of the modulations which is k = 2π a0 e x due to the initial condition. We can extract the complex amplitude of a k mode by a demodulation of ψ, through the convolution with a Gaussian kernel denoted by · , which filters out the small-scale variations, Equation (40). The modulation length scale a 0 and the equilibrium value η 0 of |η k | are found numerically to be a 0 = 10.6 and η 0 = 0.20 for the given parameters. From the order parameter η k , we can construct the defect density ρ (η 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. 4(c)), we clearly notice that each dislocation detected by ρ (η 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 Equation (11) associated to the charged defect density ρ (Q) which is localised in well-defined blobs of opposite signs around a dislocation as illustrated in Fig. Figure 4(b-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 Figure 4(e). Notice that the −1/2 defect slows down while the +1/2 acquires a net speed related to its selfpropulsion.
To summarise 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.

D. 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 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 Sec. II A, 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) [45,46], the order parameter is a weakly distorted periodic scalar field ψ(r), and can be approximated as whereψ and {η n } N n=1 are slowly varying (on the lattice unit length scale) amplitude fields, and N is the number of reciprocal lattice vectors {q (n) } N n=1 taken into consideration. Disordered or liquid phases are described by η n (r) = 0. For a perfect lattice,ψ(r) = ψ 0 and η n (r) = η 0 are constant, and an affine displacement r → r − u amounts to a phase change η n = η 0 e −q (n) ·u . The displacement field u supports dislocations, which are line topological defects. For a path ∂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, (minus sign by convention). The corresponding dislocation density tensor α ij is defined through the integral of some 2D surface M bounded by ∂M where n is the normal vector to the surface element dS. By multiplying Equation (24) with a reciprocal lattice vector q (n) of the structure, we get where s n is an integer by definition of the reciprocal lattice vector. This shows that the phase of an amplitude θ n ≡ (−q (n) ·u) is a topological order parameter that has integer winding numbers, i.e., θ n ∈ S 1 . The amplitude η n acts as an order parameter in D 2 , i.e., Ψ (n) 1 = (η n ) and Ψ (n) 2 = (η n ). A topological description of dislocations using the HM-framework has been provided in two and three dimensions in Refs. [35,37]. Here, we adopt an alternative and convenient description using the charge density from Equation (12), which is a vector field for 3D lattices where (28) By contracting Equation (25) with q j , we can relate the dislocation density tensor with the defect charge density in a given amplitude [37] where d is the spatial dimension. The amplitudes η n used to calculate D (n) are extracted from the phase-field ψ as in Equation 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 [47,48] F sq ψ = d 2 r where L X = X + ∇ 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 where {q (n) } = {(1, 0), (0, 1), (1, 1), (1, −1)} are the reciprocal lattice vectors of the square lattice with lengths 1 and √ 2. This sets the characteristic length a 0 = 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 µ = 16B 2 sq [48]. The dislocation density tensor can be factorized as α ij = t i B j , where 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 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. [49] 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 0 e x dislocation dipole. Once formed, these dislocations move under the action of the Peach-Koehler force [50], 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 Sec. IV. the nucleation is singled by a precursory localised 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 [51] 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 B prior to nucleation are due to non-linear elastic strains. From the signal profile, Fig. 5(c), we observe that these large non-linear elastic strains can be connected to a bound dislocation dipole.
From the defect density corresponding to η 1 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. 5(e). Like for the nucleation of de-fects 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) [52,53]. This approach gives access to phases and lattice deformation directly rather than through the demodulation of Equation (22). It builds on the definition of a free energy functional F η derived from the PFC free energy F sq ψ under the approximation of slowlyvarying 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 (n) , and details of the simulation setup are given in Sec. IV. For the given lattice structure, the extension of its core depends on the parameters r and s in the free energy F η . The parameter r corresponds to a phenomenological temperature controlling a first-order order-disorder phase transition at r = r 0 with r 0 the critical point and ordered (disordered) phase for r < r 0 (r > r 0 ), and s is a constant scaling the elastic moduli [54,55]. ∆r = r 0 − r is referred to as the quenching depth. These parameters affect the competition among gradient terms and the bulk energy terms in F η . Fig. 6(a)-(b) illustrate two different core sizes for the same dislocation obtained with different values for r and s. They show the reconstructed densities obtained by computing Equation (23) with the numerical solution for the amplitudes (first column), the Burgers vector density component B x (second column), a plot of B x (x, 0) and B x (0, y) (third column, empty symbols) with Gaussian fits (solid lines). The data fitting is obtained via G exp(−x 2 /2σ 2 x − y 2 /2σ 2 y ) 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 [56] or within strain-gradient elasticity theories [57,58].
The amplitude expansion defined in Equation (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 Equation (1) to the phase of the amplitudes one gets − q (n) · u = 2πq (n) · b. At least for perfect dislocations, those having a translation vector of the lattice as Burgers vector, we have that q (n) · 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. 6(c), the fields η n e iq (n) ·r entering the sum in Equation (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, , V the interface velocity along its normal and g a parameter in the free energy which multiplies gradient terms and scales the elastic constants [55,59,60] (see also Sec. IV). For a given set of parameters, we determine the specific amplitude pro-file 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 Sec. IV, Measuring the size of the dislocation core through σ x and σ y from a Gaussian fit as in Fig. 6(a)-(b), we find that it scales linearly with w when varying g, while a different scaling is observed when varying r , c.f. Fig. 6(d).
Here g is an energy scale associated with amplitudes gradients, similar to theories based on Ginzburg-Landau energy functionals [60]. r , 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. 6(c)). Also, for r = r 0 , interfaces move, which affects the width w [61]. 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. 6(d). σ y /σ x ≈ 0.75 throughout the whole range of parameters investigated here as also illustrated in Fig. 6(e). 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 phasefield 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 where R (1) bcc are the reciprocal lattice vectors of the bcc Bravais lattice with unit length [48]. This sets the length of the bcc unit cell as a 0 = 2π √ 2. We consider spherical inclusions with radius 17a 0 , rotated at an angle θ rot about the [1, 1, 1]-axis. The initial condition is relaxed by dissipative dynamics with an appropriate symmetryconserving free energy; see further details in Sec. IV. We choose three representative angles θ rot and calculate the Frobenius norm |α| = α ij α ij of α at each angle. Since |α| > 0, we plot its isosurface at half its maximum value |α| M = max r (|α|(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 nonlinear 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 [62,63]. 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 [64]. Such a regime shift echoes the onset of active turbulence in the nematic liquid, where the description in terms of the order parameter ρ (ηq) also breaks down.

III. DISCUSSION
In-depth understanding and tailoring of collective behaviors require a unified description of defects associ-ated 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 nonsingular 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 [1] 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 [28] 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 [65]. 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  Figure 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 |α| = α ij αij at 50% of its maximal value |α|M = maxr(|α|(r)), which is given in the panels.
for melting of 2D crystals via the hexatic phase [66].
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 [64], 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 solidliquid 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.

A. Bose-Einstein condesates
The damped Gross Pitaevskii equation, Equation (16), is solved by using a Fourier pseudo-spectral integration scheme which is described in detail in Ref. [31]. 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 ψ = Π 2 α=1 χ(|r − r α |)e iqαθα , where r α is the position of the vortex labeled α, θ α = arctan [(y − y α )/(x − x α )], and 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.

B. Active nematic liquid crystals
The evolution of the Q-tensor follows dissipative dynamics coupled with an incompressible Stokes flow [44] 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 = (∂ i v j − ∂ j v i ) rotates the nematic structure, λ is the flow alignment parameter which aligns the nematic orientation in the direction of shear with the trace less strain rate 2E ij = ( 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 = √ 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 = √ 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, Equation (36), in Fourier space and evolve the equation for the Q tensor, Equation (35), using the same scheme as for the BEC.

C. 2D square lattice PFC
To simulate the PFC dynamics, we use the sHPFC model proposed in Ref. [49], namely coupled to a momentum equation · is a convolution with a Gaussian kernel given by which filters out variations on length scales smaller than w. The quench depth in Equation (30) is set to r = −0.3 and the average density toψ = −0.3. Parameters are set to Γ = 1, ρ 0 = Γ S = 2 −6 , 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 0 /7. Further details can be found in Ref. [49].
In the simulation reported in Figure 5, the perfect lattice is indented by an applied external force density given by a Gaussian profile Above a critical strength f 0 = 3.5µ/a 0 and width w = a 0 , this force causes the nucleation of a dislocation dipole.
We simulate a stationary system hosting dislocations with the APFC model exploiting the (FEM) numerical approach with adaptive grid refinement outlined in Refs. [67,68]. The semi-implicit integration scheme adopted for numerical simulations can be found therein. We consider dislocations with spacing L = 50a 0 arranged in a periodic, 2D matrix with alternating Burgers vectors ±a 0x . The system is initialized by setting the displacement field of dislocation known from classical continuum mechanics [50] in the phase of amplitudes, −q (n) · u, and let relaxed according to the amplitudes evolution law (41). We can consider a system 2L × 2L by exploiting periodic boundary conditions. In Sec. II D, we characterize the extension of the core of dislocations through the field D (n) as entering the definition of the dislocation density tensor α, Equation (29). We compare the size of the defects extracted with the aid of Gaussian fits (see Figure 6(a)-(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 Equation (41) until reaching a steady state. Then, a fit of each amplitude with a function φ i =Ā i [1−tanh(x−x i )/w i ], representing a travelling wave solution for a solid-liquid interface [55,59,61], is performed withĀ i ,x i andw i parameters and the solidliquid interface thickness extracted as w = where F bcc ψ is a free energy functional that produces a stable bcc lattice, given by As parameters, we use r = −0.3 and ψ 0 = −0.325 with spatial discretization ∆x = ∆y = ∆z = a 0 /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 0 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 Equation (43) for 300 time steps with ∆t = 0.1.

DATA AVAILABILITY STATEMENT
The data will be provided by the corresponding author upon reasonable request.

CODE AVAILABILITY STATEMENT
The code will be provided by the corresponding author upon reasonable request.

FIGURE CAPTIONS
• Figure 1: Different types of excitations in a 2D vector field theory.
• Figure 3: Annihilation of a vortex dipole in a Bose-Einstein condensate.
• Figure 4: Onset of active turbulence in a nematic liquid crystal mediated by the nucleation of topological defects.
• Figure 5: Nucleation of a dislocation dipole in a square PFC model.
• Figure 6: Dislocation core size near melting by APFC modeling.
• Figure 7: Rotated inclusions in the bcc PFC model. Figure 8. Real 3D space and order parameter space for a D 3 order parameter. ∂M is the boundary of a 3D subvolume M, on which variations along the surface (dx µ 2 , dx µ 3 ) lead to variations of the order parameter (dΨ (2) , dΨ (3) ).

SUPPLEMENTARY NOTES
For the proofs in this section, we follow the notation conventions of Ref. [69]. We consider topological defects for S k order parameters, which is the space consisting of (k + 1)-dimensional unit vectors. It is known that the i-th homotopy group of S k is trivial for i < k. In particular, every loop (i = 1) on the two-sphere (k = 2) can be reduced to a point by a continuous deformation. Thus, to get a description of the topological defects for S k order parameters, we need to consider the k-th homotopy groups π k (S k ) Z corresponding to topological defects with integer charges. The dimension of the defect is given by d top = d − (k + 1), where d is the physical space dimension. The homotopy classification of loops in S k is useful beyond the direct application to models from this group because, in many systems, their order parameter space can be mapped or decomposed into products of S k spaces. We define an order parameter Ψ = (Ψ 1 , ..., Ψ n ) which resides in the order parameter space D n . The subvolume of D n swept by Ψ and n − 1 independent variations {dΨ (k) } n k=2 is given by the (signed) volume of the n-dimensional cone 1 n˜ ν1...νn Ψ ν1 dΨ (2) ν2 ...dΨ (n) νn , Here, ∂M is a (n − 1)-dimensional submanifold of R d , the boundary of some n-dimensional submanifold M, and {dΨ (k) } n k=2 are changes in Ψ due to displacements dx µ on ∂M. Formally, we write the integrand in terms of the coordinates {x µ } d µ=1 of R d as follows Since ω µ2...µn is completely anti-symmetric under interchange of indices, ω is a (n − 1)-form and we can apply Stokes' generalized theorem where dω is the exterior derivative of ω, whose components are given by where the notation [...] is the antisymmetrization over free indices (µ 1 ...µ n ). Thus, we get In the case of n = 1, for which V 1 = 2, integrals over M are typically evaluated in this way, i.e., as line integrals over the one-dimensional manifold M. We then immediately recover the defect density for n = 1, given by For higher values of n, however, one evaluates the integral in coordinates on the manifold M. Thus, at each point, we choose d − n orthogonal unit vectors {n (k) } d−n k=1 normal to the manifold M and introduce local coordinates on M. Expressed in these coordinates, we have where τ -indices iterate from 1 to n. Now, we invoke the identity Identifying N i1...i d−nˆ as the oriented volume element of M, we identify the defect density as where We now turn to finding the general equation for the velocity of the defect density. To differentiate Equation (57) where, in going from line 2 to 3, we have used that due to the contraction with both the anti-symmetric Levi-Civitas, the terms are invariant under simultaneously interchanging µ k ↔ µ k and ν k ↔ ν k , so that we can write every term in the parenthesis like the first. In going from line 4 to 5, we have used that the contraction with µ1...µn i1...i d−n ensures that only the first term survives when applying the product rule. Thus, we find ∂ t ρ i1..
Like for the d = n = 2 case, we want to identify this expression with the density current v µ1 ρ i1...i d−n . They are related up to a divergence free contribution In the d = n = 2 case, the charge density on the left-hand side had no free indices and so we could simply divide by the charge density to solve for v. In the general case, however, we project the equation by contracting with In order to fix the gauge K µ1 i1...i d−n , we look at the evolution of the order parameter Ψ as advected by a velocity field v (Ψ) These are n linearly independent linear equations to determine d components of the velocity v (Ψ) . If n < d, it is under-determined and therefore d − n additional equations are needed to determine v µ1 uniquely. We define ν n (∂ t Ψ ν1 )(∂ µ1 Ψ ν 1 ) n l=2 (∂ µ l Ψ ν l )(∂ µ l Ψ ν l ) δ [ν1 ν 1 ...δ νn] ν n n l=1 (∂ µ l Ψ ν l )(∂ µ l Ψ ν l ) (66) where it is understood that the repeated indices are summed over independently in the numerator and denominator. By inserting this expression in the LHS of Equation (63) (Venus) (67) We split the term (Venus) into ν 1 = k and ν 1 = k as follows (Mars) is identically zero, which the following argument shows: Because of the antisymmetrization over ν 1 ...ν n , in every term in (Mars) ν 1 , ..., ν n will take every value 1, ..., n. Since ν 1 = k, it means that there is some m > 1 such that ν m = k. Isolating the corresponding factor from the product, we get because the factor (∂ µ1 Ψ ν 1 )(∂ µ1 Ψ k )(∂ µm Ψ k )(∂ µm Ψ ν m ), this is symmetric under the interchange ν 1 ↔ ν m , but the Kronicker-delta product is antisymmetric. Now consider (Mercury). As before, in every term, ν 1 , ..., ν n will take every value 1, ..., n. Thus, in each term of (Mercury), there will be an m such that ν m = k, so we write Now, renaming ν m → ν 1 , ν 1 → ν 2 , ν 1 → ν 2 , and so on up to ν m−1 → ν m , and µ m → µ 1 , µ 1 → µ 2 , and so on up to µ m−1 → µ m , we get This, in turn, means (Mercury) = (Venus), which shows that v µ1 candidate is a solution to Equation (63). We have verified this calculation explicitly up to n = d = 5, using symbolic mathematical software. In addition, it is straight-forward to show that v µ1 candidate is orthogonal to ρ i1...i d−n in the sense that v i k candidate ρ i1...i k ...i d−n = 0 for all k. Identifying these as the (d − n) necessary conditions to determine v (Ψ) , we get v (Ψ) = v candidate , and fix the gauge on v by v = v (Ψ) . which gives, finally, the closed expression for the velocity v µ = −n δ ν n (∂ t Ψ ν1 )(∂ µ1 Ψ ν 1 ) n l=2 (∂ µ l Ψ ν l )(∂ µ l Ψ ν l ) δ While this derivation holds in general, we note that in the case of n = d, there is no contraction in getting to Equation (62), so the velocity can be equivalently written as Special case n = d : v µ = J µ ρ = −n µ1...µn˜ ν1...νn (∂ t Ψ ν1 ) n l=2 (∂ µ l Ψ ν l ) µ1...µn˜ ν1...νn n  Figure 9. Examples of reduced defect field corresponding to stable topological defects in O(n) models for n = 1, 2, 3 in d = 1, 2, 3.
While the manuscript has centered on systems where n = 2, the methodology can be applied to other cases, such as interfacial systems with n = 1 with wall defects [39] or systems with n = 3 such as the 3D Heisenberg model which features emergent magnetic monopoles as topological defects [40].