Spectral Damping without Quasiparticle Decay: The Dynamic Structure Factor of Two-Dimensional Quantum Heisenberg Antiferromagnets

Two-dimensional Heisenberg antiferromagnets play a central role in quantum magnetism, and yet the nature of dynamic correlations in these systems at finite temperature has remained poorly understood for decades. We solve this problem by using a novel quantum-classical duality to calculate the dynamic structure factor analytically and find a broad frequency spectrum despite the very long quasiparticle lifetime. The solution reveals new multi-scale physics to which the conventional renormalization group approach is blind. We also challenge the common wisdom on static correlations and perform Monte Carlo simulations which demonstrate excellent agreement with our theory.

The dynamic structure factor encodes the fundamental physical processes involved in the response of a system to an external probe, and is the most common experimental observable in studies of magnetic systems, determined, for example, using inelastic neutron [1], or resonant X-ray scattering spectroscopy [2]. However, theoretical analyses of these processes which are both quantitatively accurate and physically insightful can be elusive. The two-dimensional quantum Heisenberg antiferromagnet (2DQHA) plays an important role in the field of quantum magnetism precisely because of the theoretical challenges it poses in addition to its descriptive power: First, the model describes the parent compounds of cuprate high temperature superconductors [3]. Second, while the model supports long-range order at zero temperature, order is destroyed at any finite temperature [4]. Because of the importance of thermal fluctuations, 2DQHAs manifest highly non-trivial classical and quantum long-range dynamics which are not fully understood [5,6]. The nature of quantum critical points to and from quantum spin liquid phases is also a problem of intense theoretical interest (see Ref. [7] for a review). Somewhat surprisingly, the physics of thermal fluctuations in isotropic 2DQHAs is closely related to the zero temperature quantum Lifshitz phase transition between antiferromagnetically ordered states and a spin liquid phase in systems with long-range frustrated interactions (e.g., the J 1 -J 3 model) [8].
In their seminal work, Chakravarty, Halperin and Nelson used the O(3)-symmetric nonlinear σ model (NLSM) to describe the long wavelength physics of 2DQHAs at low temperature and argued that the spin-spin correlations in the so-called "renormalized classical" regime are essentially classical in nature [9]. Crucially, their analysis relied on a quantum-classical mapping which integrates out all dynamics of the quantum model. Consequently, this approach allowed the authors to derive a scaling form for the static structure factor but not for the dynamic structure factor. Later studies of the dynamic structure factor raised surprising questions. First, a direct per-turbative calculation of the magnon decay rate due to scattering from the thermal bath predicted the dynamic structure factor should have a very narrow linewidth [10]. Similarly, a 1/N expansion of the O(N ) NLSM predicted a narrow quasi-Lorentzian frequency distribution [11]. In contrast, classical time-dependent numerical simulations showed a broad dynamic structure factor [12], and it has so far remained unclear how to rigorously reconcile this apparent contradiction.
We resolve the long-standing discrepancy in this paper with a novel analytical calculation of the dynamic spin structure factor of the isotropic O(N ≥ 3) NLSM at finite temperature. In recent works [8,13], we demonstrated that infrared-divergent fluctuations-either thermal or quantum-lead to the emergence of a new quantumclassical duality: when an external probe interacts with the system, it creates a classical field which contains an infinite number of quanta with finite total energy. This concept actually originates from particle physics where it was first developed by Bloch and Nordsieck to solve the problem of the radiation field of accelerating electrons [14]. Since the O(2) NLSM is exactly solvable, we were able to rigorously show that despite the infinite quasiparticle lifetime, the dynamic structure factor at nonzero temperature is broad and non-Lorentzian [13].
The O(N ≥ 3) models are not exactly solvable, and hence, the diagrammatic expansion we derived in Ref.
[13] is not applicable. However, we leverage the concept of the infrared catastrophe to develop a new analytical technique and use it to show for the first time that the dynamic spin structure factor of the O(N ) quantum NLSM at finite temperature is very broad and non-Lorentzian. Our analysis demonstrates that this broadening is not due to short-lived quasiparticles but instead is due to the radiation of multiple spin waves by the external probe. With this result, we also obtain the static structure factor by integrating over frequency and find similarities with the scaling form known in the literature [6,9,11]. However, our static structure factor has a different tem-arXiv:2007.11827v2 [cond-mat.str-el] 26 Jul 2020 The dominant contribution to the dynamic structure factor is the emission by the probe (dashed) of a quasiparticle with energy ω k (dark wavy) and a second "soft" particle with energy |ω − ω k | ω k (light wavy).
perature dependence which originates from the underlying quantized nature of the highly-classical radiation field. Fortunately, unlike the dynamic factor, the static structure factor can be calculated numerically using path integral quantum Monte Carlo-referred to hereafter as Monte Carlo (MC). Therefore, to confirm our result for the static structure factor we also perform extensive MC simulations of the O(3)-symmetric NLSM and find excellent agreement.
Formalism.-The long-range dynamics of 2DQHAs at low temperature can be described by the O(3) NLSM with Lagrangian L = (ρ 0 /2)[c −2 (∂ t n 0 ) 2 − (∇ n 0 ) 2 ], n 2 0 = 1, where ρ 0 and n 0 are the spin stiffness and staggered magnetization order parameter, respectively, defined at the ultraviolet scale Λ 0 ∼ π/b, and b is the lattice spacing [9,15]. Quantum fluctuations are ultraviolet-divergent as a power of the momentum scale, and at a scale Λ 1 Λ 0 corresponding to several lattice spacings, reduce the order parameter down to n 0 = | n 0 | < 1. To describe physics at the scale Λ 1 , the quantum fluctuations can be integrated out, leading to the low energy Lagrangian where ρ and n are the spin stiffness and order parameter normalized at Λ 1 . Both ρ and n 0 as a function of the dimensionless coupling constant g = c/ρ 0 b can be calculated numerically using MC, and those of the O(3) NLSM are shown in Fig. 1(a). In this paper, we address the regime g < g c ≈ 1.46 which describes 2DQHAs [9]. For the sake of generality, we consider from hereon the O(N )-symmetric model in terms of n = (n 1 , n 2 , . . . , n N ), with N ≥ 3, and use units of = c = b = 1. At zero temperature, the ground state of the model has longrange collinear antiferromagnetic order-n = constwhich spontaneously breaks the O(N ) symmetry. How-ever, the Hohenberg-Mermin-Wagner theorem guarantees the destruction of long-range order at any finite temperature [4,16]. Despite this, at sufficiently low temperatures T ρ, the system remains ordered on scales up to the exponentially-large correlation length ξ ∝ exp[2πρ/(N − 2)T ] [6,9,11,17]. This separation of scales implies a notion of quasi-long-range order and allows for a perturbative treatment of the NLSM on momentum scales Λ satisfying ξ −1 Λ ≤ Λ 1 . For momentum scales on the order of temperature to fall within this range, it suffices for T ρ. The effects of fluctuations on scales Λ ∼ T can be determined within the leading order of perturbation theory. However, there are two types of contributions governing the physics of fluctuations on scales Λ < T : (i) Renormalization group (RG) "running" of physical parameters due to interactions occurring at the same scale. (ii) Beyond RG contributions originating from multi-scale interactions.
The RG contributions are well-understood [5,6,9], so we only summarize the general principles. The unit vector constraint of Eq. (1) generates interactions between the components of n, leading to renormalization of the spin stiffness-ρ → ρ Λ -and fields-n → n Λ = Z 1/2 Λ n, where Z is the quasiparticle residue. To one-loop accuracy at the momentum scale Λ < T < Λ 1 , The ultraviolet cutoff for the fluctuations in (2a) is the temperature T rather than Λ 1 due to the bosonic statistics of the quasiparticles; this is an important quantum correction to classical thermodynamics [18]. In Supplemental Material (SM) Sec. I, we give a derivation of (2) and show that higher-loop contributions are negligible when T ρ [19]. Dynamic structure factor.-The dynamic structure factor (DSF) is the Fourier transform of the order parameter correlation function n i (r, t)n i (0) [20], and is independent of the polarization index i due to the absence of long-range order at finite temperature. Expanding in a spectral representation in the basis of excited quasiparticle Fock states |α and |β yields [21], where Z is the quantum partition function, and ω α and k α are the energy and momentum of the state |α . In this paper we always work with ω > 0, since (3) implies that S(k, −ω) = e −ω/T S(k, ω). First, we account for RG contributions to the DSF by renormalizing the fields in (3) at the scale of the incoming momentum k, so that n i k = n i /Z 1/2 k ; since ξ −1 k < T , local order exists at this scale. The absence of long-range order means that the single magnon intermediate state does not contribute to (3) [13]. Therefore, we eliminate the n 2 = 1 constraint by writing the order parameter field as n k = ( π k , √ 1 − π 2 k ), where π = (π 1 , . . . , π N −1 ) are small transverse fluctuations, and use the component of n in the direction of local order-n N k = √ 1 − π 2 k -to compute the DSF. However, this approach assigns all dynamics to the directions transverse to the local moment, and hence, does not respect the O(N ) symmetry which must remain unbroken at finite temperature. To restore symmetry, we rotationally average the DSF over all N polarizations by multiplying (3) by (N − 1)/N [6,11]. Therefore, suppressing Boltzmann factors and δ functions for notational clarity, the DSF is The leading contribution is then obtained by expanding √ 1 − π 2 k 1 − π 2 k /2 → − π 2 k /2 and using Fermi's golden rule to calculate the probability of two magnon radiation. More precisely, if ω > ω k the external probe excites two quasiparticles [see Fig. 1(b)], and if ω < ω k one quasiparticle is emitted and a second is absorbed. When |ω − ω k | . = |∆| ω k both processes have the same contribution to the sum over initial and final states (see SM Sec. II [19]): However, by examining the structure of the phase space integral yielding (5), we find that one emitted quasiparticle will have energy ∼ ω k and the other will have energy ∼ |∆| ω k . Hence, the two magnon intermediate state is an inherently multi-scale process and contributions at the "soft" scale Λ ∼ |∆| are not properly accounted for; conventional RG is not sufficient to describe the process accurately. We understand from our exact solution of the O(2) NLSM that the physics of the soft scale is characterized by an interplay between thermal fluctuations and the radiation of arbitrarily low energy quasiparticles, the "probabilities" of both of which are logarithmically infrared-divergent; this divergence implies that no finite number of quasiparticles can be excited by the probe [13]. Therefore, the "second quasiparticle" with energy |∆| emitted/absorbed by the probe is actually accompanied by a classical radiation field containing infinitely-many quanta; the total energy of the soft magnon and the classical field is |∆|.
To account for beyond RG contributions to the DSF, it is important to allow for the running of the quasiparticle residue and spin stiffness down to the soft scale. This is simple since: (i) Soft radiation factorizes-to one-loop order-from the emission of the "hard" quasiparticle. (ii) The classical radiation field cuts off the divergent static thermal fluctuations below the soft scale. This leads to a correction to (5)-denoted without a tilde-with one factor of the spin stiffness modified ρ k → ρ ∆ and an additional field strength factor Z k /Z ∆ -since the running of Z Λ is multiplicative. From this, we obtain the DSF for the O(N ) NLSM in the regime ξ −1 |∆| ω k T : It is common to express the structure factor in terms of appropriate length/time scales. In the present case, the only length scale is so that in terms of the one-loop expressions for Z ∆ , ρ k and ρ ∆ given by (2), the full form of (6) is Of course, this result assumes λ|ω − ω k | 1, and hence, represents a very broad frequency distribution decaying slower than 1/|∆|. The limit N → 2 reproduces the solution obtained in Ref.
[13] using a direct summation of diagrams. While we used a very different technique in Ref. [13], the hierarchy of multi-particle contributions was compatible with the present approach. Since the case N > 2 is not exactly solvable, in this work we used the running of parameters to correctly account for the beyond RG contributions.
We have so far neglected the lifetime of quasiparticles in the O(N ≥ 3) NLSM. The dominant decay process for an on-shell quasiparticle with energy ω k T is 2 → 2 Raman scattering from a particle in the thermal bath, which leads to the well known inverse lifetime [10,11,22], Importantly, in our regime of interest (ξ −1 ω k T ρ) Γ k /ω k is an O(T 2 /ρ 2 ) small quantity. It is then clear from the analysis in this section that when Γ k < |∆|, radiative broadening of the DSF due to multiple emissions/absorptions dominates over 1/|∆| 2 Lorentzian lifetime broadening. As a side note, since Γ k ∝ N , radiative broadening may be hidden in a 1/N expansion around N = ∞. However, for N not much larger than O(ρ/T ), the region |∆| < Γ k remains very narrow. Regardless, the decay of the hard particle cannot be neglected near resonance. In particular, it serves to regularize the nonintegrable singularity at ω = ω k in (8). Since the calculation does not offer any new insights, it is presented in SM Sec. II [19]. In Fig. 2(a) we compare the DSF calculated in this paper [SM Eq. (S29)], to a Lorentzian lineshape with the same ω integrated spectral weight. Clearly, the resonant response of the DSF is greatly suppressed compared to the Lorentzian, with significant spectral weight shifted to the tails of the frequency distribution.
Equal-time correlations.-The static structure factor can be calculated directly from the DSF by integrating over frequency in the interval −ω k < ∆ < ω k (see SM Sec. III [19]): The equal-time order parameter correlation function, which is N times the Fourier transform of (10), reads The static structure factor (10) has the same functional k dependence as the well-known scaling form [6,9,11]. However, (10) contains log(λk) instead of log(ξk) in those Refs., where ξ is the correlation length 11,17]. The replacement ξ → λ leads to a particularly drastic difference for the case N = 3, where the pre-exponential factor of the correlation length is temperature-independent.
To confirm our results (10) and (11), we performed MC simulations of the O(3) NLSM and measured the equaltime order parameter correlation function. The zero temperature spin stiffness ρ and staggered magnetization n 0 presented in Fig. 1(a) have been calculated on a 64 3 size lattice. To measure the correlation function we used lattices with L x = L y = 512 and L β = 4, 6, 8 imaginary time slices which correspond to different temperatures T = gρ 0 /L β . In Fig. 2(b) we present the MC correlation function for dimensionless coupling g = 1, corresponding to ρ/ρ 0 ≈ 0.504 and n 0 ≈ 0.673. The solid lines show the theoretical prediction (11) for N = 3; note that the theory has no adjustable fitting parameters. At r 2-in units of the lattice spacing-deviations from theory are due to the dominance of ultraviolet quantum fluctuations on short length scales, and at r 128, finite-size effects from the periodic boundary conditions become important. Since the parameters used lie within the domain of validity of the theory, we find excellent agreement between our theoretical predictions and the MC simulation data. Further data and analysis can be found in SM Sec. IV [19]. The dashed lines show the correlation function (11) with λ replaced by ξ and disagree very clearly with the MC simulations.
To avoid misunderstanding we note the following: (i) The correlation length is defined in terms of the exponential decay of correlations on large length scales n(r) · n(0) ∼ e −r/ξ when r ξ. (ii) In this work, we are operating in the opposite limit r ξ. We are not claiming that the well known expression (12) for the correlation length is incorrect. However, we claim that correlations on shorter scales are characterized by the parameter λ, and not ξ.
Summary.-We have calculated for the first time the finite temperature dynamic structure factor of the 2D O(N ) quantum nonlinear σ model in the regime describing a Heisenberg antiferromagnet. The dynamic structure factor displays a very broad frequency distribution which decays slower than the first power of the detuning from resonance. Since the quasiparticle lifetime remains very long, it is irrelevant to broad tails of the spectrum. Instead, the broadening is driven by the emission and absorption of multiple soft excitations by the probe. To perform this calculation, we developed a new analytical technique which accounts for both conventional single scale renormalization group contributions and "beyond RG" effects from multi-scale physics. We expect this method to be applicable to studying the dynamics of a wide range of finite temperature interacting quantum field theories.
Using our new result for the dynamic structure factor we also calculated the static structure factor and found agreement of the functional momentum dependence with the previously known result. However, we predicted a significant modification of the characteristic length scale of correlations in the so-called scaling regime. This result implies an important correction to the temperature dependence of static correlations from the bosonic statistics of the quasiparticles. To confirm this prediction we performed extensive path integral quantum Monte Carlo simulations and demonstrated perfect agreement between the numerical data and our analytical formula. To the best of our knowledge, this is also the first Monte Carlo study of correlations in the scaling regime.
Acknowledgments In the main text section Formalism we discuss how quantum and thermal fluctuations are taken into account via renormalization. Here we provide more details on how the RG equations (2) can be derived.
The O(N ) nonlinear σ model (NLSM) in 2 + 1 dimensions, normalized at the scale Λ 0 ∼ π/b where b is the lattice spacing, is given by the Lagrangian where ∂ µ = (c −1 ∂ t , ∂ x , ∂ y ) and ρ 0 is the spin stiffness. From hereon we set c = 1. The unit vector constraint can be eliminated explicitly by writing n 0 = ( π 0 , σ 0 ) = ( π 0 , 1 − π 2 0 ). The Lagrangian in terms of the transverse fluctuations π 0 is Expanding around the zero temperature state of spontaneous symmetry breaking σ 0 = 1 and π i 0 = 0 yields where the ellipsis denotes terms of O(π 6 0 ). The interactions have the effect of renormalizing the spin stiffnessρ 0 → ρ Λ -and the fields-n 0 → n Λ . = Z 1/2 Λ n 0 . To see this, we calculate the self-energy by performing a one-loop decoupling of the quartic term where π 2 Λ are the fluctuations of the π fields with momenta in the interval (Λ, Λ 0 ). At the scale Λ 1 Λ 0 corresponding to several lattice spacings, the fluctuations reduce the effective length of the σ 0 component down to n 0 = | n 0 | < 1. Therefore, we require that the renormalized fields satisfy σ Λ1 = 1; the renormalized ground state should have the same form σ Λ1 = 1, π i Λ1 = 0. Far away from the quantum critical point, the first order perturbative calculation gives and hence, In practice, ρ and n 0 are calculated numerically, which we do using path integral Monte Carlo (see Section IV). By integrating out the ultraviolet quantum fluctuations, we can obtain the low energy Lagrangian normalized at Λ 1 L = 1 2 ρ(∂ µ n) 2 , n 2 = 1.
Turning to the case of finite temperature, we note that long-range order is destroyed by thermal fluctuations and no state of spontaneously broken symmetry exists [4,16]. However, the correlation length ξ ∝ exp[2πρ/(N − 2)T ] remains exponentially large in the low temperature regime T ρ [9,11,17]. Therefore, on momentum scales ξ −1 Λ < Λ 1 , we can apply the same analysis as above by expanding the low energy Lagrangian (S7) around a locally ordered state. The thermal fluctuations π 2 Λ can be evaluated directly when Λ < T < Λ 1 : Within the Matsubara imaginary time formalism, this same result can be obtained by noting that the dominant contribution comes from the zero Matsubara frequency in the low temperature T ρ regime. This statement is equivalent to the common wisdom that the low temperature regime of the square lattice Heisenberg antiferromagnet is characterized by classical static (zero Matsubara frequency) thermal fluctuations [9]. However, we note that the ultraviolet cutoff of the logarithm in (S8) is imposed by the Bose occupation factor in the momentum integral. If the thermal fluctuations were purely classical the ultraviolet cutoff would be Λ 1 [18]. Instead, the quantization of the field leads to a correction in 2 + 1 dimensions.
The renormalization group (RG) equations governing the flow of ρ and Z from the scale Λ 1 down to Λ then follow directly from Eqs. (S5), (S6) and (S8) by considering infinitesimally small l = log(T /Λ): When N ≥ 3, there is a non-trivial renormalization group flow. Integrating the system of equations (S9) yields These results are very well known, and follow directly from expressions given in Refs. [9,11].

B. Two-loop contributions to renormalization
In the main text, at the end of the Formalism section, we claimed that higher-loop order corrections to the spin stiffness and order parameter could be neglected at scales ξ −1 Λ < Λ 1 . Here, we consider the well-known RG flow equations for the 2D NLSM [5, 23], which describe the interactions of classical fluctuations in the (2 + 1)D model [9]. Since all calculations should be consistent to leading order with perturbation theory, we again use the temperature T instead of Λ 1 as the ultraviolet normalization point of the RG flow. The two-loop equations in terms of the dimensionless temperature t = T /2πρ are [23], Re-writing these differential equations in terms of ρ, we find the exact solution where W (x) is the inverse function of W e W = x-otherwise known as the Lambert W function or product logarithm-which has two branches for real x. The −1 branch satisfies W −1 (xe x ) = x if x ≤ −1, which guarantees that the initial condition of the differential equation is satisfied.
In the regime we are interested in-T ρ and Λ ξ −1 -the argument of the W function in (S12a) is close to zero. Using the asymptotic expansion W −1 (x) = log(−x) − log(− log(−x)) + O(1), we find Therefore, neither the quasiparticle residue nor the spin stiffness are modified-to a good degree of accuracy-by two-loop contributions in the low temperature regime.
This result may surprise, given that it is well known that the correlation length ξ is heavily modified by two-loop corrections [6,9]. However, these differences only appear in the far infrared limit. First, observe that the one-loop spin stiffness (S10a) vanishes at the scale Λ = λ −1 , where However, since W −1 (−1/e) = −1, the two-loop spin stiffness vanishes at the scale Λ = Ξ, where For comparison, the exact (inverse) correlation length is [6], where Γ(x) is the Gamma function. Most importantly, for the case N = 3 both ξ −1 and Ξ have a temperature independent pre-exponential factor-to leading order in T /ρ. This is a considerable difference compared to the temperature dependence of λ. Obviously, any perturbative calculation like RG is not valid at and beyond the strong coupling scale. However, it is clear from the above analysis that the behavior of the spin stiffness and order parameter are heavily modified by two-loop contributions when approaching that scale.
Finally, we also acknowledge that there are two-loop corrections to the speed c. However, the renormalized speed as reported in Ref. [11] is only modified at O(T 2 /ρ 2 ) when Λ ξ −1 .

II. DYNAMIC STRUCTURE FACTOR
In the main text section Dynamic structure factor, we start from the definition of the dynamic structure factor (DSF) as the Fourier transform of the correlation function and derive an expression which includes both conventional RG contributions as well as multi-scale physics. Here, we provide more of the details of this derivation. First, we consider the radiation-dominated region, and then the lifetime-suppressed resonant peak.
A. Derivation of main text equation (8) Given the absence of long-range order, we can write where the average is taken over the thermal ensemble. The matrix element can be expanded in a spectral representation which turns the integral into a sum over initial and final states [21], where Z is the partition function, |α , |β are excited Fock states, ω α , ω β and k α , k β are the energy and momentum of those states, respectively, and ω αβ = ω α − ω β , and similarly for k αβ . Before proceeding, we recall that in the absence of long-range order, the DSF will contain no elastic spectral weight; the spectrum will not have a Bragg peak at k = 0. Then, as in the main text, we account for conventional RG contributions to the spectrum by re-writing the spectral expansion in terms of the fields renormalized at the momentum transfer from the external probe k FIG. S1. The transverse response of the structure factor also contains infinitely-many multiparticle states. These are obtained from the different ways of cutting diagrams (vertical dashed lines). The bold line represents the exact π i propagator, normal lines represent the bare π i propagator, and horizontal dashed lines represent the source.
If we have ξ −1 k T , then local order exists on these momentum scales. This allows us to eliminate the unit vector constraint n 2 = 1 by artificially breaking the O(N ) symmetry and writing the order parameter as n k = ( π k , √ 1 − π 2 k ). However, the DSF must remain rotationally invariant, and in particular, spectral sum rules must be satisfied. Therefore, we account for the fact that we have "broken" the symmetry by rotationally averaging the DSF over all N polarizations; this amounts to multiplying by a factor of (N − 1)/N [6, 11]: Naively, this expression appears to contain a Bragg peak contribution, seemingly in contradiction to our earlier point. However, we re-emphasize that order only locally exists on momentum scales 0 < ξ −1 < k; it is incorrect to use (S20) at 0 ≤ k < ξ −1 .
On the one hand, the DSF must be rotationally invariant. On the other hand, the n N = √ 1 − π 2 component contains all even powers of π, while the transverse components are linear in π i . There is no actual contradiction here. To see this, consider Fig. S1, where we illustrate the simplest loop corrections to the interaction of the source with a transverse π i component. The first two diagrams correspond to the one-and two-loop contributions to the self energy of the emitted quasiparticle. However, the third diagram shows that the interactions (and self interactions) between the π i components allow the probe to create three real and on-shell particles via an intermediate virtual state. The Feynman rules for the interaction vertex in (S3) are given in textbooks (e.g., Ref. [24]). In particular, the propagator is G(q) = iρ −1 /q 2 , and the amplitude for an off-shell particle with three-momentum q to decay into three on-shell particles is A(q) = −iρq 2 = 1/G(q). Therefore, the total quantum amplitude for the intermediate state is G(q)A(q) = 1. Essentially, the virtual particle "contracts" the interaction to a single point. It is straightforward to see that all higher-order interactions in the expansion of (S2) will lead to the same amplitude for similar multiparticle emissions. We emphasize that this means that the transverse components also lead to the emission of infinitely-many quasiparticles; whether this infinity is "odd" or "even" is irrelevant, implying the preservation of O(N ) symmetry [13].
It is, however, much simpler to perform calculations using the n N longitudinal component. We obtain the first non-trivial contributions by expanding the square root inside the matrix element √ 1 − π 2 k 1− π 2 k /2: (i) The emission of two particles (ω > ω k ). (ii) The emission of one and absorption of the other (ω < ω k ). Since these two processes are mathematically similar, we focus on emission (of any polarization), where the total contribution to the sum over initial and final states of this form is given by the two particle phase space integral where is the effective finite temperature two particle emission matrix element (squared) obtained after performing the Gibbs averaging in (S20); the factor of 1/2 from the expansion coefficient of the square root is canceled by the number of Wick contractions of the matrix element. Note that due to our earlier choice of renormalization scale, the spin stiffness is evaluated at k. The integral can be evaluated exactly after expanding the Bose occupation factors to leading order in T /ω,Ĩ 2 (k, ω) and it can be observed that when |∆| ω k , the dominant contribution comes from regions of phase space where one particle has energy ∼ ω k and the other has energy ∼ |∆|; the process involving absorption which occurs for ω < ω k reduces to the same expression in the limit |∆| ω k .
As discussed in the main text, this large difference in the momentum scales of the two quasiparticles means that contributions at the soft scale are not accounted for. First, observe that we can perform a post-hoc simplification of (S21) using our knowledge of the momentum distribution, and find that the integral factorizes as into high and low energy processes. This implies that all corrections to the external particle legs will also factorize, allowing us to account for the running of the parameters of the two particles independently. First, we allow the spin stiffness to run with the momenta of the particles ρ k → ρ ki and account for the running of the quasiparticle residue from the normalization point k to the momenta of the particles, which gives a factor of Z k /Z ki for each particle. We note that there are also two-loop corrections to the source vertex which do not factorize. However, this will be a higher-order effect, so we neglect it. Therefore, the "beyond RG" version of (S23)-denoted with no tilde-will be We know from the exact solution of the O(2) NLSM that the classical radiation field created by the source cuts off infrared-divergent static thermal fluctuations with momenta smaller than that of the soft particle [13]. Therefore, all leading-order beyond RG contributions to the DSF are accounted for by the running of the quasiparticle residue of the soft leg. Hence, the DSF is given by which can also be written in terms of the length scale λ = (1/T ) exp[2πρ/(N − 2)T ], yielding Eq. (8) of the main text. It is simple to see using this expression that the exact solution of the O(2) DSF is reproduced by taking the limit N → 2; the first square bracket → 1, while the second square bracket exponentiates → (|∆|/T ) T /2πρ .

B. Lifetime broadening near resonance
At the end of the main text section Dynamic structure factor, we point out that the finite lifetime of quasiparticles becomes relevant in the narrow region of the frequency spectrum |∆| < Γ k , where Γ k is the Raman scattering rate for a particle with energy ω k . The leading contributions come from the imaginary part of the two-loop self energy diagram shown in Fig. S2. The well-known expression for the scattering rate is [10,11], where we prefer to use the convention that Γ k is the full width at half maximum of the imaginary part of the single quasiparticle Green's function; this gives a factor 2 difference compared to the expressions reported in Refs. [10,11]. In our regime of interest ξ −1 |∆| ω k , it is clear that the scattering rate of the soft quasiparticle is negligible compared to that of the higher energy particle. Therefore, the lifetime broadening of the DSF will be dominated by the decay of the hard particle.
The lifetime has the effect of "broadening" the energy conserving δ function in (S25a) However, there is an important subtlety in accounting for contributions from different pieces of phase space. Without lifetime broadening, the following cases are possible: (i) If ω > ω k > 0, two particles are emitted.
(ii) If ω k > ω > 0, one particle is emitted with energy ω k and one is absorbed with energy |∆|.
(iii) If 0 > ω > −ω k , one particle is absorbed with energy ω k and one is emitted with energy |∆|.
In principle, with account of lifetime broadening, any of these processes can occur for any values of energy and momentum transfer from the source. However, since we assume |∆| ω k , we can safely assume no mixing between the positive and negative frequency branches of the spectrum. However, for ω > 0, we must allow for mixing between processes (i) and (ii). Therefore, we generalize the integral (S25a) to give us the full form of the DSF which is plotted in Fig. 2(a) of the main text. Note that we must retain an infrared momentum cutoff for this expression. Given that we have already established that the characteristic momentum scale of the DSF in this regime is λ −1 , we use this as the cutoff. For process (i), the two emitted particles are indistinguishable bosons, but we distinguish between them, so we must impose the ultraviolet cutoff ω k to avoid double counting states. For process (ii), the dominant contribution comes from the absorption of particles with energy < ω k . Finally, we note that in the limit |∆| Γ k , (S29) reduces to (S26).

III. EQUAL-TIME CORRELATIONS
In the main text we present the static structure factor and equal-time correlation function. Here, we simply provide more detail on how the calculations are performed. The static structure factor is defined in terms of the dynamic structure factor by the relation Since the main text Eq. (8) has a non-integrable singularity at ω = ω k , we must use (S29) to compute the static structure factor. Additionally, the spectral representation (S18) implies S(k, −ω) = e −ω/T S(k, ω) S(k, ω) when ω T , which allows us to correctly include the negative frequency portion of the spectrum. Therefore, we have We can also check the total sum rule. Since the dynamic structure factor we derived was valid for ω T , we should integrate (S31b) up to T : Therefore, summing up over the N polarisations, we recover the correct normalization of the order parameter. We note that this sum rule is not satisfied if all parameters are normalized at the same scale-either ω k or |∆|. This observation validates our approach to including multi-scale physics.
The equal-time correlation function then follows directly from (S31b) by taking the Fourier transform. When r λ, it is straightforward to find, to logarithmic accuracy, since the complex exponential will oscillate rapidly and average to zero when k 1/r.

IV. PATH INTEGRAL QUANTUM MONTE CARLO SIMULATIONS
In the main text, we present a subset of our measurements of the equal-time correlation function using path integral quantum Monte Carlo simulations. Here we summarize the details of our simulations-the Monte Carlo update algorithm we implement and how we measure physical observables-and also present a larger selection of data.

A. Heat bath algorithm for O(3) NLSM
The quantum partition function for the O(3) NLSM in imaginary time is given by the path integral [9,25], where ρ 0 is the bare, un-renormalized spin stiffness defined at the lattice spacing b, n = (n (x) , n (y) , n (z) ), and the δ function in the integration measure enforces the unit vector constraint at every point in space. Discretizing the action over a uniform simple cubic lattice with spacing b yields [25], where x i = (cτ i , x i ), g = c/ρ 0 b = L β T /ρ 0 , L β is the size (in number of lattice spacings) of the imaginary time dimension, and the summation is over pairs of nearest neighbors. From hereon, we set = ρ 0 = b = 1. In these units, the bare coupling constant g = c = L β T .
We performed path integral quantum Monte Carlo simulations of the O(3) model by implementing a heat bath algorithm following Ref. [25]. To summarize: (1) Initialize the lattice in a uniformly magnetized grid.
(2) To update a lattice site at position x i , calculate the local action so that the probability density for the vector n(x i ) to lie in some element of solid angle is where C is the normalization constant of the distribution, and θ is measured from the axis directed along ω.
(3) Generate a new configuration for n(x i ) by picking θ and ϕ from this distribution, convert from local to crystal axis coordinates, and then update.
Defining a "sweep" of the lattice to be an update of every lattice site once, we allowed 2500 sweeps for the system to thermalize before starting measurements. We then performed 50,000 sweeps, measuring once every 10 sweeps to minimize correlations between measured configurations; we estimated a correlation time from measurements of the average action per site to be ≈ 2 -3 sweeps.

B. Measurement methods
To measure the zero temperature staggered magnetization n 0 , we used the standard Monte Carlo estimator where L 3 is the total number of lattice sites (at zero temperature all dimensions are of equal size), and the ensemble average · is estimated by an average over measurements.
The equal-time correlation function measurements were obtained using the formula where e µ is a unit vector along the µ direction, and we averaged over positive and negative displacements along the two equivalent spatial dimensions to improve our measurement statistics; at finite temperature, the imaginary time direction is not equivalent, so is not included in the sum over directions µ. Summing over all lattice sites and division by L β approximates the integral over the imaginary time dimension used to obtain the equal-time correlation function.
To measure the zero temperature renormalized spin stiffness, we adapted the approach described in Ref.
[26], which we summarize here. The spin stiffness measures the response of the system to a twist of the boundary conditions of dimension µ by a relative angle Φ = QL µ . At zero temperature, where E(Q) = −(g/L β ) log Z(Q) is the ground state energy functional in the presence of the twist and Z is the quantum partition function. The twisted boundary conditions can be eliminated by transforming to a "rotating" frame of reference where the twist instead modifies the local interaction: where R is a 3 × 3 rotation matrix in spin space. For a rotation about thex axis in spin space, along direction µ in real space, expanding the action to second order in Q leads to a modification of the energy where we have defined where summation is over lattice bonds directed in the µ direction. However, it is necessary to account for the fact that the direction of the twist in spin space is not generally perpendicular to the local magnetization. Therefore, the spin stiffness is obtained by averaging over the other two twist axes, weighted by 3/2-see the discussion in the main text regarding rotational averaging. At zero temperature, all three Euclidean dimensions are equivalent (L x = L y = L β . = L), so we also averaged over all bond directions to obtain a more accurate Monte Carlo estimator of the spin stiffness: where S is the average of the action (S35).

C. Results & further analysis
In Table SI we present a subset of measurements of the zero temperature spin stiffness and average staggered magnetization on a L x = L y = L β = 64 size lattice. These results are practically identical to measurements on a 32 3 lattice showing that finite size scaling effects are negligible (away from the O(3) quantum critical point g = g c ≈ 1.46). We also present the temperature in units of the renormalized spin stiffness and the length scale λ. Evidently, as the coupling g is increased, reducing the spin stiffness ρ, the relative importance of thermal fluctuations increases.
In Figs. S3 and S4 we present measurements of the equal-time correlation function on L x = L y = 512 and L β = 4, 6, 8 size lattices, for a range of values for the coupling g. The solid lines show the theoretical prediction (S33) for  Table SI. Symbols are Monte Carlo data, solid lines are theory (S33) for N = 3, and dashed lines are theory replacing λ → ξ where ξ is the correlation length (S16). Panel (b) is the same data as Fig. 2(b) in the main text, reproduced here for comparison. All vertical scales are identical. The L β = 4 data is omitted from panel (c) since T /ρ > 1 is outside the domain of validity of the theory. N = 3 with the zero temperature spin stiffness and magnetization measured on the 64 3 lattice. We emphasize that the theory has no adjustable fitting parameters. Evidently, the agreement between the data and theoretical curves is excellent. As stated in the main text, disagreement on short length scales (r 2) is to be expected due to the dominance of ultraviolet quantum fluctuations, and on larger length scales (r 128) finite-size effects originating from our choice of periodic boundary conditions become important. Therefore, we see that the Monte Carlo data agrees perfectly with the characteristic length scale λ, but not at all with the correlation length ξ given by (S16).
Strictly speaking, our theory is valid in the regime T ρ. However, the exponentially-large length scales (both λ and ξ) mean that it is still possible to study correlations at r λ when T ∼ 0.5ρ. For example, consider the case g = 1.25 and L β = 6, shown in Fig. S3(c) in red. Here T /ρ ≈ 0.692 and the theory still agrees quite well with the data. This is because λ ≈ 53,000 remains more than two orders of magnitude larger than the lattice. In contrast, when L β = 4, T /ρ > 1 and λ is just over twice the (linear) size of the lattice, and as expected the theory (S33) did not agree at all with the data [omitted from Fig. S3(c) for clarity] since the temperature is outside the domain of validity.  Table SI. Symbols are Monte Carlo data, solid lines are theory (S33) for N = 3, and dashed lines are theory replacing λ → ξ where ξ is the correlation length (S16). All vertical scales are identical.