Symmetry breaking at a topological phase transition

Spontaneous symmetry breaking is a foundational concept in physics. In condensed matter, it characterizes conventional continuous phase transitions but is absent at topological phase transitions such as the Berezinskii-Kosterlitz-Thouless (BKT) transition - as in the BKT case the expected norm (i.e., the magnitude) of the $U(1)$ order parameter vanishes in the thermodynamic limit at all nonzero temperatures. Phenomena consistent with low-temperature broken symmetry have been observed, however, in many different BKT experiments. Examples include recent experiments on superconducting films and the seminal work on two-dimensional arrays of Josephson junctions. While the inaccessibility of the above thermodynamic limit partially explains this paradox in finite systems, the full dynamical framework of symmetry breaking at the BKT transition remains unresolved. Here we provide this by introducing the broader concept of general symmetry breaking. This encompasses both spontaneous symmetry breaking and the BKT case by allowing the expected norm of the order parameter to go to zero in the thermodynamic limit, provided its directional phase fluctuations are asymptotically smaller. We demonstrate this asymptotically slow directional mixing in the low-temperature BKT phase. This explicitly shows that the order parameter arbitrarily chooses some well-defined direction in the thermodynamic limit, predicting negligible phase fluctuations compared to the expected norm in arbitrarily large experimental BKT systems. Our results provide a model for directional mixing timescales across the diverse array of experimental BKT systems. We suggest various experiments.

Topology and symmetry are two of the most fundamental concepts in physics.The former classifies objects via properties that are preserved under continuous deformation.It provides a framework for a multitude of physical phenomena, from topological insulators [1] and the quantum Hall effect [2] to fault-tolerant quantum computation [3] and knots in light fields [4] and liquid crystals [5].It was also fundamental to the defect-mediated description of topological phase transitions in condensed matter [6][7][8], where the concept of symmetry -the preservation of system properties under symmetry transformations -has been possibly most powerful in the characterization of conventional continuous phase transitions.For example, at low temperature, the two-dimensional (2D) Ising model restricted to single spin-flip dynamics breaks its global Z 2 symmetry by arbitrarily choosing either a positive or negative magnetization -the global Z 2 order parameter -and keeping this sign on a timescale that diverges with system size [9].However, it is the simplest case of continuous symmetry -the U (1) group of planar rotations -that showcases some of the most interesting and subtle symmetry properties, but at a topological phase transition.
In an ergodic U (1) system, the directional phase ϕ m ∈ (−π, π] of the global U (1) order parameter m = (∥m∥, ϕ m ) ergodically explores (−π, π] on some finite directional mixing timescale: the mean phase converges to its expected value Eϕ m = 0 on this timescale, where * michael.faulkner@warwick.ac.uk it is independent of global rotations due to its fluctuations also converging to their expected value Eϕ 2 m > 0 [the expected value Ef := f (x)π(x)dx of some observable f (x) is that predicted by the Boltzmann distribution π(x) ∝ e −βU (x) , with β the inverse temperature and U the potential].If symmetry is broken under the chosen dynamics, however, the directional mixing timescale diverges with system size.This asymptotically slow directional mixing reflects the order parameter arbitrarily choosing some well-defined direction (with zero phase fluctuations) in the thermodynamic limit.In most cases, this can be expressed mathematically by calculating the expected order parameter Em under the influence of a fixed-direction symmetry-breaking field, before taking the thermodynamic and then zero-field limits.The resultant nonzero vector then aligns with the direction of the field, while the thermodynamic limit is singular [10] because exchanging the order of the two limits returns zero.The singular limit reflects measurable zero-field discrepancies between experimental observations and predictions of the Boltzmann model (i.e., small experimental phase fluctuations compared to their expected value).This elegantly characterizes the dynamical zero-field phenomenon of spontaneous symmetry breaking [11] in terms of thermodynamic expectations, generalizing to any symmetry group and occurring at all conventional continuous phase transitions.Topological phase transitions are, however, an exception, described instead in terms of a topological ordering that typically breaks ergodicity more generally [12].For example, spontaneous symmetry breaking is absent at the Berezinskii-Kosterlitz-Thouless (BKT) transition [6,7,13] because (in zero field) the expected norm E∥m∥ of the U (1) order parameter goes to zero in the thermodynamic limit at all nonzero temperatures [14,15], but a suppression (under Brownian spin dynamics) of global topological defects breaks ergodicity at this paradigmatic topological phase transition [16].
Phenomena consistent with low-temperature broken symmetry have been measured, however, on experimental timescales in a diverse array of BKT systems [17][18][19][20][21][22][23][24][25][26][27][28][29], including on very long timescales [18] and in systems that approach idealized U (1) symmetry [19,20].This suggests that the rather elegant mathematical formalism of spontaneous symmetry breaking may be too restrictive to account for all physical observations.For example, recent measurements of the electrical resistances of films of lanthanum strontium copper oxide (LSCO) exhibited strongly nonergodic probability density functions (PDFs) between the BKT and mean-field superconducting transition temperatures, in contrast with Gaussian Aslamazov-Larkin-type [30][31][32] fluctuations involving both the amplitude and phase of the condensate wavefunction above the mean-field transition [18] (the BKT transition was significantly lower than the mean-field transition, resulting in a broad temperature range dominated by BKT phase fluctuations, i.e., where amplitude fluctuations are negligible [17]).As large condensate-phase differences (over long distances) induce resistance and the experimental timescale was on the order of ten hours, this is likely to be due to increasingly large regions of symmetry-broken (and therefore persistent on some significant timescale) condensate-phase coherence during the approach to the transition -before an onset of distinct behavior in which a single symmetry-broken region spans the entire zero-resistance system at low temperature.(Here, condensate-phase coherence describes positional coherence of the local phase of the condensate wavefunction.)Similarly, zero-resistance measurements on 2D arrays of Josephson junctions also provided direct experimental evidence of system-spanning condensate-phase coherence at low temperature [21,22].Moreover, results consistent with the phenomenon have been measured in superfluid [19,20], magnetic [23][24][25][26] and cold-atom [27][28][29] films, including very recent experiments on cold atoms [29] and a monolayer magnet [26].Elsewhere, experiments on cylindrical arrays of superconducting qubits measured a nonzero order-parameter norm [33], consistent with broken symmetry on some timescale.Indeed, we are not aware of a BKT experiment that contradicts broken symmetry at low temperature.
A partial explanation for the paradox in finite systems was provided by the expected low-temperature norm going to zero very slowly and at the same rate as its fluctuations [34,35].This led to much success in describing magnetic-film experiments [23-26, 36, 37], but the thermodynamic limit was not addressed, and a rigorous dynamical framework for the U (1) order parameter arbitrarily choosing some well-defined direction in the thermodynamic limit (i.e., for broken symmetry) re-mains unresolved.The latter is particularly pertinent to the superconducting film and Josephson-junction array, as the electrical resistance is a directly measurable quantity (on very long timescales [18]) that is conjugate to the directional condensate phases.Moreover, measurements of the magnetization vector in BKT magnetic films should provide direct experimental evidence of system-spanning symmetry-broken spin-phase coherence.In contrast, condensate-phase coherence cannot be directly measured in the superfluid helium film as there is no conjugate field.This is despite the theory establishing its BKT transition [19,38] implying a significant expected low-temperature norm in macroscopic systems [20], consistent with the accompanying systemspanning condensate-phase coherence required for macroscopic superflow.
Here we show that topological nonergodicity in the low-temperature BKT phase [16] induces asymptotically slow directional mixing of the U (1) order parameter, reflecting its phase fluctuations going to zero in the thermodynamic limit (for a model system with Brownian spin dynamics).Moreover, the phase fluctuations are asymptotically smaller than the expected norm.This explicitly demonstrates that the U (1) order parameter arbitrarily chooses some well-defined direction in the thermodynamic limit, providing a theoretical framework for negligible U (1) phase fluctuations compared to the strictly decreasing expected norm in arbitrarily large experimental BKT systems, and thus for broken symmetry to the thermodynamic limit throughout the low-temperature BKT phase.This is an example of general symmetry breaking, which broadens the definition to encompass all phenomena that result in directional (phase) fluctuations going to zero in the thermodynamic limit while being asymptotically smaller than the expected norm (these are the only requirements for broken symmetry, and are fulfilled by cases of spontaneous symmetry breaking as the expected norm is nonzero in the thermodynamic limit).This framework is also described by a singular thermodynamic limit because taking the long-time limit of the phase fluctuations before the thermodynamic limit returns their expected value, while inverting the order of the limits returns zero.These results demonstrate an interplay between topology and symmetry, parting with orthodox theory in the form of a topological phase transition that does in fact break symmetry, but outside the elegant yet restrictive definition of spontaneous symmetry breaking (we stress that symmetry breaking is typically absent at topological phase transitions, but that these results provide the community with a single counterexample).This provides a model for directional mixing timescales across the wide and diverse array of experimental BKT systems.We also introduce the concept of long-time directional stability and infer from earlier work [34,35,39] that the expected norm does not reach its thermodynamic value of zero at arbitrarily large system size.A very recent renormalization-group analysis of U (1) phase fluctuations in BKT magnetic films [37] further enhances the timeliness of this article.
The paper is organized as follows.In Section I, we introduce the U (1) symmetry-breaking order parameter and its phase fluctuations.We then demonstrate that the choice of system dynamics can result in asymmetric lowtemperature simulations on significant timescales, consistent with the existence of a singular thermodynamic limit of the phase fluctuations.In Section II, we define globaltwist dynamics, allowing us to demonstrate that topological ergodicity ensures U (1) symmetry.In Section III, we demonstrate that Brownian spin dynamics break U (1) symmetry in the low-temperature BKT phase, but that an alternative choice of local spin dynamics ensures U (1)symmetric simulations on non-divergent timescales at all nonzero temperatures -in analogy with Swendsen-Wang/Wolff simulations [40,41] of the 2D Ising model.We discuss our results and suggest various experiments in Section IV.
Spontaneous symmetry breaking is therefore absent at the BKT transition.The entire low-temperature phase is, however, critical, as the spin-spin (fluctuational) correlation length diverges for all finite β > β BKT , with β BKT the phase transition.The expected low-temperature norm does not therefore reach its thermodynamic value of zero at arbitrarily large system size, as the correlations are cut off on long length scales.More precisely, with g(x) the PDF of m and σ ∥m∥ the standard deviation of ∥m∥, the PDF σ ∥m∥ g(x) of the fluctuation-normalized order parameter m := m/σ ∥m∥ maintains a well-defined low-temperature sombrero form in the thermodynamic limit because E∥ m∥ is system-size independent [34,35].In general, for any given scalar observable, only if the ratio of its fluctuations with its expectation can be made arbitrarily small with increasing system size does there exist some finite system size at which the expectation can be considered to have reached the thermodynamic limit, which is not the case for E∥m∥ given this system-size independence.All analysis of the global U (1) phase ϕ m therefore holds to the thermodynamic limit.This is also implied by the equivalence of the directional phases of m and m/σ ∥m∥ .We note also that the expected low-temperature norm itself goes to zero very slowly: E∥m∥ ∼ N −1/(8πβJ) as β → ∞ [23] and E∥m∥ ∼ N −1/16 at the finite-size transition [51].The former was even shown to hold in an N ∼ 10 10 superfluid film [20], while the latter implies that E∥m∥ would be ∼ 10 −2 at the finite-size transition in a magnetic film "the size of Texas" [51].
Fig. 1 shows evolutions of the order parameter m over the course of single simulations at various systems sizes and temperatures.The data reflect a PDF with sombrero form at low temperature and the central well of the symmetric phase at high temperature (1/β BKT ≃ 0.887J [52]).In Figs 1(a)-(c), we investigate the 2DXY model with Brownian (spin) dynamics by simulating the model with the local Metropolis algorithm, where each Metropolis iteration proposes a (local) single-spin perturbation and we draw observations after every N such iterations, defining the Metropolis Monte Carlo time step as τ /n with n the number of observations (Metropolis dynamics converge on Brownian dynamics [53] with an N -independent physical time step that is proportional to the Monte Carlo time step, as outlined in Appendix C).  3).The local Metropolis simulations comprise 10 5 observations with acceptance rate aMetrop ≃ 0.6.The event-chain simulations comprise 10 3 observations.Straight lines connect adjacent observations.The Metropolis data suggest a low/high-temperature directional mixing timescale that does/does not diverge with system size, whereas the event-chain data suggest a non-divergent directional mixing timescale for all finite β.In particular, the low-temperature Metropolis simulations become less symmetric with increasing system size, i.e., the (unbiased) simulation variance s 2 ϕm of the global U (1) phase ϕm deviates further from its expected value Var[ϕm] = π 2 /3 [the variance of the uniform distribution U(−π, π)].This is consistent with the U (1) order parameter arbitrarily choosing some well-defined direction in the thermodynamic limit.
ics against the event-chain Monte Carlo algorithm [54] as the latter induces rotations of the order parameter [and therefore U (1) symmetry] on short timescales.As outlined in detail in Appendix C, this is achieved by continuously advancing some active spin at fixed (anticlockwise) velocity until a Metropolis rejection would have occurred.This defines a particle event and induces a new active spin.We draw observations at every N th particle event and define the event-chain Monte Carlo time step as τ /n.The number of observations n is fixed at 10 5 and 10 3 in the Metropolis and event-chain simulations.
In Figs 1(a)-(c), the low-temperature Metropolis simulations are asymmetric for N ≥ 64×64, while the hightemperature Metropolis simulations appear to be symmetric for a broad range of system sizes (with symmetric simulations defined by a simulation variance s 2 ϕm that has converged to its expected value Var[ϕ m ] within some small ϵ > 0).Moreover, the low-temperature directional mixing timescale appears to increase with system size, suggesting a singular limit.Upon defining the long-time directional stability the Metropolis results are consistent with This is due to vanishing low-temperature phase fluctuations in the thermodynamic limit, as presented in detail below.This thermodynamic limit is singular because exchanging the order of the limits in Eq. ( 2) returns zero at all finite β (we note that singular semiclassical limits analogously involving long times are commonplace in quantum chaos [55]).All event-chain results in Figs 1(d)-(f) suggest, by contrast, non-divergent directional mixing timescales, consistent with zero long-time directional stability for all finite β.We assume local Metropolis/Brownian spin dynamics below unless otherwise stated.
We additionally note that the low-temperature Metropolis outputs in Figs 1(a)-(b) contain small fluctuations towards m = 0 at random values of the directional phase ϕ m , while the event-chain outputs in Figs 1(d)-(f) appear to exhibit well-converged simulation variances in ∥m∥.This may be a result of the asymmetry of the lowtemperature ∥m∥ distributions in Fig. 2 of Ref. [34], with the Metropolis dynamics mixing poorly in the heavytailed regions towards m = 0.

II. GLOBAL-TWIST DYNAMICS
We show below that Eq. (3) holds for local Metropolis simulations, but we first demonstrate that zero long-time directional stability is guaranteed at all finite β by supplementing these local dynamics with externally applied global spin twists for all spins k ∈ {(1, 1), (2, 1), . . ., ( √ N , √ N )} along the x/y dimension (q ∈ Z 2 ).Global-twist dynamics are then defined by one such Metropolis proposal with q µ = ±1 along each Cartesian dimension µ at each Monte Carlo time step.Global-twist events occur in pairs that form the compound tunnelling events seen in Fig. 2(a), each due to a global-twist event taking the system to small ∥m∥ before another global-twist event returns the system to the well of the sombrero potential.Fig. 2(b) shows estimates of the probability of 2DXY global-twist events as a function of temperature and system size.The data demonstrate that, as for the case of the topological-sector events that ensure topologically ergodic simulations (defined alongside topological ergodicity in Appendix B) of the 2D electrolyte [16], this probability is system-size independent sufficiently far from the transition.Moreover, the probability is non-negligible at all system sizes and nonzero temperatures (despite being small at low temperature) reflecting it scaling like exp(−2π 2 βJ) as N → ∞ in the absence of other excitations.Assuming Eq. (3), it then follows that, for low-temperature Metropolis simulations supplemented with global-twist dynamics, any sequence of histograms of the global U (1) phase ϕ m at any fixed simulation timescale tends to some normalized sum over randomly distributed (around the well of m sombrero potential) Dirac distributions in the thermodynamic limit, where a large number of Dirac distributions becomes more likely with increased (fixed) simulation timescale and we have assumed suitably chosen his- togram bin sizes.Simulations supplemented with globaltwist dynamics are therefore symmetric on non-divergent timescales, with zero long-time directional stability at all finite β.This is reflected in Fig. 2(c), which shows estimates of the squared phase fluctuations as a function of system size (at βJ = 10) for fixed-timescale Metropolis simulations both with and without supplemental globaltwist dynamics.The data are consistent with the phase fluctuations going to zero in the thermodynamic limit for purely local Metropolis simulations, while supplemental global-twist dynamics increase the phase fluctuations by an amount that increases with system size.As outlined in Appendix B, supplementing local Metropolis simulations with the global-twist dynamics also ensures (and is required for) topologically ergodic low-temperature simulations on non-divergent timescales of the same order (in analogy with topological-sector dynamics in the 2D electrolyte) [16,49].Topological ergodicity therefore ensures U (1) symmetry.

III. SYMMETRY BREAKING
The low-temperature Metropolis data in Figs 1(a)-(c) are consistent with vanishing phase fluctuations (at any fixed simulation timescale and under local dynamics) in the thermodynamic limit, as described by the hypothesized long-time directional stability in Eq. (3).Conversely, at any fixed system size and under local dynamics, the squared phase fluctuations ⟨s 2 ϕm ⟩ will eventually converge to their expected value of Var[ϕ m ] = π 2 /3 [see Fig. 2(c) at small N ] on the directional mixing timescale τ mix > 0. Assuming this is described by ⟨s 2 ϕm ⟩ ∝ τ /τ mix for large enough τ < τ mix under the diffusive Metropolis dynamics, the directional mixing timescale then quantifies the scaling of the fixed-timescale phase fluctuations with system size.
To characterize the directional mixing timescale, we plot the empirical cumulative distribution functions (ECDFs) (a) Cramér-von Mises statistic [defined below Eq. ( 6)] vs reduced temperature βBKT/β and system size N for Metropolis (circles) and event-chain (stars) simulations with supplemental global-twist dynamics, with βBKT := 1/(0.887J).Results indicate a directional mixing timescale τmix ∼ N z/2 with z = 2 and z = 0 for local Metropolis and event-chain dynamics at low temperature (and z = 0 at high temperature in both cases).Metropolis data sets intersect near the transition, marked by vertical grey lines in the inset.Data is averaged over 560 realizations with n = 10 6 at βBKT/β > 1.2 and n = 3×10 7 /n = 10 7 at βBKT/β < 1.2 for N ≷ 40×40.Local Metropolis acceptance rates a ≃ 0.6.(b) Estimated reduced intercept temperatures vs 1/ (ln N ) 2 with a straight-line fit (black dashed line) that extrapolates to one (i.e., the phase transition) as N → ∞, within the estimated error.Inset: Estimated intercept values vs N with a power-law fit indicating approximate ∼ N scaling.
Their high-temperature counterparts suggest, by contrast, symmetric convergence on the target CDF F (x) on some non-divergent directional mixing timescale.This also holds for the low-and high-temperature event-chain simulations in the inset of Fig. 3(a).The Cramér-von Mises mean square distance [56,57] between F ϕm,n (x) and the target CDF F (x) measures the deviations from the target CDF F (x), providing access to the (directional mixing) timescale on which symmetry is achieved via local dynamics.Fig. 4(a) shows estimates of the Cramér-von Mises statistic ⟨nω 2 ϕm,n ⟩ as a function of temperature and system size.We use supplemental global-twist dynamics to improve the statistics at high temperature [58].For large enough n and within the simulation error, the data converge on the displayed values and indicate that ⟨τ ω 2 ϕm,n ⟩ ∼ N for low-temperature Metropolis dynamics and that ⟨τ ω 2 ϕm,n ⟩ is system-size independent for both event-chain dynamics (sufficiently far from the transition) and high-temperature Metropolis dynamics.Since ⟨ω 2 ϕm,n ⟩ → 0 as τ → ∞ and ⟨τ ω 2 ϕm,n ⟩ converges on the directional mixing timescale τ mix , it follows that ⟨ω 2 ϕm,n ⟩ ∝ τ mix /τ for all τ > τ mix .The data therefore imply that τ mix ∼ N z/2 , where at low temperature z = 2 for the physical Metropolis dynamics and z = 0 for the symmetric event-chain dynamics (reminiscent of previous investigations [59]).The data also imply that z = 0 at high temperature in both cases.All data sets meet at a high-temperature plateau, reflecting the flattening of the Boltzmann distribution with increasing temperature.The event-chain results decrease with temperature, reflecting improved directional mixing (see Appendix C).
Topological nonergodicity therefore induces broken U (1) symmetry in the low-temperature BKT phase.Symmetry is broken because the algebraic correlations combine with the diffusive Brownian spin dynamics to provoke a divergence (with system size) of the (directional mixing) timescale on which the directional phase of the U (1) order parameter ergodically explores (−π, π], but symmetry can be restored by non-physical globaltwist dynamics that tunnel through the U (1) sombrero potential and also restore topological ergodicity [16,49].The scaling of the directional mixing timescale implies that the low-temperature U (1) phase fluctuations go to zero in the thermodynamic limit while being asymptotically smaller than the expected norm of the order parameter.This case is distinct from spontaneous symmetry breaking as it cannot be identified via a singular limit of the expected U (1) order parameter, leading us to define general symmetry breaking as directional (phase) fluctuations going to zero in the thermodynamic limit while being asymptotically smaller than the expected norm.This includes spontaneous symmetry breaking and corresponds to the order parameter arbitrarily choosing some well-defined direction in the thermodynamic limit, reflecting all cases of strictly decreasing directional (phase) fluctuations being negligible compared to the expected norm in arbitrarily large systems.
In contrast with the above Metropolis conclusions, event-chain simulations are U (1)-symmetric on nondivergent timescales at all nonzero temperatures.

IV. DISCUSSION
Previous static BKT studies focused on the expected norm of the U (1) order parameter in both thermodynamic [14,15] and finite [34,35] systems, where the latter led to much success in modelling experimental quantities that are functions of the expected norm [19, 20, 23-26, 36, 37].Our dynamical study, by contrast, explicitly demonstrates that the low-temperature order parameter arbitrarily chooses some well-defined direction in the thermodynamic limit.This asymptotically slow directional mixing predicts negligible U (1) phase fluctuations compared to the strictly decreasing expected norm in arbitrarily large experimental BKT systems, thus providing a theoretical framework for broken symmetry to the thermodynamic limit throughout the low-temperature BKT phase.This constitutes a model for the directional mixing (or memory) timescale τ mix ∼ N .We suggest experiments on a single Josephson junction formed from two nodes of 2D superconducting film as an ideal showcase of the phenomenon.Moreover, measurements of the magnetization vector in XY magnetic films (with a sixfold crystal field [47]) should provide further direct experimental evidence, as should the orientational order parameter in the hexatic phase of colloidal films [44,45].Indeed, this paper provides the full dynamical framework for system-spanning symmetry-broken spin/condensatephase coherence in the wide array of critical BKT experimental systems [17][18][19][20][21][22][23][24][25][26][27][28][29]45], motivating experiments on BKT systems in general.
Our results also imply that the strongly nonergodic electrical-resistance PDFs measured in the LSCO films [18] are likely to be due to increasingly large regions of symmetry-broken condensate-phase coherence as the BKT transition is approached from high temperature, described by a critical slowing down of the symmetry-breaking order parameter analogous to that of the 2D Ising model.This is because increasingly large symmetry-broken regions persisting on significant nonergodic timescales are a necessary precursor to the system-spanning symmetry-broken region that we have demonstrated at low temperature (for any system with Brownian dynamics, nonergodic timescales/slow dynamics accompany the flattening of its configurational freeenergy landscape upon transitioning between symmetric and symmetry-broken phases).We leave the remainder of this hypothesis to a separate publication and suggest that this second effect should also be detectable near the BKT transition in 2D Josephson-junction arrays and XY magnetic films, and similarly near the solid-hexatic transition in colloidal films.Indeed, our hypothesis also applies to 3D superconductors, Josephson-junction arrays and XY magnets, but the 2D case is particularly significant to experiment as the BKT transition occurs over a broad temperature range in finite systems [23,36].
We additionally showed that event-chain 2DXY simulations are U (1)-symmetric on non-divergent timescales at all nonzero temperatures.This is analogous to the accelerated mixing (relative to Metropolis) of Swendsen-Wang/Wolff simulations [40,41] of the 2D Ising model, but due in this case to the deterministic event-chain dynamics efficiently exploring the configurational freeenergy landscape along low-gradient directions, with our event-chain simulations also suggesting significantly improved exploration of the heavy tail of the asymmetric ∥m∥ distribution [34] relative to Metropolis.This fundamental characteristic of event-chain Monte Carlo suggests that it should also circumvent critical slowing down by advancing deterministically through the flat regions of the configurational free-energy landscape found at continuous phase transitions -possibly related to the superdiffusive dynamics of the location of the 2DXY active spin [60], and again in analogy with the Swendsen-Wang/Wolff simulations.This hypothesis is also suggested by the √ N speed-up of lifted (relative to standard) Metropolis-Hastings simulations of the Curie-Weiss model at its phase transition [61], as its limiting behavior [61] is a similar piecewise deterministic Markov process that originated in Bayesian computation.This motivates additional studies into the power of piecewise deterministic Markov processes at phase transitions across statistical physics, Bayesian computation and applied probability.Indeed, the foundational viewpoint of broken symmetry presented here -asymptotically slow mixing between equilibrium states of equal probability density -could be viewed as the canonical choice in Bayesian computation and applied probability.This cross-pollination of knowledge may motivate further innovations across all three fields.We also introduced general symmetry breaking.While both fall within this general concept, the present result is distinct from spontaneous symmetry breaking because it cannot be identified by taking the thermodynamic and then zero-symmetry-breaking-field limits of the expected order parameter, though this elegant mathematical formalism [14,15] of singular limits (of statistical expectations) paved the way to the subtleties of the thermodynamic limit in a critical system [34,35].Indeed, the vanguard of these first and second waves of BKT theory correctly identified both an absence of spontaneous symmetry breaking and an experimentally relevant order-parameter norm, both of which are united by our general concept as it allows the expected norm to go to zero in the thermodynamic limit provided the phase fluctuations are asymptotically smaller.We note that, as E∥ m∥ is system-size independent, it is tempting to define the broader concept with respect to m and without stating "while being asymptotically smaller than the expected norm".However, cases of phase fluctuations being asymptotically greater than or equal to E∥m∥ must be excluded from general symmetry breaking, which might not be satisfied by this adaptation.We add that the present work suggests a symmetry-breaking singular limit of E m at low temperature, but preliminary simulation results suggest that E∥ m∥ is nonzero for all finite β in zero field.It would be interesting to explore whether 1/ lim N →∞ E∥ m(h = h 0 )∥ (or potentially lim N →∞ σ ∥m(h=h0)∥ / lim N →∞ E∥m(h = h 0 )∥) can be made arbitrarily small with decreasing ∥h 0 ∥ at low temperature.to the electrolyte-field components is then only approximate.One may circumvent this by using the 2DHXY potential in the decomposition recipe.This defines the local and global topological defects, but the global topological defects will not always correspond to the global twist-relaxation field.
In Figs 2 and 4, we supplemented local 2DXY dynamics (described in Appendix C) with externally applied (to the non-decomposed spin field) global spin twists.Due to the coupling between each component of even the 2DHXY spin field, the global twist-relaxation field t cannot be isolated and uniquely manipulated by these globaltwist dynamics, as such dynamics may alter all three spin excitations.In contrast, topological-sector dynamics alter only the topological sector of the 2D electrolyte, because the auxiliary gauge field is independent of the remaining electric field [16,63].However, Fig. 2(b) shows that, as for the case of topological-sector events in the 2D electrolyte [16], the probability of global-twist events is system-size independent sufficiently far from the transition, and non-negligible at all system sizes and nonzero temperatures, despite being small at low temperature.Defining topological order /nonergodicity under some dynamics by vanishing t fluctuations ⟨s 2 t ⟩ 1/2 in the thermodynamic limit, topological ergodicity then corresponds to and global-twist dynamics ensure topologically ergodic simulations on non-divergent timescales (with topologically ergodic simulations defined by a simulation variance s 2 t that has converged to its expected value Var[ t] within some small ε > 0, reflecting ergodic exploration of the global twist-relaxation field t ∈ Z 2 ).This is analogous to topological-sector dynamics ensuring topologically ergodic simulations (on non-divergent timescales) with respect to w in the 2D electrolyte [16].Since, to leading order, topologically ergodic simulations require ergodic exploration of t only over the five-element set {(0, 0), ±(1, 0), ±(0, 1)} ⊂ Z 2 , it is reasonable to assume that these timescales are on the order of the nondivergent timescales on which symmetric Metropolis simulations are ensured by supplemental global-twist dynamics (r ∈ N of the tunnelling events described in Section II result in r of the Dirac distributions also described there).Global-twist dynamics also ensure ergodic simulations on non-divergent timescales with respect to a non-annealed analogue of the global twist-relaxation field.This will be set out in detail in a future review article [64].
Appendix C: Simulation and fitting algorithms

Local Metropolis Monte Carlo
We used the local Metropolis algorithm with Gaussian noise to investigate the 2DXY model with Brownian spin dynamics.Each Metropolis proposal attempts to perturb a single spin by an amount ε ∼ N (0, σ noise ), with σ 2 noise the variance of the Gaussiannoise distribution.Such proposals are accepted with probability min [1, exp(−β∆U )], with ∆U the potential of the proposed configuration minus its current value.For a set of linearly coupled 1D harmonic oscillators, Metropolis dynamics were proven rigorously [53] to converge on Brownian dynamics (in the thermodynamic limit) with emergent physical time step ∆t phys := a Metrop σ 2 noise ∆t Metrop /2 per unit (Brownian) diffusivity.Here, a Metrop is the Metropolis acceptance rate and ∆t Metrop is the Metropolis Monte Carlo time step, which we define as the elapsed Monte Carlo time between N attempted single-spin moves.We tune σ noise such that a Metrop ≃ 0.6 throughout and draw observations after each Metropolis Monte Carlo time step.

Event-chain Monte Carlo
We benchmarked the diffusive local Metropolis dynamics against the event-chain Monte Carlo algorithm [54].Event-chain Monte Carlo inverts Metropolis: rather than proposing a single discrete spin translation before drawing a random number to decide whether to accept the move, event-chain Monte Carlo draws the random number and then continuously advances some active spin at fixed (anticlockwise) velocity v > 0 until a Metropolis rejection would have occurred at some particle-event time t η > 0.
To sample the particle-event times, we consider a sequence of m proposed Metropolis translations of length δ > 0 in the positive (i.e., anticlockwise) direction, starting from some initial time t 0 ≥ 0.
We draw observations at every N th particle event and define the event-chain Monte Carlo time step as τ /n.This Monte Carlo time step then results in directional mixing timescales in units of CPU time, to allow comparison with the directional mixing timescales of the Metropolis simulations.When estimating quantities related to the norm of the magnetisation ∥m∥, however, one should draw observations every N units of event-chain time or via some homogeneous Poisson process with intensity function proportional to N (this generates a sample of ∥m∥ that reflects the time spent at all values of ∥m∥, rather than a sample of ∥m∥ at the event times).
Finally, in the description of Fig. 4(a) in Section III, we stated that the event-chain Cramér-von Mises statistic ⟨nω 2 ϕm,n ⟩ decreasing with temperature reflects improved directional mixing at lower temperatures.We hypothesize that this is due to increased long-range spin-spin correlations with decreasing temperature.

Polynomial-fitting algorithms
For the local fittings in Fig. 4(a) (inset), we applied natural logarithms to each data set (corresponding to each system size) and then performed second-order polynomial fittings (within each data set) to the three data points nearest to each intercept temperature (defined in Section III).Estimated Monte Carlo errors were used in the fittings.We then performed first-order fittings to a) the resultant intercept temperatures vs 1/ (ln N ) 2 , and b) the natural logarithm of the resultant intercept values vs ln N .We used the NumPy [66] package polyfit for all fittings.
FIG. 1. Evolutions of the U (1) order parameter m := i (cos φi, sin φi)/N over the course of single Metropolis [(a)-(c)] and event-chain [(d)-(f)] simulations suggest that local Metropolis/event-chain dynamics do/do not result in the singular thermodynamic limit described below Eq. (3).The local Metropolis simulations comprise 10 5 observations with acceptance rate aMetrop ≃ 0.6.The event-chain simulations comprise 10 3 observations.Straight lines connect adjacent observations.The Metropolis data suggest a low/high-temperature directional mixing timescale that does/does not diverge with system size, whereas the event-chain data suggest a non-divergent directional mixing timescale for all finite β.In particular, the low-temperature Metropolis simulations become less symmetric with increasing system size, i.e., the (unbiased) simulation variance s 2 ϕm of the global U (1) phase ϕm deviates further from its expected value Var[ϕm] = π 2 /3 [the variance of the uniform distribution U(−π, π)].This is consistent with the U (1) order parameter arbitrarily choosing some well-defined direction in the thermodynamic limit.

FIG. 2 .
FIG.2.Local 2DXY dynamics may be supplemented with global-twist dynamics [defined below Eq. (4)] that ensure U (1)symmetric simulations on non-divergent timescales at all nonzero temperatures.Global-twist events occur in pairs that form compound tunnelling events [through the sombrero potential, as in (a)] and are system-size independent sufficiently far from the transition [see (b)].This results in zero long-time directional stability [defined in Eq. (2)] at all nonzero temperatures, reflected in supplemental global-twist dynamics increasing the squared phase fluctuations ⟨s 2 ϕm ⟩ [defined below Eq. (1)] by an amount that increases with system size [see (c)].Simulations use local Metropolis dynamics with acceptance rate aMetrop ≃ 0.6.(a) Evolution of the U (1) order parameter m over the course of a single Metropolis simulation (of an N = 256×256 system) comprising 5 × 10 5 observations with supplemental global-twist dynamics.(b) Probability of global-twist events vs reduced temperature βBKT/β and system size N based on 560n attempts, with n = 10 6 at βBKT/β > 1.2 and n = 3×10 7 /n = 10 7 at βBKT/β < 1.2 for N ≷ 40×40 [ βBKT := 1/(0.887J)].(c) Low-temperature (βJ = 10) squared phase fluctuations vs 1/ ln N with and without supplemental global-twist dynamics, averaged over 5600 simulations, each of n = 10 6 observations.

FIG. 3 .
FIG. 3. ECDFs [defined in Eq. (5)] of the global U (1) phase ϕm under local Metropolis dynamics [(a)-(b)] reflect the strikingly different symmetry properties of the low-and high-temperature phases [red and black schematics in (c), respectively].Different line styles represent different realizations.(a)-(b) Eight Metropolis simulations at low and high temperature, each of n = 10 6observations with acceptance rates a ≃ 0.6.Results suggest that any sequence of low-temperature ECDFs at any fixed simulation timescale tends to some Heaviside step function as system size N → ∞.The high-temperature simulations suggest, by contrast, symmetric convergence on the target CDF F (x) of the uniform distribution U(−π, π) on some non-divergent directional mixing timescale.Inset of (a): Four event-chain simulations at low and high temperature, each of n = 10 4 observations.Results are consistent with symmetric convergence on the target CDF F (x) on non-divergent directional mixing timescales at all nonzero temperatures.

FIG. 4 .
FIG.4.Brownian spin dynamics break symmetry throughout the low-temperature phase, in contrast with event-chain dynamics.(a) Cramér-von Mises statistic [defined below Eq. (6)] vs reduced temperature βBKT/β and system size N for Metropolis (circles) and event-chain (stars) simulations with supplemental global-twist dynamics, with βBKT := 1/(0.887J).Results indicate a directional mixing timescale τmix ∼ N z/2 with z = 2 and z = 0 for local Metropolis and event-chain dynamics at low temperature (and z = 0 at high temperature in both cases).Metropolis data sets intersect near the transition, marked by vertical grey lines in the inset.Data is averaged over 560 realizations with n = 10 6 at βBKT/β > 1.2 and n = 3×10 7 /n = 10 7 at βBKT/β < 1.2 for N ≷ 40×40.Local Metropolis acceptance rates a ≃ 0.6.(b) Estimated reduced intercept temperatures vs 1/ (ln N ) 2 with a straight-line fit (black dashed line) that extrapolates to one (i.e., the phase transition) as N → ∞, within the estimated error.Inset: Estimated intercept values vs N with a power-law fit indicating approximate ∼ N scaling.

FIG. 5 .
FIG. 5.Each spin configuration is composed of local topological defects (spin vortices), global topological defects (internal global spin twists) and continuous fluctuations around these topological defects (spin waves).Red arrows represent spins.(a) Typical 2DHXY configuration at βJ = 1.5 with fixed topological defects [i.e., for illustrative purposes, we fixed two local defects (and one global defect) in position and allowed spin waves to fluctuate around this constrained configuration -this is distinct from the non-constrained simulations presented in the main body of the paper].(b) Zero-temperature minimization of (a), i.e., with spin waves annealed away.(c)-(d) Configuration in (b) split into its (c) vortex and (d) internal global-twist components.

3 .
Other simulation details Simulations for Figs 1, 2(a) and 3 started from randomized initial configurations.Those for Figs 2(b)-(c) and 4 started from ordered initial configurations.10 4 and 10 5 initial equilibration observations were discarded (respectively) in each event-chain and Metropolis simulation.All non-visible error bars are smaller than the marker size.We consider simulations to have directionally mixed (i.e., ⟨nω 2 ϕm,n ⟩ to have converged) if ⟨s 2 ϕm ⟩ > 0.98Var[ϕ m ].