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

Adrien Florio<EMAIL_ADDRESS>Department of Physics, Brookhaven National
Laboratory, Upton, New York 11973-5000, USA  Eduardo Grossi
<EMAIL_ADDRESS>Dipartimento di Fisica, Università di Firenze and INFN
Sezione di Firenze, via G. Sansone 1, 50019 Sesto Fiorentino, Italy  Derek
Teaney<EMAIL_ADDRESS>Center for Nuclear Theory, Department of
Physics and Astronomy, Stony Brook University, New York 11794-3800, USA

###### Abstract

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)\times SU_{R}(2)\simeq O(4)$ symmetry group with corresponding
conserved charges. Below $T_{c}$ the symmetry is broken to $SU_{V}(2)\simeq
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 ($\phi_{a}=(\sigma,\vec{\pi})$), 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\,{\rm 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
$\epsilon$-regime of QCD at finite temperature111 The $\epsilon$ regime is
when the quark mass $m_{q}\propto H$ is small and the volume $V$ is large, but
$H\bar{\sigma}(T)V/T\sim 1$. Here $\bar{\sigma}(T)$ is the chiral condensate
as a function of temperature – see [45] for an overview. . 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 $\zeta=d/2$ across multiple observables, ranging from the pion
velocity and damping rate, to the diffusion of the $O(4)$ charges222Here $d=3$
is the dimensionality of space.. 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

###  II.1 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 $\phi_{a}\equiv(\sigma,\vec{\pi})$, a proxy for the quark condensate
of real QCD, dynamically coupled to vector and axial charges333 Here $a$ and
$b$ denote $O(4)$ indices; $s$, $s_{1}$, $s_{2}$, etc. denote the three
isospin indices, i.e. the components of $\vec{\pi}=(\pi_{1},\pi_{2},\pi_{3})$;
finally, spatial indices are notated $i$, $j$ and $k$. The dot product
indicates an appropriate contraction of indices when clear from context, e.g.
$\phi\cdot\phi=\phi_{a}\phi_{a}$, $\vec{\pi}\cdot\vec{\pi}=\pi_{s}\pi_{s}$,
and $\nabla\cdot\nabla=\partial_{i}\partial^{i}$. , $n_{V}^{s}$ and
$n_{A}^{s}$ respectively. They correspond to the original iso-vector and iso-
axial vector currents, $\vec{n}_{V}\sim\bar{q}\gamma^{0}\vec{t}_{I}q$ and
$\vec{n}_{A}\sim\bar{q}\gamma^{0}\gamma^{5}\vec{t}_{I}q$, respectively. These
can be conveniently combined into an antisymmetric tensor of charge densities
$n_{ab}\in\mathfrak{o}(4)$:

$\displaystyle(n_{V})_{s}$
$\displaystyle=\frac{1}{2}\epsilon_{0ss_{1}s_{2}}n_{s_{1}s_{2}}\,,$ (1)
$\displaystyle(n_{A})_{s}$ $\displaystyle=n_{0s}\,.$ (2)

The free energy (or effective Hamiltonian) in the presence of an explicit
symmetry breaking $H_{a}=(H,\vec{0})$ can be studied using an $O(4)$ Landau-
Ginzburg action

$\mathcal{H}\equiv\int
d^{3}x\left[\frac{n^{2}}{4\chi_{0}}+\frac{1}{2}\partial_{i}\phi_{a}\partial^{i}\phi_{a}+V(\phi)-H\cdot\phi\right]\,,$
(3)

where

$V(\phi)=\frac{1}{2}m_{0}^{2}\,\phi^{2}+\frac{\lambda}{4}(\phi\cdot\phi)^{2}\,,$
(4)

with $m_{0}^{2}$ negative. In equilibrium, the charges have Gaussian
fluctuations set by a susceptibility $\chi_{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

$\displaystyle\partial_{t}\phi_{a}+\frac{g}{\chi_{0}}\,n_{ab}\phi_{b}$
$\displaystyle=-\Gamma_{0}\frac{\delta\mathcal{H}}{\delta\phi_{a}}+\theta_{a}\,,$
(5a)
$\displaystyle=\Gamma_{0}\nabla^{2}\phi_{a}-\Gamma_{0}(m_{0}^{2}+\lambda\phi^{2})\phi_{a}+\Gamma_{0}H_{a}+\theta_{a}\,,$
(5b)
$\displaystyle\partial_{t}n_{ab}+g\,\nabla\cdot(\nabla\phi_{[a}\phi_{b]})+H_{[a}\phi_{b]}$
$\displaystyle=\sigma_{0}\nabla^{2}\frac{\delta\mathcal{H}}{\delta
n_{ab}}+\partial_{i}\Xi_{ab}^{i}\,,$ (5c)
$\displaystyle=D_{0}\nabla^{2}n_{ab}+\partial_{i}\Xi_{ab}^{i}\,.$ (5d)

Here, for example, $H_{[a}\phi_{b]}$ denotes the anti-symmetrization,
$H_{a}\phi_{b}-H_{b}\phi_{a}$. The coefficients $\Gamma_{0}$ and $\sigma_{0}$
are the bare kinetic coefficients associated with the order parameter and the
charges. The bare diffusion coefficient of the charges is
$D_{0}=\sigma_{0}/\chi_{0}$. The constant $g$ is a coupling of the field
$\phi$, and has the units of $({\rm action})^{-1}$ in our conventions.
Finally, $\theta_{a}$ and $\Xi_{ab}$ are the appropriate noises, which are
defined through their two-point correlations [2]

$\displaystyle\langle\theta_{a}(t,x)\theta_{b}(t^{\prime},x^{\prime})\rangle$
$\displaystyle=2T_{c}\Gamma_{0}\,\delta_{ab}\,\delta(t-t^{\prime})\delta^{3}(x-x^{\prime})\,,$
(6a)
$\displaystyle\langle\Xi^{i}_{ab}(t,x)\Xi^{j}_{cd}(t^{\prime},x^{\prime})\rangle$
$\displaystyle=2T_{c}\sigma_{0}\,\delta^{ij}\left(\delta_{ac}\delta_{bd}-\delta_{ad}\delta_{bc}\right)\,\delta(t-t^{\prime})\delta^{3}(x-x^{\prime})\
.$ (6b)

The dissipative Model G dynamics relaxes to the equilibrium distribution for
the fields $\phi_{a}$ and $n_{ab}$ [2]

$P[\phi,n]=\frac{1}{Z}e^{-\mathcal{H}[\phi,n]/T_{c}}\,,$ (7)

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\delta_{a0}$ the flavor index can always be
decomposed in a longitudinal $\sigma$ and transverse $\vec{\pi}$ direction,
therefore it is natural to define:

$\displaystyle
G_{\sigma\sigma}(t,{\bm{k}})=\frac{1}{V}\langle\phi_{0}(t,{\bm{k}})\phi_{0}(0,-{\bm{k}})\rangle_{c}\,,$
(8)

for the scalar channel, and

$\displaystyle
G_{\pi\pi}(t,{\bm{k}})=\frac{1}{3V}\sum_{s=1}^{3}\langle\pi_{s}(t,{\bm{k}})\pi_{s}(0,-{\bm{k}})\rangle_{c}\,,$
(9)

for the pseudo-scalar channel. The statistical correlators for the charges are
defined analogously

$\displaystyle G_{AA}(t,{\bm{k}})$
$\displaystyle=\frac{1}{3V}\sum_{s=1}^{3}\langle
n_{As}(t,{\bm{k}})n_{As}(0,-{\bm{k}})\rangle_{c}\,,$ (10) $\displaystyle
G_{VV}(t,{\bm{k}})$ $\displaystyle=\frac{1}{3V}\sum_{s=1}^{3}\langle
n_{Vs}(t,{\bm{k}})n_{Vs}(0,-{\bm{k}})\rangle_{c}\ .$ (11)

Note that all these correlators can be studied in mean field, see [27].

At zero explicit breaking $H=0$, we will also consider isotropized $O(4)$
charge and field correlators

$\displaystyle G_{\phi\phi}(t,{\bm{k}})$
$\displaystyle=\frac{1}{4V}\sum_{a}\langle\phi_{a}(t,{\bm{k}})\phi_{a}(0,-{\bm{k}})\rangle_{c},$
(12) $\displaystyle G_{nn}(t,{\bm{k}})$
$\displaystyle=\frac{1}{6V}\sum_{a>b}\langle
n_{ab}(t,{\bm{k}})n_{ab}(0,-{\bm{k}})\rangle_{c}\ .$ (13)

###  II.2 Simulations

This paper is a continuation of [1] and presents results obtained by
numerically solving equations (5a)-(5d), using an algorithm detailed in
Appendix A of [1]. The stepper combines a symplectic integrator for the
Hamiltonian evolution and Metropolis updates for the dissipative step.
Throughout we will present our results in lattice units, setting $g=T_{c}=1$,
and setting the lattice spacing to unity, $a=1$. This amounts to setting
$\hbar=1$ and $T_{c}=1$ and choosing the microscopic lattice spacing so that
the model reproduces the correlation length of the physical system. The pole
frequency of the pion near the critical point is reproduced by this choice of
units. The quartic coupling is set to $\lambda=4$ and the model is simulated
close to its critical mass, $m_{c}^{2}(\lambda)=-4.8110(4)$. The reduced
temperature and magnetic field are defined as in [12] :

$t_{r}\equiv\frac{m_{0}^{2}-m_{c}^{2}}{|m_{c}^{2}|}\,,\qquad
h\equiv\frac{H}{H_{0}}\,,$ (14)

where $H_{0}=5.15(5)$ is chosen so that the magnetization scales as

$\bar{\sigma}=h^{1/\delta}\,,$ (15)

on the critical line, $t_{r}=0$. For $H=0$ and in the broken phase the
magnetization scales as

$\bar{\sigma}=B^{-}(-t_{r})^{\beta}\ ,$ (16)

with amplitude444Note an unfortunate misprint in [1]. Equation (29) in the
published version should read
$\Sigma=b_{1}(m_{c}^{2}-m^{2})^{\beta}(1+C_{T}(m_{c}^{2}-m^{2})^{\omega\nu})$.
The then quoted parameters $b_{1}=0.544\pm 0.004$ and $C_{T}=0.20\pm 0.02$
lead to a parametrization
$\bar{\sigma}=B^{-}(-t_{r})^{\beta}(1+B_{1}(-t_{r})^{\nu\omega})$ , with
$B^{-}=|m_{c}^{2}|^{\beta}b_{1}=0.988\pm 0.007$ and
$B_{1}=C_{T}|m_{c}^{2}|^{\omega\nu}=0.49\pm 0.05$. , $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
$\beta=0.380(2)$ and $\delta=4.824(9)$ from [12]. For the correction-to-
scaling exponent we take $\omega=0.755(5)$ [46] and the for the dynamical
critical exponent we take $\zeta=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 2 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_{c}^{2}$) and
restored ($m^{2}_{0}>m_{c}^{2}$) phases. The dynamical parameters were chosen
as in our previous work, $\chi_{0}=5$, $\Gamma_{0}=1$ and $D_{0}=1/3$. Our
full dataset is summarized in Table 1.

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 [47, 5]. 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 $\tfrac{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.

Phase | $m^{2}_{0}$ | $t_{r}$ | $L$ | $H$
---|---|---|---|---
Restored |  | -4.6300 |  | -4.7100 |  | -4.7337
---|---|---|---|---
-4.6800 |  | -4.7200 |  | -4.7600
-4.7005 |  | -4.7280 |  | -4.7800
| 0.03762 |  | 0.02099 |  | 0.01608
---|---|---|---|---
0.02723 |  | 0.01891 |  | 0.01060
0.02296 |  | 0.01725 |  | 0.006444
$96$ | $0$
Critical line | $m^{2}_{c}=-4.8110$ | 0 | 80 |  | 0.002 | 0.006
---|---
0.003 | 0.01
0.004 |
Broken |  | -4.8236 |  | -4.9118 |
---|---|---|---
-4.8362 |  | -5.0127 |
-4.8614 |  | -5.2143 |
| -0.002620 |  | -0.02096 |
---|---|---|---
-0.005239 |  | -0.04191 |
-0.01048 |  | -0.08383 |
| 48 | 128
---|---
64 |
96 |
0
Table 1: Dataset analyzed in the present work. For all simulations, we fix
$\lambda=4$, $\chi_{0}=5$, $\Gamma_{0}=1$ and $D_{0}=1/3$. $L$ is the size of
our lattices. Figure 1: The correlation length versus temperature in our
simulations close to the $O(4)$ critical point. The correlation length in the
restored phase and on the critical line is obtained from fits described in
App. B.1 and B.2. We use $T/2f^{2}$ as a proxy for the correlation length in
the broken phase, with $T=1$ in our adopted units (see text).

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\neq 0$).

###  II.3 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\xi(T)\ll 1$, the
diffusion equation remains a valid hydrodynamic description of the system:

$\partial_{t}n_{ab}=D\nabla^{2}n_{ab}\,.$ (17)

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):

$G_{AA}(t,{\bm{k}})=G_{VV}(t,{\bm{k}})=T\chi_{0}e^{-Dk^{2}t}\,.$ (18)

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 $\tau\propto\xi^{\zeta}$
with $\zeta=d/2$. The diffusive modes have a relaxation time of $1/Dk^{2}$,
which is of order $\xi^{2}/D$ at the boundary of applicability of the
diffusion equation, $k\sim\xi^{-1}$. In the spirit universality, all modes of
with wavelengths of order $\xi$ should have a similar relaxation time,
dictating the expected scaling of the diffusion coefficient,
$D\propto\xi^{2-\zeta}$.

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\sim t_{r}^{(\zeta-2)\nu}$. 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 non-universal critical amplitude. Nevertheless, we
will extract the critical contribution to $D$ in Sect. III.1.

Figure 2: Left: Axial and vector charge correlators in the restored phase, for
two different momenta. As expected, they decay exponentially (note the log
scale) and the decay rate depends on $k^{2}$. Right: Vector correlator at a
given finite $\bm{k}$ for different masses in the restored phase. A weak
dependence of the diffusion coefficient on the mass is observed, pointing at a
weak critical dependence for the diffusion coefficient; see Fig. 5 for a
quantitative extraction.

###  II.4 The critical line

We next study the critical line, $m^{2}_{0}=m_{c}^{2}$, and study the effect
of explicit symmetry breaking, $H\neq 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 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.

(a)  (b)

Figure 3: Left: Axial correlators as a function of time at fixed momentum on
the critical line ($t_{r}=0$) for different values of the magnetic field $H$.
The axial charge is propagating but strongly damped and we see typical
oscillations. Right: Vector correlator as a function of time at fixed momentum
on the critical line for different magnetic fields. The vector charge remains
diffusive. 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.

###  II.5 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
condensate555For clarity we will restore the temperature in this section..
Following [5, 6], these fluctuations are parametrized by an angle
$\vec{\varphi}(t,{\bm{x}})$ with
$\phi_{a}=(\bar{\sigma},\bar{\sigma}\vec{\varphi}(t,{\bm{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

$\frac{1}{V}\int_{{\bm{x}}}d^{3}x\,\phi_{a}(t,{\bm{x}})=\bar{\sigma}\,n_{a}\,,$
(19)

where $n_{a}$ is a unit vector that is approximately parametrized by
$n_{a}\simeq(1,\vec{\varphi}_{0})$ to first order in the angular deviation.
When the wavelength of $\varphi$ is long compared to the correlation length
$\xi\equiv m_{\sigma}^{-1}$, a hydrodynamic description of $\varphi$ is valid,
and the fluctuations of modes with wavelength order $m_{\sigma}^{-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

$\mathcal{H}[\varphi,n]=\int{\rm
d}^{3}x\,\frac{\vec{n}_{V}^{2}}{2\chi_{0}}+\frac{\vec{n}_{A}^{2}}{2\chi_{0}}+\tfrac{1}{2}f^{2}(T)\,\partial_{i}\vec{\varphi}\cdot\partial^{i}\vec{\varphi}+\tfrac{1}{2}f^{2}m^{2}\varphi^{2}\,,$
(20)

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

$f^{2}m^{2}=H\bar{\sigma}\,.$ (21)

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
form666 In [26] the dissipative coefficient $\Gamma/f^{2}$ was called
$\zeta^{(2)}$. In mean field theory the hydrodynamic parameter $\Gamma$ equals
the bare parameter $\Gamma_{0}$ [27]. A notable feature of these equations
when the symmetry is explicitly broken is that the same coefficient
$\Gamma/f^{2}$ controls both the damping rate at finite momentum and the mass
term, see eq. (65). This is a general result dictated by entropy
considerations [3, 26, 48] and related locality arguments [49]. This became
widely recognized through a sequence of explicit holographic computations in
different contexts with explicitly broken symmetries [50, 51, 52, 53] – for a
review see [54]. The pion dynamics in a holographic model of the phase
transition is described in [43] and can be profitably compared to mean field
theory results of [27].

$\displaystyle\partial_{t}\vec{\varphi}-\frac{\vec{n}_{A}}{\chi}=$
$\displaystyle\frac{\Gamma}{f^{2}}\,\nabla\cdot(f^{2}\nabla\vec{\varphi})\,,$
(22a) $\displaystyle\partial_{t}\vec{n}_{A}-\nabla(f^{2}\nabla\vec{\varphi})=$
$\displaystyle D\nabla^{2}\vec{n}_{A}\,,$ (22b) which reflects from the
Poisson bracket between the phase and the charge,
$\\{\varphi_{s},n_{As^{\prime}}\\}=\delta_{ss^{\prime}}$. The vector charge
remains diffusive $\partial_{t}\vec{n}_{V}=D\nabla^{2}\vec{n}_{V}\,.$ (22c)

The hydrodynamic system is solved in App. A.1 and the resulting pion eigen-
waves have dispersion relations

$\omega(k)-\frac{i}{2}\Gamma(k)=\pm vk-\frac{i}{2}D_{A}k^{2}\,,$ (23)

where

$v^{2}\equiv\frac{f^{2}}{\chi_{0}}\quad\mbox{and}\quad D_{A}\equiv\Gamma+D\,.$
(24)

The form of the correlation function from the hydrodynamic theory is
determined in App. A.1

$G_{\varphi\varphi}(t,{\bm{k}})=\frac{T}{f^{2}k^{2}}e^{-\frac{1}{2}D_{A}k^{2}t}\left[\cos(vkt)+\frac{1}{2vk}(D-\Gamma)k^{2}\sin(vkt)\right]\,.$
(25)

At large wavelengths the pion effective theory determines the static and
dynamic correlation functions. The static correlation function for $k\ll
m_{\sigma}$ can be read off from the free energy in (20) and is singular in
the massless limit:

$G_{\pi\pi}(0,{\bm{k}})=\bar{\sigma}^{2}G_{\varphi\varphi}(0,{\bm{k}})=\bar{\sigma}^{2}\,\frac{T}{f^{2}k^{2}}\,,\qquad
k\ll m_{\sigma}\,.$ (26)

By contrast, the static fluctuations of the $\sigma$ are bounded (if
diverging) and determined by the static susceptibility of the $O(4)$ critical
point, which scales as

$G_{\sigma\sigma}(0,{\bm{k}})\propto m_{\sigma}^{\eta-2}\qquad k\ll
m_{\sigma}\,.$ (27)

The $\sigma$ and $\vec{\pi}$ 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\sim m_{\sigma}$. Since the vev scales as
$\bar{\sigma}^{2}\propto m_{\sigma}^{(d-2+\eta)}$, the pion constant must
scale as $f^{2}\propto m_{\sigma}^{(d-2)}$ [47]. 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 $\phi_{a}$ are dominated by the phase
fluctuations. Specifically, the static correlation function reads

$\displaystyle
G_{\phi\phi}(0,{\bm{k}})=\frac{1}{4}\left(G_{\sigma\sigma}(0,{\bm{k}})+3G_{\pi\pi}(0,{\bm{k}})\right)\simeq\frac{3}{4}\,\bar{\sigma}^{2}\,\frac{T}{f^{2}k^{2}}\,,\qquad
k\ll m_{\sigma}\,,$ (28)

up to corrections of order $(k/m_{\sigma})^{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

$G_{\phi\phi}(t,{\bm{k}})\simeq\frac{3}{4}\,\bar{\sigma}^{2}G_{\varphi\varphi}(t,{\bm{k}})\,,$
(29)

and Fig. 4 (Left) exhibits the $\phi\phi$ 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 $\xi_{\varphi}$ with variance

$\left\langle\xi_{\varphi s}(t,{\bm{x}})\xi_{\varphi
s^{\prime}}(t^{\prime},{\bm{x}}^{\prime})\right\rangle=\frac{2T\Gamma}{f^{2}}\delta_{ss^{\prime}}\delta(t-t^{\prime})\delta({\bm{x}}-{\bm{x}}^{\prime})\,,$
(30)

should be added to (22a). Integrating this equation over the volume shows that
the angular zero mode
$\vec{\varphi}_{0}(t)\equiv\int_{\bm{x}}\vec{\varphi}(t,{\bm{x}})/V$ has a
variance which increases in time

$2\,\mathfrak{D}\,t=\frac{1}{3}\left\langle\vec{\varphi}_{0}(t)\cdot\vec{\varphi}_{0}(t)\right\rangle=2\left(\frac{T\Gamma}{f^{2}V}\right)t\,.$
(31)

This result (which is presented in detail in App. A.2) relates the vev
diffusion coefficient to the pion kinetic coefficients,
$\mathfrak{D}=T\Gamma/f^{2}V$. The diffusion rate decreases with increasing
volume and is suppressed relative to the pion damping rate, $\Gamma
k^{2}\sim\Gamma/L^{2}$, by one power of length, $T/f^{2}L\propto
1/m_{\sigma}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)

$G_{AA}(t,{\bm{k}})=T\chi_{0}e^{-\tfrac{1}{2}D_{A}k^{2}t}\left[\cos(vkt)-\frac{1}{2vk}(D-\Gamma)k^{2}\sin(vkt)\right]\,.$
(32)

The vector charge continues to obey the diffusion equation with corresponding
response functions

$G_{VV}(t,{\bm{k}})=T\chi_{0}e^{-Dk^{2}t}\,.$ (33)

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,{\bm{k}})$

$G_{nn}(t,{\bm{k}})=\frac{1}{6}\left[3G_{AA}(t,{\bm{k}})+3G_{VV}(t,{\bm{k}})\right]\,,$
(34)

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.

  1. (i)

The characteristic frequency of the pion waves below the critical point is of
order $\omega(k)\sim vk$ where $v=\sqrt{f^{2}/\chi_{0}}\propto
m_{\sigma}^{(d-2)/2}$. The scaling of the velocity will be extracted from the
real time correlations of $\phi$, eq. (25), and from its static behavior, eq.
(28), and the results are summarized in Fig. 7 (Right).

  2. (ii)

At the boundary of applicability of the pion EFT $k\sim m_{\sigma}$, the pion
frequency is $\omega(k)\sim vk\sim m_{\sigma}^{d/2}$. But, since the pions are
part of an $O(4)$ multiplet, the order parameter must inherit this dynamical
timescale, $\tau\sim m_{\sigma}^{-d/2}$, setting the dynamical critical
exponent of Model G, $\zeta=d/2$. $\tau$ will be extracted in the restored
phase, Fig. 5 (Left), and on the critical line, Fig. 6 (Left).

  3. (iii)

Close to the critical point and for $k\lesssim m_{\sigma}$, the real and
imaginary parts of $\omega(k)$ should be commensurate if the critical point is
characterized by a single time scale. This reasoning dictates the scaling of
the pion damping coefficient, $D_{A}\propto m_{\sigma}^{d/2-2}$. The scaling
of $D_{A}$ is exhibited in Fig. 7 (Left).

  4. (iv)

Since the axial and vector charges must become degenerate above $T_{c}$, the
diffusion coefficient must have the same scaling as $D_{A}$, leading to the
scaling, $D\sim m_{\sigma}^{d/2-2}$. The scaling of $D$ is exhibited in the
restored phase, Fig. 5 (Right), and on the critical line, Fig. 6 (Right).

Finally, the vev diffusion coefficient $\mathfrak{D}$ is also consistent with
these scalings. Taking the scaling of $f^{2}\propto m_{\sigma}^{d-2}$ and
$\Gamma\propto m_{\sigma}^{d/2-2}$ we see the angular variance given in (31)
grows as

$2\,\mathfrak{D}\,t\propto\frac{t}{(m_{\sigma}L^{2})^{d/2}}\,.$ (35)

At the critical point where $m_{\sigma}L\rightarrow 1$, the angular variance
increases as

$2\,\mathfrak{D}\,t\rightarrow\frac{t}{L^{d/2}}\,,$ (36)

which is again consistent with a universal dynamical critical timescale of
$\tau\propto L^{\zeta}$ with a common critical exponent of $\zeta=d/2$.

(a)  (b)

Figure 4:  Left: Normalized $G_{\phi\phi}$ correlator at $k=2\pi/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\pi/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

###  III.1 $H=0$, restored phase: scaling of the correlation time

We first determine the order parameter relaxation time $\tau$ from the two-
point function $G_{\phi\phi}(t,{\bm{k}})$. Specifically, we define $\tau$ as
the autocorrelation time of $G_{\phi\phi}$ at zero momentum
$G_{\phi\phi}(t)\equiv G_{\phi\phi}(t,{\bm{0}})$

$\tau\equiv\int_{0}^{\infty}\mathrm{d}t\,\frac{G_{\phi\phi}(t)}{G_{\phi\phi}(0)}\
,$ (37)

which is motivated by an exponential decay, $G_{\phi\phi}\propto e^{-t/\tau}$
(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_{\phi\phi}(t)$ is characterized by a single timescale
$\tau$ without significant structure.

(a)  (b)

Figure 5: Left: Dynamical correlation time as a function of $t_{r}$ in the
restored phase. The fit parameters are $\tau^{+}=1.570\pm 0.037,\ \
\tau_{1}=-1.49\pm 0.15$ with $\chi^{2}/\mathrm{dof}=5.74/7$. Right: Vector
diffusion coefficient as a function of $t_{r}$. Note that despite its
smallness, we clearly observe a critical dependence of $D$. The resulting fit
parameters are $D^{+}=0.0190\pm 0.0013,D_{1}=-1.7\pm 0.55$ with
$\chi^{2}/\mathrm{dof}=5.67/7$.

From scaling, we expect the correlation time to behave as

$\tau(t_{r})=\tau^{+}t_{r}^{-\nu\zeta}(1+\tau_{1}t^{\nu\omega})\,,$ (38)

with the critical exponents $\nu$, $\omega$, and $\zeta$ given in Table 2 in
App. B. The resulting fit (with the exponents fixed) is shown in Fig. 5 and we
find

$\displaystyle\tau^{+}$ $\displaystyle=1.570\pm 0.037,$ (39)
$\displaystyle\tau_{1}$ $\displaystyle=-1.49\pm 0.15\ .$ (40)

The growth of $\tau$ near the $T_{c}$ is nicely consistent with $\zeta=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

$G_{VV}(t,{\bm{k}})=T\chi_{0}e^{-Dk^{2}t}\ ,$ (41)

as discussed in the previous section. Similar to the order parameter case, we
have the relation

$D\equiv\frac{1}{k^{2}}\left(\int_{0}^{\infty}\mathrm{d}t\,\frac{G_{VV}(t,{\bm{k}})}{G_{VV}(0,\bm{0})}\right)^{-1}\
,$ (42)

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/{\tau}\sim\xi^{\zeta}$ at the boundary of applicability of the
diffusive description, $k\lesssim\xi$. This means that the diffusion
coefficient should scale with the correlation length as,
$D\propto\xi^{2-\zeta}$.

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

$D(t_{r})=D^{+}t_{r}^{-\nu(2-\zeta)}(1+D_{1}t_{r}^{\omega\nu})+D_{0}\ .$ (43)

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}\approx 0.2425(5)\ .$ (44)

This parameter is different from the bare lattice parameter (also called
$D_{0}=1/3$) which controls the Metropolis updates of the charge at the scale
of the lattice spacing (see [1] in Appendix A.c for further details). Fixing
$D_{0}$ and fitting $D^{+}$ and $D_{1}$ in (43), we get the following
determination of the non-universal amplitudes

$\displaystyle D^{+}$ $\displaystyle=0.0190\pm 0.0013\,,$ (45) $\displaystyle
D_{1}$ $\displaystyle=-1.7\pm 0.55\ .$ (46)

$D^{+}$ is small as was anticipated based on a one loop calculation using mean
field propagators [27].

Using the non-universal amplitude $\xi^{+}$ from the correlation length
$\xi\sim\xi^{+}t_{r}^{-\nu}$ determined in App. B.1, we can combine these
results into a universal dynamical amplitude ratio

$\displaystyle\mathfrak{Q}_{D}^{+}\equiv\frac{D^{+}\tau^{+}}{\xi^{+2}}=0.148\pm
0.011\ ,$ (47)

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.

###  III.2 The critical line

(a)  (b)

Figure 6: Left: Relaxation time for the $\sigma$ field, $\tau(H)$, as function
of the magnetic field $H$ on the critical line, $m^{2}_{0}=m^{2}_{c}$. The
data points are fit with
$\tau(H)=\tau^{H}H^{-\nu_{c}\zeta}(1+\tau_{1}H^{\nu_{c}\omega})$. The results
are $\tau^{H}=1.312\pm 0.022$ and $\tau_{1}^{H}=-1.45\pm 0.06$ with
$\chi^{2}/\mathrm{dof}=6.07/3$. Right: Diffusion coefficient as function of
the magnetic field $H$ on the critical line, $m^{2}_{0}=m^{2}_{c}$. The data
points are fit with
$D(H)=D^{H}H^{-\nu_{c}\zeta}(1+D_{1}^{H}H^{\nu_{c}\omega})+D_{0}$. The results
are $D^{H}=0.023\pm 0.002$ and $D_{1}^{H}=-1.02\pm 0.42$ and $D_{0}=0.2425$ is
kept fixed with $\chi^{2}/\mathrm{dof}=0.72/3$.

We next study the dynamics on the critical line, $m^{2}_{0}=m_{c}^{2}$ and
$H\neq 0$. As in the previous section we define the order parameter relaxation
time $\tau$ from the correlation of the $\sigma$ field

$\tau\equiv\int_{0}^{\infty}\mathrm{d}t\frac{G_{\sigma\sigma}(t,\bm{0})}{G_{\sigma\sigma}(0,\bm{0})}\
.$ (48)

The results for $\tau$ as a function of the magnetic field are shown in Fig.
6, and are fit with the expected scaling form

$\tau(H)=\tau^{H}H^{-\nu_{c}\zeta}(1+\tau_{1}H^{\nu_{c}\omega})\ ,$ (49)

with fixed exponents $\nu_{c}$, $\omega$, and $\zeta$ given in Table 2 of App.
B, and fitted parameters

$\displaystyle\tau^{H}$ $\displaystyle=1.312\pm 0.022\,,$ (50)
$\displaystyle\tau_{1}^{H}$ $\displaystyle=-1.45\pm 0.06\ .$ (51)

As in the restored phase, the subleading term with exponent $\omega$ is
necessary to have a reasonable fit compatible with the expected dynamical
exponent, $\zeta=d/2$.

Similarly, the diffusion coefficient is defined in the same way as in the
restored phase

$D\equiv\frac{1}{k^{2}}\left(\int_{0}^{\infty}\mathrm{d}t\,\frac{G_{VV}(t,{\bm{k}})}{G_{VV}(0,{\bm{k}})}\right)^{-1},\text{
with }|{\bm{k}}|=\frac{2\pi}{L},$ (52)

and we adopt an analogous strategy and functional form to fit the dependence
of $D$ on the correlation length

$D(H)=D^{H}H^{-\nu_{c}(2-\zeta)}(1+D_{1}^{H}H^{\nu_{c}\omega})+D_{0}\,.$ (53)

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

$\displaystyle D^{H}$ $\displaystyle=0.023\pm 0.002\,,$ (54) $\displaystyle
D_{1}^{H}$ $\displaystyle=-1.02\pm 0.42\ .$ (55)

The singular part of the diffusion coefficient diverges next to the critical
point with the expected exponent $H^{-\nu_{c}(2-\zeta)}$ 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

$\displaystyle\mathfrak{Q}_{D}^{H}\equiv\frac{D^{H}\tau^{H}}{\xi^{H2}_{L}}=0.151\pm
0.022\,,$ (56)

where $\xi_{L}^{H}$ 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).

###  III.3 $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_{\phi\phi}$ for
various $t_{r}$ and $L$ and analyze the two lowest non-zero momentum modes,
$k=2\pi/L$ and $k=4\pi/L$. We then perform some fits using the hydro
prediction (25) and extract the parameters $\omega(k)$ and $\Gamma(k)$.
$\omega(k)$ and $\Gamma(k)$ are subsequently extrapolated to infinite volume
using the forms

$\displaystyle\Gamma(k)$ $\displaystyle=D_{A}k^{2}(1+d_{1}/L),$ (57a)
$\displaystyle\omega(k)$ $\displaystyle=vk+v_{2}k^{2}\ ,$ (57b)

with $D_{A}$, $v$, $d_{1}$, and $v_{2}$ as parameters. As discussed in Sect.
II.5, finite volume corrections to the dynamics of the Goldstone modes scale
inversely with the linear extent of the lattice $\sim 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}/\chi_{0}$ and, as we show
below, the result agrees with the dynamic fit using (57). In Sect. II.5 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.

(a)  (b)

Figure 7:  Left: $D_{A}=\Gamma+D$ axial diffusion coefficient as a function of
the reduced temperature $t_{r}$. The dashed line is a fit to the scaling form
$D_{A}=D_{A}^{-}t_{r}^{-\nu(2-\zeta)}(1+D_{A1}t_{r}^{\omega\nu})+D_{A0}$. The
fit results are $D_{A}^{-}=0.044\pm 0.009$, $D_{A1}=1.19\pm 0.13$ and
$D_{A0}=-11\pm 6$. The $\chi^{2}/\mathrm{dof}=6.9/3$. Right: Velocity of the
dispersion relation as a function of $t_{r}$. The data points are the value
extracted from the time dependence of the correlation functions. The dashed
(green) line is the fit to the data with the scaling form
$v=v^{-}\>t_{r}^{\beta/2}(1+v_{1}t_{r}^{\nu\omega})$. The fit results are
$v^{-}=0.462\pm 0.003$ and $v_{1}=0.325\pm 0.03$. The $\chi^{2}/4=15.9/4$. The
continuous (blue) line is described in App. B.3 and in the text and is
obtained by determining $v^{2}=f^{2}/\chi_{0}$ from static correlation
measurements.

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

$D_{A}(t_{r})=D_{A}^{-}|t_{r}|^{-\nu(2-d)}(1+D_{A1}^{-}|t_{r}|^{\omega\nu})+D_{A0}\
,$ (58)

and obtain the following fit parameters

$\displaystyle D_{A}^{-}$ $\displaystyle=0.044\pm 0.009,$ (59) $\displaystyle
D_{A1}^{-}$ $\displaystyle=1.19\pm 0.13,$ (60) $\displaystyle D_{A0}$
$\displaystyle=-11\pm 6.$ (61)

with a $\chi^{2}/\text{dof}$ of $6.3/3$, which is mostly driven by the large
fluctuation around $t_{r}\sim-0.004$. The exponents $\nu$, $\zeta$ and the
subleading $\omega$ 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,

$v(t_{r})=v^{-}\>t_{r}^{\nu/2}(1+v_{1}t_{r}^{\nu\omega})\ ,$ (62)

with $\nu$ and $\omega$ fixed, yielding

$\displaystyle v^{-}$ $\displaystyle=0.462\pm 0.003\,,$ (63) $\displaystyle
v_{1}$ $\displaystyle=0.325\pm 0.03\ .$ (64)

The $\chi^{2}$ is $15.9$ with $4$ degrees of freedom. The relatively large
value of the $\chi^{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=\sqrt{f^{2}/\chi_{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 $\tau$ of the order parameter (which also scales
appropriately), we recast these new results as universal dynamical amplitude
ratios $\mathfrak{Q}_{D}^{+}$, $\mathfrak{Q}_{D}^{H}$, 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.5. 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 $\epsilon$-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.

### Acknowledgements

The authors are grateful to J. Pawlowski, F. Rennecke, S. Schlichting and L.
von Smekal for interesting discussions and A. Soloviev and J. Bhambure for
their collaboration on this work. This work was supported by the U.S.
Department of Energy, Office of Science, Office of Nuclear Physics, Grants
Nos. DE-SC0012704 (AF) and DE-FG88ER41450 (DT).

##  Appendix A Hydrodynamics

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

###  A.1 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.5) for the phase and axial charge
read

$\displaystyle\partial_{t}\vec{\varphi}-\frac{\vec{n}_{A}}{\chi}=$
$\displaystyle\frac{\Gamma}{f^{2}}\,\left[\nabla\cdot(f^{2}\nabla\vec{\varphi})-f^{2}m^{2}\vec{\varphi})\right]+\vec{\xi}_{\varphi}\,,$
(65a)
$\displaystyle\partial_{t}\vec{n}_{A}-\nabla(f^{2}\nabla\vec{\varphi})+f^{2}m^{2}\vec{\varphi}=$
$\displaystyle D\nabla^{2}\vec{n}_{A}+\vec{\xi}_{A}\,,$ (65b) which follows
from the Poisson bracket between the phase and the charge,
$\\{\varphi_{s},n_{As^{\prime}}\\}=\delta_{ss^{\prime}}$. The vector charge
remains diffusive
$\partial_{t}\vec{n}_{V}=D\nabla^{2}\vec{n}_{V}+\vec{\xi}_{V}\,.$ (65c)

The noise added to the right-hand-side of Eq. (65) has variances given by the
Fluctuation-Dissipation Theorem (FDT)

$\displaystyle\left\langle\xi_{\varphi,s}(t,{\bm{x}})\xi_{{\varphi},s^{\prime}}(t^{\prime},{\bm{x}}^{\prime})\right\rangle=$
$\displaystyle\frac{2T\Gamma}{f^{2}}\delta_{ss^{\prime}}\delta(t-t^{\prime})\delta({\bm{x}}-{\bm{x}}^{\prime})\,,$
(66a)
$\displaystyle\left\langle\xi_{As}(t,{\bm{x}})\xi_{As^{\prime}}(t^{\prime},{\bm{x}}^{\prime})\right\rangle=$
$\displaystyle
2T\sigma\delta_{ss^{\prime}}\delta(t-t^{\prime})\left[-\nabla^{2}\delta({\bm{x}}-{\bm{x}}^{\prime})\right]\,,$
(66b)
$\displaystyle\left\langle\xi_{Vs}(t,{\bm{x}})\xi_{Vs^{\prime}}(t^{\prime},{\bm{x}}^{\prime})\right\rangle=$
$\displaystyle
2T\sigma\delta_{ss^{\prime}}\delta(t-t^{\prime})\left[-\nabla^{2}\delta({\bm{x}}-{\bm{x}}^{\prime})\right]\,.$
(66c)

where $\sigma=\chi_{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$)

$\vec{n}_{V}(t,{\bm{k}})=\vec{n}(0,{\bm{k}})e^{-Dk^{2}t}+\int^{t}_{0}{\rm
d}t_{1}\,\vec{\xi}_{V}(t_{1},{\bm{k}})\,e^{Dk^{2}(t_{1}-t)}\,.$ (67)

The variance in the initial condition $\vec{n}(0,{\bm{k}})$ is determined from
the free energy functional in (20) and reads

$\displaystyle\frac{1}{3V}\left\langle\vec{n}_{V}(0,{\bm{k}})\cdot\vec{n}_{V}(0,-{\bm{k}})\right\rangle=$
$\displaystyle T\chi_{0}\,.$ (68)

Alternatively one can integrate from past infinity to determine the initial
condition:

$\vec{n}(0,{\bm{k}})=\int^{0}_{-\infty}{\rm
d}t_{1}\,\vec{\xi}_{V}(t_{1},{\bm{k}})\,e^{Dk^{2}t_{1}}\,.$ (69)

Computing the variance of $\vec{n}(0,{\bm{k}})$ using the variance of the
noise in (66c) reproduces (68), which demonstrates how the noise and
dissipation work together in establishing the equilibrium state. By either
method

$\frac{1}{3V}\left\langle\vec{n}(t,{\bm{k}})\cdot\vec{n}(0,-{\bm{k}})\right\rangle=T\chi_{0}e^{-Dk^{2}t}\,.$
(70)

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
variables777For simplicity, $Y$ will denote a single isospin component of
these fields. $Y=(y_{1},y_{2})=(\omega_{k}\varphi,\mu_{A})$ where
$\mu_{A}=n_{A}/\chi$ and $\omega_{k}^{2}=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)

$\frac{1}{V}\left\langle
y_{a}(0,{\bm{k}})\,y_{b}(0,-{\bm{k}})\right\rangle=\frac{T}{\chi_{0}}\delta_{ab}\,.$
(71)

As in the diffusion example, the static susceptibility sets the initial
conditions for the subsequent evolution. The equations of motion take the form
$\partial_{t}Y+\mathcal{M}Y=\xi$ and eigen-values of $\mathcal{M}$ are

$\lambda_{\pm}=\pm
i\omega_{k}+\frac{1}{2}\Gamma_{k}\,,\qquad\omega_{k}^{2}\equiv
v^{2}(k^{2}+m^{2})\,,\qquad\Gamma_{k}\equiv\Gamma(k^{2}+m^{2})+Dk^{2}\,.$ (72)

The homogeneous solutions are $Y\propto e^{-\lambda_{\pm}t}$. The full matrix
of correlation functions

$[G_{\rm
sym}(t,{\bm{k}})]=\frac{1}{3V}\begin{pmatrix}\omega_{k}^{2}\langle\varphi_{s}(t,{\bm{k}})\varphi_{s}(0,-{\bm{k}})\rangle&\omega_{k}\langle\varphi_{s}(t,{\bm{k}})\mu_{As}(0,-{\bm{k}})\rangle\\\
\omega_{k}\langle\mu_{As}(t,{\bm{k}})\varphi_{s}(0,-{\bm{k}})\rangle&\langle\mu_{As}(t,{\bm{k}})\mu_{As}(0,-{\bm{k}})\rangle\end{pmatrix}\,,$
(73)

reads

$[G_{\rm
sym}(t,{\bm{k}})]=\frac{T}{\chi_{0}}e^{-\frac{1}{2}\Gamma_{k}t}\begin{pmatrix}\cos(\omega_{k}t)+\frac{\Delta_{k}}{\omega_{k}}\sin(\omega_{k}t)&\sin(\omega_{k}t)\\\
-\sin(\omega_{k}t)&\cos(\omega_{k}t)-\frac{\Delta_{k}}{\omega_{k}}\sin(\omega_{k}t)\end{pmatrix}\
,$ (74)

with $\Delta_{k}\equiv(Dk^{2}-\Gamma(k^{2}+m^{2}))/2$. In the body of the text
this result is used in the massless limit in (25) and (32).

###  A.2 Diffusion of the chiral condensate on the coset manifold

####  A.2.1  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 (65), but limiting the
discussion to $H=0$. In finite volume the fields take the form

$\displaystyle\vec{\varphi}(t,{\bm{x}})=$ $\displaystyle
V\vec{\varphi}_{0}(t)+\sum_{{\bm{k}}\neq
0}e^{i{\bm{k}}\cdot{\bm{x}}}\vec{\varphi}(t,{\bm{k}})\,,$ (75)
$\displaystyle\vec{n}_{A}(t,{\bm{x}})=$ $\displaystyle\sum_{{\bm{k}}\neq
0}e^{i{\bm{k}}\cdot{\bm{x}}}\vec{n}(t,{\bm{k}})\,.$ (76)

The axial charge is exactly conserved and there is no net charge in the box,
and consequently the ${\bm{k}}=0$ mode is zero for $\vec{n}_{A}$. Then,
integrating (22a) over volume, the zero mode of $\varphi$ satisfies a random
walk

$\partial_{t}\vec{\varphi}_{0}=\vec{\xi}_{0}\,,\qquad\vec{\xi}_{0}(t)\equiv\int_{\bm{x}}\vec{\xi}_{\varphi}(t,{\bm{x}})/V\,.$
(77)

Solving for $\varphi_{0}$ and computing its variance leads to (31), which is
reproduced here for convenience

$2\,\mathfrak{D}\,t\equiv\frac{1}{3}\left\langle\vec{\varphi}_{0}(t)\cdot\vec{\varphi}_{0}(t)\right\rangle=2\left(\frac{T\Gamma}{f^{2}V}\right)t\,,$
(78)

i.e. $\mathfrak{D}\equiv T\Gamma/f^{2}V$. The noise satisfies

$\left\langle\xi_{\varphi 0s}(t)\xi_{\varphi
0s^{\prime}}(t^{\prime})\right\rangle=2\mathfrak{D}\delta_{ss^{\prime}}\delta(t-t^{\prime})\,.$
(79)

The probability density $\mathcal{P}(\vec{\varphi}_{0})$ satisfies the
diffusion equation

$\partial_{t}\mathcal{P}=\mathfrak{D}\frac{\partial}{\partial\varphi_{0s}}\left(\delta^{ss^{\prime}}\frac{\partial\mathcal{P}}{\partial\varphi_{0s^{\prime}}}\right)\,.$
(80)

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, $\varphi_{0s}$, assuming they are small. The coordinates
$\vec{\varphi}$ provide a local Cartesian coordinate patch for the three
sphere close to the pole. In a general coordinate system the diffusion
equation reads

$\partial_{t}\mathcal{P}=\mathfrak{D}\nabla_{I}\nabla^{I}\mathcal{P}\,,$ (81)

where $\nabla_{I}\nabla^{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})=(\varphi_{S},\theta,\Phi)$

$ds^{2}=\mathrm{d}\varphi_{S}^{2}+\sin^{2}\varphi_{S}\,(\mathrm{d}\theta^{2}+\sin^{2}\theta\,\mathrm{d}\Phi^{2})\,,$
(82)

where $\varphi_{S}$ is approximately
$\varphi_{S}\simeq\sqrt{\vec{\varphi}^{2}}$ close to the north pole.

####  A.2.2 Explicit symmetry breaking, forced diffusion, and the
$\epsilon$-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

$\displaystyle\partial_{t}\vec{\varphi}_{0}-\frac{\vec{n}_{A0}}{\chi}=$
$\displaystyle-\Gamma m^{2}\vec{\varphi}_{0}+\vec{\xi}_{\varphi_{0}}\,,$ (83)
$\displaystyle\partial_{t}\vec{n}_{A0}+f^{2}m^{2}\vec{\varphi}_{0}=$
$\displaystyle 0\,.$ (84)

where the statistics of $\vec{\xi}_{\varphi_{0}}$ are given in (79). The
variables $\vec{\varphi}_{0}$ and $\vec{n}_{A0}$ are canonical conjugates
$\\{\varphi_{0s},n_{A0s^{\prime}}\\}=\delta_{ss^{\prime}}$. To illustrate the
structure of these equations we introduce a traditional notation

$(Q^{s},P_{s})\equiv(\varphi_{0s},n_{A0s})\,,$ (85)

and the equations of motion are written

$\displaystyle\partial_{t}Q^{s}+\\{\mathcal{H}_{0},Q^{s}\\}=$
$\displaystyle-\beta V\mathfrak{D}\frac{\partial\mathcal{H}_{0}}{\partial
Q^{s}}+\xi_{0}^{s}\,,$ (86)
$\displaystyle\partial_{t}P_{s}+\\{\mathcal{H}_{0},P_{s}\\}=$ $\displaystyle
0\,.$ (87)

Here $\beta$ and $V$ are the inverse temperature and volume respectively, and
$\mathfrak{D}=T\Gamma/f^{2}V$ is the diffusion coefficient introduce earlier.
Finally the Hamiltonian for the zero modes is

$\mathcal{H}_{0}(Q,P)=\frac{{\vec{P}}^{\,2}}{2\chi}+\frac{1}{2}f^{2}m^{2}\vec{Q}^{2}\,.$
(88)

The Fokker-Planck equation for the phase space probability density
$\mathcal{P}(Q,P)$ follows from the Langevin dynamics and reads [56, 57, 58]

$\partial_{t}\mathcal{P}+\left\\{\mathcal{P},\mathcal{H}_{0}\right\\}=\beta
V\mathfrak{D}\,\frac{\partial}{\partial
Q^{s}}\left(\delta^{ss^{\prime}}\frac{\partial\mathcal{H}_{0}}{\partial
Q^{s^{\prime}}}\,\mathcal{P}\right)+\mathfrak{D}\frac{\partial}{\partial
Q^{s}}\left(\delta^{ss^{\prime}}\frac{\partial\mathcal{P}}{\partial
Q^{s^{\prime}}}\right)\,.$ (89)

We emphasize that the canonical measure for the probability density
$\mathcal{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 (82).
The Hamiltonian is generalized to include larger angles and reads

${\mathcal{H}}_{0}(q,p)=\frac{g^{IJ}(q)}{2\chi}\,p_{I}p_{J}-f^{2}m^{2}\cos(\varphi_{S})\,,$
(90)

where $g^{IJ}(q)$ is the metric of the sphere in (82). Here we used the
original free energy in (3) and the Gell-Mann, Oakes, Renner relation to
deduce the potential for the zero mode

$\frac{-1}{V}\int_{\bm{x}}\,H_{a}\phi_{a}(t,{\bm{x}})=-H\bar{\sigma}\cos\varphi_{S}=-f^{2}m^{2}\cos(\varphi_{S})\,.$
(91)

The change of coordinates leads to the Fokker-Planck equation in its final
form

$\partial_{t}\mathcal{P}+\left\\{\mathcal{P},\mathcal{H}_{0}\right\\}=\beta
V\mathfrak{D}\,\nabla_{I}\left(\nabla^{I}\mathcal{H}_{0}\,\mathcal{P}\right)+\mathfrak{D}\,\nabla_{I}\nabla^{I}\mathcal{P}\,,$
(92)

where $\nabla_{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

$\mathcal{P}(q,p)\,{\mathrm{d}}^{3}q\,\mathrm{d}^{3}p=e^{-\beta
V\mathcal{H}_{0}}\mathrm{d}^{3}q\,\mathrm{d}^{3}p\,.$ (93)

If the momenta are of no interest they can be integrated over yielding

$\mathcal{P}(q)\,\sqrt{g}\,{\mathrm{d}}^{3}q\,=e^{\beta
Vf^{2}m^{2}\cos\varphi_{S}}\sqrt{g}\,\mathrm{d}^{3}q\,,$ (94)

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 $\beta Vf^{2}m^{2}\sim 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.

##  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 2.

exponent | definition | value | Ref.
---|---|---|---
$\beta$ | $\bar{\sigma}\propto(-t_{r})^{\beta}$ | 0.380(2) | [12]
$\nu$ | $\xi\propto t_{r}^{-\nu}$ | 0.7377(41) | [12]
$\gamma$ | $\chi\propto t_{r}^{-\gamma}$ | 1.4531(104) | [12]
$\delta$ | $\bar{\sigma}\propto H^{1/\delta}$ | 4.824(9) | [12]
$\nu_{c}$ | $\xi_{L,T}\propto H^{-\nu_{c}}$ | 0.4024(2) | [12]
$1-1/\delta$ | $\chi_{L,T}\propto H^{-1+1/\delta}$ | 0.7927(4) | [12]
$\omega$ | $Y\propto\xi^{\theta(Y)}(1+c(Y)\xi^{-\omega}+\dots)$ | 0.755(5) | [46]
$\zeta$ | $\tau\propto\xi^{\zeta}$ | $d/2$ | [3]
Table 2: Critical exponents of the $O(4)$ model used in this work. The top of
the table is for zero field, $H=0$ and $t_{r}=(m_{0}^{2}-m_{c}^{2})/m_{c}^{2}$
is the reduced temperature. The middle of the table is on the critical line,
$t_{r}=0$ and $H\neq 0$. The bottom of the table shows the correction-to-
scaling exponent $\omega$, which is also universal. In the table we denote
“any scaling quantity” by $Y$, and $\theta(Y)$ and $c(Y)$ are its associated
critical exponent and non-universal subleading amplitude respectively. $\zeta$
is the dynamical scaling exponent of “Model G” where $d=3$ is the number of
spatial dimensions.

###  B.1 Restored phase

The magnetic susceptibility in the restored phase is defined as

$\chi({\bm{k}})=\frac{1}{4V}\sum_{a=0}^{3}\left\langle|\phi_{a}({\bm{k}})|^{2}\right\rangle\
,$ (95)

and the zero momentum limit is notated with $\chi$

$\chi\equiv\chi(\bm{0})\ .$ (96)

A standard estimator of the correlation length is the so-called “second-moment
correlation length” [46], computed as

$\xi=\sqrt{\frac{\chi/F-1}{4\sin^{2}{\left(\frac{\pi}{L}\right)}}}\ ,$ (97)

with $F$ the susceptibility of the first Fourier mode

$F=\chi\left({\bm{k}}\right)|_{k=\frac{2\pi}{L}}.$ (98)

We fit their dependence on the reduced temperature
$t_{r}=(m_{0}^{2}-m^{2}_{c})/{m_{c}^{2}}$ to

$\displaystyle\xi(t_{r})$
$\displaystyle=\xi^{+}t_{r}^{\nu}(1+\xi_{1}t_{r}^{\omega\nu})\ ,$ (99)
$\displaystyle\chi(t_{r})$
$\displaystyle=C^{+}t_{r}^{\gamma}(1+C_{1}t_{r}^{\omega\nu})\ .$ (100)

The exponent $\omega$ is the first subleading scaling exponent. To perform the
fit, we fix $\beta$, $\nu$ and $\omega$ to the values presented in Tab. 2 and
fit $\xi^{+}$, $\xi_{1}$, $C^{+}$ and $C_{1}$. The fits are presented in Fig.
8. We obtain

$\xi^{+}=0.4492\pm 0.0036,\ \ \xi_{1}=0.332\pm 0.085,\ \ C^{+}=0.2077\pm
0.0032,\ \ C_{1}=1.11\pm 0.15\ .$ (101)

The chi-square for the fit of the correlation length is
$\chi^{2}/\mathrm{dof}=9.1/7$ and the one for the susceptibility is
$\chi^{2}/\mathrm{dof}=4.83/7$. These values, together with the value the non-
universal value for the magnetization in the restored phase

$\bar{\sigma}(t_{r})=B^{-}(-t_{r})^{\beta}(1+B_{1}(-t_{r})^{\nu\omega})\ ,$
(102)

that was extracted in Ref. [1]888Note an unfortunate misprint. Equation (29)
in the published version should read
$\Sigma=b_{1}(m_{c}^{2}-m_{0}^{2})^{\beta}(1+C_{T}(m_{c}^{2}-m_{0}^{2})^{\omega\nu})$.
The then quoted parameters $b_{1}=0.544\pm 0.004$ and $C_{T}=0.20\pm 0.02$
lead to $B^{-}=|m_{c}^{2}|^{\beta}b_{1}=0.988\pm 0.007$ and
$B_{1}=C_{T}|m_{c}^{2}|^{\omega\nu}=0.49\pm 0.05$ .

$B^{-}=0.988\pm 0.007\,,$ (103)

allow us to compute the universal ratio [59]

$Q_{c}=\frac{(B^{-})^{2}\xi^{3}}{C^{+}}=0.4266\pm 0.0126\,.$ (104)

This value is indeed compatible with the $Q_{c}=0.4231(9)$ of [59].

(a)  (b)

Figure 8: Left: Correlation length as a function of the reduced temperature.
The fitted non universal amplitude is $\xi^{+}=0.4492\pm 0.0036$ and the
subleading amplitude is $\xi_{1}=0.33\pm 0.09$. The $\chi^{2}/\mathrm{dof}$ is
$9.1/7$. Right: Magnetic susceptibility as a function of the reduced
temperature. The fitted non universal amplitude is $C^{+}=0.208\pm 0.003$ and
the subleading correction is $C_{1}=1.11\pm 0.15$. The $\chi^{2}/\mathrm{dof}$
is $4.83/7$. In both cases, the error bars are too small to be visible on the
plot.

###  B.2 Critical line $m_{0}^{2}=m_{c}^{2}$

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 $\chi_{L}$ susceptibilities:

$\displaystyle\chi_{L}({\bm{k}})$
$\displaystyle=\frac{1}{V}\langle|\phi_{0}({\bm{k}})|^{2}\rangle_{c}\,,$ (105)
$\displaystyle\chi_{T}({\bm{k}})$
$\displaystyle=\frac{1}{3V}\sum_{s=1}^{3}\langle|\phi_{s}({\bm{k}})|^{2}\rangle_{c}\,.$
(106)

The static susceptibilities are shown in Fig. 9. As a fitting function, we
chose the scaling prediction with the first subleading correction

$\chi(H)=\chi H^{-1+\frac{1}{\delta}}(1+\chi_{1}H^{\nu_{c}\omega})\,.$ (107)

(a)  (b)

Figure 9: Left: Longitudinal susceptibility $\chi_{L}$ as a function of the
magnetic field $H$. The data points are fit with
$\chi_{L}(H)=\chi_{L}^{H}H^{-1+1/\delta}(1+\chi^{H}_{L1}H^{\nu_{c}\omega})$.
The results are $\chi_{L}^{H}=0.1564\pm 0.0016$ and $\chi^{H}_{L1}=0.45\pm
0.06$ with $\chi^{2}/\mathrm{dof}=0.26/3$. Right: Perpendicular susceptibility
$\chi_{T}$ as function of the magnetic field $H$. The data points are fit with
$\chi_{T}(H)=\chi^{H}_{T}H^{-1+1/\delta}(1+\chi^{H}_{T1}H^{\nu_{c}\omega})$.
The results are $\chi_{T}^{H}=0.241\pm 0.004$ and $\chi^{H}_{T1}=0.25\pm 0.08$
with $\chi^{2}/\mathrm{dof}=6.17/3$.

To extract the correlation length in the longitudinal and the transverse
directions we measure the second-moment correlation length as in (97) with

$\displaystyle F_{L}$
$\displaystyle=\chi_{L}\left({\bm{k}}\right)|_{k=\frac{2\pi}{L}}\,,$ (108)
$\displaystyle F_{T}$
$\displaystyle=\chi_{T}\left({\bm{k}}\right)|_{k=\frac{2\pi}{L}}\,.$ (109)

The scaling prediction fixes the fitting function as

$\displaystyle\xi_{L}(H)$
$\displaystyle=\xi_{L}^{H}H^{-\nu_{c}}(1+\xi_{L1}^{H}H^{\nu_{c}\omega})\,,$
(110) $\displaystyle\xi_{T}(H)$
$\displaystyle=\xi_{T}^{H}H^{-\nu_{c}}(1+\xi_{T1}^{H}H^{\nu_{c}\omega})\,,$
(111)

and the result of the fit is shown in Fig. 10.

(a)  (b)

Figure 10: Left: Longitudinal correlation length $\xi_{L}$ as function of the
of the magnetic field $H$ . The data points are fit with
$\xi_{L}=\xi_{L+}H^{-\nu_{c}}(1+\xi_{L1}H^{\nu_{c}\omega})$. The results are
$\xi_{L+}=0.447\pm 0.026$ and $\xi_{L1}=-0.14\pm 0.32$ with
$\chi^{2}/\mathrm{dof}=0.19/3$. Right: Perpendicular correlation length
$\xi_{T}$ as function of the magnetic field $H$. The data points are fit with
$\xi_{T}=\xi_{T+}H^{-\nu_{c}}(1+\xi_{T1}H^{\nu_{c}\omega})$. The results are
$\xi_{T+}=0.826\pm 0.027$ and $\xi_{T1}=0.24\pm 0.19$ with
$\chi^{2}/\mathrm{dof}=6.76/3$.

The universal ratio $\xi_{T+}/\xi_{L+}=1.85\pm 0.12$ is consistent with the
more precise estimate done in [59].

###  B.3 Broken phase and the pion velocity

Our goal in this appendix is to extract the infinite volume chiral condensate
$\bar{\sigma}(T)$ and pion velocity $v^{2}(T)=f^{2}(T)/\chi_{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

$M_{a}(t)\equiv\frac{1}{V}\sum_{{\bm{x}}}\phi_{a}(t,{\bm{x}})\,,$ (112)

the static correlation function $G_{\phi\phi}(0,{\bm{k}})$, and their their
dependence on volume to extract $\bar{\sigma}(T)$ and
$(\bar{\sigma}^{2}/f^{2})(T)$ in the infinite volume limit.

In infinite volume the magnetization remains constant in time and determines
the chiral condensate

$\lim_{V\rightarrow\infty}M_{a}=\bar{\sigma}(T)\,n_{a}\,.$ (113)

Here $n_{a}$ is an arbitrary unit vector on the three-sphere characterizing
the orientation of the condensate. However, in any finite volume

$\left\langle M_{a}\right\rangle=0\,,$ (114)

since the condensate orientation vector $n_{a}$ stochastically explores the
$S_{3}$ in time – see App. A.2.

One way to extract $\bar{\sigma}$ is to look at the fluctuations of $M_{a}$,
evaluating $\left\langle M^{2}\right\rangle=\left\langle
M_{a}M_{a}\right\rangle$, which is approximately $\bar{\sigma}^{2}$ at large
volumes. The leading deviation of $\left\langle M^{2}\right\rangle$ from
$\bar{\sigma}^{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 $\bar{\sigma}$ is also a parameter.

Apart from a slightly different fitting procedure, the current determination
of ${\bar{\sigma}(T)}$ simply repeats the analysis done in our earlier work
for the larger lattices and datasets used here [1]. The expansion relating
$\left\langle M^{2}\right\rangle$ and $\bar{\sigma}^{2}$ takes the form [47,
1]

$\left\langle
M^{2}\right\rangle=\bar{\sigma}^{2}\left[1+\frac{0.677355}{f^{2}L}+\frac{0.156028}{f^{4}L^{2}}+O\left((f^{2}L)^{-3}\right)\right]\,,$
(115)

and is valid for $f^{2}(T)L\gg 1$. The numerical coefficients are known
analytically in terms of “shape coefficients” $\beta_{1}$ and $\beta_{2}$, but
are of no interest here.

To determine the chiral condensate, we plotted $\left\langle
M^{2}\right\rangle$ versus $L$ at fixed temperature, and made various fits
with (115) to determine $\bar{\sigma}^{2}(T)$ and
$(\bar{\sigma}^{2}/f^{2})(T)$. We estimated the systematic uncertainty in the
extracted values of $\bar{\sigma}(T)$ by adding a $C/L^{3}$ term to the fit
form and comparing to a fit based on (115). Except for the point closest to
$T_{c}$ the systematic uncertainty $\bar{\sigma}(T)$ is smaller than our
statistical uncertainty, which would not have been the case if only the first
term $\propto 1/f^{2}L$ in the expansion had been used.

Our results for $\bar{\sigma}(T)$ are shown in the right panel of Fig. 11,
This is fit with the functional form

$\bar{\sigma}(t_{r})=B^{-}(-t_{r})^{\beta}\left(1+B_{1}(-t_{r})^{\omega\nu}\right)\,.$
(116)

with fit parameters $B^{-}=0.988\pm 0.007$, $B_{1}=0.51\pm 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})^{\beta}$. 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

$k\sim\frac{1}{L}\ll m_{\sigma}\ll\frac{1}{a}\,,$ (117)

is dominated by the massless Goldstone mode. Recalling eq. (28) with $T=1$,
the correlation function takes the form

$\frac{4}{3}G_{\phi\phi}(0,{\bm{k}})=\frac{\bar{\sigma}^{2}}{f^{2}k^{2}}\ ,$
(118)

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_{\sigma}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

$\frac{4}{3}G_{\phi\phi}(0,{\bm{k}})=\frac{\bar{\sigma}^{2}}{f^{2}k^{2}}\left[1+\frac{\beta_{1}}{f^{2}L}-\frac{N-2}{f^{2}k^{2}V}+\mathcal{O}\left(\frac{1}{(m_{\sigma}^{2}L)^{2}}\right)\right]\,.$
(119)

where $\beta_{1}=0.225785$ is a “shape coefficient” reflecting the
fluctuations of pions in a cubic box. Numerically for $k=2\pi/L$ and $N=4$:

$\displaystyle\frac{4}{3}k^{2}G_{\phi\phi}(0,{\bm{k}})\Big{|}_{k=2\pi/L}=$
$\displaystyle\frac{\bar{\sigma}^{2}}{f^{2}}\left(1+\frac{0.175124}{f^{2}L}\right)\,.$
(120)

To extract $\bar{\sigma}^{2}/f^{2}$ we plotted $k^{2}G_{\phi\phi}(0,{\bm{k}})$
for the first Fourier mode as a function of system size and fitted the result
for $(\bar{\sigma}^{2}/f^{2})(T)$ and $1/f^{2}$ using (120). The subleading
correction $\sim 1/f^{2}L$ in (120) was treated as a free parameter, although
it is broadly consistent with the $\bar{\sigma}^{2}(T)$ extracted above and
the $(\bar{\sigma}^{2}/f^{2})(T)$ determined here.

The data on $(\bar{\sigma}^{2}/f^{2})(T)$ are shown in Fig. 11.
Asymptotically, $\bar{\sigma}^{2}/f^{2}$ should behave as

$\frac{\bar{\sigma}^{2}}{f^{2}}\sim(-t_{\rm r})^{\eta}\,,$ (121)

where $\eta=0.030$ is a critical exponent. We are unable to extract $\eta$
from our measurements. Instead, we have simply fit a form

$\frac{\bar{\sigma}^{2}}{f^{2}}=C_{0}(-t_{\rm r})^{C_{1}}\,,$ (122)

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 $\eta$ is non-zero.

Figure 11: Left: Magnetization in the broken phase for zero magnetic field
after an extrapolation to infinite volume developed in [47]. The dashed line
is a fit to the data with
$\bar{\sigma}(t_{r})=B^{-}(-t_{r})^{\beta}(1+B_{1}(-t_{r})^{\omega\nu})$
including the subleading correction. The dotted line shows the leading term of
the fit. Right: An estimate of the pion constant $f^{2}$ (or more precisely
$\bar{\sigma}^{2}/f^{2}$) based on static correlations, together with a
phenomenological fit ansatz.

Given the fit results for $\bar{\sigma}^{2}(T)$ in (116) and (122) for
$(\bar{\sigma}^{2}/f^{2})(T)$ we can extract the velocity

$v^{2}=\frac{f^{2}}{\chi_{0}}=\frac{1}{\chi_{0}}\frac{\bar{\sigma}^{2}(T)}{(\bar{\sigma}^{2}/f^{2})(T)}\,.$
(123)

The uncertainty is dominated by the contribution from $\bar{\sigma}^{2}(T)$ as
$\bar{\sigma}^{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.

##  Appendix C Axial channel systematics in the broken phase

We record in this appendix our systematic analysis leading to the extraction
of $\Gamma(k)$ and $\omega(k)$ in the broken phase.

To validate the fitting form (57)

$\displaystyle\Gamma(k)$ $\displaystyle=D_{A}k^{2}(1+d_{1}/L)\,,$ (124)
$\displaystyle\omega(k)$ $\displaystyle=vk+v_{2}k^{2}\,,$ (125)

and the specific dependence on the size $L$ of the axial diffusion coefficient
$D_{A}$, we fit the damping coefficient $\Gamma$ as function of the reduce
temperature $t_{r}$ using three different datasets: $k=2\pi/L$, $k=4\pi/L$ and
the combined sample. Figure 12 shows the results obtained for $D_{A}$ and the
leading volume correction $d_{1}$ across the three different datasets. The
consistency of the fit among the different datasets and the values of $d_{1}$
exemplify the importance of the size-dependent $d_{1}/L$ term in (57) in the
axial diffusion coefficient. Conversely, the remaining systematic dispersion
between the different datasets can be considered as an evaluation of the
remaining systematic errors.

The same strategy is adopted for the extraction of the parameter for the
dispersion curves, $\omega(k)$ that is fitted as

$\omega(k)=vk+v_{2}k^{2}\ ,$ (126)

where the quadratic term $v_{2}k^{2}$ was needed to improve the quality of the
fit. In Figure 13 we show the results for the coefficient $v$ and $v_{2}$. We
see that in this case, the systematic errors are smaller and that the fit
agrees very well across all different datasets.

Let us conclude this appendix by commenting on our assessment of statistical
errors for this analysis. We subdivided each simulation in the time direction
into blocks and fit each block independently. We then used the mean and the
standard error of this set of block fits as an estimator for the best fit
parameters and their statistical error.

(a)  (b)

Figure 12: Result of the fit of the damping rate
$\Gamma(k)=D_{A}k^{2}(1+\frac{d_{1}}{L})$. Left: $D_{A}$ coefficient as a
function of $t_{r}$ extracted for $k=2\pi/L$ (black dots), $k=4\pi/L$ (purple
dots) and the two datasets together (green dots). Right: $d_{1}$ coefficient
as a function of $t_{r}$ extracted for $k=2\pi/L$ (black dots), $k=4\pi/L$
(purple dots) and the two datasets together (green dots).

(a)  (b)

Figure 13:  Result of the fit of the dispersion curves
$\omega(k)=vk+v_{2}k^{2}$. Left: $v$ coefficient as a function of $t_{r}$
extracted for $k=2\pi/L$ (black dots), $k=4\pi/L$ (purple dots) and the two
datasets together (green dots). Right: $v_{2}$ coefficient as a function of
$t_{r}$ extracted for $k=2\pi/L$ (black dots), $k=4\pi/L$ (purple dots) and
the two datasets together (green dots).

## References

  * Florio  _et al._ [2022] A. Florio, E. Grossi, A. Soloviev, and D. Teaney, Dynamics of the $O(4)$ critical point in QCD, Phys. Rev. D 105, 054512 (2022), arXiv:2111.03640 [hep-lat] .
  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin, Theory of Dynamic Critical Phenomena, Rev. Mod. Phys. 49, 435 (1977).
  * Rajagopal and Wilczek [1993] K. Rajagopal and F. Wilczek, Static and dynamic critical phenomena at a second order QCD phase transition, Nucl. Phys. B399, 395 (1993), arXiv:hep-ph/9210253 [hep-ph] .
  * Son [2000] D. T. Son, Hydrodynamics of nuclear matter in the chiral limit, Phys. Rev. Lett. 84, 3771 (2000), arXiv:hep-ph/9912267 [hep-ph] .
  * Son and Stephanov [2002a] D. T. Son and M. A. Stephanov, Real time pion propagation in finite temperature QCD, Phys. Rev. D66, 076011 (2002a), arXiv:hep-ph/0204226 [hep-ph] .
  * Son and Stephanov [2002b] D. T. Son and M. A. Stephanov, Pion propagation near the QCD chiral phase transition, Phys. Rev. Lett. 88, 202302 (2002b), arXiv:hep-ph/0111100 [hep-ph] .
  * Pisarski and Wilczek [1984] R. D. Pisarski and F. Wilczek, Remarks on the Chiral Phase Transition in Chromodynamics, Phys. Rev. D 29, 338 (1984).
  * Wilczek [1992] F. Wilczek, Application of the renormalization group to a second order QCD phase transition, Int. J. Mod. Phys. A 7, 3911 (1992), [Erratum: Int.J.Mod.Phys.A 7, 6951 (1992)].
  * Toussaint [1997] D. Toussaint, Scaling functions for O(4) in three-dimensions, Phys. Rev. D 55, 362 (1997), arXiv:hep-lat/9607084 .
  * Engels and Vogt [2010] J. Engels and O. Vogt, Longitudinal and transverse spectral functions in the three-dimensional O(4) model, Nucl. Phys. B 832, 538 (2010), arXiv:0911.1939 [hep-lat] .
  * Engels and Karsch [2012] J. Engels and F. Karsch, The scaling functions of the free energy density and its derivatives for the 3d O(4) model, Phys. Rev. D 85, 094506 (2012), arXiv:1105.0584 [hep-lat] .
  * Engels and Karsch [2014] J. Engels and F. Karsch, Finite size dependence of scaling functions of the three-dimensional O(4) model in an external field, Phys. Rev. D 90, 014501 (2014), arXiv:1402.5302 [hep-lat] .
  * Butti  _et al._ [2003] A. Butti, A. Pelissetto, and E. Vicari, On the nature of the finite temperature transition in QCD, JHEP 08, 029, arXiv:hep-ph/0307036 .
  * Braun and Klein [2009] J. Braun and B. Klein, Finite-Size Scaling behavior in the O(4)-Model, Eur. Phys. J. C 63, 443 (2009), arXiv:0810.0857 [hep-ph] .
  * Braun  _et al._ [2011] J. Braun, B. Klein, and P. Piasecki, On the scaling behavior of the chiral phase transition in QCD in finite and infinite volume, Eur. Phys. J. C 71, 1576 (2011), arXiv:1008.2155 [hep-ph] .
  * Braun  _et al._ [2012] J. Braun, B. Klein, and B.-J. Schaefer, On the Phase Structure of QCD in a Finite Volume, Phys. Lett. B 713, 216 (2012), arXiv:1110.0849 [hep-ph] .
  * Braun  _et al._ [2020] J. Braun, W.-j. Fu, J. M. Pawlowski, F. Rennecke, D. Rosenblüh, and S. Yin, Chiral Susceptibility in (2+1)-flavour QCD, (2020), arXiv:2003.13112 [hep-ph] .
  * Ding  _et al._ [2019] H. T. Ding _et al._ (HotQCD), Chiral Phase Transition Temperature in ( 2+1 )-Flavor QCD, Phys. Rev. Lett. 123, 062002 (2019), arXiv:1903.04801 [hep-lat] .
  * Kotov  _et al._ [2021] A. Y. Kotov, M. P. Lombardo, and A. Trunin, QCD transition at the physical point, and its scaling window from twisted mass Wilson fermions, Phys. Lett. B 823, 136749 (2021), arXiv:2105.09842 [hep-lat] .
  * Cuteri  _et al._ [2021] F. Cuteri, O. Philipsen, and A. Sciarra, On the order of the QCD chiral phase transition for different numbers of quark flavours, JHEP 11, 141, arXiv:2107.12739 [hep-lat] .
  * Kaczmarek  _et al._ [2020] O. Kaczmarek, F. Karsch, A. Lahiri, L. Mazur, and C. Schmidt, QCD phase transition in the chiral limit (2020) arXiv:2003.07920 [hep-lat] .
  * Nijs  _et al._ [2021] G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Transverse Momentum Differential Global Analysis of Heavy-Ion Collisions, Phys. Rev. Lett. 126, 202301 (2021), arXiv:2010.15130 [nucl-th] .
  * Mazeliauskas and Vislavicius [2020] A. Mazeliauskas and V. Vislavicius, Temperature and fluid velocity on the freeze-out surface from $\pi$, $K$, $p$ spectra in pp, p-Pb and Pb-Pb collisions, Phys. Rev. C 101, 014910 (2020), arXiv:1907.11059 [hep-ph] .
  * Devetak  _et al._ [2020] D. Devetak, A. Dubla, S. Floerchinger, E. Grossi, S. Masciocchi, A. Mazeliauskas, and I. Selyuzhenkov, Global fluid fits to identified particle transverse momentum spectra from heavy-ion collisions at the Large Hadron Collider, JHEP 06, 044, arXiv:1909.10485 [hep-ph] .
  * ALI [2018] Expression of Interest for an ALICE ITS Upgrade in LS3 (2018).
  * Grossi  _et al._ [2020] E. Grossi, A. Soloviev, D. Teaney, and F. Yan, Transport and hydrodynamics in the chiral limit, Phys. Rev. D 102, 014042 (2020), arXiv:2005.02885 [hep-th] .
  * Grossi  _et al._ [2021] E. Grossi, A. Soloviev, D. Teaney, and F. Yan, Soft pions and transport near the chiral critical point, Phys. Rev. D 104, 034025 (2021), arXiv:2101.10847 [nucl-th] .
  * Berges  _et al._ [2010] J. Berges, S. Schlichting, and D. Sexty, Dynamic critical phenomena from spectral functions on the lattice, Nucl. Phys. B 832, 228 (2010), arXiv:0912.3135 [hep-lat] .
  * Täuber [2017] U. C. Täuber, Phase Transitions and Scaling in Systems Far from Equilibrium, Ann. Rev. Condensed Matter Phys. 8, 185 (2017), arXiv:1604.04487 [cond-mat.stat-mech] .
  * Schlichting  _et al._ [2020] S. Schlichting, D. Smith, and L. von Smekal, Spectral functions and critical dynamics of the $O(4)$ model from classical-statistical lattice simulations, Nucl. Phys. B 950, 114868 (2020), arXiv:1908.00912 [hep-lat] .
  * Schweitzer  _et al._ [2020] D. Schweitzer, S. Schlichting, and L. von Smekal, Spectral functions and dynamic critical behavior of relativistic $Z_{2}$ theories, Nucl. Phys. B 960, 115165 (2020), arXiv:2007.03374 [hep-lat] .
  * Schweitzer  _et al._ [2021] D. Schweitzer, S. Schlichting, and L. von Smekal, Critical dynamics of relativistic diffusion, (2021), arXiv:2110.01696 [hep-lat] .
  * Schaefer and Skokov [2022] T. Schaefer and V. Skokov, Dynamics of non-Gaussian fluctuations in model A, Phys. Rev. D 106, 014006 (2022), arXiv:2204.02433 [nucl-th] .
  * Chattopadhyay  _et al._ [2023] C. Chattopadhyay, J. Ott, T. Schaefer, and V. Skokov, Dynamic scaling of order parameter fluctuations in model B, (2023), arXiv:2304.07279 [nucl-th] .
  * Son and Stephanov [2004] D. T. Son and M. A. Stephanov, Dynamic universality class of the QCD critical point, Phys. Rev. D70, 056001 (2004), arXiv:hep-ph/0401052 [hep-ph] .
  * Canet  _et al._ [2003] L. Canet, B. Delamotte, D. Mouhanna, and J. Vidal, Nonperturbative renormalization group approach to the Ising model: A Derivative expansion at order partial**4, Phys. Rev. B 68, 064421 (2003), arXiv:hep-th/0302227 .
  * Benitez  _et al._ [2009] F. Benitez, J. P. Blaizot, H. Chate, B. Delamotte, R. Mendez-Galain, and N. Wschebor, Solutions of renormalization group flow equations with full momentum dependence, Phys. Rev. E 80, 030103 (2009), arXiv:0901.0128 [cond-mat.stat-mech] .
  * De Polsi  _et al._ [2021] G. De Polsi, G. Hernández-Chifflet, and N. Wschebor, Precision calculation of universal amplitude ratios in O(N) universality classes: Derivative expansion results at order O($\partial$4), Phys. Rev. E 104, 064101 (2021), arXiv:2109.14731 [cond-mat.stat-mech] .
  * Roth and von Smekal [2023] J. V. Roth and L. von Smekal, Critical dynamics in a real-time formulation of the functional renormalization group, (2023), arXiv:2303.11817 [hep-ph] .
  * Floerchinger [2012] S. Floerchinger, Analytic Continuation of Functional Renormalization Group Equations, JHEP 05, 021, arXiv:1112.4374 [hep-th] .
  * Floerchinger [2016] S. Floerchinger, Variational principle for theories with dissipation from analytic continuation, JHEP 09, 099, arXiv:1603.07148 [hep-th] .
  * Dupuis  _et al._ [2021] N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. M. Pawlowski, M. Tissier, and N. Wschebor, The nonperturbative functional renormalization group and its applications, Phys. Rept. 910, 1 (2021), arXiv:2006.04853 [cond-mat.stat-mech] .
  * Cao  _et al._ [2022] X. Cao, M. Baggioli, H. Liu, and D. Li, Pion dynamics in a soft-wall AdS-QCD model, JHEP 12, 113, arXiv:2210.09088 [hep-ph] .
  * Flory  _et al._ [2022] M. Flory, S. Grieninger, and S. Morales-Tejera, Critical and near-critical relaxation of holographic superfluids, (2022), arXiv:2209.09251 [hep-th] .
  * Damgaard [2011] P. H. Damgaard, Chiral Random Matrix Theory and Chiral Perturbation Theory, J. Phys. Conf. Ser. 287, 012004 (2011), arXiv:1102.1295 [hep-ph] .
  * Hasenbusch [2022] M. Hasenbusch, Three-dimensional $O(N)$-invariant $\phi^{4}$ models at criticality for $N\geq 4$, Phys. Rev. B 105, 054428 (2022), arXiv:2112.03783 [hep-lat] .
  * Hasenfratz and Leutwyler [1990] P. Hasenfratz and H. Leutwyler, Goldstone Boson Related Finite Size Effects in Field Theory and Critical Phenomena With O($N$) Symmetry, Nucl. Phys. B 343, 241 (1990).
  * Armas  _et al._ [2021] J. Armas, A. Jain, and R. Lier, Approximate symmetries, pseudo-Goldstones, and the second law of thermodynamics, (2021), arXiv:2112.14373 [hep-th] .
  * Delacrétaz  _et al._ [2022] L. V. Delacrétaz, B. Goutéraux, and V. Ziogas, Damping of Pseudo-Goldstone Fields, Phys. Rev. Lett. 128, 141601 (2022), arXiv:2111.13459 [hep-th] .
  * Ammon  _et al._ [2019] M. Ammon, M. Baggioli, and A. Jiménez-Alba, A Unified Description of Translational Symmetry Breaking in Holography, JHEP 09, 124, arXiv:1904.05785 [hep-th] .
  * Amoretti  _et al._ [2019] A. Amoretti, D. Areán, B. Goutéraux, and D. Musso, Universal relaxation in a holographic metallic density wave phase, Phys. Rev. Lett. 123, 211602 (2019), arXiv:1812.08118 [hep-th] .
  * Donos  _et al._ [2021] A. Donos, P. Kailidis, and C. Pantelidou, Dissipation in holographic superfluids, JHEP 09, 134, arXiv:2107.03680 [hep-th] .
  * Ammon  _et al._ [2022] M. Ammon, D. Arean, M. Baggioli, S. Gray, and S. Grieninger, Pseudo-spontaneous $U(1)$ symmetry breaking in hydrodynamics and holography, JHEP 03, 015, arXiv:2111.10305 [hep-th] .
  * Baggioli and Goutéraux [2023] M. Baggioli and B. Goutéraux, Colloquium: Hydrodynamics and holography of charge density wave phases, Rev. Mod. Phys. 95, 011001 (2023), arXiv:2203.03298 [hep-th] .
  * Yao and Täuber [2022] L. H. Yao and U. C. Täuber, Critical dynamics of the antiferromagnetic O(3) nonlinear sigma model with conserved magnetization, Phys. Rev. E 105, 064128 (2022), arXiv:2204.11145 [cond-mat.stat-mech] .
  * Kittel [1961] C. Kittel,  _Elementary Statistical Physics_ (Wiley, 1961).
  * Arnold [2000a] P. B. Arnold, Symmetric path integrals for stochastic equations with multiplicative noise, Phys. Rev. E 61, 6099 (2000a), arXiv:hep-ph/9912209 .
  * Arnold [2000b] P. B. Arnold, Langevin equations with multiplicative noise: Resolution of time discretization ambiguities for equilibrium systems, Phys. Rev. E 61, 6091 (2000b), arXiv:hep-ph/9912208 .
  * Engels  _et al._ [2003] J. Engels, L. Fromme, and M. Seniuch, Correlation lengths and scaling functions in the three-dimensional O(4) model, Nucl. Phys. B 675, 533 (2003), arXiv:hep-lat/0307032 .
  * Dimitrovic  _et al._ [1991] I. Dimitrovic, P. Hasenfratz, J. Nager, and F. Niedermayer, Finite size effects, Goldstone bosons and critical exponents in the d = 3 Heisenberg model, Nucl. Phys. B 350, 893 (1991).