{'fraction_non_alphanumeric': 0.03433977690825206, 'fraction_numerical': 0.009986567774338609, 'mean_word_length': 4.725008358408559, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Topological-sector fluctuations and ergodicity breaking at the Berezinskii-Kosterlitz-Thouless transition

The Berezinskii-Kosterlitz-Thouless (BKT) phase transition drives the unbinding of topological defects in many two-dimensional systems. In the two-dimensional Coulomb gas, it corresponds to an insulator-conductor transition driven by charge deconfinement. We investigate the global topological properties of this transition, both analytically and by numerical simulation, using a lattice-field description of the two-dimensional Coulomb gas on a torus. The BKT transition is shown to be an ergodicity breaking between the topological sectors of the electric field, which implies a definition of topological order in terms of broken ergodicity. The breakdown of local topological order at the BKT transition leads to the excitation of global topological defects in the electric field, corresponding to different topological sectors. The quantized nature of these classical excitations, and their strict suppression by ergodicity breaking in the low-temperature phase, afford striking global signatures of topological-sector fluctuations at the BKT transition. We discuss how these signatures could be detected in experiments on, for example, magnetic films and cold-atom systems.

An important discovery of BKT and later authors [11][12][13]26] is that the defect-mediated transition of the plane rotator (or 2D-XY) model and its analogues can be mapped to the insulator-conductor transition of a twodimensional Coulomb gas [29]. The long-range Coulomb interactions emerge from a purely local Hamiltonian, so that the mapping at the microscopic level, although complete [30], is far from transparent. However, as Maggs and co-workers [31][32][33][34][35][36] have shown in three dimensions, a Coulomb fluid can be transformed into a local problem by using an electric-field representation and introducing a freely fluctuating auxiliary gauge field. Following this * peter.holdsworth@ens-lyon.fr work, it is straightforward to show that the XY Hamiltonians that admit a BKT transition map onto this generalized electrostatic problem in two dimensions. A practical consequence of the phase-space extension to a fluctuating auxiliary gauge field is the development of purely local algorithms for the simulation of Coulomb fluids in both three [31][32][33][34][35][36] and two dimensions [37], which circumvent the technical difficulties associated with longrange interactions. In particular, the logarithmic potential that governs charge-charge interactions in the twodimensional Coulomb gas is dealt with locally, allowing a new approach to the efficient simulation of twodimensional Coulombic systems.
In this paper we exploit these developments to formulate and simulate a lattice-field description of the twodimensional Coulomb gas for the purpose of investigating the topological properties of the BKT transition. The BKT transition is topological in the sense that it separates a topologically ordered phase from a disordered one. Topological order in this context means that the local topological defects (charges in the two-dimensional Coulomb gas) are confined. Vallat and Beck [30] considered the two-dimensional XY model on a torus, and showed how a winding field can be associated with the global topology of the system. In the high-temperature phase, where charge is deconfined, non-zero values of this winding field define global topological defects that are distinct from the local topological defects driving the BKT transition. Here we show that the lattice-field description naturally lends itself to classifying and investigating this property. In this paper, we treat only the twodimensional Coulomb gas, but in a further publication we will extend our analysis to the case of two-dimensional XY models on the torus.
Our key observation is that the topology of the torus on which the Coulomb gas is placed generates a multiplicity of states in the lattice electric-field representation that are equivalent for charge configurations but not energetically degenerate. Given an arbitrary charge distribution, one is at liberty to add an integer multiple of some constant to each component of the harmonic mode of the electric field while leaving the charge distribution unchanged. This global topology associated with the BKT transition describes the winding of charges around the torus. In the high-temperature phase, charge deconfinement allows for fluctuations in the winding component of the harmonic mode, which can be classified as different topological sectors. Below the transition, however, the binding of charge pairs causes the winding component to be zero. Topological-sector fluctuations in the electric field therefore mark the appearance of the hightemperature, topologically disordered phase at the BKT transition.
The present study of topological-sector fluctuations in the two-dimensional Coulomb gas may be compared to previous studies on the three-dimensional Coulomb phase of spin-ice materials and models [38][39][40][41]. In spin ice, the onset of topological-sector fluctuations is shown to signal a Curie law crossover [38] for the zero-field susceptibility and a Kasteleyn transition in the presence of an applied field [39,40]. Our study of the BKT transition reveals aspects of topological-sector fluctuations that are not found in either of these established cases. For example, by our analysis, the two-dimensional Coulomb gas should be considered to present an ergodicity-breaking transition to a topologically ordered phase in the absence of an applied field, whereas spin ice has no equivalent phase.
The paper is organised as follows. In Section II, we introduce the lattice-field representation of the twodimensional Coulomb gas on a torus and use this to define the partition function and the topological sectors of the electric field. We use numerical simulations to demonstrate that topological-sector fluctuations appear in the high-temperature (conducting) phase but not in the lowtemperature (insulating) phase. In Section III, we show that the reason for the strict suppression of topologicalsector fluctuations in the low-temperature phase is ergodicity breaking at the transition. A finite-size scaling analysis, given in Section IV, confirms that in the thermodynamic limit, ergodicity is broken precisely at T BKT . Conclusions, and comparisons with experimental systems, are discussed in Section V.

II. TOPOLOGICAL-SECTOR FLUCTUATIONS
Using the unit system defined in Appendix A, we formulate the two-dimensional Coulomb gas using discrete vector calculus on a square lattice with periodic boundary conditions (PBCs) applied. The PBCs enforce the toroidal topology but not the curvature of a true torus. All functions are defined to be the discrete counterparts of smooth vector fields [42], and any lattice vector field F is defined component-wise via [42] F(x) where x is any lattice site and e x/y is the unit vector in the x/y direction. The operators∇ and∇ are the forwards and backwards finite-difference operators on a lattice, respectively, and the lattice Laplacian is defined by ∇ 2 :=∇ ·∇ [42]. The most general electric field E may be Helmholtz decomposed into the sum of a Poisson (divergence-full) component −∇φ, a rotational compo-nentẼ and a harmonic componentĒ: This electric field is the most general solution to Gauss' law on a lattice:∇ Here ρ(x) := qm(x)/a 2 is the charge density at each lattice site x, q is the elementary charge, the integer m denotes the charge species, a is the lattice spacing and 0 is the electric permittivity of free space (see Appendix A). Using the field E adds an auxiliary fieldẼ to the usual solution of electrostatics, as in the electrostatic model of Maggs and Rossetto (MR) [31]. This allows us to simulate the physics of Coulombic interactions on a lattice via local electric-field updates, avoiding the need to treat computationally intensive long-range interactions, as outlined in Appendix B. The validity of introducing the auxiliary field is seen in the context of the separability of the partition function into its Coulombic and auxiliary components: the auxiliary field contributes to the internal energy of the system, but it is statistically independent of the Coulombic element. In Appendix B, we give a full description of the algorithm and a derivation of the partition function for the Coulomb gas of multi-valued charges.
The internal energy of the electric fields corresponding to a given charge and auxiliary-field configuration is given by where D is the set of all lattice points. To represent the Coulomb gas in the grand canonical ensemble, we add a core-energy term U Core , given by where c (m) is the core-energy constant of each charge mq, and c (m) = c (−m) since charges are excited in neutral pairs. The grand-canonical energy of the system U = U 0 + U Core may be expanded by combining Eqs.
(2), (4) and (5) to give a sum of terms arising from the different field components, which add to the core energy: Here (respectively), U Rot := 0 a 2 x∈D |Ẽ(x)| 2 /2 and U Harm := 0 L 2 |Ē| 2 /2 are the auxiliary-field and harmonic-mode components of the grand-canonical energy, and U Self and U Int are the self-energy and Coulombic charge-charge interaction components. As outlined in detail in Appendix B, the latter two components may be expressed in terms of the lattice Green's function G, according to . Note that, while U Int can be negative, the sum U Self + U Int is necessarily ≥ 0, as it arises from the term in |∇φ| 2 .
Using the above results we may define the chemical potential for the introduction of a neutral pair of charges ±mq: In the following, we specialise to a Coulomb gas of n pairs of elementary charges of chemical potential µ := µ 1 , by setting c (m = 0, ±1) = 0 and c (m = 0, ±1) = ∞ [43]. The harmonic modeĒ is a uniform field found by averaging the total electric field E(x) over x. In a simply connected space, the average field may be conveniently related to the average polarization P arising from an effective surface charge distribution byĒ = −P/ 0 . For a charge-neutral system in a simply connected space, P := x∈D xρ(x)/N is invariant with respect to the origin shift x → x + x 0 , and is therefore origin-independent. The situation is more complicated on a toroidal surface asĒ can also depend on a harmonic-field component that corresponds to a charge winding around the torus, and it is necessary to adopt a convention to define distances between points (the concepts 'close together' and 'far apart' are ambiguous on a torus). In Appendix C, we show in detail how it is possible to define origin-independent po-larizationĒ p and windingĒ w components of the harmonic mode such that where w is an integer-valued winding field chosen such thatĒ and L is the lattice length. This decomposition of the harmonic field has the following interpretation. A charge pair may unbind and wind around the torus in opposing directions before assuming its original configuration. When a single charge winds around the torus in the x/y direction, the x/y component of the harmonic mode of the electric field E x/y increases by ±q/L 0 . As shown in Appendix C, the lowest-energy harmonic mode that describes an arbitrary charge distribution is therefore an element of the set (−q/2L 0 , q/2L 0 ] and is defined as the polarization component in Eq. (10) by applying modular arithmetic toĒ. The remainder is the winding component. The modulo operation removes any need for a 'distances' convention to define the polarization component, as well as any origin dependence of the field components (see Appendix C for further details).
With these results we may use the integer-valued winding field w to define the topological sector of the system as the number of times charges wind around the torus in the x and y directions, with all non-trivial topological sectors given by w = 0. The topological sector of the system changes any time a charge pair unbinds and winds around the torus and hence thermal fluctuations of the topological sector are closely related to the unbinding of charge pairs, as elucidated further below.
The statistical mechanics of the topological-sector fluctuations may now be formulated by considering how the polarization and winding components of the harmonic mode enter the lattice partition function. As shown in Appendix B, the partition function splits into two statistically independent components. One component is the Coulombic partition function Z Coul and contains all information about the charge-charge correlations, while the other is the auxiliary-field partition function Z Rot and contains all information about the auxiliary field: the auxiliary field can freely fluctuate without affecting charge-charge correlations (see Appendix B). The partition function is written as where Z Coul is given by Here, β := 1/k B T is the inverse temperature and the sum The first exponential of Eq. (12) describes the anharmonic charge-charge interactions, the second describes the sum of the self-energies associated with each elementary charge, and the third describes the polarization and winding state of the system. The sum over the winding field w is necessitated by the degeneracy of the harmonic mode of the electric field: infinitely many topological sectors describe any given charge configuration. In addition to the local updates of the MR algorithm (see Appendix B), we also consider global updates, which correspond to independently sampling this winding field.
Henceforth we set the elementary charge q = 2π, the lattice spacing a = 1, the electric permittivity 0 = 1, and Boltzmann's constant k B = 1. The choice q = 2π recognises the standard BKT theory, where a charge emerges as a local 2π winding in an associated lattice field, such as the spin differences in the 2D-XY model [26].
The BKT transition drives the deconfinement of charge pairs in the two-dimensional lattice Coulomb gas, which generates topological-sector fluctuations. The transition occurs at T BKT = 1.35 (to 3 significant figures) [44] in the thermodynamic limit (a value specific to a gas of elementary charges with c (m = 1) = 0), which is scaled to higher temperatures in finite-size systems (see below). Fig. 1 shows the evolution of the (normalized) x-component of the harmonic mode of a system of linear size L = 16, simulated using local moves only (numerical simulation details are described in Appendix D). No topological-sector fluctuations are visible just below the BKT transition temperature T BKT = 1.35, but they become important at temperatures above T BKT .

III. ERGODICITY BREAKING
A convenient measure of the topological-sector fluctuations is the winding-field susceptibility χ w , defined by In a fully ergodic system, χ w is nonzero, even in the absence of charge fluctuations, as can be seen by limiting the Gibbs ensemble contributing to Z Coul to field configurations of zero charge. In this case it is straightforward to show, using Eqs. (9) and (12), that the constrained susceptibility is given by for k B T q 2 /2 0 . The system-size dependence falls out of this expression so that a fully ergodic system would show small but finite topological-sector fluctuations in the low-temperature phase.
Assuming local charge dynamics, a topological-sector fluctuation requires the separation of a pair of charges over a distance greater than L/2 in either the x or the y direction (see Eq. (8) and the subsequent discussion). As the charge concentration falls to zero at low temperature, screening becomes negligible and the energy barrier against such configurations diverges logarithmically with the linear system size L [11,12,29]. As the charge concentration increases with temperature, however, entropy and charge screening break down the free-energy barrier, making it finite at the BKT transition. Above the transition, charge pairs are free to unbind and trace closed paths around the torus, giving finite-valued winding fields, as observed in Fig. 1. In contrast, in the low-temperature phase, the probability of separation becomes strictly zero in the thermodynamic limit.
The BKT transition is therefore an ergodicity breaking: a change in the phase space explored by a system with local dynamics. In detail, it is an ergodicity breaking between topological sectors, signalled by the strict suppression of topological-sector fluctuations in the electric field at T < T BKT . If the dynamics were non-local (including global updates of the winding component of the harmonic mode [31][32][33][34][35][36]), χ w would remain finite at all temperatures. To explore this ergodicity breaking, we have simulated the two-dimensional Coulomb gas, first with local field updates only, and second with both local and global field updates [31][32][33][34][35][36]. Corresponding to each case, we define the winding-field susceptibilities χ local w and χ all w , respectively. Differences between χ local w and χ all w reflect the inability of local moves to explore a fully representative phase space on the time scale of the simulation. To quantify this, we introduce the susceptibility quotient χ local w /χ all w , which may be used to analyse the ergodicity of the system. Fig. 2 clearly shows that ergodicity is broken in the vicinity of the BKT transition. For T > 1.6, χ local w = χ all w , indicating that the free-energy barrier for a topologicalsector fluctuation via local moves is small. For T < 1.2, the quotient is zero, indicating that the energy barrier prevents topological-sector fluctuations via local charge moves. In between these low-and high-temperature regions there are strong fluctuations in the quotient because charge deconfinement via local updates represents increasingly rare events, an inevitable precursor to loss of ergodicity. In Section IV, this ergodicity breaking is shown to occur precisely at T BKT in the thermodynamic limit.
Our analysis thus leads to a precise definition of topological order for the two-dimensional Coulomb gas through the ergodic freezing of the topological sector to its lowest absolute value. Two-dimensional systems with U (1) symmetry are often associated with an absence of an ordering field at finite temperature [45]. Here we explicitly show that, in the case of the BKT transition, the ordering of a conventional order parameter is replaced by topological ordering through an ergodicity breaking between the topological sectors. The topological order is directly related to the confinement-deconfinement transition of the charges, the local topological defects of the electric field. This type of ergodicity breaking is distinct from either the symmetry breaking that characterizes a standard phase transition, or that due to the rough free-energy landscape that develops at a spin-glass transition [46].

IV. FINITE-SIZE SCALING
In order to explore the approach to the thermodynamic limit, the two-dimensional Coulomb gas was simulated by the Monte Carlo method as a function of system size, using the MR algorithm. The global update was employed in order to improve the statistics (numerical simulation details are described in Appendix D). Fig. 3 shows the simulated winding-field susceptibility χ w as a function of temperature for L × L Coulomb gases of linear sizes between L = 8 and L = 64. There is a marked increase in the winding-field susceptibility χ w as the system passes through the BKT transition temperature T BKT = 1.35 [44] for all system sizes. Susceptibility curves for successive values of L intersect at temperatures above T = 1.8 and below T = 1.5. Between these two temperatures, the winding-field susceptibility increases for a given temperature as the linear system size L increases. These results are consistent with the finite-size scaling of the BKT transition temperature [22]: as the system size decreases the effective transition temperature T * (L) increases.
Below T BKT , the probability of a charge pair separating over a distance greater than L/2 increases with decreasing system size. This, combined with the finitesize transition temperature T * (L) also increasing with decreasing system size, results in the winding-field susceptibility curves for successive values of L intersecting in the vicinity of T BKT . The inset in Fig. 3 shows that the low-temperature crossover points of the susceptibility curves are at T = 1.45, T = 1.40 and T = 1.37. To extrapolate the trend of the data shown in Fig. 3 to the thermodynamic limit, we define the crossover temperature T Cross (L) to be the lower temperature at which χ w (L) = χ w (L/2). Fig. 4 shows the crossover temperature T Cross as a function of inverse linear system size 1/L, along with straight-line fits to the data. In the thermodynamic limit, T Cross extrapolates to the value T Cross = 1.35 to within the estimated error of the extrapolation, that is, it extrapolates to the BKT transition temperature [44]: The 1/L scaling of T Cross is unusual for the BKT transition, for which the finite-size BKT transition temperature typically scales as a simple function of the logarithm of L [22,47]. However, Minnhagen and Kim [48]   found that a fourth-order cumulant that measures fluctuations of the helicity modulus in the 2D-XY model also scales as 1/L: as this closely relates to fluctuations in the harmonic-mode susceptibility [30], it seems likely that we are observing the same finite-size scaling here. The magnitude of the winding-field susceptibility at the crossover points χ Cross w (L → ∞) similarly extrapolates to ∼ 4 × 10 −4 in the thermodynamic limit, with an estimated error of the same order. This small number is not measurably different to the winding-field susceptibility due to global moves only, which, at T BKT , evaluates to ∼ 5 × 10 −5 for all system sizes (see Eq. (14)). The inference is that topological-sector fluctuations due to local moves only turn on precisely at the universal point T Cross (L → ∞) = T BKT in the thermodynamic limit. This confirms that topological-sector fluctuations signal charge deconfinement and the high-temperature phase of the BKT transition.
Given that the topological-sector fluctuations turn on at the temperature at which the system experiences the famous universal jump in the helicity modulus [13,26,30,48], it is interesting to estimate the contribution that topological-sector fluctuations make to the universal jump. To do this, we define the harmonicmode susceptibility χĒ and the polarization susceptiblity χ p by replacingĒ w in Eq. (13) withĒ andĒ p , respectively. The helicity modulus is then given by Υ = −1 0 (1 − χĒ/2) [30], so that χĒ makes a jump of order unity at T BKT . We find that the ratio (χĒ−χ p )/χĒ is less than 5×10 −2 for all T ≤ 1. small. This reflects the near-cancellation of Ē 2 w with the coupling term 2 Ē p ·Ē w in the evaluation of χĒ, reflecting strong correlations between the polarization and winding fields at the transition.

V. CONCLUSIONS
In conclusion, the BKT transition has long been a paradigm for the importance of topological defects in condensed-matter physics [1]. Vallat and Beck showed that XY-type systems on tori generate global topological defects at the BKT transition that reflect the torus topology [30]. Here we have used lattice-field simulations to reveal topological-sector fluctuations in the electric field of a two-dimensional lattice Coulomb gas on a torus. We have shown how these provide a striking and sensitive measure of the topological, and ergodicitybreaking, character of the BKT transition, allowing a precise definition of topological order in terms of this broken ergodicity.
The topological-sector fluctuations at the BKT transition are very clearly revealed in the lattice electric field description of the two-dimensional Coulomb gas, but we expect them to be equally relevant to any system that has a BKT transition. In suitable systems, the winding-field susceptibility that signals the onset of topological-sector fluctuations will contribute to experimentally measurable responses of the system. For example, in a cylindrical or toroidal magnetic film with XY symmetry, winding-field fluctuations in the Coulomb gas representation correspond to measurable spin configurations in the magnetic representation. As we will show in future work [49], fluctuations of an appropriate topological sector accompany the destruction of the finitesize magnetization of an XY spin system through vortex deconfinement. They could therefore be observable in ultrathin ferromagnetic metallic films [50] or magnetic Langmuir-Blodgett films [51,52].
Another promising system on which to measure these topological-sector fluctuations is the one-dimensional quantum lattice Bose gas. When the system is placed on a ring, its angular momentum is no longer a good quantum number. The angular momentum can therefore fluctuate quantum mechanically, and the system should undergo a dramatic increase in these fluctuations as it passes through the superfluid -Mott insulator quantum phase transition [53,54]. This dramatic increase in the fluctuations corresponds to finitevalued global topological defects in the quantum system, and therefore, via the Feynman path-integral mapping, to topological-sector fluctuations in the two-dimensional classical lattice Coulomb gas on a torus. Murray et al. measured the angular momentum of ring-shaped Bose-Einstein condensates via the vortex-density profile of the system [55]. Our measure of the BKT transition could therefore correspond to equivalent, experimentally measurable topological-sector fluctuations in the cold-atom system [56].
Finally, it is worth noting that it is natural to associate a conducting phase with the excitation of winding fields, as may be seen by considering a loop of wire in a changing magnetic field. Recalling that the magnetic field does no work on a test charge, the induced electromotive force must arise from a divergence-free electric field running round the loop. The curl of this field obeys the Maxwell-Faraday law, ∇ × E = ∂B/∂t. Hence, in three dimensions, electromagnetic induction provides a practical method of exciting topological winding fields analogous to those discussed here.

ACKNOWLEDGMENTS
It is a pleasure to thank A. C. Maggs, S. T. Banks, V. Kaiser and G. B. Davies for valuable discussions, and A. Gormanly for help with automating repeated simulations. We are also grateful to T. Roscilde for pointing out the possible application to the one-dimensional lattice Bose gas. M.F.F. is grateful for financial support from the CNRS and University College London. P.C.W.H. acknowledges financial support from the Institut Universitaire de France. In the following, [ . . . ] denotes the dimensions of some quantity, L denotes the dimensions of length, d is the spatial dimensionality of the system, and 0 is the vacuum permittivity in d−dimensional space. Consider Gauss' law for the MR algorithm, and the dimensions of electric charge density, which generates From a consideration of the exponent of the Boltzmann factor (with β = 1/k B T ) we find a dimensionless group, Setting the charge to be dimensionless, it follows, in d = 2: The MR electrostatic model is a lattice-field model from which it is possible to form the lattice partition function of electrostatics. To show this, we describe the MR algorithm in terms of microscopic variables that represent the local field updates. A conjugate lattice D is defined such that each of its sites is at the centre of each plaquette of D. Each site in D is associated with a realvalued variable ϕ whose adjustment corresponds to an update of the auxiliary field, while each pair of nearestneighbour sites is associated with an integer-valued variable s whose adjustment corresponds to a charge-hop update. Both sets of variables are subject to PBCs.
Component-wise, we now define the field and identify A charge hop in the positive x/y direction corresponds to a decrease/increase in the relevant s variable by an amount q, as shown in Fig. 5 (where s ij represents the s variable between sites i and j of the conjugate lattice).  With the internal energy of the electric fields corresponding to a given charge and auxiliary-field configuration given by U 0 in Eq. (4), it is possible to write the partition function in the microscopic-variable representation. For ease of manipulation, we allow charge-hop updates to create charge pairs out of the vacuum and include the possibility of all integer-valued multiples of the elementary charge. Combining Eqs. (4) and (B2), the partition function in the microscopic-variable representation is given by where and Here, the grand-canonical energy of the Coulombic system U = U 0 + U Core is used.
This representation reproduces Gauss' law: where Q Γ is the charge enclosed within some subset of the lattice Γ ⊆ D, ∂Γ ⊂ D is the boundary enclosing Γ, and l traces an anticlockwise path along ∂Γ and has dimensions of length. This equation results from the ϕ variables cancelling and the s variables being integer valued. It follows that recovering Eq. (3), as required. The constraints imposed upon the electric field (Gauss' law and the form of the harmonic mode, the latter of which is derived in detail in Appendix C) are combined with the grand-canonical energy of the system to write the partition function in terms of the electric field. We define the set X := qZ/a 2 , such that the partition function is given by where the functional integral for any vector field F.
This partition function may be separated into two components by defining the new rotational field The partition function is then given by where and are the Coulombic and auxiliary-field components of the partition function, respectively. Here Y := qZ/ 0 a 2 and we have used the fact that all coupling terms in the grand-canonical energy sum to zero. The MR algorithm therefore reproduces Coulombic physics since the separation of the auxiliary-field partition function from the Coulombic partition function ensures that the chargecharge correlations are independent of the auxiliary field. The lattice Green's function G(x, x ) between two charge-lattice sites x and x is defined such that where the subscript x denotes with respect to which coordinate system the lattice Laplacian is applied. We define the k-space lattice Green's functionG, and the set k∈B := i∈{x,y} ki∈Bi , where where the k = 0 part of the lattice Green's function is set to zero since the harmonic component of E is attributed toĒ. It follows that The internal energy of the Poisson component of the electric field is given by hence, the Coulombic partition function can be written as

Appendix C: Polarization
We consider the sum of each component of the electric field over the entire lattice in order to analyse the harmonic mode. The sum of the x/y-component is split into separate sums over all x/y-components that enter a particular strip of plaquettes of width a that wrap around the torus in the y/x direction. Each component of the harmonic modeĒ x/y is then expressed in terms of the charge enclosed along each of the strips of plaquettes: which follows from applying Gauss' law to each strip of plaquettes that wrap around the torus in the y direction.
The same argument holds for the y-component, hence, the harmonic mode is given bȳ where P := x∈D xρ(x)/N is the origin-dependent polarization vector of the system and w 0,x := 0 a L y=a E x (a/2, y)/q is the x-component of the origindependent winding field, with the y-component defined analogously. Here, P and w 0 are measured from a specific origin. Note that the above applies to systems composed of either single-or multi-valued charges.
We have thus shown thatĒ, which is originindependent, is given by the sum of two origin-dependent terms. One of these is attributed to the polarization of the system, while the other describes the winding of charges around the torus given that the polarization is measured with respect to the chosen origin.
Restricting our attention to the gas of elementary charges, we now devise an origin-independent measure of the topological sector of the system. First, we note that adding ω windings to either component of the harmonic modeĒ corresponds tō and that this results in a change in the grand-canonical energy of the system given by Hence, given an arbitrary charge distribution, the lowestenergy harmonic mode that describes the charge distribution is an element of the set in Eq. (10). We therefore define a convention in which the harmonic mode is given by Eq. (8), where the polarization component of the harmonic mode is an element of the set in Eq. (10) and the winding component of the harmonic mode is given by Eq. (9).
Appendix D: Simulation details The system was simulated using the MR algorithm on an L × L lattice of lattice spacing a = 1. One charge-hop sweep corresponded to picking a charge site at random, picking the x or y direction at random, then proposing a charge hop in the positive or negative direction (at random), repeating this 2N times (replacing each site / field bond after each proposal). One auxiliary-field sweep corresponded to picking a charge site at random and proposing a field rotation around the site, repeating this N times. One global sweep corresponded to proposing a winding update in the positive or negative (at random) x and y directions. For all simulations, we performed five auxiliary-field sweeps per charge-hop sweep, and, for those simulations that also employed the global update, we performed one global update per charge-hop sweep. One charge-hop sweep corresponds to one Monte Carlo time step.
The data sets in Sections III and IV were averaged over multiple runs of 10 6 charge-hop sweeps per lattice site. The data set in Fig. 2 was averaged over 608 and 446 runs between T = 1.15 and 1.45 with the global update off and on, respectively, over 384 runs between T = 1.5 and 1.6, and over 256 runs between T = 1.65 and 1.75.
The L = 8 data set in Fig. 3  We also simulated the L = 10, L = 20, and L = 40 systems over small temperature ranges to calculate additional crossover points for Fig. 4: all data sets were averaged over 512 runs.