Topological excitations in statistical field theory at the upper critical dimension

We present a high-precision Monte Carlo study of the classical Heisenberg model in four dimensions, showing that in the broken-symmetry phase it supports topological, monopole-like excitations, whose properties confirm previous analytical predictions derived in quantum field theory. We discuss the relevance of these findings and their possible experimental applications in condensed-matter physics.


Introduction
Lattice models of interacting vector spins with global O(N ) symmetry have many important theoretical as well as experimental applications. Considering only nearest-neighbor interactions, their Hamiltonian can be written as [1] where s(x) is an N -component real vector of unit length, defined on the site x of a regular Euclidean lattice in D dimensions, the summation runs over all distinct pairs of nearest-neighbor sites, and J > 0 corresponds to ferromagnetic coupling. As particular cases, eq. (1) includes the self-avoiding random-walk model (for N = 0), the Ising model (for N = 1), the XY model (for N = 2), the Heisenberg model (for N = 3), a toy model for the Higgs sector in the Standard Model of elementary particle physics (for N = 4), and the spherical model (for N → ∞). In the present work, we study the Heisenberg model in D = 4, where a long-range-order phase is known to exist at sufficiently low temperatures [2] and, D = 4 being the upper critical dimension, the critical exponents at the phase transition that separates the low-temperature phase, from the disordered, high-temperature one, are equal to their mean-field values, up to logarithmic corrections [3]. Nevertheless, the model encodes non-trivial dynamics: in particular, in this work we study a class of non-local, finite-energy excitations (topological defects) in the low-temperature, brokensymmetry phase and show, by numerical Monte Carlo simulations, that their properties can be successfully predicted using quantum-field-theoretical tools in a continuum formulation of the model [4]. For lower-dimensional systems, the theoretical expectations derived in that work are consistent with exact results for the field theory describing the continuum limit of the Ising model in two dimensions [5] and with numerical results for the XY model in three dimensions [6]. Here, for the first time, we provide evidence supporting these predictions, at the quantitative level, also in a four-dimensional Euclidean spacetime, where these excitations have some analogies with monopole-like states (even though this terminology should be taken with a grain of salt, since the model that we are considering has no gauge symmetry nor gauge fields). As will be shown below, however, we also find that the scaling properties of the parameters of the quantum-field-theory model are different from what was originally conjectured in ref. [4], and challenge a possible interpretation of the topological configurations discussed in this work as physical particles.

Computation setup
We consider the system defined by the Hamiltonian introduced in eq. (1), working on a fourdimensional isotropic, hypercubic lattice of spacing a. We denote the linear extent of the system in each of the three spatial directions as L (with spatial coordinates ranging from −L/2 to L/2), while R is the size of the system in the Euclidean-time direction. We also denote the system temperature as T and define β = J/T ; moreover, we also introduce the reduced temperature t = (T − T c )/T c , where T c is the critical temperature. In this work, we use the most recent estimate of the critical temperature, T c /J = 2.19879 (2), which was obtained in ref. [7] through a sophisticated finite-size-scaling analysis of the Binder cumulant associated with the magnetization. In our Monte Carlo simulations, Markov chains of vector-field configurations are generated by a combination of local heat-bath [8] and overrelaxation [9] updates. For a subset of our runs, we also use non-local, single-cluster updates [10]. Our production runs were executed in perfectly parallel workloads on machines equipped with Intel Xeon Skylake processors.
We study the system both in the high-temperature and in the low-temperature phases. In the high-temperature phase, periodic boundary conditions are assumed in the four directions. We consider the zero-spatial-momentum spin operators and extract the mass m of the lightest physical state by fitting the two-point correlation function to the functional form for sufficiently large values of τ . In the symmetric, high-temperature phase, we run numerical simulations on lattices of sizes (L/a) 4 ranging from 40 4 to 104 4 and for reduced temperatures in the range 0 t 0.06. Conversely, in the low-temperature phase we impose boundary conditions enforcing the existence of a "monopole-like" spin configuration. This is done by taking advantage of the topological nature that this type of excitations have in the continuum: they can be constructed by mapping the spatial boundary of the continuum system at fixed Euclidean time, which has the topology of the S 2 sphere, to the manifold of (classical) vacua of the theory in the broken-symmetry phase, which is also the S 2 sphere. The simplest non-trivial mapping of this type, enforcing the existence of a single, isolated, monopole, is the one that, for the points at the spatial boundary of the system, identifies the direction of s(x) with the direction of the spatial component of x (with respect to the center of the system). We also impose the boundary conditions identifying the direction of s(x) with the direction of the spatial component of x also for all points at the "initial" (x 0 = −R/2) and "final" (x 0 = R/2) Euclidean times. Then, up to an overall normalization, the partition function of the system Z can be identified with the probability amplitude for the monopole propagation from x 0 = −R/2 to x 0 = R/2. We determine the i-th component of the magnetization s i (x i ) and the energy density profile ε(x i ) along the i-th spatial axis through the center of the system (so that ε(x i ) is proportional to the scalar product of s with the sum of the spins on the nearest-neighbor sites, and is normalized to 1 for a uniform field configuration). To increase statistics, we average over the three spatial axes.

Results
The results of our calculations in the symmetric phase confirm the estimate of the critical temperature reported in ref. [7]. As an example, figure 1 shows our results for the two-point, zeromomentum correlation function defined in eq. (3) on the smallest (L 4 = (40a) 4 , left-hand-side panel) and on the largest (L 4 = (104a) 4 , right-hand-side panel) lattices. In particular, the results of our Monte Carlo calculations close to T c (shown as black symbols) indicate that the mass of the lightest physical state that propagates in the theory, vanishes in the thermodynamic limit. This statement is made more quantitative by the fit results for ma for T ≈ T c , which are reported in table 1. The estimated uncertainties on all results are always of the order of 1% or less. We note that, while essentially all of our data for L/a 50 appear to be consistent with a simple, 1/L-decay and with the expected exactly zero asymptotic value for L/a → ∞, it is worth mentioning that the precise form of finite-size corrections in this model is not yet a completely settled issue (see, e.g., the discussion in the recent ref. [7] and in the references therein). Moreover, it should be emphasized that a fully systematic study of the value of ma in the thermodynamic limit and very close to criticality would also require taking into account subleading corrections, which are neglected in eq. (4): these include, for instance, short-distance corrections to r that depend on the lattice geometry [11], contributions to G(τ, R) from additional periodic images of the zero-momentum operator, terms related to heavier excitations (which are exponentially suppressed with respect to the lightest mode), etc. While their discussion lies beyond the scope of this work, we remark that our present analysis of the G(τ, R) results provides a reliable and robust way to extract the mass of the lightest physical excitation of the theory for small but nonvanishing values of t. Following this strategy, we then evaluate the mass in units of the inverse  lattice spacing ma at fixed t and for different lattice sizes, obtaining the results extrapolated to the thermodynamic limit through a fit to We note, in particular, that for t ≈ 0 the thermodynamic-limit value of am = −0.0058(37) is compatible with 0 within less that two standard deviations. In addition, we also note that the systematic uncertainty associated with the choice of the functional form in eq. (5) can be (roughly) estimated by studying how the thermodynamic-limit value of am varies, when a different functional form is chosen; if one includes an additional term a 2 k 2 /L 2 on the right-hand side of eq. (5), the extrapolated value of am in the large-volume limit changes to am = −0.023 (15), indicating that the compatibility of am with zero is robust, and that the systematic uncertainty may be of a size similar to the statistical one. Neglecting logarithmic corrections, the masses extrapolated to the thermodynamic limit are expected to depend on t as where Λ + is a constant with the dimensions of an energy, and a renormalization-group analysis predicts the critical exponent near the Gaußian fixed point to be ν = 1/2. Fitting our results to eq. (6), we obtain the amplitude value m/Λ + = 1.995 (24) and the results shown in figure 2.
We now turn to the results of our Monte Carlo simulations in the broken-symmetry phase at T < T c . In this case, we investigated the spin profile and energy density in the presence of a topological defect induced by the boundary conditions of the system, as described above. The left-hand-side panel of figure 3 shows the spin profile along lines parallel to the three main spatial axes, and touching the edges of the cube at the center of the lattice: more precisely, the plot displays the average value of the i-th component of the spin with the boundary conditions enforcing a topological defect, as a function of the coordinate of the i-th axis (in units of the lattice spacing). This quantity is averaged over the three directions, and is shown for a lattice of spatial sizes L = 90a and extent R = 20a in the Euclidean-time direction, and for different values of the reduced temperature t (denoted by symbols of different colors). Our Monte Carlo results are compared with the analytical predictions, denoted by solid curves, derived in ref. [4]: where z = x i 2M/R, with M , which, according to ref. [4], would represent the mass associated with the topological defect, and the asymptotic value of s i at large distances from the defect core v as fit parameters. All fits are done in the range −25 ≤ x i /a ≤ 25 range to avoid systematic effects due to the boundaries of the system. We find excellent agreement between the theoretical curves and the numerical results, as shown in table 2.  Table 2: Results for the fits of our numerical results for the spin profile, from simulations with L/a = 90 and R/a = 20, to eq. (7).  7), that was derived in ref. [4] from quantum-field-theory arguments. Right-hand-side panel: the average spin value v far from the defect core, obtained from fits to eq. (7), plotted against the temperature in units of the coupling, from different sets of simulations on lattices with L = 90a and R = 20a (including those displayed in the plot on the right-hand-side panel), and their fit to eq. (8).
In addition, we also observe that the values for v extracted from these fits exhibit the expected scaling when the temperature approaches T c : in the infinite-volume limit, the modulus of the magnetization is expected to scale as v ∝ (−t) β , where (neglecting logarithmic corrections) the critical exponent is expected to take its Gaußian value β = 1/2. This is indeed confirmed by our results for v, which can be successfully fitted by with A v = 1.376 (5) and T c /J = 2.2195 (7), as shown on the right-hand side of figure 3. The slight mismatch between the fitted and the actual value of critical temperature provides an indication of the impact of finite-volume effects and logarithmic corrections on the critical exponent (that were neglected in this fit): their combined effect is below 1%. We also note that the precision of our results for v is sufficient to rule out different values of the critical exponent. For example, fitting our data to linear (in T − T c ) form instead of eq. (8) yields a reduced χ 2 larger than 200. It is interesting to study the scaling of the fitted value for M a, as a function of the parameters of the theory, L, R and t. The analytical calculations presented in ref. [4] are done in the thermodynamic limit, hence we restrict our analysis to lattices whose spatial volume L 3 is sufficiently large, in order to suppress finite-volume effects. More quantitatively, L has to be much larger than the other length scales of the theory, i.e. the inverse of M (implying 1 M L) and the lattice spacing (1 L/a). Both inequalities are satisfied in our data samples. The dependence of M a on R, the lattice extent in the Euclidean-time direction, is more subtle. According to the interpretation of the topological field configuration as a particle excitation [4], M should be interpreted as the mass of the particle, and, as such, should be independent from R (possibly up to discretization effects, for values of R comparable with a). Our fit results, however, indicate that the dimensionless parameter M a extracted from the fits scales approximately proportionally with R, which questions the interpretation of the topological excitation as a physical particle. This is clearly revealed by the values of the spin profile s i (x i ) computed numerically at fixed temperature and for different values of R, an example of which (for T /J = 2.062725) is shown in the plot on the left-hand-side panel of figure 4, which only shows results at x i /a ≥ 0: the data obtained on lattices with Euclidean-time extent R/a = 16, 20, 28, and 36, and spatial sizes L 3R collapse on the same curve (except for the points close to the boundaries of the lattice with L = 60a), which can be described well by eq. (7). In turn, since the latter depends on M only through the M/R ratio, it follows that M is proportional to R, or, equivalently, that the topological excitation is characterized by an approximately constant µ = M/R. Note that, if this is the case, i.e. if µ should be thought of a physical quantity, then in a set of simulations at fixed R/a and at different temperatures, the aM parameter extracted from the fits should scale as aµR = µ(R/a)a 2 ∝ |t| 2ν . Remarkably, this behavior is indeed seen in our numerical data: as an example, the plot on the right-hand side of figure 4 shows the results for aM obtained from lattices with R/a = 36 (red diamonds), R/a = 28 (violet squares), and R/a = 20 (black circles), at various temperatures (and for L/a = 90). The data set corresponding to each of the three different values of R/a are fitted to the form (with the critical temperature fixed to its value computed in ref. [7]), and, for the three cases, we always find values of the E exponent very close to 1, i.e. to twice the value of ν predicted in the Gaußian approximation. The fitted curves are shown as dashed lines on the right-hand-side plot of figure 4.  Table 3: Results of the two-parameter fits of the data shown in the right-hand-side panel of figure 4 to eq. (9).
Finally, in figure 5 we present our results for ε, defined as the average energy density (in units of J) along the main spatial axes through the edges of the cube at the center of the lattice. This quantity is plotted as a function of the coordinate x i along that axis (in units of a), averaging over the three spatial axes to enhance statistical precision. At any fixed temperature T < T c , in the infinite-volume limit ε(x i ) is expected to tend to a spatially uniform value C at sufficiently large distances. As for the average spin profile, we observe that, for sufficiently large values of L and R, also the average energy-density profile depends only on the temperature (for all points, except those close to the boundaries of the lattice), and, in particular, the numerical results obtained from Monte Carlo simulations at the same T , but for different values of R, collapse on the same curve. An example of this behavior is shown in the plot on the left-hand side of figure 5, displaying the results for ε(x i ) from simulations at T /J = 2.07557, for R/a = 16 (red triangles), R/a = 20 (yellow squares), R/a = 28 (green circles), and R/a = 36 (blue diamonds); the spatial volume is L 3 = (90a) 3 , except for R/a = 16, for which L 3 = (60a) 3 .
The plot on the right-hand-side of figure 5 shows a comparison of our numerical results for the average energy-density profile to the theoretical prediction derived in ref. [4], which is The  3 . The solid curves are obtained from fits to eq. (10), that was derived in ref. [4].
profile, the numerical results exhibit the expected qualitative features (for instance, the energy density has a peak at the center of the system, whose width becomes larger when T → T − c ) and are in excellent quantitative agreement with the theoretical model.

Discussion and conclusions
To summarize, we studied the Heisenberg model on a four-dimensional Euclidean lattice, through massively parallel numerical simulations, based on state-of-the-art algorithms and run on highperformance-computing clusters. Even though D = 4 is the upper critical dimension, the model reveals interesting dynamical features. As a benchmark, our study of correlators of zeromomentum spin operators in the disordered, high-temperature phase confirms the results recently reported in ref. [7] (and based on the analysis of a different observable), yielding, in particular, full consistency with the critical temperature determined in that work. This may also be interpreted as indirect evidence for the non-trivial subleading finite-size correction terms discussed in that reference.
In the low-temperature phase, we implemented boundary conditions enforcing the existence of topologically non-trivial field configurations (which, in a continuum description, can be characterized by a non-zero "winding number" under the second homotopy group π 2 (S 2 )) and studied their propagation from an initial to a final Euclidean time. As discussed in ref. [4], in the continuum limit the dynamical properties of this type of configurations can be studied with analytical tools from quantum field theory, using only general assumptions about the form factors relevant for different observables of the model as input. The predictions derived in this approach are supported by exact analytical solutions in two dimensions [5] and by simulation results in three dimensions [6]. In the present work, for the first time, they have also been quantitatively confirmed in four Euclidean dimensions, to the high level of numerical precision that was possible to achieve through a large set of dedicated simulations, in a range of volumes spanning two orders of magnitude (the number of degrees of freedom on the largest lattices being approximately 2.34 · 10 8 ), and for a fine scan of temperatures close to T c . It is remarkable that the analytical approach advocated in ref. [4] has very strong predictive power, yielding quantitatively accurate expectations for various observables, from a very limited number of unknown parameters only, and for quite different models, based on continuous or discrete degrees of freedom, invariant under Abelian or non-Abelian symmetries, and defined in different spacetime dimensions.
Our simulation results also reveal that the scaling of the parameters (which could be interpreted as "low-energy constants") that appear in the quantum-field-theoretical description of the spin model for T → T − c is non-trivial, and suggest that a straightforward interpretation of these topological configurations as particles is problematic. In particular, the data lead us to the conclusion that the parameter playing the rôle of the particle mass in ref. [4] is actually a quantity proportional to the "duration" (in Euclidean time, in this case) of its propagation. Our numerical evidence for this scaling is twofold: one the one hand, the results of simulations at fixed temperature and at different values of R reveal that the characteristic widths of average spin and energy-density profiles do not grow with √ R. Rather, their squared width appears to saturate at a given value of R/M , so that, even when the topological configuration is allowed to propagate for a longer Euclidean-time interval, it remains spatially localized. On the other hand, we showed that the aM parameter, that we extracted from the fits, scales with the reduced temperature t proportionally to |t| 2ν , which indicates that the appropriate "physical" quantity that characterizes these objects is µ = M/R, rather than M .
These results have interesting implications. In particular, the scaling properties observed in four dimensions can be compared and contrasted with those reported in an analogous numerical study in three dimensions [6], which investigated the vortices in the XY model. In that work, it was found that the vortex has a well-defined mass, which is approximately 2.1 times larger than the mass of the lightest physical excitation in the disordered phase, and which, when expressed in units of the inverse lattice spacing, scales like |t| ν . Moreover, it was also pointed out that those findings provided indirect evidence that Derrick's argument [12] (a theorem in classical field theory implying that non-trivial, regular, static, soliton-like configurations in scalar field theory cannot exist in more than two spacetime dimensions, as their energy is unstable under scale transformations) may be violated at the quantum level. The possibility of an anomalous violation of Derrick's theorem is particularly intriguing, especially in the light of some arguments that have been recently put forward in axiomatic quantum field theory [13].
It should be noted that an anomalous violation of Derrick's theorem would entail far-reaching consequences, in many different areas of physics. As an example with particularly striking phenomenological implications, it is worth remarking that in the context of relativistic astrophysics, Derrick's theorem implies that bosonic stars consisting of particles that are the excitations of real scalar fields do not exist [14]. In fact, typical ways to evade Derrick's theorem consist in relaxing any of its assumptions, e.g. by introducing a local internal symmetry and gauge fields [15] or replacing the real scalar field with a complex one, which may be charged under an unbroken, global, continuous internal symmetry and sustain stationary (i.e. time-dependent and oscillating), rather than static, configurations [16].
Alternatively, one can construct topologically stable, monopole-like field configurations of finite energy in physical systems characterized by an intrinsic cutoff scale: these include, in particular, those relevant for condensed-matter theory. It is, in fact, in this setting that the findings of our study may be particularly relevant: it is also worth pointing out that, while our calculations are carried out in a classical statistical mechanics setting, one could alternatively interpret this model as the lattice regularization of the corresponding quantum theory in three spatial dimensions [17]. The possibility that this theory supports artificial monopole-like objects [18] may have important applications, given that synthetic magnetism is expected to provide a route towards quantum simulation [19] and information communication [20]. For a recent example of experimental work in this area, see, for instance, the study of manganese germanide reported in ref. [21].
Finally, while the present numerical study has been carried out in a conventional, equilibrium Monte Carlo setting, it would be interesting to generalize it to address aspects related to outof-equilibrium dynamics. A possible way to do this, harnessing well-established mathematical theorems in non-equilibrium statistical mechanics [22], was recently discussed in ref. [23], which presented a non-perturbative study of gauge theories subject to boundary conditions enforcing a non-zero minimum action field configuration, and determined the value of the physical running coupling, in a well-defined renormalization scheme [24], by measuring the response of the system under a sequence of quantum quenches that "deform" the field configurations at the boundary. Repeating a similar type of computation for the setup considered in the present work could shed further light onto the properties of the topological configurations that we studied here, and onto their potential applications in condensed-matter systems.