Dynamics of the $O(4)$ critical point in QCD

We perform a real-time simulation of the $O(4)$ critical point of QCD, which lies in the dynamic universality class of Model G. The axial charge and the order parameter $\phi_a =(\sigma, \vec{\pi})$ exhibit a rich dynamical interplay, which reflects the qualitative differences in the hydrodynamic effective theories above and below $T_c$. From the axial charge correlators on the critical line we extract a dynamical critical exponent of $\zeta=1.47 \pm 0.01 ({\rm stat})$, which is compatible with the theoretical expectation of $\zeta = d/2$ (with $d=3$) when systematic errors are taken into account. At low temperatures, we quantitatively match the $O(4)$ simulations to the superfluid effective theory of soft pions.


I. INTRODUCTION
Chiral symmetry breaking and the chiral phase transition play a prominent role in QCD at finite temperature. In the limit of two massless flavors, the transition from a chirally restored phase T > T c , to a chirally broken phase T < T c , is second order and is in the O(4) universality class [1,2]. Although the static properties of the O(4) critical point have been studied in detail both numerically and theoretically [3][4][5][6][7][8][9][10][11][12], the dynamic scaling properties of the critical region demand additional study, which is the goal of this work.
This study may appear academic: in the real world the finite quark mass explicitly breaks chiral symmetry, reducing the influence of the O(4) critical point on the static and dynamic correlators of QCD at finite temperature. However, lattice QCD computations of the chiral condensate as a function of quark mass show that the qualitative and, to some extent, even quantitative properties of the chiral crossover can be understood using static O(4) scaling functions [13,14]. These scaling functions predict the singular behavior of the chiral condensate near the pseudocritical temperature T pc and other static observables. Motivated by the lattice effort, we will perform the real-time simulations of the O(4) critical region, which we hope can provide an analogous understanding of the scaling of dynamic correlators in QCD in the crossover region.
The current study is also motivated by experimental results on the momentum spectra of particles produced during the collisions of heavy ions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). In much of the accessible momentum range these spectra are remarkably well described by ordinary viscous hydrodynamics [15]. From a theoretical perspective, the relevant symmetry group of QCD close to the chiral limit is approximately SU L (2) × SU R (2), leading to the conservation of iso-vector charge, and the approximate conservation of the iso-axial vector charge. The corresponding densities should be included as additional fields in the hydrodynamic description. When chiral symmetry is spontaneously broken, the pions π (which are the associated Goldstone bosons) should also be added to the hydrodynamic fields, and the appropriate hydrodynamics resembles a non-abelian superfluid [16,17]. Finally, close to the O(4) critical point where the σ meson is also light, the O(4) order parameter φ a ∼ (σ, π) also should be included as an additional hydrodynamic degree of freedom [2,18]. Current hydrodynamic simulations do not include the iso-axial charge, the pions, or the order parameter as explicit hydrodynamic variables. By including these variables as explicit degrees of freedom we hope to increase the predictive power of hydrodynamic simulations in the crossover region.
In fact, there is an excess of soft pions relative to the predictions of current ordinary hydrodynamic simulations [19][20][21]. We have previously suggested that this excess may reflect the cavalier treatment of chiral symmetry breaking and the O(4) transition in almost all hydrodynamic simulations of heavy ion collisions to date [18]. To corroborate this suggestion, we will need to simulate the real-time dynamics of the O(4) phase transition for an expanding fluid. As discussed in [18], the dynamics of the order parameter matches smoothly onto a pion-hydro effective field theory (EFT) for T T c ; this pion EFT subsequently matches onto a kinetic description for soft pion particles coupled to the background fluid flow [17]; and finally, the kinetics can be used to propagate pions to freezeout with definite predictions for soft pions yields and their correlations. In addition there is an experimental proposal to measure soft pions and their correlations over a wide range in rapidity [22], which is ideally suited to unravel this physics and to probe the (in)applicability of ordinary hydrodynamics in this regime.
As a first step, in this paper we will compute the real-time correlation functions of an O(4) critical system and study their scaling properties. There is considerable theoretical interest in the critical correlators themselves. Many years ago Rajagopal and Wilczek determined that the dynamic universality class of QCD is similar to "Model G" of [23], where the order parameter φ a = (σ, π) is not conserved, but has a non-trivial Poisson bracket with vector and axial vector charges [2]. They also determined dynamical critical exponent to be ζ = d 2 , which we find in the simulations presented here. Because the critical theory must transition between ordinary hydro at high temperatures and a non-abelian superfluid hydro at low temperatures, the expected structure of the hydrodynamic correlations functions is rich [18]. It would be nice to see this structure in a simulation.
Earlier numerical studies on the critical dynamics of field theories (including O(4) symmetric ones) have been performed in the "classical-statistical" framework [24][25][26]. Given some relativistic quantum field theory, the high-temperature spectral functions are saturated by their classical counterparts close to the critical point. Since the non-anomalous symmetries and conservation laws of the classical field theory are shared with the quantum one, the classical dynamics belongs to the same dynamic universality class as the full quantum theory. Of particular relevance to our work was the study done in [25], which studied a classical relativistic O(4) model, and determined the spectral functions of the order parameter. The spectral functions were shown to display the appropriate behavior as a function of temperature, and pion quasiparticle poles were observed in the broken phase. Because the classical model has O(4) Noether charge densities, n ab ∼ φ [a ∂ t φ b] , which have a non-trivial Poisson bracket with the order parameter, the dynamics of this model should lie in the universality class of "Model G". However, within this set-up studying the "Model G" dynamics is difficult, since rapidly oscillating UV modes (which build up the charge densities microscopically) must be carefully evolved. Consequently, a conclusive extraction of the dynamical critical exponent was not possible, and the interplay between the order parameter and the axial charge was not studied. Very recently [27], the same group adopted an approach for "Model B" and "Model D", which is somewhat closer in spirit to the one taken here for "Model G", where the charge densities are treated as additional slow variables. In their recent work, Israel-Stewart-like diffusion models belonging to the specified dynamical class ("B" or "D") were simulated, and a careful study of the dynamical scaling function and of their momentum dependence was performed. In particular, this work constitutes an important stepping stone towards the study of "Model H", which is believed to describe the universality class of the speculated QCD critical point [28].
An outline of the paper is as follows. In Sect. II we discuss the model equations we will solve. Of particular interest is the numerical strategy presented in Sect. II C, which may be useful for other model systems. In Sect. III A we will review the thermodynamics of the model and fix the non-universal (thermodynamic) parameters of the model. Finally in Sect. IV we turn to the dynamical properties of the model presenting the principal results. In Sect. IV B we present a qualitative overview of the phase transition, and examine the dynamics in the chirally restored limit. Then in Sect. IV C, we examine the low temperature limit where the O(4) dynamics should match with the pion EFT. We examine the Gell-Mann-Oakes-Renner relation, and the dissipative pion dynamics proposed by Son and Stephanov [29,30]. In the last section, Sect. IV D, we examine the scaling of correlation functions along the critical line. We extract the dynamical critical exponent ζ and find ζ 1.47 ± 0.01 (stat.), which is very close to the predictions of Rajagopal and Wilczek of ζ = d/2. Finally, a short outlook is presented in Sect. V.

II. MODEL
A. Model equations QCD with two degenerate massless quarks is well known to have a second order phase transition and is in the universality class of the O(4)-critical point. Dynamical properties of a theory near a continuous phase transition are also universal, but theories with the same static properties can lead to different dynamical universality classes. Different dynamics arise because of the existence or non-existence of conserved charges in the theory [23]. As pioneered in [2] (see [18] for a recent review) the dynamics of the QCD O(4)-critical point is the one of an O(4) antiferromagnet, "Model G" of [23].
Model G consists of an O(4) order parameter 1 φ a = (σ, π) field, and adjoint charge densities n ab . The field φ is a proxy for the quark condensate q R q L , and as a result is not a conserved quantity. The antisymmetric tensor of charge densities n ab can be decomposed into a vector part, n s V = 1 2 ss 1 s 2 n s 1 s 2 , and an axial part, n s A = n 0s . They represent the original iso-vector n V ∼qγ 0 t I q and iso-axial n A ∼qγ 0 γ 5 t I q charge densities. The vector current is exactly conserved for equal quark masses, while the axial current is only approximately conserved, since the finite quark mass explicitly breaks the chiral SU L (2) × SU R (2) ∼ O(4) symmetry. The equilibrium action (or effective Hamiltonian) in the presence of an external magnetic field H a = (H, 0) parametrizing the explicit symmetry breaking, takes a Landau-Ginzburg form Here n 2 = n ab n ab and with m 2 0 negative. The relevant equations of motion for these fields are Here, for example, H [a φ b] denotes the anti-symmetrization, H a φ b − H b φ a . χ 0 is the isovector and the iso-axial-vector charge susceptibility; these susceptibilities are equal and 1 Here a and b denote O(4) indices; s, s 1 , s 2 , etc. denote the isospin indices, i.e. the components of π; finally, spatial indices are notated i, j and k. The dot product indicates an appropriate contraction of indices when clear from context, e.g. φ · φ = φ a φ a , π · π = π s π s , and ∇ · ∇ = ∂ i ∂ i . approximately constant near the critical point. µ ab is the chemical potential, n ab /χ 0 . The coefficients Γ 0 and σ 0 are the bare kinetic coefficients associated to the order parameter and the charges. The bare diffusion coefficient of the charges is D 0 = σ 0 /χ 0 . The constant g 0 is a coupling of the field φ, and has the units of (action) −1 in our conventions. Finally, θ a and Ξ ab are the appropriate noises, which are defined through their two-point correlations [23] The equations of motion naturally break up into an ideal hydrodynamic evolution (the left hand side of the equations) with viscous damping (the right hand side of the equations). If the right hand side is set to zero, it is easy to show that the ideal evolution leaves H = const. More generally one can write down the Fokker-Plank equation associated with the stochastic process and straightforwardly show that the equilibrium probability distribution is The thermodynamics of this model is recalled in [18] and we will determine some static properties of relevance in Sec. III A. The real-time correlation functions we will study here are where . . . c refers to a connected two point function. We will limit this study to k=0. Close to the critical point, the dynamics is expected to be controlled by "dynamic scaling" [23]. In particular, we expect the time development of our two-point functions to scale with the correlation length ξ as Here the functions χ (k), χ ⊥ (k) are the static order parameter susceptibilities, and depend on k and ξ; χ 0 is corresponding charge susceptibility which lacks these dependencies; Ω is a non-universal constant normalizing the time; finally, z is the familiar static scaling variable involving the reduced temperature and magnetic field (see below). Y σ , Y π , and Y A are universal dynamical scaling functions and ζ is the corresponding dynamical critical exponent of the theory. The expected dynamical critical exponent for "Model G" is ζ = d/2 [2]. The scaling form (with ζ = d/2) implies that if the correlation length increases by a factor of two, then the characteristic relaxation time increases by a factor of 2 3/2 , thereby exhibiting a "critical slowing down".

B. Lattice units and matching the model to QCD
To simulate the model, we begin by taking g 0 and T c as our microscopic units of (action) −1 and energy, respectively, setting g 0 = T c = 1 in the computer code. Similarly, we will choose a microscopic length a as the cutoff in our problem, setting the lattice spacing to unity in the code.
As a result of these choices, the quantities we measure directly from our simulations are expressed in lattice units, and they are dimensionless numbers. To convert these quantities to physical predictions, we need to assign a physical value to g 0 , T c , and a. The critical temperature T c can be matched directly to the QCD critical temperature. Once T c is fixed, g 0 T c is adjusted so that the model reproduces the pole frequency of the pion. Lastly, the cutoff a can be adjusted so that our system reproduces the correlation length of QCD. The aim of this section is to explain this procedure in greater detail. Before doing so, let us note that our set of units, g 0 = T c = a = 1, will be implicit both in the figures and the text. However in this section, and if necessary for clarity, we will adopt a "hat" notation for variables in lattice units, e.g.n = n a 3 andχ = T c χ a 3 are the dimensionless charge density and charge susceptibility, respectively.
The model has an O(4) critical point at a critical mass parameterm 2 c (λ). At infinite volume and close to the critical point, the dependence of the model condensate (in units of T c /a) on the mass parameter and magnetic field takes the conventional scaling form [6] where δ is the critical exponent, and f G (z) is a universal function with f G (0) = 1. Here h is the reduced magnetic field and z is the scaling variable, whilem 2 c (λ),Ĥ 0 (λ), and m 2 (λ) are order one non-universal constants that are fit to our numerical data on thermodynamics (see Sect. III A and eq. (26)). In physical unitŝ In QCD, the chiral condensate close to the critical point takes the same form qq /B QCD = h 1/δ f G (z), but with scaling variables Evidently, to match the two systems we are to equate the scaling variables, h and z, and equate the order parameters: qq For H small and T < T c , the universal function f G (z) behaves as z β , and the model condensate takes the form:σ while chiral condensate takes an analogous form at the corresponding z providing an explicit map between (T − T c )/T c and the mass parameter of the model. The constant B QCD has units of (meters) −3 and H QCD has unit of energy. They can be chosen arbitrarily, but not independently, as the parameter fixes a microscopic unit of length. The diverging correlation length of QCD near the critical point is a universal function times this length [31]. As we show in App. A 1, by choosing the model will reproduce both the correlation length and pole frequency of the pion in QCD.
Having set three of our parameters to unity to fix our units of space, time, and energy, we are still left with three more dimensionless parameters which must be specified, namelŷ The susceptibility χ 0 sets the magnitude of charge fluctuations relative to the fluctuations of the order parameter, while Γ 0 and D 0 determine the relaxation of the order parameter and the charge diffusion, respectively. Switching to the conventional = c = 1 units for this paragraph, for the system under study there really is only one scale T c ∼ Λ QCD . We expect that the microscopic (i.e. cutoff) length and time are both of order ∼ 1/T c . The susceptibility in units of T c is also of order unity. Indeed, we expect that all dimensionless constants are of order unity, and therefore, in this study we will takeχ for definiteness. It may be worthwhile to explore the dependencies on these parameters further, but we have not done so here.

C. Numerical strategy
To simulate the real-time dynamics, we will discretize the stochastic evolution equations in (3a) and (3c), placing the system on a spatial lattice of size L and volume V = L 3 . We briefly present our algorithm in this section; the interested reader can find detailed explanations in App. A.
As the equations naturally separate into an ideal and a dissipative part, we use an "operator splitting" approach. In spirit, we first evolve our fields for a short time, neglecting the dissipative part where spatial derivatives are discretized appropriately. We then neglect the ideal part and solve for the dissipative dynamics Decoupling the equations in such a way allows us to use methods specifically tailored to the two different dynamics. In particular, we use a symplectic integrator to evolve the ideal part, preserving in this way the underlying Poisson bracket structure. To simulate the dissipative Langevin dynamics, we use a Metropolis algorithm. A similar strategy to simulate the Langevin dynamics was used previously to calculate the sphaleron transition rate in hot non-abelian plasmas 2 [32]. At every lattice site x, the order parameter is updated as where for each flavor index a the increment is Here ξ 0 is a random number with unit variance ξ 2 0 = 1. The update proposal is accepted with probability min(1, e −∆H ), where ∆H is the change in the discretized Hamiltonian. If the proposal is rejected, then φ(t + ∆t, which can be used to show straightforwardly that the mean and variance of the accepted proposals reproduce the dissipative and stochastic terms of the Langevin process (see also App. A for more details) The charges are updated in a similar way, with the extra difficulty that the noise term generated by the updates must be a total divergence. This is tackled by updating the lattice 2 In the sphaleron case the time scales between the metropolis and Langevin times must be carefully

cells in pairs, making a Metropolis proposal for the charge transfer between two cells -see App. A.
Using a Metropolis update to solve for the non-ideal part of the dynamics has several advantages over a direct time evolution of the Langevin process. For instance, it allows us to design a scheme whose equilibrium properties are independent of the time-stepping. It also allowed us to use larger time steps compared to a naive discretization of the equations of motion. The numerical code is implemented with PETSc and MPI [33,34].

A. Thermodynamics
Our goal in this section is to fix the static non-universal parameters of the model from its thermodynamics. The magnetization of the system is an average over the volume at a given time moment: and its time average, denoted with . . . , determines the condensateσ At infinite volume, the dependence of the condensate on the temperature and magnetic field takes the scaling form given in (8). The non-universal constants m 2 c , H 0 , m 2 are fit to our numerical data onσ. We first determine m 2 c , then we simulate on the critical line with m 2 0 = m 2 c to determine H 0 , and finally, we simulate at H = 0 to find m 2 . Anticipating the results of this section, we obtain with λ = 4 Following standard technique [35], we determined the critical coupling of the model m 2 c by measuring Binder cumulants and determining when they cross a nominal value, which was taken from previous simulations [36]. Further details are given in App. B 1.
To determine H 0 we made a scan on the critical line, with details presented in App. B 2. The data forσ on 32 3 and 64 3 lattices on the critical line are shown in the left panel of Fig. 1. They were fit to a a finite size functional form given by Engels and Karsch [6], which fixes the value of H 0 given in (26), and quantifies finite size corrections. The fit is reasonable and has χ 2 /dof=2. The magnetization at infinite volume from the results of this fit are shown by the dashed line. We see that already at L = 64 we are essentially at infinite volume for the range of H considered in this work. Our dynamical simulations in Sect. IV are all done with L 3 = 80 3 . This analysis on the critical line suggests that finite volume corrections are modest.
In the next step we performed simulations at H = 0 with T < T c , in order to fix the nonuniversal constant m 2 . Details are presented in App. B 3. The infinite volume magnetization Σ at zero field is defined as The fits and extraction procedure are discussed in the text. Also shown is the fit result without the subleading correction.
Extracting the magnetization Σ is difficult as, in any finite volume, This is because when HΣV ∼ 1, the orientation of magnetization vector M a begins to wander on the group manifold, averaging to zero in the limit of zero external magnetic field. One way to extract Σ is to look at the fluctuations of M a , evaluating M 2 = M a M a , which is approximately Σ 2 at large volume. The leading deviation of M 2 and Σ 2 at finite volume comes from the fluctuations of long wavelength Goldstone modes, and can be neatly analyzed with a Euclidean pion EFT [37]. We detail these corrections in App. B 3, which were essential to a reliable extraction of Σ(T ).
Our results for Σ(T ) are shown in the right panel of Fig. 1, and are fit with the functional form with critical exponents β and δ from [6] and ω from [36]. Here we are using instead oft r , and we defined b 1 ≡ (|m 2 c |/m 2 ) β . The second term in (29) captures the first subleading correction to scaling.
Our fit to Σ(T ) is shown in the right panel of Fig. 1 and yields b 1 = 0.544(4) and C T = 0.20(2) with a χ 2 /dof = 1.4. We have excluded the largest value of (−t r ) from the fit. For comparison, we also show the fit results for the first term b 1 (−t r ) β . Clearly, for precision work the subleading corrections are important in the temperature range we are considering.
The parameter b 1 determines the scale m 2 described earlier (i.e. m 2 = |m 2 ) yielding the results presented in (26).
To summarize, in this section we have established the non-universal parameters m 2 c , H 0 , and m 2 which determines the map between the model and the conventionally parameterized O(4) critical point. The results are given in (26).

B. The Static Pion EFT and Gell-Mann-Oakes-Renner
Before turning to the dynamics we will determine the validity of the Euclidean pion EFT discussed above, relegating all details to App. B 4. At all temperatures, O(4) symmetry guarantees that the transverse susceptibility is determined by the condensateσ where G ππ (k) is the static correlation function. At low temperatures the magnitude of the condensate φ 2 is approximately frozen toσ, and the long wavelength order parameter fluctuations are determined by the fluctuations in the phase ϕ, π s (x) σϕ s (x). The static action for the Gaussian effective theory describing the phase fluctuations takes the form [29,38] and makes a definite prediction for the static correlator where f 2 is the decay constant and m is the screening mass. Comparing the predicted correlator to the susceptibility yields the Gell-Mann-Oakes-Renner (GOR) relation At a finite negative z, the GOR relation is only approximate, receiving corrections due to fluctuations of the σ field. We have fit the static ππ correlator to find the decay constant f 2 and the screening mass m 2 at a nominal point in the broken phase, z = −2.2011 and H = 0.003 (see Fig. 2). Comparing f 2 m 2 to Hσ yields Evidently, already at z = −2.2, the Euclidean pion EFT works to better than a percent.
Having studied the statics of the pions, in the next section we will turn to the dynamics, making use of these results in Sect. IV C.
A parametrization of the longitudinal susceptibility χ ∝ f χ (z) and mean magnetization σ ∝ f G (z) taken from the simulations of Engels and Karsch [6]. The black points are the z values which will be simulated in this work. Further simulation details are given in Fig. 3.

A. The simulations
To get an overview of the phase diagram, in Fig. 2 we show the scaling function of the magnetization f G (z) and the corresponding function for longitudinal susceptibility f χ [6] The susceptibility shows a prominent maximum at the pseudocritical point with z value of z pc 1.35. In order to scan the dynamics of the transition, we have performed real-time simulations at the black points. We also made a scan on the critical line z = 0 for various values of the magnetic field. The dynamical parameters as well as the run times and other information is gathered in Fig. 3.

B. Overview
We will start by presenting an overview of the critical dynamics as the temperature is scanned across the phase transition. At high temperatures, the order parameter is small and simply dissipates through the damping term in the equations of motion. Since there is no preferred direction, the longitudinal and transverse order parameters excitations, δσ and π, are nearly degenerate. In the vector channel, the total charge is constant in time and the dissipation affects only non-zero Fourier modes, which are not studied here. The situation is different in the axial channel, since the axial charge is not conserved. However, the explicit symmetry breaking term in the action, Hσ, is tiny, since it is proportional to the magnetic field H (or quark mass) and the order parameter, which is small at high temperatures,σ ∝ H. As a result, the axial charge will dissipate rather slowly over a time scale of order Hσ/χ 0 ∝ H 2 . In this regime, the dynamics of the axial charge is unrelated to the pions. However, as we lower the temperature, the order parameter acquires an Hindependent expectation value, and the axial channel gets modified; the order parameter field and the axial charge are now entangled. In the deeply broken phase at low temperatures, the axial charge and the transverse part of the field will no longer just dissipate. Indeed, their dynamics become intrinsically locked, and they acquire the quasiparticle characteristics of the Goldstone modes associated with the broken symmetry. By contrast, in this regime the longitudinal excitation of the order parameter (the σ) has a large mass and its dynamics remains purely dissipative.
These qualitative behaviors are precisely observed in our data. In Fig. 4, we start by showing the results of a simulation performed in the unbroken phase, z = 3.87. In the left plot we show the statistical correlator for the σ, π, and axial channels as a function of time. Noting that the x-axis is on a logarithmic scale, the slow dissipation in the axial channel is apparent. It is also apparent that the σ and π channels are almost degenerate and dissipate on a much shorter timescale. This is also clearly seen in the corresponding Fourier transforms (right), where σ and π correlators appear as a single dissipative peak, which is much broader than corresponding peak in the axial correlator.
In Fig. 5, we show the behavior of the axial charge correlator at the pseudocritical and critical temperatures, and at lower temperatures, in the broken phase. In the left panel, we show the correlation functions as a function of time, while in the right panel we show their Fourier transforms. At the pseudocritical temperature (the red curves), the axial charge correlator is still purely dissipative, but the peak is much broader than in Fig. 4, indicating that the charge is no longer approximately conserved. As we lower the temperature to T c (the blue curves), we start seeing the emergence of propagating pions, which appear as oscillations in the correlator as a function of time, or equivalently, as quasi-particle peaks in the Fourier transform. At the critical temperature there are no drastic changes (this is expected in a finite magnetic field), and the correlator behaves as it does in the broken phase, with propagating pions which are clearly visible in the axial channel. As one moves further into the broken phase (the purple curves), the pion peaks become increasingly separated, and the real-time pion EFT discussed below becomes valid (see Sect. IV C). It seems that around the pseudocritical temperature z pc (the red curve) the axial charge propagator starts changing its behavior from purely dissipative to quasiparticle-like. Indeed, at z = 0.7z pc (the green curve), i.e. slightly below the pseudocritical point, the dissipative peak is already quite deformed, which reflects the nascent formation of the two quasiparticle peaks.
In the left and right panels of Fig. 6, we show the corresponding statistical correlators for the π and σ fields as a functions of frequency, with z spanning the phase transition. In the deeply unbroken phase the two channels are mostly indistinguishable (the grey bands), as pointed out before. Lowering the temperature to the pseudocritical point, the pseudoscalar channel acquires a double peak structure, while the scalar channel remains purely dissipative. Going further down in temperature, the quasi-particle peaks in the pseudoscalar channel separate. Interestingly at z pc , the pion correlator already has a quasiparticle peak, while the axial charge correlator is still dissipative (Fig. 5); only past the pseudocritical point do their correlation functions become closely related.

C. Broken phase: pion EFT
Deep in the broken phase, the fluctuations of the order parameter are dominated by the phase fluctuations π s (t, x) σϕ s (t, x), which are tightly correlated to the axial charge fluctuations through the Josephson constraint, ∂ t ϕ µ A . The dissipative hydrodynamic theory for the phase fluctuations has been worked out in [16,17,29], and provides a real time analog of the static Gaussian effective theory described in Sect. III B.
The linear response of hydrodynamic theory has been analyzed in [18,29], and the hydrodynamic prediction for the dynamical correlators in the k = 0 case is Here m 2 p is the pole mass of the pion excitation, m is the transverse static screening mass, Γ is a dissipative coefficient correcting the Josephson constraint, and finally χ 0 and χ ⊥ are the appropriate static susceptibilities, which are required to normalize these expressions The fact that the pions are pseudo-Goldstone bosons, and correspondingly that the axial current is partially conserved (PCAC), leads to the well-known and remarkable property that the dynamical pole mass m p can purely be computed from the static properties discussed in Sect. III B. In particular, at low-enough temperatures, we have a finite temperature Gell-Mann-Oakes-Renner (GOR) relation [29,30,39] where v 2 ≡ f 2 /χ 0 is the pion velocity. Already in Fig. 5 we saw the appearance of pion excitations. We will now try to assess the validity of the pion EFT. To do so, we attempt to fit expressions (37)-(38) from our statistical correlators. To perform these fits, we first fix the normalizations by extracting from our data the susceptibilities, χ 0 and χ ⊥ . We then use a two parameter model, involving m p and Γ p = Γm 2 , and simultaneously fit the statistical correlators in the π and axial channels.
This extraction of the pole mass allow us to verify the dynamical part of the GOR relation. Referring the reader again to Sect. III B for the corresponding extraction of the static quantities, we find We see again that already at z = −2.2, the deviations from GOR are remarkably small, of order 1%, which could be due to corrections of order ∼ (Γ p /2m p ) 2 . Note also that part of this 1% deviation could be due to some remaining systematic errors in our time evolution, see App. A 2 a for more details.

D. Critical line: dynamical scaling
Moving on to the critical line z = 0, we consider the scaling of the critical dynamics. Focusing on k = 0 modes and recalling that the (e.g. longitudinal) correlation length scales as on the critical line, the "dynamic scaling hypothesis" (7a)-(7c) gives us the following scaling forms for the correlators on the critical line with for example, Y c A H ζνc t = Y A Ω ξ −ζ t, 0, 0 . To verify the validity of the hypothesis and to determine the dynamical exponent ζ, we studied a set of simulations at m 2 0 = m 2 c for H = 0.002, 0.003, 0.004, 0.006, 0.01. We present the results obtained for the axial-axial channel in Fig. 8. The left plot shows the timedependent correlator G AA (t) for the different magnetic fields, while the right panel displays its corresponding Fourier transform G AA (ω). Qualitatively at least the curves show a scaling behavior.
To quantitatively assess the scaling ansatz (49) and to extract the exponent ζ from our data, we located the time when G AA (t, H) reaches its first minimum, t min (H), which can be determined with reasonable accuracy. From the scaling ansatz, we see that, given two magnetic fields H 1 , H 2 , we expect We show this ratio as a function of H 1 /H 2 in the left panel of Fig. 9. The data are well described by the power law form, and we obtain a nominal value for the dynamical exponent of taking ν c = 0.4024 from [6].
With an estimate of the critical exponent in hand, we can verify the ansatz (49). Indeed by appropriately rescaling times and frequencies, we expect to see our correlators G AA (t, H) and spectral function collapse to a single curve. The scaling of G AA (t, H) is shown in the right panel of Fig. 9, while the scaling ρ AA (ω, H) is shown in Fig. 10. To obtain this data collapse, we have rescaled the time and inverse frequency by (H ref /H) ζ fit νc with H ref = 0.002. Dynamical scaling is also expected to hold in the other channels and in particular we expect the same ζ fit to govern the dynamics in the σ channel. This indeed happens, which we illustrate in Fig. 11 by showing the σσ correlator (left) and the corresponding collapsed spectral function (right).
Before moving on, let us emphasize that our numerical estimate of the critical dynamical exponent is close to the critical scaling prediction [2,23], ζ = d/2. Considering, for example, the small violations of scaling seen in Fig. 9, we do not consider the deviation of ζ fit from d/2 to be significant.

V. DISCUSSION
In this work, we numerically studied the universal critical dynamics relevant to two-flavor QCD close to the chiral phase transition. More precisely, we simulated the dynamics of an O(4) antiferromagnet, "Model G" of [2,23]. After reviewing the model and explaining our conventions in Sect. II, we performed some "scale setting" in Sect. III, where we studied the thermodynamic properties of the model and extracted the relevant non-universal constants. We also determined some of the static properties of pions in the broken phase such as their screening masses and decay constants.
With this data in hand, we moved on to the main section of this work, Sec. IV, which studied the dynamics. Focusing on correlators at zero spatial momentum, we first performed a scan in temperature across the phase transition. We qualitatively confirmed that the dynamics takes place as expected, by studying the real-time correlation functions in the σσ, ππ and axial-axial channels. At high temperature, the σ and π are degenerate and the axial charge is almost conserved. In the broken phase, the σ remains purely dissipative, while the π propagates and carries axial charge. In particular, we were able to observe that the coupling of the π to the axial charge precisely happens in the vicinity of the pseudocritical point, z pc , defined as the line in the phase diagram where the static susceptibility peaks. This observation is yet another link between the static and dynamical properties of this critical model. We also performed a quantitative study of the pion properties in the broken phase. We were able to fit the dynamical correlator to a particle resonance ansatz predicted by the chiral hydrodynamic effective theory, and extract the pole mass and decay width. Furthermore, we verified that the Gell-Mann-Oakes-Renner relation, which relates the dynamical pole mass of the pions to their static screening mass, holds at the sub-percent level. Last but not least, we performed a set of simulations along the critical line and extracted the dynamical critical exponent ζ = 1.47 ± 0.01 (stat), very close to the critical scaling prediction ζ = 1.5 [2].
The numerical determination of ζ can be considered as a first step towards a complete quantitative characterization of the dynamics of the O(4) antiferromagnet. Such a characterization would include additional studies at finite spatial momentum as in [27], and a more complete investigation of the dynamics in the chiral limit at finite volume with an appropriate real-time EFT. (The corresponding finite volume static EFT was written down long ago [37], and was helpful in the thermodynamic analysis in Sect. B 3). In order to use the model to analyze heavy-ion data as discussed in [17,18], it will be important to analyze the critical O(4) dynamics for an expanding fluid, which introduces a rich hierarchy of scales. Finally, it will be interesting to apply the algorithm presented in App. A to other stochastic and critical systems.

ACKNOWLEDGMENTS
The authors would like to thank Stony Brook Research Computing and Cyberinfrastructure , and the Institute for Advanced Computational Science at Stony Brook University for access to the SeaWulf computing system, which was made possible by a National Science Foundation grant (#1531492).
In addition, this research used resources resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award 001-ERCAP6306. This work is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, grants Nos. DE-FG-02-08ER41450. AS is supported by the Austrian Science Fund (FWF), project no. J4406. E.G. is supported in part by the GLUODY-NAMICS project funded by the "P2IO LabEx (ANR-10-LABX-0038)" in the framework "Investissements d'Avenir" (ANR-11-IDEX-0003-01) managed by the Agence Nationale de la Recherche (ANR), France.
Appendix A: The model on the Lattice

Relating lattice units to QCD
Continuing the discussion of Sect. II B, our goal is to fix the spatial cutoff a in units of meters, so that the model will reproduce the physical correlation length. In the computer code the cutoff is the lattice spacing, and is set to unity. Similarly the times are measured in units of 1/(g 0 T c ), and we will set the g 0 in physical units so that the model reproduces the pion pole frequency. The results of this section are summarized by eq. (16) For any critical system with canonically normalized magnetic Hamiltonian ∆H = d 3 xHσ, and mean order parameter of the formσ = Bh 1/δ f G (z) with h = H/H 0 , the combination of parameters H 0 B/T c has dimensions (length) −d and defines a non-universal length: The longitudinal correlation length of the generic critical system takes the form where f ξ (z) is a universal function, including its normalization 3  (A5) Next we discuss the dynamics. There is a time scale set by the frequency of the pion pole 4 In the O(4) model it is easy that the corresponding frequency in physical units is Comparing the two expressions, using qq = B QCDσ and the identification leads to the result Eqs. (A5) and (A12) are presented in the body of the text in eq. (16).

Overview of the algorithm
In this section we will describe the update algorithm in detail. We first discretize the fields on a spatial lattice (with lattice spacing a = 1) writing the effective Hamiltonian as Here x = (x 1 , x 2 , x 3 ) labels the lattice sites, and e i = e 1 , e 2 , e 3 is a unit vector in the corresponding direction. Variational derivatives in the equations of motion get replaced by ordinary derivatives, δH/δφ → ∂H/∂φ, etc. The equations of motion for U ∈ [φ a , n ab ] can be written schematically 4 This formula assume that the total charge operators Q ab are unitless and satisfy the O(4) commutation The susceptibility is defined by the averages The first operator describes the evolution under the ideal equations of motion, while the second two operators describe the dissipative dynamics of the order parameter φ and the charges n ab respectively. We will use operator splitting to solve for the total time evolution. The most straightforward procedure is to update the fields sequentially for a small period time ∆t U More explicitly we have: A Stage: B Stage: C Stage: We will view each step as part of a Markov Chain. A technical complication is that the C step takes approximately six times longer than the B step, because there are many more random numbers to generate. The ideal step A also is about twice slower than the B step. So as a practical matter, for a complete step over a time ∆t we will take the following updates where the time increment for B is ∆t B = ∆t/6 while the time step for A is ∆t A = ∆t/3. An optimal time step thermalizes modes of order the lattice spacing in a short period of walltime. We have found ∆t = 0.24/Γ 0 is approximately optimal (see below), for the algorithm discussed here.

a. Ideal step
In order to perform our ideal step, let us first rewrite the ideal part of our continuous equation as follows Eq. (A23) makes it apparent that the ideal evolution of the order parameter is simply an O(4) rotation by the currents. More explicitly, in the O(4)-algebra matrix notation, we have with and K, J the generators of so(4) In particular, (A26) can be solved as With this in mind, before describing our time evolution, we need to discretize (A24) in space. For f and g functions evaluated on a discrete lattice, we discretize terms of the sort ∂ x (g∂ x f ) in a straightforward way by integrating over finite volume cells, e.g. in one dimension 1 a where we have approximated the value g| x i+ 1 2 at each interface as the mean of the central one, g| x i+ 1 2 = 1 2 (g| x i+1 + g| x i ). Using the shorthand notation f ±i ≡ f (t, x ± e i ), leads us to define the following discrete evolution kernels σπ s,+i − π s σ +i − σ −i π s + π s,−i σ .
To evolve this system, we use a "position Verlet"-like symplectic integration. We start by computing Φ at half-integer steps, use it to evolve the currents by a time step δt and finish updating Φ by an extra half-time step, which gives A practical way to perform the rotations in (A34) and (A37) is to rewrite the O(4) rotation as a direct product of SU L (2) × SU R (2) and to use the explicit form of the SU (2) matrices for a given set of angles. The ideal evolution is associated to the conservation of the the discretized energy, H, for δt → 0. Our symplectic evolution leads to a violation ∆E = H(t + δt) − H(t) ∼ O(δt 2 ). One approach to this violation would be to just ignore it. Then the equilibrium action will be modified slightly by terms of order O(δt 2 ) from (A13), shifting T c by a small amount. We have seen indications of these shifts but have not explored this in detail. Instead we have added a Metropolis "accept-reject" step to the ideal evolution using min (1, exp (−∆E)) as the accept-reject probability. For our 80 3 lattices (which represent the majority of the simulations presented here) the reject probabilities are presented in Table I. The differences imply that the relative size of the dissipative and real parameters of the pions will weakly depend on δt.
The downside of having an acceptance probability p different from one is that it introduces a non-trivial renormalization of our time. Effectively, when the ideal step is rejected, the next dissipative step should be thought of as a way to generate a new candidate configuration for the ideal step; the clock freezes. The leading effect of a non-zero rejection probability for the ideal step can then be absorbed by rescaling ∆t by the acceptance probability p. As a result, for all our simulations, the time variable we use is defined as with the corresponding p read from Table I and ∆t is the global timestep of our algorithm.
As the results presented through this work support, this procedure allows us to faithfully correct our time variable. It is nonetheless true that it introduces some uncontrolled subleading systematic errors which may impede us from performing precision measurements in the future. This, together with the fact that the acceptance rate degrades for larger lattices, will lead us to use smaller ideal time steps for future simulations. It may also be worth investigating higher order symplectic integrators, which would help to keep the reject probability small even for large volumes.

b. Viscous steps for φ
The spatially discretized equation to be solved is where the noise correlator is given by a discretized eq. (4a). We will realize the Langevin process with Metropolis updates. Briefly, an update proposal is made for a lattice site x where for each flavor index a the increment is Here ξ 0 is a random number with unit variance ξ 2 0 = 1. In practice ξ 0 is generated from a flat distribution between [− √ 3, √ 3], since this is faster than generating Gaussian random numbers. The update proposal is accepted with probability min(1, e −∆H ), where ∆H is the change in the discretized Hamiltonian. If the proposal is rejected φ(t + δt, x) = φ(t, x). For ∆φ small, and then the mean and variance of the accepted proposals reproduce the dissipative and stochastic terms of the Langevin process: For the sake of clarity, let us rederive this result. Consider the Markov process generated by the Metropolis algorithm: ∆φ is accepted if e −∆H − 1 is positive, otherwise it is accepted only with probability e −∆H . Employing the step function θ(x), the update rules for each lattice site can be written as = φ a (t, x) + ∆φ a + θ(1 − e −∆H ) e −∆H − 1 ∆φ a .
In the limit of small δt, ∆φ a is small, and we can Taylor expand the energy and the probability, obtaining Taking averages and noting that the θ-function vanishes for half of the realizations, one immediately reproduces Eqs. (A42a)-(A42b).
To iterate over the sites, we loop over the lattice in a checkerboard pattern, first updating all of the even sites, and then updating the odd sites. Since the interactions are nearest neighbors only, the even site updates are independent of each other and can be done in any order. In addition, the checkerboard Metropolis updates maintain the lattice translational invariance and are easy to implement with PETSc and MPI [33,34].
Finally we turn to the step size δt. We would like the computer time required to thermalize modes of wavelength ∼ a to be as short as possible. If δt is small, then the steps are always accepted, but lead only to a small change in φ; equilibration then requires many steps. If δt is large, then ∆φ is large, but the updates are always rejected, again requiring many steps. We have found that choosing δt = 0.04/Γ 0 leads to an accept-reject probability of approximately 0.5, optimizing these considerations.
c. Viscous steps for charges n A and n V We are considering the evolution equation ∂ t U = O C (U ). Since each charge in the tensor n ab is independent we will dispense with the flavor indices in the rest of this section. All the updates described here will be applied in sequence to the three axial charges n s A and the three vector charges n s V . The continuum equation to be solved is the stochastic diffusion equation and equilibrium effective Hamiltonian is 5 To generate the Langevin dynamics in (A45) we will again use Metropolis steps. In order to get the correct diffusive dynamics at long wavelengths the charge must be exactly conserved by the update proposals. We therefore update the cells in pairs by making a Metropolis proposal for the charge transferred between two cells over a time δt.
The figure below shows a few sites of the lattice, with the even sites painted grey. Integrating (A45) over the spatial volume of lattice cell A and time δt, the discretized equation of motion for the charge takes the form 5 Again this action describes only one isospin component of the iso-axial or iso-vector charge. In general and n ab n ab is written n 2 in the majority of the text.
where, for example, Q x + is the charge transfer between A and B over a time δt. (For clarity below we have restored the hats to indicate quantities in lattice units, e.g.n = na 3 )

A B
The proposed Metropolis flux through the interface is where again ξ 0 is a uniform random number with unit variance. Thus the proposed update for cells A and B isn and change in action by the proposed change is The proposed updated is accepted with probability min(1, exp(−∆Ĥ)). Then it is easy to see that mean charge transfer is which is the expected charge transfer for a diffusive step. Finally, is easy to show that the flux Ξ x q/(δta 2 ) has the expected variance. Thus for small δt the Markov updates produce an equivalent update to the Langevin step.
To iterate over the faces of the lattice we again divide the cells into a checkerboard pattern. We first do the Metropolis updates for all of the x + interfaces for all of the even cells, i.e. cell A is even and cell B is odd as shown in the figure above. These updates are independent of each other and can be done in any order. This step is followed by Metropolis updates of the x − interfaces of the even cells, i.e. now cell A is odd and cell B is even. Then we proceed to update the y and z directions in a similar manner. To eliminate potential bias, the order of the (x, y, z) iterations and the (+, −) iterations are each randomly shuffled for each iteration of the C stage of the Markov chain.  approaches a universal value U c 4 at the critical temperature for large L [41]. From a precision study by Hasenbusch [36], we have taken this asymptotic value to be U c 4 = 1.0945(1). Now, for each L, we determined a value of m 2 0 , denoted m 2 × (L), where the Binder cumulant reaches U c 4 (see figure). For large L the expected scaling of m 2 × (L) − m 2 c is and this scaling can be used to determine m 2 c . Then we fit m 2 × (L), with the functional form in (B4) to find the value of m 2 c quoted in the body of the text in (26). A plot of our fit is shown in the right panel of Fig. 12 and has χ 2 /dof = 0.5, suggesting that the error bars have been slightly overestimated.

Thermodynamics on the critical line
In this section we continue the discussion in Sect. III A, and determine the non-universal parameter H 0 by making a scan on the critical line.
In practice, we made a scan only approximately on the critical line at m 2 = −4.8130, and then used reweighting to determineσ at our nominal value of m 2 c = −4.8110. In the scaling theory the condensate takes the form where in addition to the scaling variables h and z, we have included an additional dependence on the system size L through the scaling variable z L = L 0 /Lh νc , and its associated constant L 0 [6]. The scaling function f G (z, z L ) has been parameterized for z L = [0, 1.2] by Engels and Karsch [6] (see their Eq. 29). We have also included a subleading scaling function, f G (z, z L ), which provides a correction to the leading asymptotics close to T c . We are working on the critical line where z = 0, and the dependence on z L is weak for L large. Thus we will neglect the dependence on z L in the subleading term and describe our data with the form σ = h 1/δ (f G (0, z L ) + h ωνc C H ) . (B6) For the critical exponents here and below we use the results from [6] β = 0.380(2) , δ = 0.4824(9) , and then used the hyperscaling relations to determine all others, e.g. dν = β(1 + δ) 2.213. We have taken ω = 0.77 for the subleading exponent from [36]. The data forσ on 32 3 with χ 2 /dof = 2. H 0 is recorded in Eq. (26). Also shown in the left panel of Fig. 1 is the predicted magnetization from the fit at infinite volume (the dashed line). We see that already at L = 64, we are essentially at infinite volume for the range of H considered in this work. Indeed a simple two parameter fit to our L = 64 results (not shown) with a simple form,σ = (H/H 0 ) 1/δ (1 + C H (H/H 0 ) ωνc ), yields compatible results for H 0 .
To reduce the remaining effects of higher states, we reduce the range of our fit to [x min , L − x min − 1] and study the dependence of the parameters on x min . Their nominal value is then extracted by fitting the resulting plateaus. This procedure is illustrated in Fig. 13 and lead to the following determination leading to f 2 m 2 = (1.053 ± 0.007 (stat.)) · 10 −3 .