Dynamics of the $O(4)$ critical point in QCD: critical pions and diffusion in Model G

We present a detailed study of the finite momentum dynamics of the $O(4)$ critical point of QCD, which lies in the dynamic universality class of Model G. The critical scaling of the model is analyzed in multiple dynamical channels. For instance, the finite momentum analysis allows us to precisely extract the pion dispersion curve below the critical point. The pion velocity is in striking agreement with the predictions relation and static universality. The pion damping rate and velocity are both consistent with the dynamical critical exponent $\zeta = 3/2$ of Model G. Similarly, although the critical amplitude for the diffusion coefficient of the conserved $O(4)$ charges is small, it is clearly visible both in the restored phase and with finite explicit symmetry breaking, and its dynamical scaling is again consistent with $\zeta=3/2$. We determine a new set of universal dynamical critical amplitude ratios relating the diffusion coefficient to a suitably defined order parameter relaxation time. We also show that in a finite volume simulation, the chiral condensate diffuses on the coset manifold in a manner consistent with dynamical scaling, and with a diffusion coefficient that is determined by the transport coefficients of hydrodynamic pions. Finally, the amplitude ratios (together with other non-universal amplitudes also reported here) compile all relevant information for further studies of Model G both in and out of equilibrium.


I. INTRODUCTION
High-multiplicity reactions in colliders are a unique tool to explore the phase diagram of QCD. Their ability to probe different phases and phase transitions crucially relies on a precise understanding of its dynamics in a hot and dense state. Arguably one of the most impactful results of the scientific program of the Relativistic Heavy-Ion Collider (RHIC) was the discovery of a very short thermalization time after which the quark-gluon plasma created in the collision exhibits collective behavior, faithfully reproduced by hydrodynamics. This phenomenon is key to the predictive power of heavy-ion collisions.
This work is a direct continuation of [1] and touches on the nature of the appropriate hydrodynamic theory required by heavy-ion collisions. Away from a phase transition, the hydrodynamics of a system is fully characterized by symmetries and conserved charges. At criticality, the dynamics of the order parameter and its coupling to the conserved charges also needs to be considered [2]. This is crucial to the dynamics of two-flavor QCD in the chiral limit, close to its second order phase transition [3]. In the chiral limit, the hydrodynamic theory above T c is specified by an unbroken SU L (2) × SU R (2) O(4) symmetry group with corresponding conserved charges. Below T c the symmetry is broken to SU V (2) O(3) and the associated massless Goldstone modes (the pions) must be incorporated into the hydrodynamic effective theory [4][5][6]. Finally, in the critical region the appropriate macroscopic description includes an O(4) order parameter (φ a = (σ, π)), and the effective theory interpolates between the two hydrodynamic limits above and below the critical point.
Beyond the theoretical appeal of the chiral limit, the real-world up and down quark masses are small compared to QCD scales, to the point where it is legitimate to ask whether aspects of real-world QCD close to its phase transition can be recovered by deforming the massless limit. For static quantities this question was asked in [7,8] and led to further studies of the chiral critical point with explicit symmetry breaking [9][10][11][12][13][14][15][16][17]. These studies, when combined with important results from Lattice QCD at various quark masses [18][19][20][21], indicate that the static properties of the chiral crossover in real world QCD can be reasonably understood using universality considerations and the existence of a critical point in the massless two-flavor limit.
It remains to be seen whether the dynamics of the chiral crossover can be similarly understood using the proximity of QCD to the O(4) critical point. The relevant hydrodynamics of the phase transition lies in the dynamical universality class of "Model G" in the Halperin-Hohenberg classification scheme [3], which predicts the existence of long-range chiral modes that describe the evolution of soft pions in the crossover region. When the O(4) symmetry is explicitly broken, the pions develop a mass, tempering the impact of these modes on the dynamics on the system. Nevertheless, although these modes contribute negligibly to the energy momentum tensor, they survive over time scales much longer than the typical relaxation time and most definitely affect the long-range dynamics of chiral observables in QCD in the crossover region.
For an experimental perspective, although it is challenging to systematically measure soft pions of transverse momenta smaller than 400 MeV, recent analyses uncovered some discrepancies between state-of-the-art hydrodynamic models and available data in this momentum range. The results are indicative that conventional hydrodynamics does not fully describe the long wavelength dynamics of these pseudo-Goldstone modes [22][23][24].
Taking these indications into consideration, and further encouraged by the proposed upgrade to the ALICE detector [25] which will be more sensitive to low p T pions, a natural question, already raised in [3,6] and subsequently studied further [26,27], is the following: what is the effect of the critical chiral dynamics on the hydrodynamics of real-world QCD? To answer this question a more precise and detailed understanding of the dynamics of "Model G" itself is required.
Over the years, a series of classical-statistical simulations have been presented in the literature [28][29][30][31][32][33][34], building up towards studies of the dynamics of "Model G" and "Model H", the critical model of the conjectured second order phase transition at the finite baryon chemical potential [35]. In a similar spirit, functional Renormalization Group (fRG) studies have been carried out in [36][37][38][39][40][41][42]. Although limitted to mean field, holographic models can provide additional analytical insight into the dynamics of the chiral phase transition [43].
This work expands on [1] (where the first successful numerical simulations of "Model G" were presented) by extending the zero momentum analysis to finite momenta. With this extension, we are able to quantify the dispersion curve of critical pions in detail and study the critical diffusion of the conserved charges. We show that in a finite volume simulation the orientation of the chiral condensate diffuses on the three sphere (the coset manifold), with a diffusion coefficient determined by the transport coefficients of hydrodynamic pions in infinite volume. When chiral symmetry is explicitly broken, the diffusive dynamics is coupled to a Hamiltonian structure determined by the pion mass, and the combined evolution dynamically realizes the -regime of QCD at finite temperature 1 . We also compute new universal dynamical ratios relating the diffusion coefficient to the (appropriately defined) order parameter relaxation time, which we also determine both above T c and on the critical line. In total, these measurements corroborate the dynamical scaling of "Model G" with critical exponent ζ = d/2 across multiple observables, ranging from the pion velocity and damping rate, to the diffusion of the O(4) charges 2 . These results will prove essential to the next step of this series of work, which will quantitatively analyze the non-equilibrium dynamics of the expanding critical system.

II. OVERVIEW OF MODEL G AND ITS PHASES
A. Model G The theory we are considering is a generalization of "Model G" of [2] and was first introduced in [3]. It is the relevant hydrodynamic theory constructed around the O(4) critical point of "2-flavor chiral QCD", namely QCD with massless up and down quarks (m u = m d = 0). It contains a O(4) order parameter φ a ≡ (σ, π), a proxy for the quark condensate of real QCD, dynamically coupled to vector and axial charges 3 , n s V and n s A respectively. They correspond to the original iso-vector and iso-axial vector currents, n V ∼qγ 0 t I q and n A ∼qγ 0 γ 5 t I q, respectively. These can be conveniently combined into an antisymmetric tensor of charge densities n ab ∈ o(4): The free energy (or effective Hamiltonian) in the presence of an explicit symmetry breaking H a = (H, 0) can be studied using an O(4) Landau-Ginzburg action where with m 2 0 negative. In equilibrium, the charges have Gaussian fluctuations set by a susceptibility χ 0 , which is the same coefficient for both the iso-vector and the iso-axial-vector charges near the critical point. The interactions between the charges and the order parameter enters the dynamics through the equations of motion and are determined by the non-trivial Poisson brackets between these fields -see [3,27] for explicit expressions and a derivation. This construction leads to the following evolution equations Here, for example, H [a φ b] denotes the anti-symmetrization, H a φ b − H b φ a . The coefficients Γ 0 and σ 0 are the bare kinetic coefficients associated with the order parameter and the charges. The bare diffusion coefficient of the charges is D 0 = σ 0 /χ 0 . The constant g 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 [2] θ a (t, The dissipative Model G dynamics relaxes to the equilibrium distribution for the fields φ a and n ab [2] which reproduces the static critical behavior of the O(4) model. The basic observable that will be analyzed are statistical two-point function between fields as a function of relative momentum and time. In the presence of the explicit breaking H a = Hδ a0 the flavor index can always be decomposed in a longitudinal σ and transverse π direction, therefore it is natural to define: on the critical line, t r = 0. For H = 0 and in the broken phase the magnetization scales as with amplitude 4 , B − = 0.988 (7). The non-universal parameters B − and H 0 completely fix the properties of the magnetic equation of state close to the critical point. For the static critical exponents quoted here and below we take the values β = 0.380(2) and δ = 4.824(9) from [12]. For the correction-to-scaling exponent we take ω = 0.755(5) [46] and the for the dynamical critical exponent we take ζ = d/2 [3], where d = 3 is the number of spatial dimensions. A summary of all exponents used in this work is given in Table II of  App. B. We performed additional analyses along the critical line using the data generated in [1] for different lattice sizes L. We also generated an extensive set of data at H = 0, both in the broken (m 2 0 < m 2 c ) and restored (m 2 0 > m 2 c ) phases. The dynamical parameters were chosen as in our previous work, χ 0 = 5, Γ 0 = 1 and D 0 = 1/3. Our full dataset is summarized in Table I.
For context, we also display the correlation length (in lattice units) for our different simulations in Fig. 1. The correlation lengths in the restored phase and on the critical line correspond to the fits presented in App. B 1 and B 2. As discussed below, in the broken phase the hydrodynamic theory is characterized by a pion EFT at the longest wavelengths. The pion parameter f 2 (T ) is a matching coefficient in the EFT, which characterizes the pion's fluctuations and scales as the inverse correlation length close to the critical point [5,47]. We will use 1/(2f 2 ) as an estimate for the correlation length in the broken phase which is determined numerically in App. B 3. Here the factor of 1 2 is conventional and is motivated by the fact that finite volume corrections to correlation functions (which are computable using the pion EFT [47]) are organized in powers of 1/(2f 2 L) with order one coefficients, see App. B 3. The black data points correspond to the different temperatures we simulated. Before diving into a quantitative study, we use the rest of this section to present the qualitative behavior of the different channels in the restored and broken phases (with H = 0) and on the critical line (t r = 0 and H = 0).

C. The restored phase with H = 0
In the restored phase, the O(4) symmetry is unbroken and there is no distinction between the axial and vector channels. The order parameter simply dissipates towards equilibrium on a finite (if long) timescale. But, both the iso-vector and iso-axial vector charges are exactly conserved. As a result, the non-trivial dynamics at the longest wavelengths is fully encoded in the correlators of conserved charges at finite momenta. Indeed, for wavelengths long compared to the (diverging) correlation length kξ(T ) 1, the diffusion equation remains a valid hydrodynamic description of the system: Thus, the relaxation rate of a diffusive mode of wavenumber k is controlled by the diffusion coefficient Dk 2 leading to the correlators (see App. A 1):  The charge correlations in the simulation are illustrated on the left-hand side of Fig.2, where we plot the axial correlator (orange) and the vector correlator (blue), for the first two momenta, for a given temperature. We clearly see that the two channels are degenerate. The k dependence is consistent with the expected diffusive behavior.
In writing the diffusion equation in eq. (17) we have integrated out the contributions to the currents from the order parameter, which then leads the scaling behavior of the transport coefficient. As discussed below, the order parameter has a dynamical relaxation time of order τ ∝ ξ ζ with ζ = d/2. The diffusive modes have a relaxation time of 1/Dk 2 , which is of order ξ 2 /D at the boundary of applicability of the diffusion equation, k ∼ ξ −1 . In the spirit universality, all modes of with wavelengths of order ξ should have a similar relaxation time, dictating the expected scaling of the diffusion coefficient, D ∝ ξ 2−ζ .
On the right-hand side of Fig. 2, we show the dependence of the vector correlator (we could have equivalently chosen the axial) on t r . This very weak dependence seems at first to be in contradiction with the expected critical scaling of the diffusion coefficient D ∼ t (ζ−2)ν r . As we will see more clearly in the next section, this apparent contradiction originates in the presence of a large regular constant contribution to D, and by the smallness of its nonuniversal critical amplitude. Nevertheless, we will extract the critical contribution to D in Sect. III A.

D. The critical line
We next study the critical line, m 2 0 = m 2 c , and study the effect of explicit symmetry breaking, H = 0. One of the qualitative outcomes of our previous work [1] was to demonstrate the emergence of damped pion waves already on the critical line. This qualitative observation  The symmetry is broken and the vector channel is no longer degenerate with the axial channel. As in the restored phase (see Fig. 2) we see a weak dependence on the temperature, indicating of a small critical amplitude which is extracted in Fig. 6.
was obtained from the k = 0 modes of the axial current and the pions correlators. The same qualitative observation holds for non-zero momentum k and is observed in Fig. 3. At the level of the charge correlators, the magnetic field lifts the degeneracy between the axial and vector channels. The axial correlator then displays characteristic oscillations associated with axial charge propagation in the form of damped pion waves. The vector channel on the other hand remains purely diffusive. As in the restored phase, we already observe that the critical behavior of the vector diffusion constant is weak.

E. Broken phase with H = 0
Below the critical point the longest wavelength modes are characterized by massless Goldstone bosons, which describe phase fluctuations of the chiral condensate 5 . Following [5,6], these fluctuations are parametrized by an angle ϕ(t, x) with φ a = (σ,σ ϕ(t, x)). We are assuming that the chiral condensate is nearly aligned with the pole of the three sphere (the coset manifold) and takes the form where n a is a unit vector that is approximately parametrized by n a (1, ϕ 0 ) to first order in the angular deviation. When the wavelength of ϕ is long compared to the correlation length ξ ≡ m −1 σ , a hydrodynamic description of ϕ is valid, and the fluctuations of modes with wavelength order m −1 σ simply determine the critical behavior of the parameters of the hydrodynamic theory.
The free energy associated with the phase and charge fluctuations in the hydrodynamic approximation is where we have included a mass term, which is relevant only when the magnetic field is non-zero and is determined by the Gell-Mann, Oakes, Renner relation We will set H = m = 0 for the remainder of this section. The resulting hydrodynamic equations of motion for the pions and axial charge take the form 6 which reflects from the Poisson bracket between the phase and the charge, {ϕ s , n As } = δ ss . The vector charge remains diffusive The hydrodynamic system is solved in App. A 1 and the resulting pion eigen-waves have dispersion relations where The form of the correlation function from the hydrodynamic theory is determined in App. A 1 At large wavelengths the pion effective theory determines the static and dynamic correlation functions. The static correlation function for k m σ can be read off from the free energy in (20) and is singular in the massless limit: By contrast, the static fluctuations of the σ are bounded (if diverging) and determined by the static susceptibility of the O(4) critical point, which scales as The σ and π are part of an O(4) multiplet and their correlators must be the same order of magnitude at the boundary of applicability of the pion EFT, k ∼ m σ . Since the vev scales . We have anticipated this scaling in Fig. 3 where T /2f 2 (with T = 1 in our adopted units) served as a definition of the correlation length in the broken phase.
Because of these scalings, the long-wavelength behavior of the static and dynamic correlation functions of φ a are dominated by the phase fluctuations. Specifically, the static correlation function reads up to corrections of order (k/m σ ) 2 . Additional finite volume corrections to this formula are given in App. B 3 and are suppressed by 1/f 2 L [47]. Phase fluctuations also determine the dynamic correlations and Fig. 4 (Left) exhibits the φφ correlation function, which shows characteristic pion oscillations of (25) become increasingly damped near the critical point. The critical scaling of v and D A are extracted numerically in the next section and are ultimately summarized in Fig. 7.
In a finite volume simulation the orientation of the chiral condensate is not fixed, but diffuses in time due to thermal noise. Indeed, a stochastic variable ξ ϕ with variance should be added to (22a). Integrating this equation over the volume shows that the angular zero mode ϕ 0 (t) ≡ x ϕ(t, x)/V has a variance which increases in time This result (which is presented in detail in App. A 2) relates the vev diffusion coefficient to the pion kinetic coefficients, D = T Γ/f 2 V . The diffusion rate decreases with increasing volume and is suppressed relative to the pion damping rate, Γk 2 ∼ Γ/L 2 , by one power of length, T /f 2 L ∝ 1/m σ L. When the magnetic field is non-zero the Fokker-Planck equation describing the motion of the chiral condensate on the three sphere features a rich interplay between Hamiltonian and diffusive dynamics, which is developed in App. A 2.
Turning to the conserved charges, the correlation function of the vector and axial-vector charges in the broken phase are intrinsically different. The axial charge is coupled through in equations of motion to the pions and its correlation function reflects this coupling (see App. A 1) The vector charge continues to obey the diffusion equation with corresponding response functions Separating the vector and axial vector response in the absence of explicit symmetry breaking is challenging, since the condensate forms in an arbitrary direction. Moreover, because of finite volume, this direction is not completely fixed, but slowly wanders along the coset manifold. This can be addressed in future work. For now we have simply computed the combined correlator G nn (t, k) which is exhibited in Fig. 4 (Right). The result curves show a mix of exponential decay and oscillating pion waves. It is striking how the hydrodynamic EFT and the pattern of chiral symmetry break essentially dictates the scaling dynamical response of "Model G" [5,6]. A summary of the reasoning proves an overview of the next section.
(i) The characteristic frequency of the pion waves below the critical point is of order The scaling of the velocity will be extracted from the real time correlations of φ, eq. (25), and from its static behavior, eq. (28), and the results are summarized in Fig. 7 (Right).
(ii) At the boundary of applicability of the pion EFT k ∼ m σ , the pion frequency is we see the angular variance given in (31) grows as At the critical point where m σ L → 1, the angular variance increases as which is again consistent with a universal dynamical critical timescale of τ ∝ L ζ with a common critical exponent of ζ = d/2.  4. Left: Normalized G φφ correlator at k = 2π/128 as function of time for different reduced temperatures t r in the broken phase. We see that this correlator is sensitive to the Goldstone modes. Right: Normalized G nn correlator at k = 2π/128 as function of time for different reduced temperatures t r in the broken phase. This correlator is also sensitive to the Goldstone dynamics.

III. RELAXATION TIME AND FINITE MOMENTA DYNAMICS
A. H = 0, restored phase: scaling of the correlation time We first determine the order parameter relaxation time τ from the two-point function G φφ (t, k). Specifically, we define τ as the autocorrelation time of G φφ at zero momentum which is motivated by an exponential decay, G φφ ∝ e −t/τ (in practice, we cutoff the integral at some T max and make sure the results are independent of this technicality). Qualitatively, the decorrelation of G φφ (t) is characterized by a single timescale τ without significant structure. From scaling, we expect the correlation time to behave as with the critical exponents ν, ω, and ζ given in Table II in App. B. The resulting fit (with the exponents fixed) is shown in Fig. 5 and we find The growth of τ near the T c is nicely consistent with ζ = d/2. The diffusion coefficient D of the O(4) charges is also expected to show dynamical critical scaling, but since the charges are conserved, the coefficient has to be extracted from finite momentum correlators. In particular, we expect the vector correlator in the restored phase to behave as as discussed in the previous section. Similar to the order parameter case, we have the relation which we use in practice to extract D from our data. It is important that the correlator be well described by an exponential, which is clearly seen in Fig. 2. Again, we introduce a cutoff T max and make sure the results are independent of this choice. As discussed above, we expect Dk 2 to be comparable 1/τ ∼ ξ ζ at the boundary of applicability of the diffusive description, k ξ. This means that the diffusion coefficient should scale with the correlation length as, D ∝ ξ 2−ζ .
The critical behavior of D appears to be weak and its regular part cannot be neglected. To take this into account, we add a constant D 0 to our fit In an abuse of notation, D 0 here is the renormalized diffusion coefficient far above T c , which we determine independently by running our simulations without the coupling to the order parameter, leading to D 0 ≈ 0.2425 (5) .
D + is small as was anticipated based on a one loop calculation using mean field propagators [27].
Using the non-universal amplitude ξ + from the correlation length ξ ∼ ξ + t −ν r determined in App. B 1, we can combine these results into a universal dynamical amplitude ratio which numerically encodes the coupling between the charge diffusion and the order parameter relaxation at criticality. The universal dynamical ratio presented here is new and will prove useful to any future study which realizes the critical dynamics of Model G.

B. The critical line
We next study the dynamics on the critical line, m 2 0 = m 2 c and H = 0. As in the previous section we define the order parameter relaxation time τ from the correlation of the σ field The results for τ as a function of the magnetic field are shown in Fig. 6, and are fit with the expected scaling form with fixed exponents ν c , ω, and ζ given in Table II of App. B, and fitted parameters As in the restored phase, the subleading term with exponent ω is necessary to have a reasonable fit compatible with the expected dynamical exponent, ζ = d/2. Similarly, the diffusion coefficient is defined in the same way as in the restored phase and we adopt an analogous strategy and functional form to fit the dependence of D on the correlation length As in the previous section, D 0 = 0.2425 is the regular part of the diffusion coefficient and is extracted from our diffusion only runs. The fit results are shown in Fig. 6 and the fit parameters are The singular part of the diffusion coefficient diverges next to the critical point with the expected exponent H −νc(2−ζ) and the subleading correction improves the quality of the fit. Including the regular part of the diffusion coefficient D 0 is crucial to extracting the scaling of D with H, since the critical amplitude D H is small, as anticipated in the one loop computation in [27]. The results for the order parameter and the diffusion coefficient can be concisely recast as a universal amplitude ratio where ξ H L is the amplitude associated with the longitudinal correlation length and is determined in App. B 2. Note that within numerical errors the dynamical ratio on the critical line matches the one in the restored phase (47).

C. H = 0 and the broken phase
Our aim in this section is to extract the matching coefficients of the pion hydrodynamic EFT by studying the momentum dependence of the correlators and analyzing the pion dispersion curve. To this end, we study G φφ for various t r and L and analyze the two lowest non-zero momentum modes, k = 2π/L and k = 4π/L. We then perform some fits using the hydro prediction (25) and extract the parameters ω(k) and Γ(k). ω(k) and Γ(k) are subsequently extrapolated to infinite volume using the forms with D A , v, d 1 , and v 2 as parameters. As discussed in Sect. II E, finite volume corrections to the dynamics of the Goldstone modes scale inversely with the linear extent of the lattice ∼ 1/f 2 L and are therefore important to control. Including d 1 and v 1 as fit parameters captures this characteristic volume dependence and proved crucial to reliably extracting D A and v from this analysis. The fitting procedure, together with other systematic checks, are presented in detail in App. C. In the static case, a finite volume pion EFT was developed many years ago by Hasenfratz and Leutwyler and can be used to determine the first 1/f 2 L correction to the chiral condensate analytically [47]. We used their expansion in App. B 3 to determine velocity v 2 = f 2 /χ 0 and, as we show below, the result agrees with the dynamic fit using (57). In Sect. II E we have taken the first steps in developing a real time finite volume analog of their work by deriving the vev diffusion coefficient in (31) and a corresponding forced Fokker-Planck equation on the coset manifold in App. A 2. But, a more complete analysis is left for future work.
Knowing v and D A , we can proceed further and study their scaling with temperature. We start with a study of the critical behavior of D A , shown on the left hand side of Fig. 7. We use the following scaling ansatz and obtain the following fit parameters with a χ 2 /dof of 6.3/3, which is mostly driven by the large fluctuation around t r ∼ −0.004. The exponents ν, ζ and the subleading ω were held fixed. As in all the previous cases, including a subleading correction and a regular part is necessary to obtain a precise fit. We did not find a way to independently fix the regular part for this analysis. This leads to a degeneracy between D A0 and D − A1 , resulting in large errors for both of these parameters. This problem does not affect the extracted value of the leading amplitude, D − A . The Goldstone velocity v is shown on the right hand side of Fig. 7. The data are fitted with the expected scaling form, with ν and ω fixed, yielding v − = 0.462 ± 0.003 , The χ 2 is 15.9 with 4 degrees of freedom. The relatively large value of the χ 2 is due to the very high precision of our data. It is driven by the points farther away from the critical point and probably reflects the fact that our data are sensitive to sub-subleading corrections. We also show the determination of v = f 2 /χ 0 obtained from the static analysis in App. B 3. The analysis makes use of (28) but exploits the known finite volume analysis of Hasenfratz and Leutwyler [47]. The agreement between the static analysis and the fits to the real time correlation functions is remarkable and is a highly non-trivial confirmation of the low-energy EFT.
Ideally, at this point we would complete this study by determining the critical amplitude of the diffusion coefficient in the broken phase. However, because the amplitude is small, and because finite volume corrections are only suppressed by 1/f 2 L, we were unable to unambiguously extract this amplitude with the current dataset and our current understanding of finite volume corrections. We can clearly see the diffusive behavior, but the coefficient is constant within uncertainty.

IV. DISCUSSION
In this work, we performed an extensive study of the finite momentum dynamics of "Model G", the critical dynamical model corresponding to two-flavor QCD in the chiral limit. Across all different phases of the theory, we exhibited striking agreement between theoretical predictions and our full-fledged numerical simulations. In particular, although the critical amplitude is small, we showed that the critical contribution to the vector diffusion coefficient D scales with the correct dynamical critical exponent, both in the restored phase and along the critical line. After additionally extracting the critical relaxation time τ of the order parameter (which also scales appropriately), we recast these new results as universal dynamical amplitude ratios Q + D , Q H D , which reflect the coupling between the order parameter and the diffusive dynamics. Next, we took on the challenge of studying the finite momentum dynamics of the broken phase at zero magnetic field. Because the simulation is at finite volume, the chiral condensate is not constant but wanders through the coset manifold at a finite rate. We show that the diffusion coefficient for this process is consistent with dynamical scaling and is determined by the transport coefficients of a pion hydrodynamic EFT, which we review in Sect. II E. Despite this subtlety, we were able to clearly observe the critical pion dynamics. The culmination of this part is presented in Fig. 7, which compares the pion dispersion relation measured in this work to the scaling predictions for the pion damping rate and velocity close to the critical point. The agreement is striking and further exemplifies the precision of our numerics.
Future directions to be explored are diverse. One avenue is to refine the treatment of finite volume effects in the broken phase by developing a finite volume real time EFT for the dynamics of the chiral condensate in the spirit of [47] and [55]. The first steps in this direction are presented in App. A 2, which worked out a Fokker-Planck description for the dissipative and Hamiltonian dynamics of the condensate on the coset manifold. This provides a real time description of QCD in the -regime, suggesting connections with Random Matrix Theory [45].
A more physically motivated direction is the study of the Kibble-Zurek dynamics of this model. Physical systems are typically not tuned at their physical point, but traverse the phase transition, with their relevant coupling varying at a finite rate. This finite rate competes with the critical relaxation rate of the system and prevents the correlation length from diverging. When transiting from the restored to the broken phase, this competition results in the formation of finite-length domains with different values of the condensate. A careful examination of this phenomenon, both in the presence and absence of explicit symmetry breaking, will be crucial to any future phenomenological predictions. We plan to utilize the precision measurements presented here and in our prior work [1] to study this phenomenon in detail.

Appendix A: Hydrodynamics
This appendix extends the discussion of the hydrodynamics of the broken phase given in Sect. II E. We will keep the temperature explicit in this section, following the practice of II E.

Correlation functions from the pion EFT
Here we will review the temporal correlation functions which arise the from pion hydro equations of motion, (22), but include the explicit symmetry breaking term and noise. Although these correlators are limits of the mean field results derived in [27], they transcend mean field and follow directly from the hydrodynamic equations of motion.
The hydrodynamic equations (see Sect. II E) for the phase and axial charge read which follows from the Poisson bracket between the phase and the charge, {ϕ s , n As } = δ ss . The vector charge remains diffusive The noise added to the right-hand-side of Eq. (A1) has variances given by the Fluctuation-Dissipation Theorem (FDT) where σ = χ 0 D. The equations are easily solved. For instance, we can Fourier transform in space and write the solution to the vector equation by integrating from t = 0 up to t (with t > 0) The variance in the initial condition n(0, k) is determined from the free energy functional in (20) and reads Alternatively one can integrate from past infinity to determine the initial condition: Computing the variance of n(0, k) using the variance of the noise in (A2c) reproduces (A4), which demonstrates how the noise and dissipation work together in establishing the equilibrium state. By either method An identical strategy can be used in the axial channel by finding the eigen functions of the system. As in the vector case, we Fourier transform in space and calculate response function. It is convenient to work with the vector of variables 7 Y = (y 1 , y 2 ) = (ω k ϕ, µ A ) where µ A = n A /χ and ω 2 k = v 2 (k + m 2 ). The static susceptibility is diagonal in these variables, and can be read off from the free energy written in (20) As in the diffusion example, the static susceptibility sets the initial conditions for the subsequent evolution. The equations of motion take the form ∂ t Y + MY = ξ and eigen-values of M are The homogeneous solutions are Y ∝ e −λ ± t . The full matrix of correlation functions with ∆ k ≡ (Dk 2 − Γ(k 2 + m 2 ))/2. In the body of the text this result is used in the massless limit in (25) and (32).

Diffusion of the chiral condensate on the coset manifold
a. Free diffusion on the three sphere In a finite volume simulation the chiral condensate is not fixed but diffuses on the coset manifold. Let us examine this diffusion by analyzing the zero mode of the hydrodynamic equations of motion given in (A1), but limiting the discussion to H = 0. In finite volume the fields take the form The axial charge is exactly conserved and there is no net charge in the box, and consequently the k = 0 mode is zero for n A . Then, integrating (22a) over volume, the zero mode of ϕ satisfies a random walk Solving for ϕ 0 and computing its variance leads to (31), which is reproduced here for convenience i.e. D ≡ T Γ/f 2 V . The noise satisfies The probability density P( ϕ 0 ) satisfies the diffusion equation The direction of the chiral condensate is labeled by a unit vector on the three sphere n a n a = 1. We have parametrized this direction by three angles, ϕ 0s , assuming they are small. The coordinates ϕ provide a local Cartesian coordinate patch for the three sphere close to the pole. In a general coordinate system the diffusion equation reads where ∇ I ∇ I is the Laplacian on the sphere. More explicitly, a set of coordinates covering the three sphere is q I = (q 1 , q 2 , q 3 ) = (ϕ S , θ, Φ) where ϕ S is approximately ϕ S ϕ 2 close to the north pole.
b. Explicit symmetry breaking, forced diffusion, and the -regime of QCD When the symmetry is explicitly broken by a magnetic field, the axial charge must be reconsidered. In this case the zero modes satisfy the equations of motion where the statistics of ξ ϕ 0 are given in (A15). The variables ϕ 0 and n A0 are canonical conjugates {ϕ 0s , n A0s } = δ ss . To illustrate the structure of these equations we introduce a traditional notation (Q s , P s ) ≡ (ϕ 0s , n A0s ) , and the equations of motion are written Here β and V are the inverse temperature and volume respectively, and D = T Γ/f 2 V is the diffusion coefficient introduce earlier. Finally the Hamiltonian for the zero modes is The Fokker-Planck equation for the phase space probability density P(Q, P ) follows from the Langevin dynamics and reads [56][57][58] We emphasize that the canonical measure for the probability density P(Q, P ) is invariant under canonical change of coordinates.
We have worked with a set of spatial coordinates, which parametrize a region close to the pole of the three sphere by Cartesian coordinates. We then make a canonical change of variables to the polar coordinates q I given in (A18). The Hamiltonian is generalized to include larger angles and reads where g IJ (q) is the metric of the sphere in (A18). Here we used the original free energy in (3) and the Gell-Mann, Oakes, Renner relation to deduce the potential for the zero mode The change of coordinates leads to the Fokker-Planck equation in its final form where ∇ I is the covariant derivative on the three sphere.
It is straightforward to see that the steady state solution to this equation is the equilibrium probability distribution If the momenta are of no interest they can be integrated over yielding which is consistent with the partition functions discussed in [47]. It is satisfying that the dissipative dynamics of the finite volume chiral condensate on the three sphere is controlled by the (infinite volume) pion pole mass and damping coefficient. Indeed, the Fokker-Planck equation provides a detailed real time picture of the so called "epsilon" regime of QCD. In this regime the quark mass is very small, but βV f 2 m 2 ∼ 1 is of order unity. The thermodynamics of this regime has been extensively studied using Random Matrix Theory and chiral perturbation theory (see [45] for an overview). It would be interesting to pursue this connection further.  [12] δσ ∝ H 1/δ 4.824(9) [12]   The bottom of the table shows the correction-to-scaling exponent ω, which is also universal. In the table we denote "any scaling quantity" by Y , and θ(Y ) and c(Y ) are its associated critical exponent and non-universal subleading amplitude respectively. ζ is the dynamical scaling exponent of "Model G" where d = 3 is the number of spatial dimensions.

Appendix B: Static measurements in the O(4) model
In this appendix we determine the correlation length and the magnetic susceptibility in the restored phase and along the critical line. In particular, their non-universal critical amplitudes serve to define universal ratios, including the novel dynamical ones introduced in the main text. We also determine the pion constant f 2 using static measurements. For convenience we collect the critical exponents used in this work in Table II.

Restored phase
The magnetic susceptibility in the restored phase is defined as and the zero momentum limit is notated with χ χ ≡ χ(0) .
A standard estimator of the correlation length is the so-called "second-moment correlation length" [46], computed as with F the susceptibility of the first Fourier mode We fit their dependence on the reduced temperature t r = ( The exponent ω is the first subleading scaling exponent. To perform the fit, we fix β, ν and ω to the values presented in Tab. II and fit ξ + , ξ 1 , C + and C 1 . The fits are presented in Fig. 8. We obtain ξ + = 0.4492 ± 0.0036, ξ 1 = 0.332 ± 0.085, C + = 0.2077 ± 0.0032, C 1 = 1.11 ± 0.15 . (B7) The chi-square for the fit of the correlation length is χ 2 /dof = 9.1/7 and the one for the susceptibility is χ 2 /dof = 4.83/7. These values, together with the value the non-universal value for the magnetization in the restored phasē This value is indeed compatible with the Q c = 0.4231(9) of [59]. 8 Note an unfortunate misprint. Equation (29)  At a non-zero magnetic field the magnetization will be parallel to the magnetic field. Then it is natural to define and measure the longitudinal and transverse χ L susceptibilities: The static susceptibilities are shown in Fig. 9. As a fitting function, we chose the scaling prediction with the first subleading correction To extract the correlation length in the longitudinal and the transverse directions we measure the second-moment correlation length as in (B3) with The scaling prediction fixes the fitting function as and the result of the fit is shown in Fig. 10. The universal ratio ξ T + /ξ L+ = 1.85 ± 0.12 is consistent with the more precise estimate done in [59].

Broken phase and the pion velocity
Our goal in this appendix is to extract the infinite volume chiral condensateσ(T ) and pion velocity v 2 (T ) = f 2 (T )/χ 0 by analyzing static correlators below the critical point at H = 0. The resulting parameterization of the velocity is shown in Fig. 7 (Right). We will use the magnetization the static correlation function G φφ (0, k), and their their dependence on volume to extract σ(T ) and (σ 2 /f 2 )(T ) in the infinite volume limit.
In infinite volume the magnetization remains constant in time and determines the chiral condensate lim Here n a is an arbitrary unit vector on the three-sphere characterizing the orientation of the condensate. However, in any finite volume since the condensate orientation vector n a stochastically explores the S 3 in time -see App. A 2.
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 volumes. The leading deviation of M 2 fromσ 2 at finite volume comes from the fluctuations of long wavelength Goldstone modes and can be neatly analyzed with a Euclidean pion effective theory [47,60] In the chiral limit the only parameter of the EFT is the pion constant f 2 (T ). When a finite magnetic field is included the chiral condensateσ is also a parameter.
Apart from a slightly different fitting procedure, the current determination ofσ(T ) simply repeats the analysis done in our earlier work for the larger lattices and datasets used here [1]. The expansion relating M 2 andσ 2 takes the form [1,47] and is valid for f 2 (T )L 1. The numerical coefficients are known analytically in terms of "shape coefficients" β 1 and β 2 , but are of no interest here.
To determine the chiral condensate, we plotted M 2 versus L at fixed temperature, and made various fits with (B21) to determineσ 2 (T ) and (σ 2 /f 2 )(T ). We estimated the systematic uncertainty in the extracted values ofσ(T ) by adding a C/L 3 term to the fit form and comparing to a fit based on (B21). Except for the point closest to T c the systematic uncertaintyσ(T ) is smaller than our statistical uncertainty, which would not have been the case if only the first term ∝ 1/f 2 L in the expansion had been used.
Our results forσ(T ) are shown in the right panel of Fig. 11, This is fit with the functional formσ (t r ) = B − (−t r ) β (1 + B 1 (−t r ) ων ) .
with fit parameters B − = 0.988 ± 0.007, B 1 = 0.51 ± 0.07. The results for the vev amplitude and subleading corrections are compatible with our previous results. 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 decay constant f 2 of the pion EFT can be extracted using similar methods The correlation function of the scalar field at small wavenumbers is dominated by the massless Goldstone mode. Recalling eq. (28) with T = 1, the correlation function takes the form to leading order in the pion EFT. However, there are important corrections to this result stemming from the finite system size, which are of order 1/f 2 L and are captured by the leading order pion EFT. Short distance correlations of order (m σ L) 2 and smaller lead to higher derivative terms in the pion EFT, which are not included. Naturally, these terms become increasingly important near the critical point. Following the perturbative expansion of Hasenfratz and Leutwyler the correction takes the form where β 1 = 0.225785 is a "shape coefficient" reflecting the fluctuations of pions in a cubic box. Numerically for k = 2π/L and N = 4: To extractσ 2 /f 2 we plotted k 2 G φφ (0, k) for the first Fourier mode as a function of system size and fitted the result for (σ 2 /f 2 )(T ) and 1/f 2 using (B26). The subleading correction ∼ 1/f 2 L in (B26) was treated as a free parameter, although it is broadly consistent with theσ 2 (T ) extracted above and the (σ 2 /f 2 )(T ) determined here. The data on (σ 2 /f 2 )(T ) are shown in Fig. 11. Asymptotically,σ 2 /f 2 should behave as where η = 0.030 is a critical exponent. We are unable to extract η from our measurements. Instead, we have simply fit a formσ with C 0 = 0.910(8) and C 1 = 0.009(3), which gives a heuristic parametrization of the results shown in the figure. There is a hint in the data that η is non-zero. Given the fit results forσ 2 (T ) in (B22) and (B28) for (σ 2 /f 2 )(T ) we can extract the velocity v 2 = f 2 χ 0 = 1 χ 0σ 2 (T ) (σ 2 /f 2 )(T ) .
The uncertainty is dominated by the contribution fromσ 2 (T ) asσ 2 /f 2 (which is nearly unity) in total provides a ten percent correction. In Fig. 11 we have varied the fit parameters in the range allowed by the fits to deduce the static velocity and its uncertainty.