Spectral functions and dynamic critical behavior of relativistic $Z_2$ theories

We investigate the dynamic critical behaviour of a relativistic scalar field theory with $Z_2$ symmetry by calculating spectral functions of the order parameter at zero and non-vanishing momenta from first-principles classical-statistical lattice simulations in real-time. We find that at temperatures above the critical point $(T>T_c)$, the spectral functions are well described by relativistic quasi-particle peaks. Close to the transition temperature $(T \sim T_c)$, we observe strong infrared contributions building up. In the ordered phase at low temperatures $(T<T_c)$, in addition to the quasi-particle peak, we observe a soft mode with a dispersion relation indicative of collective excitations. Investigating the spectral functions close to $T_c$, we demonstrate that the behavior in the vicinity of the critical point is controlled by dynamic scaling functions and the dynamic critical exponent $z$, which we determine from our simulations. By considering the equations of motion for a closed system and a system coupled to a heat bath, we extract the dynamic critical behavior for two different dynamic universality classes (Models A&C) in two and three spatial dimensions.


Introduction
When a thermodynamic system comes close to a critical point, one finds anomalous behaviour in a multitude of its properties.Large fluctuations lead to scale-invariant physics, not only in static quantities, but also in dynamic observables such as multi-time correlation functions and relaxation rates [1].These observables can be measured by a range of experiments, e.g. in condensed-matter applications like ultrathin films [2], and will likely become relevant in future heavy-ion collision experiments closing in on the QCD critical endpoint (CEP) [3].
Over the last few years, the search for the QCD critical endpoint has gained a considerable amount of theoretical and experimental attention [4,5].Dedicated heavy-ion collision experiments at the Relativistic Heavy-Ion Collider (RHIC), and the future GSI-FAIR and NICA facilities are designed to probe the relevant region of the QCD phase diagram, so the task is set for theorists to predict how the QCD critical point can be located with the obtained data.Since first principle calculations based on lattice QCD are currently unable to address the interesting high-density region at large baryon chemical potential, one has to employ effective theories, making use of prior knowledge about its critical behaviour [6].
The scale-invariant nature of dynamic critical phenomena leads to universal behaviour of microscopically different systems, and to an identical value of the dynamic critical exponent z.But in addition to the well-known universality classes governing static critical phenomena, one must take into account properties of the equations of motion of the system.By factoring in conservation laws and mode couplings, Hohenberg and Halperin [1] developed a scheme of dynamic universality classes, called "Models."In some of these models, the dynamic critical exponent z follows directly from the static critical exponents via simple scaling relations.This is for example the case for Model C, a relaxational model with a conserved density.In general, however, the precise determination of z turns out to be significantly more challenging.Even for the rather simple 2D Ising model with Glauber dynamics (Model A), rather sophisticated methods [7] were needed to obtain a reasonably precise result for z.
In this study we investigate a relativistic scalar field theory in the static universality class of the QCD CEP (Z 2 Ising), yet in a different dynamic one.While the dynamic universality class at the QCD critical endpoint is believed to be that of Model H [8], we will focus here on the conceptually simpler relaxational Models A and C which do not require the explicit treatment of an additional conserved vector field and its dynamic coupling to a conserved order parameter.Nevertheless, we believe that investigating the simpler models provides a first important step towards characterizing the real-time dynamics in the vicinity of a critical point in a systematic way, ready to be extended to more complicated models such as Model H in the future, which will then be of direct relevance to QCD.
Our analysis of dynamic critical phenomena is based on the real-time correlations of the orderparameter field, specifically, its spectral function.Spectral functions are of great interest in many areas of condensed matter, nuclear and particle physics.They play an important role for our understanding of the physics of heavy-ion collisions, the Early Universe and more.Containing the spectrum of quasi-particle, multi-particle and collective excitations of the system in the given channel, they can be used to identify the relevant degrees of freedom e.g. for transport or hydrodynamic descriptions.In addition, transport coefficients are typically also obtained from particular low-energy limits of spectral functions [9].Naturally, close to a critical point, the emergent dynamic critical phenomena have a strong effect on spectral functions.Since slow modes dominate the dynamics, infrared-divergent power laws arise, whose spectral index is controlled by the dynamic critical exponents z and the static critical exponents α, β, γ, δ, and ν.Due to scale invariance in the vicinity critical point, one further expects that, similar to static quantities, the low-energy behaviour of the spectral function is fully described by universal scaling functions.
Of course, to capture such highly non-trivial infrared phenomena, one needs to calculate spectral functions non-perturbatively.Quite generally, however, this is very challenging for all equilibrium field theory methods in Euclidean space-time, as for example first-principles lattice QCD simulations, because it requires an analytic continuation of Euclidean correlations back to the real-time domain which is a numerically ill-posed inverse problem.While different reconstruction schemes exist for different purposes, ranging from Maximum Entropy methods [10][11][12], the Bachus-Gilbert method [13], Tikhonov regularization [14] or the Schlessinger-Point (or Resonances-via-Pad) method [15,16], all of these come with limited ranges of applicability (see, e.g.Ref. [17] for a recent comparison).To calculate spectral functions from non-perturbative functional methods, e.g. based on Dyson-Schwinger or Functional Renormalization Group equations [18][19][20][21][22][23][24][25][26][27][28] can provide a powerful alternative, especially when the functional equations are analytically continued or even formulated directly on the closed-time path.However, in one way or another, they all require prescriptions to obtain a closed set of equations form an originally infinite hierarchy.Controlling truncation errors then becomes important for the systematics, and any a priori knowledge of the structure of correlations in the theory is obviously beneficial as additional input or benchmark.
In this work, we calculate the spectral function of a single-component scalar field theory employing classical-statistical lattice simulations in real-time [29][30][31].We use the fluctuation-dissipation relation or Kubo-Martin-Schwinger (KMS, [32,33]) condition to obtain the spectral function from the statistical function which in the classical limit can be calculated from an elementary unequal-time correlation function of classical fields.While the method was successfully employed in a pre-cursor study [30], where the dynamic critical exponent of a scalar field theory in two spatial dimensions with conservative dynamics (Model C) was calculated, we extend this study in several regards and systematically compare the results in 2+1 and 3+1 dimensions.By extending the analysis of spectral functions to finite spatial momenta, we extract the universal scaling functions that govern the dynamics in the vicinity of the critical point.We also quantify the divergence of the characteristic time scale ξ t at the critical point and extract the critical exponent z for the different dynamic universality classes.
This paper is organized as follows: After giving a short overview of our model and simulation setup in Section 2, the classical-statistical approach to calculate spectral functions is explained in Section 3. We continue with a deeper look into the static critical behaviour of the system in Section 4, where we also extract the non-universal amplitudes of critical power laws, which allow us to present our subsequent results in terms of dimensionless scaling variables.Starting with Section 5, we present results for spectral functions at different points in the phase diagram, and analyze their behaviour away from criticality.We proceed with a detailed investigation of the dynamic critical behavior in Section 6, where we extract the dynamic scaling functions and scaling exponents.Our conclusions as well as an outlook for further applications are provided in Section 7. Some details on our analysis of the diverging time scales as well as the general structure of the non-critical spectral functions are given in the appendices.

Simulation setup
We simulate a relativistic Z 2 scalar field theory defined by the following lattice Hamiltonian: where φ x are single-component scalar field variables at lattice site x, and π x their respective conjugate momenta.By ř y"x we denote the sum over all lattice sites y adjacent to lattice site x, i.e. the nearest neighbours of x, while the sum ř x runs over all lattice sites x in the volume V .The spatial lattice spacing is denoted by a, which we set to a " 1 and drop from equations and units from now on.The parameter d P t2, 3u counts the spatial dimensions.In the absence of an explicit symmetry breaking (J " 0), the Hamiltonian in Eq. ( 1) is invariant under the Z 2 transformation (φ Ñ ´φ).However, for negative values of m 2 ă 0, the symmetry is spontaneously broken for temperatures T ă T c , and restored above T c via a second order phase-transition in the Z 2 (Ising) universality class.We study the critical behaviour of this Ising-like model in both 2+1 and 3+1 dimensions, and set our model parameters to m 2 " ´1, λ " 1, as well as J " 0 if not stated otherwise.
In contrast to quantum field theory, one can compute time-dependent observables in classical field theory directly.In particular, there is no fluctuating complex action inducing a sign problem, so one can obtain the real-time evolution by integrating the classical equations of motion of the system and averaging an ensemble over initial conditions.Clearly, there are many ways to obtain a thermal set of initial conditions; in practice, we use a Hybrid Monte-Carlo method [34].
We will consider both isolated systems and systems coupled to a heat bath.Generally, to compute a generic observable Opφptq, πptqq of the phase space variables φ and π, we solve Langevin-type equations where and η is a Gaussian white noise with zero mean and xη x ptqη y pt 1 qy " δ xy δpt ´t1 q.We subsequently evaluate the observable Opφptq, πptqq on the classical trajectory, and finally perform an averaging over the classical-statistical ensemble.We note that in Eq. ( 2) we have introduced an additional parameter γ which characterizes the strength of the coupling to the heat bath.Due to the fluctuation dissipation relation, the coefficient of the noise term is then fixed by γT and ensures a thermal Gaussian distribution of the conjugate momentum field variables π x at temperature T .Evidently, for γ Ñ 0 one recovers the Hamiltonian equations of motion, such that the total energy of the system H in Eq. ( 1) is explicitly conserved.In this case, individual realizations of the system evolve micro-canonically and the temperature of the system only enters in the characterization of the initial ensemble.
Based on the structure of the equations of motion in Eq. ( 2) and Eq. ( 3), one concludes that our system belongs to the class of relaxational models (Models A-D) in the classification scheme of Hohenberg and Halperin [1].Since our equations of motion do not conserve the order parameter, this is further limited to Model A if there are no conservation laws (γ ą 0) or Model C with conserved energy (γ " 0).

Spectral function
Spectral functions of bosonic operators Ô pt, xq in a quantum field theory are defined by a decomposition of the two-point function (with translational invariance in space and time): where T denotes the time-ordering operator, and we have introduced the spectral (ρ) and statistical (F ) correlation functions with r , s `representing the anti-commutator.While the spectral function ρ is relevant to describe the (linear) response of the system to an external perturbation, and thus provides the spectrum of possible excitations, the statistical correlation function F in general describes the quantum and thermal statistical fluctuations present in the system.Based on the Kubo-Martin-Schwinger (KMS) [32,33,35] condition in thermal equilibrium, F and ρ are related via the fluctuation-dissipation relation.
F pω, p, T q " ˆnT pωq `1 2 ˙ρpω, p, T q (7) Note that, for practical reasons, we defined the Fourier transformations of F and ρ as ρpω, p, T q " ´i ż dtd d x e ipωt´pxq ρpt, x, T q, where the additional factor of ´i in (9) ensures that the spectral function ρptq, which is an odd function under space-time reflections ρp´t, ´xq " ´ρpt, xq, is real in both the time (ρptq) and frequency domain (ρpωq).In the classical limit, the thermal distribution becomes and the spectral can be computed directly as Close to the critical point, the dynamics of the system is dominated by its slow infrared modes with frequencies ω !T very small compared to the temperature scale where quantum effects are relevant.Therefore, the classical description suffices to fully capture the dynamic critical behaviour of the model.
We are interested in the spectral function of the order parameter field φ, which in the classicalstatistical field theory is formally defined as ρpt ´t0 , x ´x0 q " ´xtφpt, xq, φpt 0 , x 0 quy cl . ( Table 1: Ising critical exponents.In 2D, they are analytically known from Onsager's solution [38].The high-precision calculations for the 3D Ising exponents were obtained from the conformal bootstrap approach [39,40].[36,37]), we follow previous works and instead exploit the classical KMS condition to calculate the spectral function in thermal equilibrium [29][30][31].By virtue of Eq. ( 11), we can compute this spectral function ρpt, x, T q via the statistical two-point function F pt, x, T q directly from the classical lattice fields, ρpt, x, T q " ´1 2T xπpt, xqφp0, 0q ´φpt, xqπp0, 0qy .

Static universality and scale setting
Before we address the real-time dynamics of the system, we briefly investigate the static critical behaviour of our model in 2 and 3 spatial dimensions. 1 Measuring the non-universal amplitudes will enable us to express our later results in terms of dimensionless quantities, for a more direct comparison as well as a plausibility check for our numerics setup.
The basic observables we want to look at are the order parameter in form of the magnetization M and the corresponding susceptibility χ, which are defined as χ " B xM y BJ " V T ´@M 2 D ´xM y Since in the absence of explicit symmetry breaking (J " 0) the magnetization defined in Eq. ( 14) vanishes identically xM y " 0 in any finite volume, we instead consider When extracting the spatial correlation length ξ, we consider the plane-correlation function Ḡpnq between the field-average Spnq over lines for d " 2 (respectively planes for d " 3): Ḡpnq which for sufficiently large separations n is expected to follow and exponential behaviour of the form Ḡpnq " A pexpp´n{ξq `exp ppn ´Lq{ξqq .
Based on this behaviour, we deduce an effective correlation length ξ eff.pnq from the data by considering the logarithmic derivative 1.5 (6) and subsequently look for a plateau in a range of separations n, which is then used to determine the spatial correlation length ξ. 2We provide a compact summary of our static results in Figs. 1 to 3, where we show the behaviour of the order parameter M and susceptibilities χ (Fig. 1), visualize the extraction of the spatial correlation length (Fig. 2), and show the behaviour of the correlation length ξ as a function of temperature T and external field J (Fig. 3).For all of the aforementioned observables, we approach the thermodynamic limit: By comparing data points of at least two different lattice volumes for any given point in the phase diagram, we can select a range of data where finite volume effects are negligible.We follow common procedure and express the temperature dependence of our results in terms of the reduced temperature τ τ " where T c is the critical temperature given in Table 2.
Based on our results in Figs. 1 and 3, we extract the non-universal amplitudes of the critical power laws either as functions of τ at vanishing external field, J " 0, or as functions of a dimensionless J " J{J 0 ě 0 at τ " 0. We follow the notation of [41] and distinguish the various non-universal amplitudes by superscripts, marking the region of the respective phase diagram that they describe.Specifically, X ˘indicates the sign of τ and marks the amplitude on the J " 0 axis, whereas X c denotes the corresponding amplitude on the τ " 0 axis.Motivated by [41] and [42], we extract the critical temperature and non-universal amplitudes, including the leading scaling corrections, by fitting the following ansätze to correspondingly selected points.
xM pτ qy " Bp´τ q β p1 `B1 p´τ q ων ´B2 τ q , τ ă 0 , φppq "    where γ c " γ{βδ " 1 ´1{δ and ν c " ν{βδ " p1 `1{δq{d from scaling and hyperscaling relations.We note that J 0 is chosen such that B c J 1{δ 0 " B, and C c J ´γc 0 " B{pδJ 0 q (as well as J ´νc 0 " pB c {Bq ν{β ); and in the absence of scaling corrections, the magnetic equation of state can be written in the Widom-Griffiths form, y " f pxq, with dimensionless magnetization Ď M " M {B, By virtue of the Z 2 symmetry of our model, we know that its critical behaviour is described within the Ising universality class.Our model and method are not optimized for computing static critical exponents; we therefore reduce the degrees of freedom in our calculations by taking the critical exponents from the literature.In d " 2 spatial dimensions, the critical exponents are exactly known from the analytic solution by Onsager [38].For the numeric values of the d " 3 exponents, we take those obtained by the conformal bootstrapping method [39,40], which are listed in Table 1 alongside their d " 2 counterparts.
The asymptotic scaling amplitudes along with the coefficients of the scaling corrections extracted from the fits to the data are listed in Table 2.The comparison between ansatz and data is shown in Figs. 1 and 3.For most quantities, we obtain a χ 2 {d.o.f.close to unity, indicating sufficient agreement between fit and data.For d " 2, our procedure for the correlation length at J " 0 does not work too well, especially for τ ă 0, leading to a relatively large resulting value of χ 2 .This is primarily due to heavily fluctuating data in the low-temperature phase, as well as probably an underestimation of the error of the correlation length ξ in our extraction.
It is a useful cross-check to validate the measured amplitudes by comparing the known universal amplitude ratios.We use a subset that relates all of the amplitudes we measure: The ratios we measured are shown alongside the amplitudes in Table 3.For most of the ratios, we find results in the right ballpark, with exception of U ξ and Q 2 in d " 2 dimensions: Those ratios are the ones including amplitudes of ξ.
Different definitions of ξ lead to slightly different U ξ , see e.g. the difference between second-moment correlation length and inverse mass gap in [41].Since our definition again differs from the aforementioned ones, we are bound to get (slightly) different amplitudes.Nevertheless, it is strongly related to the definition by the inverse mass gap, so we would expect U ξ « 2. Since we do not rely on exact estimates of the correlation length for the purpose of this study, we accept these residual ambiguities for now.
When studying dynamic critical behaviour, we will normalize our results using the non-universal amplitudes.In this process, we will also obtain a timescale by fitting the amplitude f t of the correlation time ξ t , which will be determined in a later section (see Eq. ( 69) in Sec.6.3).If not stated otherwise, we give all results in terms of dimensionless scaling variables, which are indicated by a bar and constructed as follows: p " pf `, ω " ωf t . (39)

Non-critical spectral functions
Based on our analysis of static critical phenomena, we will now investigate the behaviour of the spectral function at different points in the phase diagram.We proceed as outlined in Section 1 and calculate the classical-statistical spectral function of the order parameter field using Eq. ( 11) resp.Eq. ( 13).We prepare " 30 independent thermal configurations as initial conditions, which we then evolve for " 10 4 ´10 5 a, using and Euler-Maruyama scheme.Notably, we find that in the vicinity of the critical point, the time step ∆t in the integrator has to be chosen sufficiently small to avoid discretization errors.If not stated The spectral functions ρpω, pq for the smallest non-zero momentum mode are highlighted by a black solid line on the 3D surface.Spectral functions in the symmetric phase T " Tc are dominated by a quasi-particle peak with relativistic dispersion relation.Near the critical point T « Tc, a strongly infrared-enhanced mode emerges in addition to the quasiparticle peak that persists at higher momenta.Below the critical temperature T !Tc, a soft second mode at low frequencies emerges for finite momenta.otherwise, we employ ∆t " .00625a, and we have checked that discretization errors are negligible for the results presented here.By recording the evolution of the Fourier modes of the order parameter field, we subsequently compute the spectral function ρpt, pq in Eq. (11).Statistical errors are estimates form point-wise averages of ρpt, pq and respectively ρpω, pq over different configurations.If not stated otherwise, the spectral functions shown in this section are obtained for vanishing external field J " 0 on lattice volumes of 256 2 respectively 128 3 , where we have the smallest remaining finite volume effects.
In Fig. 4, we give an overview over the behaviour of spectral functions at different temperatures.At high temperatures (T " T c ), the spectral function is well approximated by a relativistic Breit-Wigner peak shape, with the peak position shifting with p according to the relativistic dispersion relation ω 2 ´p2 " m 2 pT q.In [29], the zero-momentum (p " 0) mode of the spectral function was investigated in a scalar theory without phase transition.The study was done at high temperatures, so that T " ω.It was found that both, masses and widths of the measured spectral function agree well with the analytical predictions from resummed two-loop perturbation theory.We find here a very similar Breit-Wigner shape, and since we also calculate the spectral function at finite spatial momenta, we can check for the correct dispersion relation of the quasi-particle mass.
When the system approaches criticality from above, the effective mass m 2 pT q decreases, and the quasi-particle peak becomes less and less pronounced.Close to the critical point, an infrared power law behaviour builds up at small frequencies and momenta which, as we will discuss shortly in Section 6, encodes the dynamic critical behaviour of the spectral function. 3With increasing spatial momentum p, the cut-off imposed by p suppresses the critical contribution.Specifically, at high momenta the spectral function retains its Breit-Wigner shape even in this near-critical regime, however with significantly smaller quasi-particle mass.By comparing the different columns of Fig. 4 one further notices that the window of reduced temperatures τ , where a critical enhancement of the spectral function can be observed, is much smaller in 3+1 than in 2+1 dimensions, as can be expected from the larger influence of infrared fluctuations in lower dimensions.
Below the critical temperature, in the ordered phase, the effective quasi-particle mass m 2 pT q increases again, and gradually approaches its mean-field value of ?´2m 2 in the limit T Ñ 0, where all thermal fluctuations are suppressed.However, in addition to the quasi-particle peak, a second low-frequency excitation arises for finite spatial momenta p, with a different spectral shape and a dispersion relation indicative of soft collective excitations such as thermally driven capillary waves, see Appendix B. While in 3+1 dimensions the contributions from this soft mode only carry a small fraction of the spectral weight, it is more pronounced in 2+1 dimensions, where by looking at the lower left panel of Fig. 4, one can easily discern the valley in the spectral function that separates it from the quasi-particle peak.
By comparing the spectral functions in the symmetry-broken phase at τ ă 0, and the vicinity of the critical point τ « 0 in Fig. 4, one is led to speculate that it is this second excitation at τ ă 0 that might turn into the critical IR divergence close to the critical point.Indeed, as we will demonstrate shortly, by tracing the maxima of both excitations, one observes an avoided-crossing behaviour near T c , before the low-frequency mode eventually disappears at τ ą 0. A similar avoided crossing was already observed in the classical-statistical lattice simulations of O(4)-model spectral functions in 3+1 dimesions in [31].Although this study was restricted to vanishing spatial momenta where in agreement with the results presented here, there is no comparable soft mode in the longitudinal σ-spectral function of the order parameter fluctuations, an avoided crossing does show up at zero momentum in the transverse π-spectral function of the O(4)-model, see Ref. [31] for details.

Dispersion relation of quasi-particle peak
Based on the results in Fig. 4, we will now investigate some properties of the spectral functions extracted from lattice data.We start by confirming the relativistic dispersion relation of the quasiparticle peak, study the temperature dependence of the peak parameters and analyse the effect of the Langevin damping γ on the spectral function.
In order to analyze the dispersion relation, we fit the peaks in the spectral function with a relativistic Breit-Wigner ansatz where the central frequency of the effective mass resonance is expected to follow the relativistic energymomentum relation Visualizations of the fits to the spectral function are given in the top panels of Fig. 5 for the symmetric phase (τ ą 0) and in Fig. 6 for the symmetry broken phase (τ ă 0), while the bottom panels of Figs. 5  and 6 show the extracted values of the peak-frequency ω c ppq and damping rates Γppq.Clearly, the fit to the Breit-Wigner function describes the spectral function very well at high temperatures, where the quasi-particle excitation is the only discernible structure in spectral functions in Fig. 5.While in the low temperature phase, shown in Fig. 6, the Breit-Wigner fit still describes the quasi-particle peak very well, one also clearly observes the additional low-energy excitation for finite spatial momenta (p ą 0).
So far we have focused on the general behaviour of the spectral function for Hamiltonian dynamics (Model C), which we will now compare to the spectral functions for Langevin dynamics (Model A).Since the key difference between Hamiltonian and Langevin dynamics lies in the introduction of an additional frequency-independent damping and noise coupling to the heat bath, one naturally expects the additional damping to contribute to the resonance decay width Γppq of the (non-critical) spectral functions.Explicit comparisons of the results for Hamiltonian (γ " 0) and Langevin (γ " 0.1) dynamics, shown in the bottom panels of Figs. 5 and 6, confirm this expectation, indicating further that differences in the numerically extracted damping rates Γ γ"0.1 ppq ´Γγ"0 ppq » γ are in fact close to the Langevin coupling γ.While away from criticality, the qualitative and quantitative features of the spectral functions are very similar between conservative and weakly dissipative systems, the change in conservation laws does of course affect the dynamic critical behaviour, as we will discuss in more detail in Section 6.
We conclude our discussion of the non-critical spectral functions by investigating the temperature dependence of the central frequencies and decay widths of the quasi-particle peak, which are summarized in Fig. 7.By defining m eff " a ω c ppq ´p2 , curves for different momenta p should coincide whenever the quasi-particle peaks satisfy the dispersion relation in Eq. (41).Starting from the low temperature phase, the resonances turn into sharp δ-peaks in the limit T Ñ 0, where m eff pT " 0q " ?´2m 2 assumes its mean-field value and the decay width ΓpT " 0q " 0 vanishes.
Even at the lowest temperatures measured, the effective masses are still rather large but then decrease as the system approaches the critical temperature, while the decay width ΓpT q increases simultaneously.Near the critical temperature, the behavior of the spectral function at low momentum can no longer be described by a simple quasi-particle structure, as already seen in Fig. 4.However, when increasing the temperature above the critical point, the quasi particle structure at low momentum is restored, with the mass increasing monotonously as a function of temperature for T ą T c .Conversely, for the modes with large spatial momenta, this process is continuous.While the effective mass reaches a finite minimum at τ " 0, the spectral function retains its Breit-Wigner form across the transition.
In summary we find that, in the non-critical regime, the measured spectral functions behave as expected.The relativistic dispersion relation is fulfilled, and we find the expected temperature dependence on both sides of the phase transition.Is is interesting to see how, for low spatial momenta p ! 1, the spectral shape changes at different temperatures and shows critical scaling, while at larger spatial momenta p ą 1, the spectral function retains its Breit-Wigner shape and moves across the transition in a completely continuous fashion.
In the non-critical regime, introducing a finite Langevin damping γ leads to broadening of resonances, as one expects from the structure of the equations of motion.Close to the critical point, we see the same effect of the Langevin damping γ on the remaining quasi-particle contribution to the spectral function, whereas the infrared-divergent critical part appears qualitatively unchanged.However, as we will show in detail in Section 6, the slope of the infrared divergent power law is modified slightly.Upper panels show the spectral functions in lattice units for vanishing Langevin coupling γ " 0 with an exemplary Breit-Wigner fit as solid line; lower panels show the resulting fit parameters in lattice units over the different spatial momenta, with a fit to a relativistic dispersion relation Eq. (41).Note that the width Γ is scaled by a factor of 5 to improve legibility.

Critical dynamics and scaling functions
In the following section we investigate the critical behaviour of spectral functions.We focus on the determination of a universal scaling function for the spectral function, as well as the extraction of the dynamic critical exponent z.
Dynamic critical phenomena have been studied since the late 60s, with the first numerical results for Ising models in the 90s [7,[43][44][45][46][47][48][49].These studies were mostly concerned with finding the dynamic critical exponent z of the Ising model with Glauber dynamics (2D Model A), where a multitude of results ranging from z " 1.7 to z " 2.7 have been published, however slowly converging towards z « 2.2.Nightingale and Blöte [7] were the first to calculate z in 2D Model A in an Ising model with high precision using a variance-reducing Monte Carlo algorithm, albeit on rather small lattices up to 15 ˆ15.They found z " 2.1665p12q, quoting a two-sigma error, in accordance with the former trend.
In 2004, Dunlavy and Venus measured critical slowing down in ferromagnetic ultra-thin films [2], governed by 2D Model A dynamics as well.The resulting critical exponent of 2.09p6q, giving the 95% confidence interval, is close to the Monte Carlo result, but the difference of nearly 2.5σ is somewhat large.
The most recent result by Zhong et al. [50] for a two-dimensional scalar φ 4 model with local Metropolis updates seems to confirm the results for 2D Model A by Nightingale and Blöte from [7], as they find z A " 2.17p3q and z A " 2.19p3q for two different values of the coupling constant.Additionally, they demonstrated that quantities derived from the two-point function follow a scaling behaviour.
In 2010, a precursor study [30] of the present one first tried to confirm the dynamic critical exponent in 2D Model C, z " 2 `α{ν " 2 for an Ising-like scalar field theory with conservative dynamics, and found z " 2.0p1q.
Naturally, numerical simulations in d " 3 are more expensive than in d " 2, and thus there are fewer and less precise results for the dynamic critical exponents.The study by Matz et al. [45]     Due to availability of data, we use τ " 0.0009p2q (d " 2) resp.τ " ´0.00008p5q (d=3) as a proxy for the critical temperature.Note that, despite τ « 0, finite size effects are not relevant here due to finite spatial momentum.However, to achieve the comparatively very low p in d " 3, a single set of data at N " 512 very close to Tc was generated.
We are not aware of any previous Monte-Carlo studies on the dynamic critical exponent for Model C in d " 3; however, based on the classification scheme by Halperin and Hohenberg the dynamic critical exponent there is known by virtue of the same scaling relation (z " 2 `α ν ) as in d " 2, which amounts to z « 2.17 for a model in the 3D Ising universality class.

Dynamic scaling functions
Since the spectral function is derived directly from the two-point correlation function, one expects the critical behavior to be governed by the following scaling form [30] ρ pω, p, τ q " s 2´η ρ ´sz ω, sp, s in the limit ω, p, τ Ñ 0, and we omit any residual dependencies on the finite volume.If not stated otherwise, we will restrict ourselves to positive frequencies pω ą 0q to compactify notation, noting that the behaviour for negative frequencies pω ă 0q is trivially obtained from the symmetry of the spectral function ρp´ω, p, τ q " ´ρpω, p, τ q.The scaling law in Eq. ( 42) allows us to define three alternative scaling functions f ω , f p and f τ according to ρ pω, p, τ q " p´p2´ηq f p ´ω{p z , τ {p 1{ν ¯, where γ is the static susceptibility exponent with γ " νp2 ´ηq from static scaling relations.Indeed, the scaling behaviour in Eqs. ( 43) to ( 45) is clearly visible in our classical-statistical simulations, as can  of data, we use τ " 0.0009p2q (d " 2) resp.τ " ´0.00008p5q (d=3) as a proxy for the critical temperature.Solid lines represent the fit to the ansatz in Eq. (57), we fit fω and ap as free parameters, and the exponent of the power law for large xp is given by p2 ´ηq{z comb. .Note that finite size effects are not relevant here due to finite spatial momentum.To achieve the comparatively very low p in d " 3, a single set of data at N " 512 very close to Tc was generated.
be seen from Figs. 8 to 13, where we present results for the scaled spectral functions in the vicinity of the critical point.We note that in order to perform the axis re-scaling in Figs. 8 to 13, one also needs the value of the dynamic critical exponent z, and if not stated otherwise, we employ the values in last row of Table 6, labeled combined, with exception of d " 3 Model C, where we used the analytic value of z « 2.17.
Since the scaling functions in Eqs. ( 43) to ( 45) are all derived from the scaling behavior of the spectral function, the scaling functions f ω ,f p and f τ are not independent.Denoting the natural arguments of the respective scaling functions as one finds the following relations between the scaling functions where the super-script in f τ is used to distinguish sgnpτ q " sgnpy p q " ˘1 above/below T c .Similarly, the Fourier transform of Eq. ( 42) yields in the time domain ρ pt, p, τ q " s p2´η´zq ρ ´s´z t, sp, s 1{ν τ ¯, from which one can derive the time-domain scaling functions analogous to Eqs. ( 43) to (45).
Based on the results of [30], we expect the scaling function f ω px ω , y ω q to be regular in the origin, with a constant value f ω " f ω p0, 0q, such that the critical zero-momentum spectral function obeys the infrared power law ρpω, 0, 0q " This expectation is verified in Figs. 8 and 9, where we present simulation results for f ω px ω , 0q " ω p2´ηq{z ρpω, p, 0q as a function of x ω at τ 9 y ω " 0 (Fig. 8) and respectively f ω p0, y ω q " ω p2´ηq{z ρpω, 0, τ q as a function of y ω at p 9 x ω " 0 (Fig. 9).Generally, one observes an excellent scaling collapse of the data for different momenta(reduced temperatures), with deviations only at very small values of x ω (y ω ) where the residual contributions from the quasi-particle peak at finite momentum(reduced temperature) effectively acts as a cut-off.Specifically, for small arguments x ω or y ω the two data sets in Figs. 8 and 9 appear to converge towards the same value f ω indicated by a solid black line, demonstrating that f ω is indeed finite and that the limits p Ñ 0 and τ Ñ 0 commute at small but finite frequency ω ‰ 0. Due to the above relations between the scaling functions f ω , f p and f τ , we can further exploit the regularity of f ω near the origin px ω , y ω q " p0, 0q to deduce the behavior of the scaling functions f p and f τ for large values of the arguments x p " 1{x ω and x τ " 1{|y ω | νz as p2´η ρpω, p, 0q " This fixes the high-frequency (ω " 1{ξ t pτ, pq) behavior of the critical scaling function f p px p , 0q at very low momenta, and that of the zero-momentum scaling function f τ px τ , 0q very close to criticality at the same time.The characteristic time ξ t pτ, pq for this asymtpotic behavior, from the definition of our scaling variables, behaves as ξ t p0, pq " 1{p z and ξ t pτ, 0q " 1{|τ | νz at small p and τ , respectively.Conversely, the low-frequency pω !1{ξ t pτ, pqq behavior of the spectral function is determined by the behavior of f ω px ω , 0q and f ω p0, y ω q at asymptotically large values of x ω and y ω , and can be determined from the behavior of f p px p , 0q and f τ px τ , 0q at small values of the arguments x p " 1{x ω and x τ " 1{|y ω | νz .Numerically, we find from Figs. 10 and 11 that the leading behavior at small x p and x τ is well described by which corresponds to Specifically, for the critical spectral function (τ " 0) at finite spatial momentum, Eq. ( 54) implies an infrared behavior, valid asymptotically for ω !pz , ρpω, p, 0q " a p p´p2´ηq ω{p z `¨¨¨.
Comparing Eq. ( 55) with Eq. ( 51), we see that the limits p Ñ 0 and ω Ñ 0 do not commute, as the critical spectral function is non-analytic in the origin.In particular, we have lim ωÑ0 ρpω, 0, 0q " 8, while lim pÑ0 ρp0, p, 0q " 0. Physically, this has the intuitive interpretation that any finite momentum p introduces an effective IR cutoff for the correlation length ξ " 1{p, which in turn is associated with a finite correlation time ξ t p0, pq " ξ z " 1{p z , that defines the characteristic frequency ω " pz where the power law changes from (51) to (55).
Similarly, for the zero-momentum spectral function off criticality we conclude that the infrared behavior for ω !|τ | νz is modified as so that for all τ " 0 we also have lim ωÑ0 ρpω, 0, τ q " 0, with the characteristic frequency where the power law changes from ( 51) to (56) given by ω " |τ | νz .Once again, this simply reflects the finiteness of the correlation time at non-vanishing τ where we now have ξ t pτ, 0q " ξ z " 1{|τ | νz .So far we have analyzed the limiting behavior of the scaling functions for very large and very small arguments.Now, in order to interpolate between these two limits, we observe that the inverse of the scaling functions f p px, 0q and f τ px, 0q is globally very well described by a sum of the reciprocal power laws ( 52) and ( 53) for large and small x.We therefore use the following parametrizations of the scaling functions to fit our data, Our results are compactly summarized in Figures 10 and 11, where we show these parametrizations for the scaling functions f p pω{p z , 0q and f τ pω{|τ | νz , 0q fitted to our data for pp2´ηq ρpω, p, 0q and |τ | γ ρpω, 0, τ q respectively.The regions where the points overlap are used to determine the universal scaling functions.
The numerical results for the parameters a p , a τ , f ω obtained from the fits are listed in Table 4. Corresponding plots for the scaling laws in the time domain are shown in Figs. 12 and 13 where we plot the Fourier transformed scaling functions fp px, 0q and f τ px, 0q, from the analogous definitions ρpt, p, τ q " p´p2´η´zq fp pp z t, τ {p 1{ν q , (59) By looking at Figures 10 and 11, one finds that scaled data sets for the spectral functions at different momenta/reduced temperatures overlap with each other to rather good accuracy, and that, apart from f τ px τ , 0q in d " 2 spatial dimensions, the ansätze in Eqs. ( 57) and ( 58) match the overlapping data points exceptionally well.Moreover, the characteristic frequencies ωc " pz and ωc " τ νz mentioned above can be read off directly from the coinciding maxima in the respective scaling variables x p and x τ .The above fits ( 57) and ( 58) in turn imply that the frequency scaling function f ω px ω , y ω q, either at criticality (y ω " τ {ω 1{νz " 0) or at zero momentum (x ω " pz {ω " 0), is equally well described by with a τ for sgnpy ω q " ˘1 above and below T c , which provides a complete description of the scaling function f ω px ω , y ω q along the two coordinate axes.By comparing our results for the scaling functions in Figures 10 and 11, we find that in general, the scaling regions are larger in d " 2 than in d " 3, and similarly larger for Model A than for Model C. Specifically, for d " 3, it seems that in some cases the standard volumes (N " 128) are still not large enough to show critical effects of reasonable strength.Hence, in order to achieve some sufficient overlap for extracting the universal scaling function at finite spatial momentum p, we have generated a single additional data set at T « T c with N " 512 in d " 3, which allows us to investigate very small spatial momenta, with p ă 0.05 in lattice units.By inspecting e.g. the upper right panel in Fig. 10, it is then clear that for very small p, the overlap region does indeed extend to the right of the maximum, where the slope in the logarithmic plot is determined by the critical power law.
One important advantage of working with finite spatial momenta p ą 0 is that finite volume effects are essentially irrelevant, as the relevant infrared cut-off is set by the momentum p rather than the system size.Conversely, at p " 0, finite volume effects inevitably appear close to criticality, i.e. for very small τ , as can be seen e.g. in Figs.11 and 13.
By closer inspection of the results in Figures 10 and 11, e.g. when comparing models A (upper row) to models C (lower row), one further notices that the spectral functions ρpω, p, τ q start to deviate from the scaling function when reaching the remnants of the quasi-particle peak, which appear as additional "shoulders" at the high-frequency end.Clearly, this effect is more pronounced in Model C, where highfrequency fluctuations do not receive the additional damping due to the coupling to the heat bath.By looking at the results for τ ă 0 in Fig. 11 we also note that the universal scaling of the spectral functions, when approaching criticality from the ordered phase, does not appear to emerge from the remnants of the quasi-particle peak, even in close vicinity of the critical point.Instead, as previously alluded to in the context of the discussion of Fig. 4, it rather seems that the universal critical behavior of the spectral function emerges from the soft collective low-frequency excitation in the ordered phase.

Extraction of the dynamic critical exponent
Having established the dynamic scaling behavior of the spectral function, we now turn to the extraction of the dynamic critical exponent z.Generally speaking, there are multiple possibilities to extract the dynamic critical exponent z from the data, and we will explore three different methods in the following, which we will refer to as the scaling method, the critical IR power law method, and the divergence of the correlation time ξ t method.
Evidently, to obtain the plots of the scaled spectral function in Figs. 8 to 13, one needs the correct value of z that maximizes the overlap of the data points.Vice versa, we can exploit this scaling behavior to extract the dynamical critical exponent z by minimizing the deviations from perfect scaling, which we quantify in terms of the L 2 -norms of pairwise distances of rescaled functions over some frequency interval rω l , ω h s.Specifically, for the p-rescaled critical spectral functions as depicted in Figure 10, based on Eq. ( 44), this amounts to minimizing the quantity where ∆ρ denotes the statistical error of the measured spectral functions, which is used to weight the deviations.Similarly, an analogous functional is minimized to optimize the scaling of the τ -rescaled spectral functions at vanishing momentum in Figure 11.
Even though this procedure is in principle very robust, we were not able to completely eliminate the dependence on the upper limit ω h of the frequency interval.Since this dependence is particularly strong for the τ -rescaled spectral functions, we have disregarded them in the final estimate.By the principle of least sensitivity, i.e. by looking for a plateau in the results for different ω h , we can then estimate reasonable values for ω h in case of the p-rescaled spectral functions, which yield a set of plausible values for the dynamic critical exponents that are shown in Table 6 in the row labeled scaling.However, since the systematic uncertainties associated with this procedure are still somewhat uncontrolled, we have also explored two alternative methods to calculate the dynamic critical exponent z.
Our second method of extracting the dynamic critical exponent z from our data exploits the large x p and x τ behavior of the scaling functions f p px p , 0q and f τ px τ , 0q.Based on the scaling form of the spectral function in Eq. ( 43), we extract the infrared power law (51) from either the critical spectral function at τ " 0 and some sufficiently small momentum p, or that at zero momentum sufficiently close to criticality.For example, Eq. (57) entails that the frequency dependence of the critical spectral function (τ " 0), at a given fixed value of p, is of the form ρpωq " 1 Likewise, cf.Eq. ( 58), the same form should also describe the frequency dependence of the zero-momentum spectral function at fixed small τ " 0, with σ " p2 ´ηq{z in each case.We note that the power-law amplitude b " f {f ω is fixed by the constant f ω governing the regular behavior of f ω px ω , y ω q near the origin at x ω " y ω " 0, whereas the parameter a is related to the amplitudes governing the leading small x p resp.x τ behaviour of f p px p , 0q resp.f τ px τ , 0q via the relations where pfit , τ fit designate the spatial momentum resp.reduced temperature where the fit was performed.By fitting this ansatz with the amplitudes a, b and the exponent σ as free parameters to the critical spectral functions (τ " 0) at fixed momentum we obtain reasonably stable results for the dynamic critical exponent z, which are listed in the second row of Table 6 with the label IR power law.In practice, we simultaneously fit the spectral functions for the smallest two (d " 3) or three (d " 2) spatial momentum indices to improve statistics.Unfortunately, in d " 3 even for our largest possible lattices we can not reach small enough spatial momenta p to obtain a reasonably well constrained signal for z, especially in Model C. .Different symbols correspond to extractions of the correlation time using fits (filled symbols) and integration (open symbols) of the spectral function; upper curves in each panel correspond to ξtpτ ą 0q, while lower curves correspond to ξtpτ ă 0q.Solid lines in each panel correspond to a power law fit according to Eq. ( 71), from which we extract the dynamic critical exponent z along with the non-universal amplitudes f t .
In order to estimate the uncertainty in z, we vary the upper limit of the frequency interval where we fit Eq. (64) to the data in a sensible range.We eliminate one third of the results with the largest χ 2 {d.o.f.. Of the remaining results, we take the highest and lowest values for z as bounds on the confidence interval, and calculate the average of z weighted with the statistical uncertainty.If the weighted average of z is not near the center of the confidence interval, we separately note the uncertainties in both directions.
If one naively repeats this process for the zero-momentum spectral functions at small τ , one arrives at implausible values of z, which drift towards the result at finite spatial momentum upon approaching the critical point at τ " 0. This is due to the fact that, as can be seen in Fig. 11, the universal scaling function f τ deviates from the asymptotic power law for intermediate x, but converges to it for large x.We take the convergence of f τ for large x to the same power law that describes f p px, 0q for x Á 1 as an indicator that this power law describes the true asymptotic behaviour of the universal scaling functions.Nevertheless, since f p converges much earlier, we believe that the results for z from fits to the spectral function at τ " 0 and small momentum p are much more reliable.

Critical behavior of the auto-correlation time ξ t
We can also extract the dynamic critical exponent z from the divergence of the auto-correlation time ξ t pτ q in the vicinity of the critical point.By use of Eq. (50) and setting s " t1{z , the spectral functions in the time domain satisfy the following scaling relation where the presence of the second scaling variable y " τ t 1{νz explicitly shows that there is a temperaturedependent characteristic time scale ξ t " |τ | ´νz associated with the auto-correlation time.For the sake of simplicity, we focus on the case of vanishing spatial momentum p " 0, where the spectral function near criticality at times larger than this characteristic correlation time pt{ξ t pτ q ą 1q, is well described by a product of a power law and an exponential, as can be seen in Fig. 13.Near the critical point, the correlation time ξ t pτ q diverges as the system approaches criticality, and the spectral function then assumes its characteristic power-law form.Nevertheless, away from criticality the correlations in the system decay exponentially on large timescales t Á ξ t , and the functional form in Eq. ( 68) can be used to define the auto-correlation time ξ t pτ q from fits to numerical data, and infer its non-universal amplitudes f t and their universal ratio U ξ,t as We further use the results for the auto-correlation time ξ t pτ q in the high-temperature phase pτ ą 0q to define t " t{f t as our dimensionless time variable when presenting spectral functions in the time domain, such that for τ ą 0 the scaling variable τ νz t employed in Fig. 13 corresponds directly to the ratio t{ξ t .Besides extracting the auto-correlation time ξ t pτ q from a fit to the ansatz in Eq. (68), it is also possible to extract ξ t pτ q by integrating the spectral functions in the time domain, as explained in detail in Appendix A. Before we present our numerical results, some further remarks are in order.Since we focus on the behavior of modes with vanishing spatial momentum p " 0, their auto-correlation time ξ t pτ q diverges at the critical point as ξ t pτ q " 1{|τ | νz an infinite system.However, in our simulations the divergence of the auto-correlation time is limited by the finite system size L, which in the immediate vicinity of the critical point τ « 0 becomes the relevant infrared cut-off.Based on the dynamical finite-size scaling hypothesis, one expects that in this regime the auto-correlation time behaves as where g ξ pxq is the finite-size scaling function of the auto-correlation time analogous to the ratio R " ξ{L used for static finite-size scaling.While for all finite values of the finite-size scaling variable x " τ L 1{ν the divergence of the auto-correlation time is effectively regulated by the finite volume, in order to recover the infinite-volume scaling in Eq. ( 69), one needs asymptotically large values of ˘x, where this finite-size scaling function satisfies We present our results for the auto-correlation time in Figs. 14 and 15, where we study the dependence of ξ t pτ, LqL ´z as a function of the finite-size scaling variable τ L 1{ν in Fig. 14 and subsequently estimate the magnitude of singular and regular contributions to the correlation length ξ t pτ q in Fig. 14.Generally, the results of the two different extraction methods (fit and integration) agree very well with each other, although for τ ă 0 the integrated ξ t are generally somewhat closer to the power law, and the slope of this power law fit produces slightly smaller results for z in Model A. Strikingly, one also observes from Fig. 14 that the data exhibits a clear finite-size scaling across different lattice sizes, which we can exploit to extract the dynamic critical exponent as explained in the following.
In order to obtain the dynamic critical exponent z from the correlation times, we first apply finite-size scaling with a plausible estimate for z, to find a region where the data for different lattice sizes shows sufficient overlap.Based on the results depicted in Fig. 14 it becomes obvious that this hardly works at τ ă 0, but gives a clear power law at large values of τ L 1{ν , for τ ą 0. We then fit the power law in Eq. ( 71) for τ ą 0 to the un-scaled data in the selected region, to get both the amplitude f t and the exponent νz.Subsequently, we estimate the amplitude ratio U ξ,t as far as possible by fitting a power law with the exponent obtained earlier to a few data points with τ ă 0. Errors are obtained in a similar way as for the power law fit of the IR divergence of the spectral function.By varying the temperature interval where we fit the power law to the correlations times ξ t pτ q and keeping two-thirds of the results with the lowest χ 2 {d.o.f. as well as eliminating outliers with χ 2 ą 2χ 2 min , we compute the averages weighted by statistical uncertainties and estimate the confidence interval by taking the highest and lowest values for ξ t , f t and U ξt .We remark that these parameters are strongly correlated, so a large uncertainty in the dynamic critical exponent z leads to large uncertainties in both f t and U ξt .
The results of this procedure for the non-universal amplitude f t and the ratio U ξ,t are given in Table 5.Those for the dynamic critical exponent z are shown in the row denoted by "ξ t power law" of Table 6.Especially the two large amplitude ratios U ξ,t in d " 2 seem quite remarkable when compared to the analytically known amplitude ratio of the spatial correlation length U ξ " 2 [41].Although the d " 2 data for ξ t below the critical temperature pτ ă 0q does not necessarily justify a power law fit all that well, by looking at Fig. 14 one is led to conclude that we might rather underestimate this ratio.
While in Fig. 14 the data at least in 2+1D perfectly fits a power law above the critical temperature τ ą 0, we find that a precise extraction of z remains difficult with the available data.In order to improve the accuracy, one could generate data closer to τ Á 0 in large volumes to minimize finite-size effects.Below the transition temperature in both 2+1D and 3+1D, we find that the data for the correlation time ξ t deviates strongly from the expected power law behaviour.One reason for this is the (much) smaller value of the non-universal amplitudes f t , which in combination with relatively large regular contributions leads to a suppression of the critical signal.We try to capture the regular contributions by fitting a regular function up to linear order in addition to the power law The comparison between the resulting fit and the data is shown in Fig. 15.The fit now also describes the data away from τ " 0 much better, which is dominated by the regular part, shown in Fig. 15 as a dashed line, especially on the low temperature side.However, by introducing these additional degrees of freedom in the fit, we lose precision in the estimate of the dynamic critical exponent z, both in terms of statistical and systematic uncertainties.
Besides providing an alternative means to extract z, one additional advantage of the auto-correlation time method is that it allows for a direct comparison of the critical dynamics of different models.In particular, to estimate the difference between the dynamic critical exponents z of Models A and C, one can look at the ratio of the correlation times at the same (reduced) temperature, which satisfy ξ t,A pτ q ξ t,C pτ q τ ą0 " f t,A f t,C ¨τ ´νpzA´zC q (73) Solid lines in the left and right panels show a power law fit of the large τ L 1{ν behavior, from which we extract the difference νpz C ´zA q of the dynamic critical exponents (cf.Eq. ( 73)).Dashed lines indicate the confidence interval of the extraction of νpz C ´zA q, which is also presented in the figure.
in the infinite volume limit.Such a direct comparison between Models A and C is presented in Figure 16, where we show the ratio of the correlation lengths in Eq. ( 73) for the symmetric phase (τ ą 0q as a function of the finite-size scaling variable τ L 1{ν .Even though this ratio reveals some tension between the two different extraction methods (exponential fit and integration), the general trends are clearly visible, where in d " 2 dimensions, z C ą z A , and the power law at large τ L 1{ν slopes downwards; while in d " 3 dimensions, the difference changes sign z C ă z A and the slope of the power law is positive.By performing a power law fit to the ratio, we can obtain a direct estimate of z C ´zA , which is also indicated in Figure 16, with the quoted errors obtained in the same way as for the power law fit of the correlation times.So far we have only considered the critical behaviour of the spectral function in absence of explicit symmetry breaking (J " 0).When introducing a non-zero explicit symmetry breaking (J ‰ 0), the magnetic scaling hypothesis states that the singular part of the free energy density can be written as f sing pτ, Jq " s ´df s ps 1{ν τ, s βδ{ν Jq (74) allowing one to express the singular behavior of the free energy in terms of scaling functions f˘, f J (note that νd " 2 ´α) Overlapping data points correspond to the universal magnetic scaling function of the auto-correlation times ξt,J pτ {| J| 1{βδ q (cf.Eq. ( 77)).Scaling breaks down when J becomes too large, as indicated by the blue data points.
which by analogy allows to define the magnetic scaling function We investigate this behavior in Fig. 17, where the auto-correlation time ξ t pτ, Jq obtained by integration for different values of J is rescaled to recover the underlying magnetic scaling function ξt,J as a function of τ {| J| 1{βδ .Data points for the two smaller values of J largely overlap confirming the magnetic scaling of the auto-correlation time.Even though the magnetic scaling starts to break down for the larger J, as one departs from the critical region, it is clear that in a certain (model dependent) range of J and τ , one can extrapolate ξ t pJ, τ q from the overlap region in Fig. 17.

Discussion of results for the dynamic critical exponent z
We provide a summary of our results for the dynamic critical exponent z in Table 6, where we give the results from the three different extractions, along with an error-weighted combined average.
For Model A in 2+1D, the overlap method and the power law fit to the correlation times ξ t give a surprisingly small result for z with a relatively large error.The power law fit to the infrared divergence of the spectral function yields a result closer to the one of [7], albeit even a bit larger.Combining the results leads to a z that is closer to, but still smaller than the result of [7].Incidently, our combined result for Model A in 2+1D is fully consistent with the experimentally measured value of z " 2.09p6q from Ref. [2].For Model C in 2+1D, we find a value that closely matches the analytic result z " 2, with all methods agreeing within their respective statistical errors.
In 3+1D, the non-critical effects on the spectral function are strong, leading to very large uncertainties in the power law fits.Since the critical window is smaller than in 2+1D, the uncertainty of overlap method increases as well.Nevertheless, the combined results in 3+1D for both Model A and Model C are compatible with earlier studies, and with the value pedicted from the scaling relation z " 2 `α{ν for Model C.
We note that, in order to obtain more precise results, it would be highly beneficial to consider an improved action and/or larger volumes at temperatures closer to the critical point.In 2+1D with larger volumes, one could probe smaller spatial momenta at τ " 0 and increase the precision of the IR power law method.When considering 3+1D, our results in Fig. 11 compared to those in Fig. 10 clearly indicate that, with our present setup, we were not able to probe low-enough momentum regimes away from the critical temperature T « T c .In both cases we are limited by small critical amplitudes of the correlation time ξ t , especially below the critical temperature pτ ă 0q.However, this study was not intended to achieve high-precision results for z, but to test the framework and lay the ground for upcoming non-equilibrium studies.

Conclusion and outlook
By using the classical field approximation close to a second order phase transition, we performed a first principles calculation of the spectral functions of the relativistic Z 2 model.By including finite spatial momenta in our analysis, we could provide a comprehensive overview of the behavior of the spectral function in all distinct parts of the phase diagram at finite temperature.We found that, in the symmetric phase, the spectral function is sufficiently well described by a Breit-Wigner quasi-particle shape, with a relativistic dispersion relation and a weakly momentum-dependent decay rate.When the Z 2 symmetry is spontaneously broken below T c , we find an additional excitation at low frequencies, with a different spectral shape and dispersion indicative of a soft collective mode, cf.Appendix B. While for sufficiently large spatial momenta, the quasi particle peak behaves continuously across the transition, we found clear indications that it is the second low-frequency mode that transforms into the dominant IR divergence at the critical point.Since our results have been obtained from first-principles numerical studies in the classical-statistical limit, it may be insightful to compare these results quantitatively to others obtained e.g. by the use of functional methods in the future.
By examining the static critical behaviour of our model, we have set the scales for the main focus of this study, which has been the quantification of dynamic critical effects.By analyzing the behavior of the spectral function of the order parameter in the vicinity of the phase transition, we explicitly verified the dynamic scaling hypothesis, and performed a detailed analysis of the scaling properties as a function of frequency, momentum and reduced temperature.We successfully extracted the corresponding universal scaling functions, which describe the infrared properties of the spectral function of the order parameter, and carefully assessed the implications for the behavior of the spectral function at small frequencies and momenta.We further developed a complete parametrization of the scaling functions in the special cases where either the spatial momentum or the reduced temperature vanishes.
We also analyzed the divergence of the correlation time ξ t in the vicinity of the critical point, by performing two different extractions of this quantity, and demonstrated its finite-size and magnetic scaling properties in the vicinity of the phase transition.By modifying the classical equations of motion to introduce a coupling to a heat bath, we were able to simultaneously study the behavior in the dynamic universality classes of Models A and C. While away from criticality, the additional coupling to the heat bath in Model A only gives rise to minor changes of the spectral function, we found that it leads to a change in the dynamic critical exponent z, as is expected due to the change in conservation laws.
We extracted the dynamic critical exponent z for a total of four different scenarios, corresponding to Hamiltonian (Model C) and Langevin (Model A) dynamics in 2+1D and 3+1D systems.Clearly, calculating the dynamic exponent z proved to be more of a challenge than expected, as e.g.deviations of the universal scaling function of the spectral function from their asymptotic behaviour give rise to sizeable systematic uncertainties in the corresponding scaling analysis, and also prohibit a precise extraction of z from fitting the asymptotic power law dependence of the scaling functions.Similarly, the extraction of the correlation times ξ t itself is difficult, and small amplitudes hinder a precise determination of the exponent z.Nevertheless, by combining the results of the three different methods to assess the systematic uncertainties, we arrive at plausible values, which are largely compatible with earlier studies and analytic results.In case of 2D Model A, we find our result closer to the one from experiment [2] than the Monte Carlo result [7], but not by a significant margin.For 2D Model C as well as 3D Models A and C, our results do not deviate significantly from scaling relations and earlier results.
By virtue of magnetic scaling, we are further able to predict the dynamic critical properties of the system in a certain radius around the critical point, which provides the baseline for studying possible signatures of criticality for non-equilibrium systems which approach the critical point in the phase diagram.Indeed, one particularly appealing feature of the classical-statistical setup is that it can be extended to the study of (non-equilibrium) dynamical phase transitions.We believe that such classical-statistical simulations of non-equilibrium phase transitions can provide an important benchmark to different macroscopic descriptions, currently being developed in the context of the QCD critical point search [6,[51][52][53][54][55], and intend to return to this problem in a future publication.
Besides applications to non-equilibrium phase transitions, it should also be possible, with rather limited modifications of the simulation setup, to extend this study to systems with a conserved order parameter (Models B and D) or systems with continuous symmetries such as the relativistic O(4)-model (Model G), to arrive at a more systematic characterization of the dynamic scaling behavior in the vicinity of a critical point.

Appendices
A. Integrating the spectral function In the ordered phase, at low temperatures τ ă 0, the least-squares fit of the ansatz in (79) does not produce satisfying results.This might be due to strong non-critical oscillations at rather early times.We try to remove them by integration using the following procedure.
If one integrates the ansatz ρptq " Ct ζ expp´t{ξ t q, one finds by definition of the upper incomplete Gamma function: The result unfortunately is not exactly proportional to the correlation time ξ t , but to a power thereof.However, if take a ratio of integrals, where we modify the power of the time variable t in one integrand, we get a result with the desired proportionality: By introducing extra exponents ∆, x ą 0, we allow for emphasizing the long-time exponential over the early, non-critical fluctuations.This is illustrated in Fig. 18: The additional power law shifts the weight of the integrand towards the tail of the exponential, to time scales of the correlation time ξ t .The proportionality factor then is a power of a ratio of incomplete Gamma functions, which rather strongly depends on the exponent ζ.Since we cannot calculate this ratio explicitly from the data, we have to limit this method to regimes where we can assume ζ to be a constant.

B. Structure of the spectral function at low temperatures
In the ordered phase, at τ ă 0, the spectral functions shows not only a quasi-particle structure, but also a second mode at lower frequencies.If one plots the spectral function over the squares of spatial momentum and frequency, it becomes apparent that the two structures are separated by the light cone, i.e. the low-frequency excitation inhabits the space-like region.Tracking the maximum ρpω max , pq " max ρpω, pq| ω 2 ăp 2 of the spectral function in that region, one finds for the dispersion relation roughly a power-law behavior of the form,  where ω 2 ´p2 ą 0 are plotted with green lines, space-like parts (where ω 2 ´p2 ă 0) with red lines.The colormap-projections at the bottom also show solid lines where ω 2 ´p2 " 0, indicating the light-cone.In the bottom panels we show the corresponding dispersion relations of the quasi-particle peaks (above) and the soft modes (below the light-cone), by tracing the positions of the respective local maxima in the spectral function.For comparison, we also plot solid lines, for an ideal sound-wave dispersion relation ω 2 " p 2 {3 in d " 2 (left), and for cappillary waves with ω 2 " p 3 in d " 3, to guide the eye.At high spatial momenta p, the soft mode dissolves in 3+1D.
in d spatial dimensions.For d " 3 this agrees with the well-known dispersion relation of thermally driven capillary waves.In d " 2 spatial dimensions the situation seems less clear.In Ref. [56] for example, the excitation specrta of two-dimensional fluid droplets have been studied with resulting dispersion relations that depend on the details of the fluid parameters.
In Fig. 19 we show exemplary low-temperature spectral functions in d " 2 and 3 spatial dimensions, as functions of frequency and momentum squared.The time-like and space-like parts are separated by different colors, and the light-cone is shown as a solid line in the colormap projection in the bottom plane.The light-cone marks the separation between the quasi-particle and the soft mode.In d " 3 dimensions this soft mode closely follows the power-law dispersion relation of thermally driven capillary waves as seen in the bottom right panel of Fig. 19 where we plot a solid line with ω " p 3{2 for comaprison.Bottom left we show the corresponding quasi-particle mode above and soft mode below the light-cone in d " 2 dimensions together with an ideal sound-wave dispersion ω " p{ ? 3 to guide the eye.Close to the critical temperature, the low-frequency part grows in magnitude and seems to merge with the quasi-particle.In the symmetric phase, at τ ą 0, there is only the quasi-particle peak left.

χ 2 χ 2 Figure 1 :
Figure1: Near-critical behaviour of the order parameter M and susceptibilities χ in d " 2 (left panels) and d " 3 (right panels) spatial dimensions.Solid lines show fits to Eq. (24) ff., with the corresponding χ 2 -values given in each figure.By taking into account the first two sub-leading corrections to scaling, we obtain excellent agreement between data and fits.

Figure 2 :
Figure 2: Extraction of the correlation length ξpτ, Jq from the two-point function Ḡpnq.Upper panels show exemplary data for the slice correlator Ḡpnq at different temperatures τ ; lower panels show the effective correlation length ξ eff.pxq defined by Eq. (20).We extract the correlation length ξpτ, Jq by looking for a constant plateau in ξ eff.pxq and fitting a constant, if possible.Black lines with vertical marks in the lower plots indicate both range and result of the fit.d" 2 d " 3

Figure 3 :
Figure3: Near-critical behaviour of the spatial correlation length ξpτ, Jq as function of τ at J " 0 (top panels), respectively as function of J at τ " 0 (lower panels).Solid lines show fits to Eq. (28) with the corresponding χ 2 -values indicated in each figure.

3 Figure 4 :
Figure4: Overview of the behaviour of the spectral function ρpω, pq for Hamiltonian dynamics (Model C, γ " 0) at different points in the phase diagram, above Tc (top panels), near Tc (central panels), and below Tc (bottom panels) in 2+1 (left panels) and 3+1 (right panels) dimensions.Heat maps at the bottom of each panel visualize support and spectral strength in the pp, ωq plane.The frequency (ω) and spectral function (ρ) axes are scaled logarithmically, the momentum (p) axis linearly.The spectral functions ρpω, pq for the smallest non-zero momentum mode are highlighted by a black solid line on the 3D surface.Spectral functions in the symmetric phase T " Tc are dominated by a quasi-particle peak with relativistic dispersion relation.Near the critical point T « Tc, a strongly infrared-enhanced mode emerges in addition to the quasiparticle peak that persists at higher momenta.Below the critical temperature T !Tc, a soft second mode at low frequencies emerges for finite momenta.

Figure 5 :
Figure5: Dispersion relation and damping rate in the symmetric phase.Upper panels show the spectral functions in lattice units for vanishing Langevin coupling γ " 0 with an exemplary Breit-Wigner fit as solid line; lower panels show the resulting fit parameters in lattice units over the different spatial momenta, with a fit to a relativistic dispersion relation Eq.(41).Note that the width Γ is scaled by a factor of 5 to improve legibility.

Figure 6 :
Figure6: Dispersion relation and damping rate in the ordered phase, illustrated analogously to Fig.5.Due to the presence of the second excitation at lower ω, the Breit-Wigner fit does no longer describe the spectral function data as well as in the symmetric phase.Hence, the fit window for the Breit-Wigner fit was narrowed to the high-frequency regime.

Figure 7 :
Figure7: Effective quasi-particle masses and damping rates as a function of reduced temperature.Near the critical point (τ « 0) the effective mass drops and the damping rate increases; while at low momentum p the spectral function is no longer described by a simple quasi-particle peak, the behaviour of high momentum modes is remarkably smooth across the phase transition.

Figure 8 :
Figure 8: Scaling function fωpxω, 0q of the critical spectral function at non-zero spatial momentum for the four different critical scenarios.Solid black lines represent the amplitude fω obtained from the fit of the universal scaling function in Fig. 10 at large xp.Due to availability of data, we use τ " 0.0009p2q (d " 2) resp.τ " ´0.00008p5q (d=3) as a proxy for the critical temperature.Note that, despite τ « 0, finite size effects are not relevant here due to finite spatial momentum.However, to achieve the comparatively very low p in d " 3, a single set of data at N " 512 very close to Tc was generated.

Figure 10 :
Figure10: Scaling function fppxp, 0q of the spectral function for the four different critical scenarios.Due to availability of data, we use τ " 0.0009p2q (d " 2) resp.τ " ´0.00008p5q (d=3) as a proxy for the critical temperature.Solid lines represent the fit to the ansatz in Eq. (57), we fit fω and ap as free parameters, and the exponent of the power law for large xp is given by p2 ´ηq{z comb. .Note that finite size effects are not relevant here due to finite spatial momentum.To achieve the comparatively very low p in d " 3, a single set of data at N " 512 very close to Tc was generated.

Figure 14 :
Figure14: Behaviour of the auto-correlation time as a function of the finite-size scaling variable |τ |L 1{ν .Different symbols correspond to extractions of the correlation time using fits (filled symbols) and integration (open symbols) of the spectral function; upper curves in each panel correspond to ξtpτ ą 0q, while lower curves correspond to ξtpτ ă 0q.Solid lines in each panel correspond to a power law fit according to Eq. (71), from which we extract the dynamic critical exponent z along with the non-universal amplitudes f t .

Figure 15 :
Figure15: Behaviour of the correlation time ξtpτ q extracted from fits and integration of the spectral function, as a function of reduced temperature τ .Note that in order to minimize finite volume effects, we only include results for |τ |L 1{ν ą 2.5 for d " 2, and |τ |L 1{ν ą 5 for d " 3 obtained on 256 2 and 128 3 lattices.Solid black lines show a fit according to Eq. (72); the regular contribution is shown separately as a black dashed line.

z 18 −Figure 16 :
Figure16: Dependence of the ratio of auto-correlation times in Models A and C, on the finite size scaling variable τ L 1{ν .Solid lines in the left and right panels show a power law fit of the large τ L 1{ν behavior, from which we extract the difference νpz C ´zA q of the dynamic critical exponents (cf.Eq. (73)).Dashed lines indicate the confidence interval of the extraction of νpz C ´zA q, which is also presented in the figure.

Figure 17 :
Figure 17: Magnetic scaling of the correlation time ξtpτ, Jq.Overlapping data points correspond to the universal magnetic scaling function of the auto-correlation times ξt,J pτ {| J| 1{βδ q (cf.Eq. (77)).Scaling breaks down when J becomes too large, as indicated by the blue data points.

Figure 18 :
Figure 18: Integrands of Eq. (81), in this case Model A spectral functions.Clearly, the late-time tail is amplified in comparison to early fluctuations.

Figure 19 :
Figure19: Spectral functions in the low-temperature phase for Hamiltonian dynamics (Model C, γ " 0) in d " 2 (left) and d " 3 (right) spatial dimensions.Time-like parts (where ω 2 ´p2 ą 0 are plotted with green lines, space-like parts (where ω 2 ´p2 ă 0) with red lines.The colormap-projections at the bottom also show solid lines where ω 2 ´p2 " 0, indicating the light-cone.In the bottom panels we show the corresponding dispersion relations of the quasi-particle peaks (above) and the soft modes (below the light-cone), by tracing the positions of the respective local maxima in the spectral function.For comparison, we also plot solid lines, for an ideal sound-wave dispersion relation ω 2 " p 2 {3 in d " 2 (left), and for cappillary waves with ω 2 " p 3 in d " 3, to guide the eye.At high spatial momenta p, the soft mode dissolves in 3+1D.

Table 2 :
Asymptotic amplitudes and corrections to scaling, as obtained from fits to Eq. (24) ff.The given uncertainties only include statistical errors.If a correction amplitude does not improve the χ 2 {d.o.f. of a fit, it is set to zero and denoted by a dash.

Table 3 :
[41]ersal amplitude ratios.The given uncertainties only include statistical errors.The ratios are shown alongside their literature values, which were taken from[41].

Table 4 :
Extracted fit parameters from fitting Eqs.(57) and (58) to the scaling functions.The first line corresponds to the results of the fits of Eq. (64) with z as a fit parameter, with an estimate of the systematic error.The other lines represent the values obtained by fitting the data in Figs. 10 and 11 with a fixed value for the dynamic critical exponent, corresponding to the last row of Table6labeled combined, see below.For the latter ones we only give the statistical fit error.

Table 5 :
Non-universal amplitude f t and universal amplitude ratio U ξ,t of the correlation time ξt, obtained by fits of the data to Eq. (69), shown in Fig.14.Since the extraction of the amplitudes is strongly correlated with the extraction of the exponent, the uncertainties are rather large.

Table 6 :
Extracted values for the dynamic critical exponent z.See Section 6.4 for detailed discussion of the different extraction methods.{| J| 1{βδ that combines the τ and J dependence.Similarly, one expects an analogous magnetic scaling behaviour for unequal-time correlation functions, such that e.g. the autocorrelation time ξ t is expected to obey the scaling form ξ t pτ, Jq " s z ξ t ps 1{ν τ, s βδ{ν Jq,