Detecting defect dynamics in relativistic field theories far from equilibrium using topological data analysis

We study nonequilibrium dynamics of relativistic $N$-component scalar field theories in Minkowski space-time in a classical-statistical regime, where typical occupation numbers of modes are much larger than unity. In this strongly correlated system far from equilibrium, the role of different phenomena such as nonlinear wave propagation and defect dynamics remains to be clarified. We employ persistent homology to infer topological features of the nonequilibrium many-body system for different numbers of field components $N$ via a hierarchy of cubical complexes. Specifically, we show that the persistent homology of local energy density fluctuations can give rise to signatures of self-similar scaling associated with topological defects, distinct from the scaling behaviour of nonlinear wave modes. This contributes to the systematic understanding of the role of topological defects for far-from-equilibrium time evolutions of nonlinear many-body systems.

Relativistic scalar fields with global O(N ) symmetry provide a paradigm for understanding far-fromequilibrium dynamics, which we focus on in this study.In particular, such a model allows for the analytic derivation of scaling exponents related to nonthermal fixed points up to anomalous dimensions based on the two-particleirreducible effective action in large-N expansions [1-4, 13, 20, 51, 54, 55].Moreover, at early times, for weak couplings and large occupations, the quantum dynamics of the model can be accurately described by means of classical-statistical simulations [51,56], which we employ.* v.noel@thphys.uni-heidelberg.de Recent experimental and theoretical studies of diverse systems indicate that for a small number of field components different initial conditions can lead to distinct nonthermal fixed points.They can be related to separate physical mechanisms: self-similar dynamics consistent with descriptions in the limit of many field components [44,47], and topological defects undergoing coarsening dynamics for few field components [14, 17, 24-26, 28, 31, 32, 35, 39, 49, 50, 57, 58].Also for the O(N ) vector model in three spatial dimensions, a number of different, N -dependent topological defects are expected to contribute to the dynamics in the infrared [19].Remarkably, for this model and typical over-occupied initial conditions, many types of equal-time correlation functions merely reveal scaling dynamics consistent with the large-N descriptions, even for low N , and in agreement with nonrelativistic complex U(N ) vector models [2-4, 19, 51, 54, 59, 60].This includes the momentumresolved distribution function, which is computed from equal-time two-point correlation functions and forms the basis for many studies of far-from-equilibrium universality.The investigation of unequal-time two-point correlation functions has revealed differences among the excitation spectra for different N [60].However, this does not contain clear links to the scaling behaviour related to the coarsening dynamics of topological defects.A thorough understanding of their role for nonthermal fixed point dynamics can be a prerequisite for the comprehensive classification of universal behaviour far from equilibrium, based on the experimentally and numerically found scaling dynamics.Hence, there is a need for topologically sensitive and computationally accessible observables.
Topological data analysis (TDA) [61,62] can provide complementary information to correlation functions.Persistent homology, which is the primary tool from the TDA toolbox, gives access to the topology of a nested hierarchy of complexes constructed from the field data.It has been successfully applied to describe nonthermal fixed points for two-dimensional nonrelativistic scalar and three-dimensional gluon fields [63,64], accompanied by a corresponding mathematical analysis [65].In the related context of phase transitions, persistent homology allowed for insights into critical phenomena and phase structures [66][67][68][69][70][71][72][73][74][75][76][77].
In this work, we consider the persistent homology of local energy density fluctuations in O(N ) vector models in the classical-statistical regime, starting from overoccupied initial conditions without imprinting defects.Local energy densities can show clear signatures of dynamically built-up defects as we demonstrate.In particular, we show that the related Betti number distributions give rise to N -dependent topological structures reminiscent of the classification of defects in condensates of relativistic O(N ) scalar fields.The time dependence of the related length scales indicates topological dynamics consistent with phase-ordering kinetics [78].
The remainder of this paper is structured as follows: Section II describes the lattice simulations, along with a description of topological defects in O(N ) theories and their presence in local energy densities.Section III introduces persistent homology and discusses the signals of topological defects in Betti number distributions.Moreover, the latter are interpreted in light of the coarsening dynamics of topological defects and energy transport.Finally, Section IV provides a conclusion.

II. O(N ) VECTOR MODEL AND TOPOLOGICAL DEFECTS
We consider a relativistic O(N )-symmetric scalar field theory with field variables ϕ a (t, x), a = 1, ..., N , in d = 3 spatial dimensions with classical action (1) where t,x ≡ dt d 3 x, summation over repeated indices is implied, m is the bare mass and λ is the coupling constant.The action (1) is invariant under global O(N ) rotations acting on the internal field components indexed by a: ϕ a → R ab ϕ b for R ∈ O(N ).

A. Lattice simulations
For highly occupied systems at not too late times and weak couplings, the quantum dynamics can be accurately mapped to a classical-statistical field theory (truncated Wigner approximation) [79][80][81][82][83][84].In classical-statistical simulations, one samples over initial conditions and each realisation is evolved according to the classical equation of motion.Expectation values of observables are obtained by averaging over their evaluations for the classical trajectories.Real-time classical-statistical simulations have been extensively used to study the dynamics of scalar fields within closed quantum systems in corresponding regimes of applicability [51,60,80,[84][85][86][87][88].
We consider over-occupied box initial conditions in momentum space for the scalar fields.More specifically, with the statistical two-point correlation function, defined as (2) for spatially homogeneous scenarios, occupation numbers can be defined as where F (t, t ′ , p) = ∆x F (t, t ′ , ∆x) exp(−ip∆x) and ∆x ≡ d 3 ∆x.Then field configurations ϕ a (t = 0, x) are sampled with large Gaussian fluctuations up to a characteristic momentum scale Q, described by with zero macroscopic fields |ϕ(t = 0, x)| = 0.The initial field configurations are time-evolved according to the equation of motion following from the action (1).For more details we refer to [51,56].We consider lattice simulations for N 3 s = 512 3 spatial lattice sites and use a leapfrog solver with Qa s = 0.8, where a s is the spatial lattice constant, temporal lattice spacing dt = 0.1 a s , coupling λ = 0.1, initial amplitude n 0 = 125 and renormalised mass squared M 2 = 4Q 2 , which is used in an iterative procedure to fix the bare mass m as described in [56].The comparably large mass suppresses fluctuations on smaller length scales, which will allow us to reveal the presence of topological features associated with defects.In persistent homology as introduced below, the topological features appear less pronounced for smaller values of the mass, see Appendix A. For later use we denote the spatial lattice by Λ s := {0, a s , . . ., (N s − 1)a s } 3 .
In this work we focus on the structures visible in local energy density fluctuations around the mean energy density.The local energy density corresponding to the action (1) is: (5) where π a = ∂ t ϕ a and space-time arguments have been suppressed on the right-hand side.Central spatial derivatives are employed on the lattice.
We consider a single classical-statistical realisation in this work, based on the self-averaging property often encountered for observables in classical-statistical simulations.The latter also holds approximately for the later introduced Betti numbers due to the large number of contributing features. 1In fact, the Betti numbers are expected to scale proportionally to the system volume for sufficiently large lattices, based on mathematical theorems [65,89].Specifically, we have verified insensitivity of Betti numbers after division by the lattice volume for N s = 128, 256, 512.We have also numerically verified their insensitivity to variations of the lattice spacing upon comparison with simulations for Qa s = 0.6 and Qa s = 1.2.
B. Topological defects in relativistic O(N ) theories

Topological defects in condensates
The chosen initial condition sets high occupation numbers up to the momentum scale Q.As discussed in previous works [15,51], the particle number redistributes towards the infrared via an inverse particle number cascade, which is a consequence of transient approximate particle number conservation.In this dynamical process, a condensate forms in the zero mode, which is initially absent and results from increasing occupancies in the deep infrared.Simultaneously, long-range order gradually builds up [19]. 2  Considering the ordering dynamics of the condensate, the phase space of spatial zero modes of the (Fourier-transformed) field variables φa (t, p) = x ϕ a (t, x) exp(−ipx), formed by ( φa , ∂ t φa ) a=1,...,N with φa ≡ φa (t, p = 0), is of relevance.Approximate particle number conservation and energy minimisation provide two constraints for the realised condensate configurations [19].Taking these into account, we denote the physically accessible condensate phase space by C N , which depends nontrivially on the number of field components N .In particular, the topology of C N can be nontrivial, so that topologically nontrivial configurations (defects) can occur.These have been classified in [19], which we review in detail in Appendix B. Specifically, the condensate can feature string defects (vortex lines) for N = 1, 2, 3, domain walls for N = 2 and monopoles for N = 4.For three spatial dimensions no other defects are expected, and condensates are defect-free for N ≥ 5.
Note that defects can be dynamically generated during the evolution of the (classical) equation of motion, and are not explicitly part of the initial conditions under consideration.As described later in this work, defects typically annihilate each other with time via related coarsening dynamics [78], such that their number decreases.On longer time scales than considered here, the condensate itself is expected to decay again due to number-changing 2 By causality, the formation of a condensate at finite evolution times requires a finite system volume.For infinite system volumes, modes in the deep infrared can mimic the dynamics of a condensate, but lack the related system-scale long-range order [19].
processes in the relativistic theory [15], suppressing topological defects as well.

Observation of defects in energy density fluctuations
We probe topological defects via their signatures in local energy density fluctuations around their mean values, given by ∆T 00 (t, x) := T 00 (t, x) − T 00 T 00 , where T 00 := (1/N3 s ) x∈Λs T 00 (t, x) is timeindependent due to energy conservation.In Figure 1(a) we display two-dimensional slices of the time-evolving energy density fluctuations for N = 1.As time elapses, the initial fluctuations in local energy densities become gradually more homogeneous, with string-like structures distinctively emerging at early times and minimal energy densities.We associate these with string defects, whose number appears to decline with time.
Comparing the two-dimensional snapshots of energy density fluctuations for N = 1, 2, 3, 4 at a fixed time (Qt = 5000), displayed in Figure 1(b), we observe similar string-like structures for N = 2 and N = 3.Moreover, for N = 2 we can also see indications of emerging domains separated by domain walls. 3These appear at locally low energy densities as closed loops in the snapshots and can be discriminated from the open elongated minima attributed to string-like structures.For N = 4 and higher N (not shown), the local energy densities become gradually more homogeneous without any distinct stringor domain-like structures forming, even when running the simulations for comparably long times up to Qt = 50000.Such defect structures agree with the classification outlined in Section II B 1, albeit with monopoles not seen for N = 4. Monopoles form point-like local minima in energy densities, similar to the many other fluctuations present in the system.The lattice configurations can therefore still contain monopoles, but might be indistinguishable from configurations without monopoles based on our methods.
The emergence of defect structures at minima in local energy density fluctuations can be heuristically understood based on energetic considerations.Defects locally minimise potential energy densities due to necessary zerocrossings in field amplitudes between oppositely-signed domains, at a string or within a monopole core, which manifest as zero local potential energy densities.Kinetic energy densities require more careful considerations.Defects move slowly in comparison to typical time scales associated with the inverse particle cascade, resulting in small contributions by a π 2 a to T 00 in Equation (5).Moreover, spatial gradients in directions tangent to the curves formed by local energy density minima due to string defects or tangent to the corresponding surfaces for domain walls are naturally small, as in these directions local energy densities do not change much (see Figure 1).In normal directions defects come with related characteristic (healing) length scales.These can be comparably large, leading to a suppression of normal gradients as well, such that kinetic energy densities of defects are locally suppressed along with potential energy densities.
We have observed this in our simulations for both kinetic and potential energy densities (not shown).

III. PERSISTENT HOMOLOGY OF ENERGY DENSITY FLUCTUATIONS
Persistent homology provides a method to calculate scale-dependent topological structures from data along with measures of their persistence.In Section III A we present the concept of persistent homology for the investigated local energy density fluctuations.In Section III B the persistent homology results are discussed in light of the previous discussion of topological defects in conden-sates.Section III C considers the time dependence of Betti number distributions in relation to coarsening dynamics, while Section III D provides an examination of its connection to energy transport towards the ultraviolet.
A. Persistent homology of energy density fluctuation sublevel sets

Cubical complexes
In this work, persistent homology is computed for cubical complexes [91], which is ideally suited for data in pixel format.Here, we focus on a short, intuitive introduction, while Appendix C provides mathematically more concise constructions.For more elaborate mathematical introductions to TDA we refer to [61,62].
A cubical complex is a collection of cubes of different dimensions, closed under taking boundaries.For instance, the boundary of a 3-cube is the union of all its six faces, i.e., boundary squares, the boundary of a 2cube (square) is the union of its four boundary edges, the boundary of a 1-cube (edge) consists of its two endpoints, and the boundary of a 0-cube (point) is empty.The full cubical complex for the spatial lattice under consideration consists of a 3-cube for each lattice site, with a lattice site located at each cube's center.For each 3-cube all boundary cubes of lower dimensions are included in the full cubical complex.Subcomplexes of the full cubical complex can be used to describe sublevel sets of functions on the spatial lattice.Intuitively, for a given function these are constructed by including a 3-cube in the subcomplex whenever the corresponding function value is below a chosen filtration parameter.While this fixes the filtration parameters when 3-cubes enter subcomplexes, the filtration parameters for lower-dimensional cubes are set inductively (the lower star filtration, see Appendix C).This way the subcomplex becomes indeed a cubical complex for every filtration parameter, which can be seen as a pixelisation of the corresponding lattice function sublevel set.In particular, a nested sequence of cubical complexes is formed upon increasing the filtration parameter.This procedure is illustrated in Figure 2 for the example of an excerpt of the two-dimensional slice of the local energy density fluctuations shown in Figure 1(b) (N = 1).While the two-dimensional slice is for illustrative purposes only, the analysis is carried out in all three spatial dimensions.For a given function on the lattice such as the image excerpt given in Figure 2(a), the cubical complexes corresponding to the sequence of sublevel sets can resemble the structure of minima, see Figure 2(b).In particular, note the persistence of the horseshoe-like accumulation of squares across filtration parameters.
More formally, the evaluation of local energy density fluctuations for a given classical-statistical realisation ϕ(t, x) at time t provides a real-valued function ∆T 00 (t, •) on the spatial lattice Λ s .Its sublevel sets are defined as The described pixelisation procedure leads to the cubical complexes C ∆T 00 (t, ν), which are the cubical complexes of interest in this work.They form a filtration of the full cubical complex, i.e., a nested sequence of subcomplexes of the full cubical complex with for all ν ≤ µ.For ν smaller than the minimum value of ∆T 00 (t, •) the cubical complex C ∆T 00 (t, ν) is empty, for ν larger or equal the corresponding maximum value the full cubical complex is recovered.

Persistent homology: holes in complexes
The full cubical complex of the three-dimensional lattice contains a 3-cube for each spatial lattice point.However, as energy density fluctuation values are swept through from the lowest to the highest, in general, the cubical complexes C ∆T 00 (t, ν) do not contain a cube for every spatial lattice point.Holes of different dimensions can appear, which are described by homology groups.For cubical complexes, such holes are illustrated in Figure 3(a).Connected components are described by dimension-0 homology classes, dimension-1 homology classes describe planar-like holes (which in three dimensions can also be viewed as tunnels), and dimension-2 homology classes describe enclosed volumes.
Within the filtration, holes are born at a birth parameter b and die again with death parameter d, possibly deforming as the (energy density fluctuation) filtration is swept through, i.e., being present in the filtration for the filtration parameter interval [b, d).The persistence p = d−b is a measure of the dominance of a feature in the filtration.This is illustrated in Figure 3(b), where two landscapes of example functions with distinct minima are displayed.The left function contains two distinct dips.As the sublevel set filtration is swept through, the cubical complex is empty as long as the filtration parameter ν is less than the minimum value of the function.As ν is increased, a dimension-0 homology class is born with birth parameter b = ν when the minimum value of the function is reached (green plane).Further increasing ν (blue plane), the single 2-cube turns into a set of 2-cubes, leaving this homology class unchanged.When the value of ν reaches the pink plane, a second dimension-0 homology class is born (on the right) corresponding to the second dip in the function.The two dimension-0 homology classes merge into one at the red plane, at which point the first homology class dies with death parameter d = ν, and only the second one survives.
Turning to the right function in Figure 3(b), a few dimension-0 homology classes are born (green plane), merging to form a dimension-1 homology class when ν is increased to the level of the blue plane, represented by the circular structure in the complex surrounding a hole.Further increasing ν, the homology is unchanged at the pink plane and the dimension-1 hole only disappears somewhere between the pink and red planes when it is fully filled by 2-cubes.
From the cubical complexes C ∆T 00 (t, ν), the different dimension-ℓ homology groups H ℓ (C ∆T 00 (t, ν)) can be computed.For three spatial dimensions, the homology groups are generally nontrivial for ℓ = 0, 1, 2, while the dimension-3 homology group only captures the toroidal lattice topology itself; all higher homology groups are trivial.Their dimensions, called Betti numbers, specify the number of independent dimension-ℓ holes (homology classes): We focus on Betti numbers as an informative persistent homology descriptor in this work.Persistent homology has a number of useful features.Notably, it is stable against small perturbations of the input [61,62] (local energy density values in this case) and well-defined large-volume asymptotics exist for suitable persistent homology descriptors such as Betti numbers, including notions of ergodicity [65,89].We compute the persistent homology of cubical complexes with Z 2 coefficients and periodic boundary conditions using the open source TDA library GUDHI [92].
To summarise, the persistent homology of sublevel sets can be particularly sensitive to the extended structures formed by local minima and saddle points.We utilise this to investigate the structure of minima in energy density fluctuations, for which nontrivial defect dynamics is expected, based on the discussion in Section II B 2.

Blockwise averaging lattice functions
By construction, the persistent homology of cubical complexes does not contain spatial metric information.However, this also means that persistent homology computed from functions on the lattice does not differentiate between different length scales and therefore lattice artefacts can enter the analysis.By blockwise averaging (coarsening) the local energy densities on the lattice over blocks of c 3 points for c > 1, both these artefacts and, partly, ultraviolet fluctuations are removed.This can pronounce contributions related to the infrared dynamics including topological defects, depending on the parameter c.We emphasise that the blockwise averaging procedure does not affect the dynamical evolution, since it is only applied a posteriori to the lattice configurations as part of the persistent homology read-out.This is demonstrated in Figure 4 for the case of N = 1 with dimension-0 Betti numbers.The number of persistent homology classes typically decreases as we average blockwise, since the number of lattice sites decreases accordingly.Betti numbers are expected to scale proportionally to system volume for sufficiently large lattices [65,89], so that we can account for this by volumenormalising to a certain number of lattice sites, here chosen to be 64 3 . 4Without blockwise averaging, we see that the number of connected components specified by dimension-0 Betti numbers increases at earlier times.Averaging local energy densities over every 4 3 lattice blocks, we notice a twofold peak structure emerging: connected components accounting for the left peak decrease in numbers, while the right peak grows with time.Increasing to averaging over every 8 3 blocks, the left peak appears more pronounced, while the right peak decreases in height.This behaviour persists when coarsening by a factor of 16 in total in each direction.
Concerning the physical interpretation, we notice that the initially large occupations give rise to dynamics towards the infrared and growing wavelengths of correspondingly dominating modes as can be seen in the occupation number distributions, which are comparable to [51] (not shown here).The height of the left peak in dimension-0 Betti numbers decreases over time for energy densities averaged over blocks of 4 3 , 8 3 or 16 3 nearby lattice sites, which implies an increase in the average distance between the related connected components.This indicates that dimension-0 Betti numbers can probe the infrared dynamics for blockwise averaged configurations.
In particular, the coarsening dynamics of topological defects merging is typically also accompanied by a growth of respective characteristic length scales [78].In the following, we show results for local energy densities blockwise averaged over cubes of 8 3 neighbouring lattice points.When quantitatively discussing the time dependence of Betti numbers below, it is investigated if the observed phenomena are stable with respect to further blockwise averaging.

B. Topological defects in Betti numbers
As shown in Figure 5(a), upon investigating the timeevolving dimension-0 Betti numbers of local energy density fluctuations from N = 1 to N = 4 with averaging over blocks of 8 3 lattice sites, we observe two distinct peaks appearing for N = 1, 2, 3, while for N = 4 and higher N (not shown) there is only a single peak.Starting from the initial conditions (displayed in grey), for all N , the single peak in dimension-0 Betti numbers first moves to larger filtration parameters at early times, i.e., local minima shift to larger energy density values.For times larger than Qt = 500, the peak position stays approximately constant for N = 1, 2, 3, while the overall number of dimension-0 structures decreases with time.For N = 1, 2, 3, the single peak splits up into two distinct peaks, the left one decreasing, the right one increasing in height.For N = 4, no split-up happens, and the peak first decreases, then increases in height, and continuously shifts to larger filtration parameters.
Since dimension-0 Betti numbers count independent connected components, a decline in peak height implies that average length scales associated with the component configurations increase, which correspond to local energy density minima in the field configurations.Heuristically, we may thus associate the left peak and its decrease as a signature of dynamics towards larger length scales.Notably, for N = 1, 2, 3, the connected components making up the left peak appear for those ∆T 00 -values where we have observed defects in the two-dimensional snapshots of ∆T 00 in Figure 1(b).This indicates that in terms of specific field configurations it is primarily the defects which dominate the declining left peaks in Betti numbers.Accordingly, the decrease in peak height can signal their coarsening dynamics, which we discuss in more detail later, see Section III C. Vice versa, Betti numbers increasing with time as for the right peak for N = 1, 2, 3, and for N = 4 for later times Qt ≳ 4000 implies refinement dynamics of the connected components, i.e., an increasing number of local minima appears.Corresponding length scales shrink in time.The peak shifting towards ∆T 00 = 0, this indicates an ongoing homogenisation process of energy densities and can be suggestive of the transport of energy towards towards smaller length scales (larger momenta), as we analyse in more detail in Section III D.
Likewise, the analysis of dimension-1 Betti numbers shows qualitatively similar results for all N , but with only a single peak structure that decreases in amplitude, except for N = 2, for which two peaks appear.This is displayed in Figure 5(b).More precisely, we notice a height decline of the peak for N = 1, 3, 4, which turns into an increase at later times.For N = 2, the behaviour is different: the height of the initial peak only decreases, and the single peak splits up into two peaks, where the right one forms at larger energy density values corresponding to ∆T 00 ≃ 0 and increases in height with time.The dimension-2 Betti numbers agree qualitatively for all N and do not come with twofold peak structures (not shown).
Still, a decline in dimension-1 Betti numbers indicates that length scales associated with the configurations of holes increase in time as for the dynamics towards the infrared.In particular, this applies to the left peak for N = 2, which constantly decreases in height and appears at energy densities, for which we have inferred the presence of defects from Figure 1(b).Similarly to the previous discussion of dimension-0 Betti numbers, at later times, the increase in height for most peaks in dimension-1 Betti numbers together with their shift towards ∆T 00 = 0 can be indicative of refinement dynamics leading to a homogenisation of local energy densities.
To summarise, we observe potentially defect-related peaks in the dimension-0 Betti numbers shown in Figure 5(a) for N = 1, 2, 3 but not for N = 4, and only for N = 2 in the dimension-1 Betti numbers displayed in Figure 5(b).This is in line with the classification of topological defects for condensates outlined in Section II B 1 and reviewed in detail in Appendix B, provided that the left peak in dimension-0 Betti numbers is primarily due to strings and the left peak in dimension-1 Betti numbers is mostly due to domain walls.Indeed, this appears well-motivated from the structure of the defects along with fluctuations on top as inferred from the local energy density snapshots displayed in Figure 1.For instance, closed strings manifest as loop-like minima in local energy densities and would thus naively appear as distinct dimension-1 persistent homology classes in the sublevel sets for correspondingly low filtration parameters (cf.also the discussion at the end of Section II B 2).Yet, smaller-scale fluctuations on top result in a landscape of local minima and maxima, which adds to this and may interrupt the clean loop-like minima in energy densities.Accordingly, a closed string would not appear as a distinct dimension-1 homological feature anymore, but as a multitude of dimension-0 structures, which in addition cannot be distinguished from open strings.Similarly, domain walls would naively show up as empty volumes in energy density sublevel sets for correspondingly small filtration parameters, which would give rise to distinct dimension-2 persistent homology classes.The addition of smaller-scale fluctuations can yield local energy density maxima, which then pierce the enclosed volumes in the sublevel sets, yielding only scaffold-like networks which surround the empty volumes.These do not give rise to distinct dimension-2 persistent homology classes anymore, but to an abundance of dimension-1 features.
The same reasoning applies to more exotic configurations of topological defects such as strings pinned to domain walls [26], which our persistent homology analysis cannot distinguish from the domain walls themselves.

C. Signatures of coarsening dynamics
On the time scales under consideration, the overoccupied initial conditions lead to nonthermal fixed point dynamics as characterised by dynamical self-similar scaling.The temporal dependence of the distribution function is restricted to spatial rescalings by time-dependent power laws while maintaining its shape in time [56].Near a nonthermal fixed point, the Betti number distributions of energy density sublevel sets can also reveal self-similarity [64].In this work, the shape of the Betti number distributions shown in Figure 5 is not globally preserved in time, in particular for dimension 0 at N = 1, 2, 3, and dimension 1 at N = 2. Yet, for these, the shape of the peaks at lower ∆T 00 -values remains approximately constant in time.We discuss this in detail in Appendix D considering potential self-similar scaling.
For clarity, here we focus on the time-dependence of the corresponding dimension-0 Betti number peak values, β topo 0,max (t) = max{β 0 (t, ν) | ν ∈ left peak}.This counts the maximum number of connected components formed by the pixelised sublevel sets as the filtration parameters corresponding to the peak are swept through, correlating with the number of defects.If a power law in time with exponent −ϑ N can be identified from β topo 0,max (t), i.e., then in three dimensional space, average length scales associated with the connected component configurations at the respective value of ν ≡ ∆T 00 grow as a power law with exponent ϑ N /3: Such a power law growth of length scales can be indicative of the self-similar scaling of corresponding Betti number distributions.Since for N = 1, 2, 3, the lower-∆T 00 peak is defect-related as noted in Section III B, this length scale can be sensitive to the dynamics of defects.
In particular, the number of strings correlating with the number of connected components of the ∆T 00 sublevel sets (cf.Section II B 2), L N (t) can serve as a proxy for the time-dependence of average length scales associated with string configurations.
In Figure 6, we display the time-dependence of β component configurations grow, as discussed previously (cf.Section III B).The decrease of β topo 0,max (t) is of a power law form for long time ranges.Fitting a power law to these curves for the indicated times 5 via standard χ 2 fits yields ϑ 1 = 0.69±0.02(N = 1), ϑ 2 = 0.64±0.02(N = 2) and ϑ 3 = 1.19 ± 0.01 (N = 3) when averaging over every 8 3 blocks.We obtain ϑ 1 = 0.67 ± 0.04, ϑ 2 = 0.61 ± 0.05 and ϑ 3 = 1.23 ± 0.03 when averaging over 16 3 blocks.Thus, within errors these exponents are insensitive to blockwise averaging in the considered regime.They describe dynamics predominantly related to string defects, which appears well-separated from the dynamics on much smaller length scales based on the described insensitivity to blockwise averaging.This is in contrast to the analogous analysis for the potentially defect-related peak at N = 2 in dimension-1 Betti number distributions, for which the power law scaling behaviour is different when averaging over blocks of 8 3 versus 16 3 lattice points.For this reason, a more detailed analysis of the defect-related peak in N = 2 dimension-1 Betti number distributions is only described in Appendix E.
The fitted exponents ϑ N correspond to the following power laws describing the time dependence of the length scales L N (t) (for averaging over every 8 3 blocks): It is instructive to compare this with the known dynamics of topological defects as captured by phaseordering kinetics [78], which describes the growth of order through coarsening dynamics when a system is quenched across a phase transition.Notably, phase-ordering kinetics captures the mutual annihilation of vortices and anti-vortices (strings) as well as the shrinking of domain walls with time.These processes can give rise to self-similar scaling which is characteristic of a nonthermal fixed point, with defect-related length scales displaying a power law in time.For instance, numerical studies of universal dynamics in one-dimensional Bose gases have revealed power law exponents of order 0.25 to 0.35 [32,39,93].Recent experiments with ultracold atoms in a quasi-one-dimensional elongated optical trap have pointed towards a growth of length scales related to vector solitons as a power law in time with exponent 0.28 ± 0.05 [49].In two spatial dimensions it is known that inter-vortex distances can grow proportional to t 0.2 for nonrelativistic systems [28,58,63,94], and similarly for a nonrelativistic projection of relativistic scalar theory [31].For three spatial dimensions, signatures of defects in scalar field theories have been detected previously, albeit without analysing their self-similar scaling 5 Note the different time intervals for 8 3 versus 16 3 blocks.This is due to initially slower dynamics in Betti number distributions for averaging over 16 3 blocks compared to 8 3 blocks, since part of the faster ultraviolet modes do not contribute to the former (cf. Figure 4).dynamics [14,24,57,90].Yet, even in this case, phaseordering kinetics predicts scaling exponents for defectrelated length scales in the range of 0.2 to 0.3, if conserved order parameters are considered [78].Our results for L 1 (t) and L 2 (t) as deduced with persistent homology are well within this range.The scaling exponent of L 3 (t) is closer to 0.5, which is considered to be the relevant exponent for length scale dynamics associated with the particle cascade in large-N expansions (up to anomalous dimensions) [51,59,60], where topological defects are absent.Still, topological defects can also give rise to such dynamics [4,78].Note that the association of persistent homology classes with defects is not a one-to-one correspondence, and not all persistent homology classes need to behave uniformly according to (12).Therefore, for our method we expect systematic errors on the deduced scaling dynamics of topological defects, which we do not consider here and are to be discussed in a future work.

D. Signatures of energy transport
In addition to the coarsening dynamics, there is an ongoing homogenisation process of local energy densities.This is apparent from the peaks with increasing amplitudes close to ∆T 00 = 0 in the Betti number distributions, which shifts towards the mean local energy density.Again, we study their maxima, in this case for dimension-1 Betti number distributions, denoted by In Figure 7, we show β energy 1,max (t) for N = 2 when averaging over every 8 3 blocks (main figure) and 16 3 blocks (bottom inset).Comparable outcomes can be obtained for other dimensions and N .We find that the number of features as counted by β energy 1,max (t) in Figure 7 grows in time, such that corresponding dimension-1 holes (tunnels) steadily decrease in size.Though slightly bent, the data can be approximately described by a power law, β energy 1,max (t) ∼ t −ϑenergy .Standard χ 2 fits yield ϑ energy = −0.54± 0.02 for averaging over every 8 3 blocks and ϑ energy = −0.53± 0.03 for 16 3 blocks, such that ϑ energy appears stable for this regime of blockwise averaging.
As before, we can estimate the time dependence of average length scales associated with the dimension-1 holes from the time dependence of β energy 1,max (t): which leads to L energy (t) ∼ t −0.18 for averaging over every 8 3 blocks.It is well known that turbulent energy transport towards the ultraviolet is accompanied by selfsimilar scaling behaviour characteristic of a nonthermal fixed point, for which related length scales dynamically shrink as ∼ t −1/7 [40,41].The power law behaviour of L energy (t) is close to this, which suggests that the corresponding structures in Betti number distributions are due to excitations in local energy densities transported towards smaller length scales.Deviations can be caused by the blockwise averaging procedure and the tentative interpretation of features in the Betti number distributions with different physical phenomena, which are in general not perfectly separated.Moreover, while the maxima of the Betti number distributions can provide estimates for the overall number of structures associated with the corresponding peaks, this number partly remains ambiguous in light of non-uniform persistences of features in the filtration.

IV. CONCLUSIONS
In this work, we have investigated the real-time dynamics of relativistic O(N ) scalar fields in over-occupied scenarios.We have applied TDA to lattice configurations, which can provide complementary information to the traditionally investigated distribution functions computed from equal-time two-point correlations.More specifically, we have considered Betti numbers computed for sublevel sets of local energy density fluctuations.These have revealed clear signals of dynamically generated topological defects.The identification of defects in the Betti numbers is based on the comparison with defect structures visible in local energy density landscapes.Crucially, the N -dependent topological features visible in the Betti numbers have been consistent with the clas-sification of defects for condensates in relativistic O(N ) scalar fields [19].
The number of connected components associated with the topological defects has decreased in time as a power law, which is indicative of self-similar scaling dynamics.This behaviour corresponds to power law growth of length scales associated with their dilution, with scaling exponents ∼ 0.2 for N = 1 and N = 2, and ∼ 0.4 for N = 3.While the values for N = 1, 2 agree well with the findings for topological dynamics in a variety of simulations, the value for N = 3 is closer to the scaling behaviour of the model at large N .Yet, topological dynamics can also lead to the latter scaling dynamics [4,78].A thorough quantitative analysis of the relation to phaseordering kinetics and a more careful study of possible dynamical contributions to the structures associated with defects in Betti numbers is beyond the scope of this paper.
In addition, we have observed signatures of energy transport and the related universal scaling behaviour in the Betti numbers.While for over-occupied scenarios the distribution function reveals a dual cascade with particle and energy transport, the Betti numbers computed for local energy density fluctuations may therefore be able to distinguish topological dynamics from energy transport.
Our work hints at the importance of topologically sensitive observables in order to uncover defect dynamics in universal regimes far from equilibrium.In particular for three-dimensional systems, defect-driven temporal scaling dynamics can be hard to access with distribution functions and appears to be suppressed in over-occupied scenarios.Based on the present work, TDA-based observables are particularly promising to investigate along with distribution functions, since they can be sensitive to extended field configurations of arbitrary size.This method is applicable without major modifications to the analysis of ultracold atom experiments.Local atom densities can play a similar role to the local energy densities considered in this work, with topological configurations also showing up as distinct minima.These are accessible for instance using commonly employed absorption imaging techniques.In this spirit, a recent work has exploited TDA to detect dark solitons in condensate density images [95].It is particularly beneficial that this method does not rely on more sophisticated experimental techniques such as response measurements to probe correlation functions at unequal times [23,59,60,[96][97][98].
Our work reinforces the potential of TDA to provide relevant information on the dynamics of strongly correlated many-body systems.The choice of cubical complexes has especially been advantageous to reveal the presence of non-local structures.TDA can facilitate the systematic study of the role of topological defects for nonequilibrium quantum many-body dynamics, with diverse regimes of applicability ranging from cold atoms to the collisions of heavy nuclei.
The simulations reported about in the main text have used a comparably large renormalised mass with M = 2Q.In this appendix, we discuss the dependence of the persistent homology results provided in Section III B on this choice.
Figure 8 provides the dimension-0 Betti numbers for simulations with N = 2 field components and renormalised mass M = Q/2.Comparison with Figure 5, which has given the corresponding results for M = 2Q, reveals that the right peak close to ν ≡ ∆T 00 = 0 increases in height for the smaller mass value, while the height of the left peak appears roughly insensitive to the choice of mass.In particular, the scaling behaviour of β topo 0,max (t) remains the same compared to the larger mass (not shown).
We can heuristically understand this behaviour as follows.For the results displayed in Figure 5 and Figure 8, we have employed averaging over blocks of 8 3 lattice sites.Removing ultraviolet fluctuations, this emphasises structures in the infrared, where the model is well-described by a nonrelativistic complex scalar field theory [31,51].The Hamiltonian of such a theory comes with a kinetic term ∼ k 2 /2M , where k denotes the spatial momentum of a Fourier mode of the complex-valued field.A larger value for M therefore suppresses the contributions of spatial gradients of local fluctuations to energy densities, which pronounces structures with small spatial gradients, for instance topological defects.
Moreover, we have argued in favour of the association of the right peak in the dimension-0 Betti numbers with energy transport towards the ultraviolet, see Section III D. In the ∆T 00 snapshots presented in Figure 1, this corresponds to the local fluctuations on small length scales.A larger mass therefore suppresses these, such that the topological defects appear pronounced relative to the small-scale fluctuations, even though their number might remain roughly invariant under the choice of mass.This explains the behaviour of the dimension-0 FIG. 8. Dimension-0 Betti numbers for ∆T 00 sublevel sets for N = 2 field components, with averaging over blocks of 8 3 lattice sites.The renormalised mass used to generate these results is Betti numbers shown in Figure 8.
For completeness, we note that the occupation number distributions computed for both mass values agree (not shown) and are similar to the results of [51].In this appendix, we derive the classification of topological defects in condensates of the O(N ) vector model in detail.The arguments are drawn from [19] and expanded here.
Intuitively, topological defects form when growing occupancies in the infrared organise into condensates along with the inverse particle cascade.To this end, the behaviour of spatial zero modes of the field variables is of relevance for the classification of topological defects.In particular, it is the topology of the phase space of zero modes of the field variables, ( φa , ∂ t φa ) a=1,...,N , where φa ≡ φa (t) ≡ φa (t, p = 0), along with dynamical constraints which determines the latter.Denoting this space for the O(N ) vector model by C N , we first motivate its general construction from first principles based on the condensate dynamics.We subsequently describe the specific structure of the C N and their topological properties, ultimately leading to the classification of defects in three spatial dimensions.

Dynamically realised condensate phase space
In our simulations, the modes in the infrared are highly occupied at early times for the considered initial condi-tions and are thus well-described classically.Here, we derive an approximate equation of motion for the condensate.First, the classical inverse propagator follows from the action (1) as with □ ≡ ∂ µ ∂ µ and where we explicitly denoted the summation over the field components.This leads to the classical self-energy matrix The action (1) yields an inhomogeneous Klein-Gordon equation as the classical equation of motion for spatial zero modes of the fields: where p ≡ d 3 p/(2π) 3 and the dependence on the classical self-energy matrix Σcl ab (t, −p) has been explicitly denoted.We assume a large fraction of the particles occupies the zero mode, so that the momentum integral in Equation (B3) is dominated by momentum zero.This finally yields the equation which governs the time evolution of the condensate.For a single field component (N = 1) φ ≡ φa=1 and the effective mass squared is M 2 = m 2 + Σcl 11 .Equation (B3) leads to the decomposition φ = φ0 cos(M t + δ), where φ0 is a real-valued peak amplitude of the field, δ a phase and M determines the oscillation frequency of the field.The amplitude | φ0 | is fixed by transient particle number conservation. 6From a topological viewpoint, the N = 1 condensate phase space C 1 is thus described by a circle, For N ≥ 2 field components the solutions φ ≡ ( φ) a=1,...,N to Equation (B4) generally have the form of oscillations sweeping out ellipses in the internal space of field components.These can interpolate between the extreme possibilities of φ and ∂ t φ parallel, so that oscillations happen along a line in this space, or φ and ∂ t φ orthogonal, so that oscillations happen along a circular orbit.It can be shown, that the energy-minimising condensate configurations are not straight lines and, therefore, have the topology of circular orbits.Indeed, repeating the argument provided in [19], we can understand this by comparing the energy cost of a circular orbit to that of a straight line in O(N ) field space.The energy density ε of the zero mode of a single-component scalar field oscillating back and forth in a quartic potential is where we have neglected the mass term.This equation can be rearranged for ∂ t φ to obtain a periodic solution, which has oscillation frequency The particle number stored in the condensate is n(ε) = ε 0 dε ′ /ω = 4ε/(3ω), and using Equation (B5), we obtain On the other hand, for a scalar field with two or more components with a circular orbital motion in the same potential, the energy density is B8) where we have used the Virial theorem (∂ t φ) 2 /2 = λ φ4 0 /4 and φ0 = φ(t 0 ) denotes the condensate configuration at some initial time t 0 .Hence, Using the same arguments as before for n(ε), we get that ε = 3 2 7/3 λ 1/3 n 4/3 .
(B10) Therefore, comparing the numerical prefactors, at fixed particle number the energy cost of a circular orbit is lower than that of a straight line.We can conclude that energyminimising orbits φ(t) are not straight lines and thus homotopy-equivalent to circles in the zero mode phase space.We expect that the addition of a non-zero mass term does not alter the topology of energy-minimising orbits.This is a posteriori reinforced by the consistency of our work with the defect classification outlined here, which relies on orbits being homotopy-equivalent to circles.
We thus proceed with the consideration of circular orbits in the zero mode phase space.For these the constancy of | φ| due to effective particle number conservation is complemented by |∂ t φ| being constant due to the circular orbit geometry, along with the orthogonality constraint φ ⊥ ∂ t φ in the internal field space.In general, a fixed-length φ is an element of the (N − 1)sphere S N −1 ⊂ R N .Any tangent vector to S N −1 at a point φ ∈ S N −1 is orthogonal to the vector φ itself within R N .The vectors ∂ t φ orthogonal to φ actually form the tangent manifold T φS N −1 .The constancy constraint for |∂ t φ| singles out tangent vectors of constant length.Since topology does not discriminate between the length of such tangent vectors, the constancy constraint for |∂ t φ| can be taken care of upon restricting to normalised tangent vectors ∂ t φ with |∂ t φ| = 1. 7Finally, for N ≥ 2, the condensate phase space C N is (homotopy-equivalent to) the unit tangent bundle of S N −1 , C N ≃ T 1 S N −1 , which as a set reads

Homotopy groups of CN
The topology of C N can be nontrivial as quantified by the low-order homotopy groups, which is the mathematical origin of topological defects.Specifically, a nontrivial zeroth homotopy group π 0 (C N ) indicates that C N comprises different connected components, i.e., the condensate is separated into different domains bounded by domain walls.If the fundamental group π 1 (C N ) is nontrivial, C N is not simply connected, and string defects can occur, i.e., vortex lines.If the second homotopy group π 2 (C N ) is nontrivial, then monopole defects can occur.A nontrivial third homotopy group π 3 (C N ) would indicate the possibility of textures.However, these have no defect core and are unstable [99], such that we exclude them from our analysis.Apart from textures, strings, domain walls and monopoles are all types of defects which can occur in three spatial dimensions.Investigating the presence of topological defects thus requires the computation of the homotopy groups π ℓ (C N ) for ℓ = 0, 1, 2 and different N .For N = 1, the condensate phase space is C 1 ≃ S 1 , for which π 1 (S 1 ) ∼ = Z and all other homotopy groups vanish, such that only string defects can occur.

b. N = 2
For N = 2, we have to consider C 2 ≃ T 1 S 1 .The unit tangent line at any point in S 1 consists of exactly two 7 For the same reason we here ignored the mass dimension of |∂t φ|.

c. N = 3
For N = 3, we have to consider C 3 ≃ T 1 S 2 .This space can be identified with SO(3), since φ describes a direction in R 3 and ∂ t φ is orthogonal to it.Together with their cross product, they single out an orthonormal coordinate frame.The fundamental group is nontrivial: 3)) ∼ = Z 2 and all other homotopy groups vanish.There are thus string defects.
e. N ≥ 5 For general N ≥ 5 one needs to consider C N ≃ T 1 (S N −1 ) ∼ = Spin(N )/Spin(N − 2).The identification with the spin group quotient comes about since the elements of the unit tangent bundle of S N −1 single out orthonormal 2-frames in R N , which form the second Stiefel manifold V 2 (R N ).Specifically, for N ≥ 3 we have that V 2 (R N ) is diffeomorphic to SO(N )/SO(N − 2) [100].This in turn is diffeomorphic to the quotient Spin(N )/Spin(N − 2) for all N ≥ 3, as can be seen with standard results for homogeneous spaces together with the consideration of relevant stabiliser subgroups of SO(N ) and Spin(N ).
The computation of the homotopy groups of this space requires more advanced methods from algebraic topology [100], which we only briefly outline here.Homotopy groups of the quotients Spin(N )/Spin(N −2) can be computed with the long exact sequence of relative homotopy groups, which reduces their computation to the homotopy groups of the spin groups themselves.The spin groups are the universal covers of the special orthogonal groups for N ≥ 3, as indicated by the related short exact sequence of groups: The spin groups are thus simply connected, so that π 0 (Spin(N )) ∼ = π 1 (Spin(N )) ∼ = 0 for all N ≥ 3.Moreover, the short exact sequence (B12) induces a long exact sequence of homotopy groups.Using π 2 (SO(N )) ∼ = 0 for all N ≥ 2 (in the stable range), we then find π 2 (Spin(N )) ∼ = 0 for all N ≥ 2. Together with the previously mentioned long exact sequence of relative homotopy groups this yields π 0 (C N ) ∼ = 0, π 1 (C N ) ∼ = 0 and π 2 (C N ) ∼ = 0 for all N ≥ 5, while higher homotopy groups can be nontrivial.Yet, it is these homotopy groups which determine the topological defects for three spatial dimensions.Topological defects are thus absent in the considered three-dimensional O(N ) vector model for all N ≥ 5.
Appendix C: Details on the lower star filtration of cubical complexes and persistent homology This appendix provides mathematical details on the construction of the lower star filtration of cubical complexes, homology groups and persistent homology.It closely follows similar expositions in [64].

Lower star filtration
We introduce the lower star filtration for a real-valued lattice function f : Λ s → R such as ∆T 00 (t, •), intuitively corresponding to a pixelisation of its lattice sublevel sets.To begin with, let C denote the full cubical complex of the lattice Λ s , consisting of one 3-cube x+[−a s /2, a s /2] 3 for each spatial lattice point x ∈ Λ s .C also includes all faces, edges and vertices of every 3-cube, such that it is closed under taking boundaries.C is equipped with the information contained in the function f by means of inductively constructing a map F : C → R. By construction, any 3-cube C ∈ C has a unique lattice point x ∈ Λ s at its center, so that we can set F (C) := f (x).For all 2-cubes D ∈ C we set Equation (C1) is then applied inductively to construct F for lower-dimensional cubes from higher-dimensional ones, until F is defined on all C.This construction is called the lower star filtration.
We define cubical complexes corresponding to lattice sublevel sets of f as Indeed, these are closed under taking boundaries.As stated in the main text, they form a filtration: whenever ν ≤ µ, we have the inclusion C f (ν) ⊆ C f (µ).

Homology groups
Let C be a cubical complex, although the construction of homology groups is the same as for simplicial complexes.More details can be found e.g. in [61] and references cited therein.In this work, we focus on chain complexes and homology groups with Z 2 -coefficients, such that the k-th chain complex C k (C) of C consists of formal sums of chains of k-cubes with Z 2 -coefficients.The boundary operator ∂ k : C k (C) → C k−1 (C) is defined to map a chain of k-cubes to its boundary, which is a (k−1)chain.Since boundaries of such chain boundaries are empty, ∂ k−1 • ∂ k = 0. We define the cycle group as Z k (C) := ker(∂ k ), which consists of all closed k-chains, i.e., k-chains without boundary.The boundary group can be defined as B k (C) := im(∂ k+1 ), consisting of all those k-chains, which are boundaries of (k + 1)-chains.As subgroups, B k (C) ⊆ Z k (C), such that their quotient groups are well-defined: These are called homology groups.
The topology of C can be studied via the homology groups H k (C).They capture similar topological information compared to homotopy groups of C, but are often not the same.Elements of H k (C) are called homology classes and form equivalence classes of k-cycles, defined modulo higher-dimensional boundary contributions.Intuitively, they can be thought of as independent holes.Their number is given by the Z 2 -dimension of H k (C): which is called the k-th Betti number.

Persistent homology groups
Let {C ν } ν∈R be a filtration of complexes, such as C ν = C f (ν) for the energy density sublevel set filtration considered in this work.Suppose we compute all their individual homology groups {H k (C ν )} ν .In addition, the filtration contains for all ν ≤ µ the inclusion maps C ν → C µ , which induce maps on the homology groups: The ι ν,µ k map a homology class in H k (C ν ) either to one in H k (C µ ), if it is still present for C µ , or to zero, if corresponding (potentially deformed) cycles appear as boundaries in H k (C µ ).Moreover, nontrivial cokernels can appear for ι ν,µ k : new homology classes can appear in C µ , which are not in C ν .The parameter µ can be chosen, so that for sufficiently small ϵ > 0: The collection {(H k (C ν ), ι ν,µ k )} ν≤µ forms a so-called persistence module, which is tame, if (C6) holds only for finitely many distinct values of µ.
By the structure theorem of persistent homology (see e.g.[61] and references therein), any tame persistence module is isomorphic to its persistence diagram, i.e., the collection of all the birth-death pairs (b, d), b < d ∈ R ∪ {∞}.Persistence diagrams are multisets, so the same birth-death pair may appear multiple times.
Appendix D: Self-similar scaling in persistent homology In this appendix, we discuss the self-similarity of Betti number distributions as introduced in [63] and similarly employed in [64].For this, we need to introduce the socalled persistence pair distribution.The persistent homology of the energy density filtration is fully described by the persistence diagram, which consists of all birthdeath pairs (b, d).We denote it for dimension-ℓ features as Dgm ℓ,i (t), computed for a classical-statistical realisation ϕ i (t, x).The dimension-ℓ persistence pair distribution is then given by [63]  In Figure 9, we show dimension-0 Betti number distributions for (a) N = 1, (b) N = 2 and (c) N = 3, rescaled in time according to Equation (D4) (averaging over every 8 3 blocks).The optimal scaling exponents η 1 and η 2 for the data to match the scaling behaviour (D4) are extracted using the self-similarity fitting protocol detailed in [51].The fits take all times in the interval from Qt = 2500 to Qt = 12500 into account.They are done for those filtration parameter ranges which correspond to the defect-related peaks in Betti number distributions: for N = 3.Indeed, rescaling the Betti number distributions with these scaling exponents consistently leads at least to approximate constancy in time, see Figure 9, in particular for N = 1 and N = 3.Yet, smaller systematic deviations remain, since the shape of the peaks is not fully independent of time.
According to the Betti number scaling (D4), the peak Betti number scales ∼ t 2η1−η2 .The numbers provided by the self-similarity fits, Equations (D5) to (D7), match the analysis results provided in Section III C. There, we also discuss the value of the exponent 2η 1 − η 2 in light of phase-ordering kinetics.
The value of η 1 is consistent with zero within errors.If connected components at these ∆T 00 values are primarily due to defects, this can indicate that defects appear locally at particular constant energy density values.We note that a zero result for η 1 is not in line with the packing relation [65], which yields η 2 = 3η 1 for energy conservation, providing a one-dimensional constraint on the filtration.This is not in contradiction with [65], since the packing relation has been proven for self-similar scaling that applies to the entire filtration range.This is not the case here, where the ongoing homogenisation of smallscale structures in energy densities provides features in Betti number distributions which cannot be rescaled with the exponents of (D5) to (D7) (cf. the right peaks in Figure 9).This appendix discusses the defect-related left peak in the dimension-1 Betti number distributions for N = 2 (cf. Figure 5).In Figure 10 we show maximum values of the peak depending on time.The main figure has been computed for averaging over every 8 3 blocks, the top inset for averaging over 16 3 blocks.As it is clearly visible, the numbers of dimension-1 features decline with time, such that the average distances associated with the structures grow.Qualitatively, the discussion of Section III C applies here in an analogous way, so that we can associate this to the coarsening dynamics, potentially for domain walls.
Again, for an intermediate time range a power law can be fitted.For averaging over every 8 3 blocks, the peak values decrease as a power law with exponent −0.85 ± 0.01, and with exponent −0.62 ± 0.03 for averaging over 16 3 blocks.While the corresponding length scales of dimension-1 features thus grow as power laws with exponents well within the regime expected for coarsening dynamics (see the discussion in Section III C), the numbers do not agree within uncertainties.This indicates that in dimension-1 Betti numbers the dynamics in the infrared is not well-separated from ultraviolet dynamics for averaging over every 8 3 or 16 3 blocks.Additional potential sources of uncertainties have been discussed at the end of Section III D.

FIG. 1 .
FIG. 1. Two-dimensional slices of local energy density fluctuations around their mean value, given by ∆T 00 = (T 00 − T 00 )/ T 00 .(a) Snapshots at N = 1 for different times, and (b) snapshots for different N at time Qt = 5000, showing signatures of string defects for N = 1, 2, 3 and of domain walls for N = 2.

FIG. 2 .
FIG. 2. Pixelisation of lattice sublevel sets leads to a filtration of cubical complexes.(a): 64 × 64 pixels excerpt from the top-left corner of the N = 1 energy density fluctuation slice shown in Figure 1(b).(b): Three nested cubical subcomplexes corresponding to the indicated sublevel sets of the image of (a), for filtration parameters ν ≡ ∆T 00 ≤ −0.5, −0.2, 0.0 from left to right.

3 FIG. 4 .
FIG. 4. Dimension-0 Betti numbers for ∆T 00 sublevel sets with local energy densities averaged over nearby blocks of every c 3 = 1 3 , 4 3 , 8 3 and 16 3 lattice sites (from left to right).Betti numbers have been normalised to 64 3 lattice size in accordance with their volume-scaling.Betti numbers at initial time Qt = 0 are shown in grey.

FIG. 5 .
FIG. 5. (a) Dimension-0 and (b) dimension-1 Betti numbers for ∆T 00 sublevel sets and N = 1 to N = 4 field components (from left to right), with averaging over blocks of 8 3 lattice sites.Betti number distributions at initial time Qt = 0 are shown in grey.For Betti number distributions with twofold peak structures we label the peak associated with topological defects by "topo.", and the peak corresponding to energy transport by "energy".

FIG. 6 .
FIG. 6. Temporal scaling of the left dimension-0 Betti number distribution peak values for (a) N = 1, (b) N = 2 and (c) N = 3, with averaging over blocks of 8 3 lattice sites.The bottom insets show the Betti number distributions, where the red circles indicate the peak positions.The top insets show peak values for averaging over every 16 3 blocks.The power law fits are based on the data in the red shaded Qt ranges.

∼ t 0. 54 ∼FIG. 7 .
FIG. 7. Temporal scaling of the right dimension-1 Betti number distribution peak values for N = 2. Averaging over blocks of 8 3 lattice sites has been employed.The top inset shows the actual Betti number distributions, where the red circles indicate the peak values.The bottom inset shows peak values for averaging over every 16 3 blocks.The power law fit is based on the data in the Qt range for which it is displayed.Note the smaller time interval shown here in comparison to Figure 6.
Appendix B: The classification of topological defects in condensates of the O(N ) vector model
Appendix E: Maxima of dimension-1 Betti number distributions at N = 2

FIG. 10 .
FIG. 10.Temporal scaling of the left dimension-1 Betti number distribution peak values for N = 3. Averaging over blocks of 8 3 lattice sites has been employed.The bottom inset shows the actual Betti number distributions, where the red circles indicate the peak values.The top inset shows peak values for averaging over every 16 3 blocks.The power law fits are based on the data in the red shaded Qt ranges.