††institutetext: Institut für Kernphysik, Technische Universität Darmstadt
Schlossgartenstraße 2, D-64289 Darmstadt, Germany

# The Sphaleron Rate from 4D Euclidean Lattices

Marc Barroso Mancha, Guy D. Moore<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We develop a new method to determine thermal activation rates, such as for
bubble nucleation, topology change, etc., using 4-dimensional Euclidean
methods. This allows nonperturbative study on the lattice. We then investigate
the strong sphaleron rate in pure-glue QCD at temperatures between 1.3
$T_{\mathrm{c}}$ and 1000 $T_{\mathrm{c}}$, making contact with previous
results but extending them down close to the critical temperature. The
extension to full QCD will be straightforward. Limitations of the proposal
(the inability to compute a certain dynamical prefactor, puzzling large-volume
behavior, and the inability to treat temperatures $T<1.3T_{\mathrm{c}}$) are
also discussed.

###### Keywords:

QCD, Quark-Gluon Plasma, topology, sphaleron, tunneling, separatrix

##  1 Introduction

Nonabelian gauge theories, including QCD and the $SU_{\mathrm{w}}(2)$ weak
sector of the Standard Model, feature gauge groups which can have
topologically nontrivial configurations in 4 spacetime dimensions
Belavin:1975fg ; tHooft:1976rip ; tHooft:1976snw . Specifically, the spacetime
integral of the contraction of the field strength tensor111Our conventions
are: we use natural units, $\hbar=1$ and $c=1$. The Euclidean metric is
positive and the Minkowski metric is mostly-positive
$g_{\mu\nu}=\mathrm{Diag}[-1,+1,+1,+1]$; the QCD/electroweak gauge fields are
$G_{\mu}=G_{\mu}^{a}T^{a}$ and $W_{\mu}=W_{\mu}^{a}T^{a}$ with $T^{a}$ the
Hermitian fundamental representation Lie algebra generators, normalized to
$\,\mathrm{Tr}\,T^{a}T^{b}=\frac{\delta_{ab}}{2}$; the fields take their
geometrical normalization such that the gauge coupling does not appear in the
covariant derivative $D_{\mu}=\partial_{\mu}-iG_{\mu}$ but as a denominator in
the gauge action
$\mathcal{L}=\ldots-\,\mathrm{Tr}\,G_{\mu\nu}G^{\mu\nu}/2g^{2}$.  $G_{\mu\nu}$
with the dual field strength tensor
$\tilde{G}^{\mu\nu}\equiv\epsilon^{\mu\nu\alpha\beta}G_{\alpha\beta}/2$,
satisfies:

$q(x)\equiv\frac{\,\mathrm{Tr}\,G_{\mu\nu}\tilde{G}^{\mu\nu}}{16\pi^{2}}\,,\qquad
Q_{\mathrm{I}}=\int\mathrm{d}^{4}x\;q(x)\;\in\mathbb{Z}\,.$ (1)

That is, the integral of $q(x)$ the “topological charge density” always
returns an integer, provided the integral is carried out over a compact domain
without boundaries or with boundaries where the field strength goes to zero.
(For the case of the $SU_{\mathrm{w}}(2)$ sector, replace $G_{\mu\nu}\to
W_{\mu\nu}$.)

These topological structures are important because they lead to anomalies in
axial quark number currents: according to the Adler-Bell-Jackiw anomaly
Adler:1969gk ; Bell:1969ts , each fundamental-representation fermionic species
suffers a nonconservation of its axial charge in the amount

$\partial_{\mu}J^{\mu}_{5}(x)=\ldots+2q(x)\,,\qquad\Delta
Q_{5}=2Q_{\mathrm{I}}\,.$ (2)

Here $J^{\mu}_{5}$ is the Nöther current associated with axial quark number
and $Q_{5}$ is the associated axial number. In the electroweak sector the
effect is dramatic: because only left-handed species couple to
$SU_{\mathrm{w}}(2)$, both baryon number and lepton number are violated,
tHooft:1976snw :
$\partial_{\mu}J^{\mu}_{B}=N_{g}\frac{\,\mathrm{Tr}\,W_{\mu\nu}\tilde{W}^{\mu\nu}}{16\pi^{2}}$
with $N_{g}=3$ the number of generations. Topological transitions in real time
can be stimulated if the fields can go over the so-called sphaleron barrier
Klinkhamer:1984di , and this process becomes efficient at high temperatures
Kuzmin:1985mm  and may explain the baryon number of the Universe (see for
instance Rubakov:1996vz ; Cline:2006ts  and references therein).

In Quantum Chromodynamics topological transitions in real time are also
relevant for interesting high-temperature dynamics. They are responsible for
the equilibration between right- and left-handed quark number in electroweak
baryogenesis scenarios McLerran:1990de ; Giudice:1993bb . In the context of
relativistic heavy ion physics, noncentral collisions contain intense magnetic
fields, which can combine with axial quark number imbalances to provide
interesting new physics signals – the so-called chiral magnetic effect
Fukushima:2008xe . This proposal has inspired a great deal of work, including
a deeper appreciation of the hydrodynamics of systems with axial quark number
Son:2009tf ; Stephanov:2012ki .

Such chiral dynamics in QCD are very dependent on the decay rate of chiral
imbalances. Close to equilibrium, a net axial charge abundance is described by
a chemical potential for axial quark number $\mu_{A}$; at leading order and in
each quark species, the axial number density $n_{A}$ is related to $\mu_{A}$
through Giudice:1993bb

$n_{A}=\mu_{A}T^{2}\,.$ (3)

The coefficient is determined by the quark number susceptibility, which is
known Bazavov:2013uja  to be within $10\%$ of the leading-order value already
at $T=2\,T_{\mathrm{pc}}$ with $T_{\mathrm{pc}}$ the pseudocritical
temperature, $T_{\mathrm{pc}}\simeq 155$ MeV Borsanyi:2013bia ; HotQCD:2014kol
. The axial chemical potential biases thermal transitions and thereby acts to
create a net $G_{\mu\nu}\tilde{G}^{\mu\nu}$ value which drives down the axial
number. The Minkowski time $t_{\mathrm{m}}$ derivative of $n_{A}$ due to this
thermal relaxation effect is:

$\frac{\mathrm{d}n_{A}}{\mathrm{d}t_{\mathrm{m}}}=-\frac{\sum_{f}\mu_{A,f}}{T}\Gamma_{\mathrm{sphal},s}=-n_{A}\frac{n_{f}\Gamma_{\mathrm{sphal},s}}{T^{3}}\,,$
(4)

where $n_{f}$ is the number of light fermionic species and
$\Gamma_{\mathrm{sphal},s}$, the so-called strong sphaleron rate, describes
the equilibrium diffusion rate of topology due to random, real-time
topological processes at temperature $T$:

$\Gamma_{\mathrm{sphal},s}\equiv\frac{1}{Vt_{\mathrm{max}}}\left\langle\left(\int_{0}^{t_{\mathrm{max}}}\mathrm{d}t_{\mathrm{m}}\int_{V}\mathrm{d}^{3}x\>q(x,t_{\mathrm{m}})\right)^{2}\right\rangle=\int\mathrm{d}^{4}x\langle
q(x)q(0)\rangle\,.$ (5)

That is, $\Gamma_{\mathrm{sphal},s}$ is the mean-squared topology per unit
Minkowski 4-volume.222Some older sources define the sphaleron rate such that
$dn_{A}/dt_{\mathrm{m}}=-\sum_{f,s}\mu_{A,fs}\Gamma/T$ where the sum also runs
over spins $s$ – that is, the sum is over all species which would be created
by a sphaleron. This definition is half as large as the one we use, that is,
$\Gamma=\Gamma_{\mathrm{sphal},s}/2$.  The definition is similar to the
topological susceptibility $\chi_{\mathrm{top}}$, except that
$\chi_{\mathrm{top}}$ refers to correlators in Euclidean spacetime, while the
time appearing in Eq. (5) is real Minkowski time. As we will discuss below,
these two rates have different physical interpretations and relevance and at
high temperatures they strongly differ from each other.

In order to understand the dynamics of chiral imbalances, it is imperative to
have at least a decent determination of the sphaleron rate
$\Gamma_{\mathrm{sphal},s}$. If it is large, chiral imbalances relax quickly
and are unlikely to lead to interesting dynamics. On the other hand, if it is
small, it may be a good approximation to treat axial quark number as
approximately conserved, at least for light species.

Unfortunately, the state of the art in determining $\Gamma_{\mathrm{sphal},s}$
is not very advanced. The equivalent $SU_{\mathrm{w}}(2)$ rate is relatively
well understood. The relevant dynamics involves soft (wave number small
compared to $T$) gauge fields which should behave like classical fields
Grigoriev:1988bd , though the relevant dynamics are subtle Bodeker:1995pp ;
Arnold:1996dy ; Bodeker:1999gx . Building on early work by Ambjørn et al
Ambjorn:1990pu , a fairly complete picture for the sphaleron rate in the
$SU_{\mathrm{w}}(2)$ sector emerged by the end of the 1990s Ambjorn:1997jz ;
Moore:1997sn ; Arnold:1997yb ; Bodeker:1998hm ; Moore:1998zk ; Bodeker:1999gx
; Arnold:1999uy . This success is based on the fact that the
$SU_{\mathrm{w}}(2)$ sector of the Standard Model has a rather small gauge
coupling, $g^{2}N/4\pi\simeq 1/15$, which ensures that the relevant dynamics
is quite infrared and the classical-field approximation is well under control.
In comparison, even at electroweak temperatures $T\sim 1000\,T_{\mathrm{pc}}$,
the strong coupling is of order $g^{2}N/4\pi\sim 0.3$, and the classical-field
approximation is not very reliable. There are only two attempts to determine
the strong sphaleron rate nonperturbatively using lattice techniques
Moore:1997im ; Moore:2010jd , and the latter concludes that the results are no
better than factor-of-2 estimates even at $T=100\,\mathrm{GeV}\sim
700\,T_{\mathrm{pc}}$. The rate below this temperature is not even estimated.

We need a new approach to compute the strong sphaleron rate – one which does
not rely on an approximation of QCD dynamics in terms of 3+1 dimensional
classical fields. The approach must be nonperturbative, because the
temperatures achieved in the early stages of heavy ion collisions are at most
a few times higher than $T_{\mathrm{pc}}$, where the theory is fully
nonperturbative. Therefore the approach must involve our only rigorous tool
for understanding QCD in this regime, which is 4-dimensional Euclidean lattice
gauge theory.

One approach proposed recently Altenkort:2020axj ; Altenkort:2021jbk  is to
measure $q(x)\,q(0)$ correlators as a function of Euclidean time and to
attempt an analytic continuation to obtain the real-time spectral function,
whose low-frequency limit determines $\Gamma_{\mathrm{sphal},s}$. This is an
interesting approach, but it is very challenging, particularly because of our
ignorance of the expected low-frequency behavior of the spectral function and
the challenge in reconstructing nontrivial low-frequency structures in such
analytical continuations.

We will provide a completely new approach for estimating
$\Gamma_{\mathrm{sphal},s}$ from Euclidean lattice path-integral data. Our
approach is based on the picture of thermal activation developed by Langer
Langer:1969bc  and further developed by Affleck Affleck:1980ac  and Linde
Linde:1980tt ; Linde:1981zj . The next section (Sec. 2) will review this
picture. Thermal activation is controlled by Euclidean configurations which
are balanced at a saddlepoint between two (topological) minima, with a rate
determined by the prevalence of saddlepoint configurations times the real-time
rate that the fields move through such saddles. In Subsection 2.2 we show that
this rate is related to the size of fluctuations in the Euclidean direction,
such that different Euclidean-time slices lie on different sides of the
saddle. This is something which can be measured purely within Euclidean
simulations. We formulate how to do this specifically for the QCD sphaleron
rate (that is, $SU(3)$ gauge theory with or without fermions) in Section 3,
and show how to really carry out such calculations in Section 4.

With the method in hand, we make the usual checks for lattice spacing, lattice
volume, and gradient-flow depth dependence, within pure-glue QCD ($SU(3)$
Yang-Mills theory) in Section 5. Having passed all tests, we explore what
range of temperatures we can investigate in Section 6, finding that the method
breaks down below about $T=1.3\,T_{\mathrm{c}}$. Our final results for
temperatures above this threshold are presented in Section 7, followed by a
discussion and outlook.

##  2 Thermal activation: a Review

The theory of thermal activation was pioneered by Langer Langer:1969bc  and
put into field-theoretical language by Affleck Affleck:1980ac  and by Linde
Linde:1981zj . The canonical example is the problem of bubble nucleation in a
system with a stable and a metastable vacuum. At low or zero temperatures,
this is controlled by a 4-dimensional bubble configuration, as argued by
Coleman and Callan Coleman:1977py ; Callan:1977pt . At higher temperatures
(roughly, when $1/T$ is of order or smaller than the diameter of the would-be
4D bubble), the relevant solution is instead a bubble which is constant in the
(periodic, $\beta=1/T$-extent, Euclidean) time direction and varies as a
function of spatial radius. The field theoretical justification is worked out
by Affleck Affleck:1980ac  and by Linde Linde:1980tt ; Linde:1981zj .

However, for our purposes this example is confusing and involves extraneous
details. Instead, we will outline the key ideas, and motivate the approach we
will take, with a much simpler toy example, proposed by Arnold and McLerran.

###  2.1 Simple example: pendulum under gravity

Consider a rigid pendulum in a gravitational potential. There is one degree of
freedom $\varphi$, and the Lagrangian is

$\mathcal{L}(\varphi,\dot{\varphi})=\frac{m}{2}\dot{\varphi}^{2}-V_{0}\left(1-\cos\varphi\right)\,.$
(6)

Arnold and McLerran used this problem to illustrate the difference between
instantons and sphalerons and the physical role each plays Arnold:1987zg , and
we are guided by their considerations.

First, note that the energy spacing between the low-lying quantum excitations
for this system is $\omega_{0}=\sqrt{V_{0}/m}$, whereas the energy required to
spin the pendulum over the top of its potential is $2V_{0}$. We will consider
the case $2V_{0}\gg\omega_{0}$ and will consider temperatures $T<\omega_{0}$,
where the vacuum dominates, and $2V_{0}\gg T\gg\omega_{0}$, where a range of
thermal states participate but states with $E>2V_{0}$ are still rare. To make
the thermal ensemble make sense and the real-time behavior be more similar to
a hot, chaotic system, we will couple $\varphi$ to a thermal bath of
oscillators. We assume that this bath causes a negligible modification of the
above Lagrangian, but that it leads to noisy, dissipative real-time dynamics.
We will not specify these couplings in detail, since this is only an
illustrative example for us and not our main item of study.

Let us understand the distinction between instantons and sphalerons, and the
role of each. First the instanton. The classical Lagrangian does not actually
specify our problem completely, because the values $\varphi=0$ and
$\varphi=2\pi$ (or any pair of values separated by $2\pi$) are physically
equivalent, and the wave function must be periodic up to a phase,
$\Psi(\varphi+2\pi)=\Psi(\varphi)e^{i\theta}$. Here $\theta$ is a parameter,
analogous to the $\Theta$-parameter of QCD, which only plays a role in the
quantum theory. Instantons are configurations which answer the question:

> Instantons: how strongly does the free energy depend on $\theta$?

We can address the $\theta$-dependence of the free energy by considering its
calculation with the Euclidean path integral. Naively we expect in the
Euclidean path integral that $\varphi$ is periodic,
$\varphi(t=0)=\varphi(t=\beta)$, but in fact we only require periodicity
modulo $2\pi$. However, configurations with
$\varphi(t=\beta)=\varphi(t=0)+2\pi n$ should be weighted with a factor
$e^{in\theta}$, such that the thermal partition function is:

$Z(T)=\sum_{n\in\mathcal{Z}}e^{in\theta}\int_{-\pi}^{\pi}\mathrm{d}\varphi_{0}\langle\varphi_{0}|e^{-\beta
H}|\varphi_{0}+2\pi n\rangle=\sum_{n}e^{in\theta}Z_{n}(T)\,.$ (7)

The $\varphi_{0}$ integral is a trace over all states, with the $n$-sum
reflecting the fact that these states are $2\pi$-periodic but allowing us to
write a path integral in terms of states which do not have this periodicity.
We name the configurations contributing to $Z_{n}$ ($n\neq 0$) $n$-instantons.

At low temperature $Z(T)=e^{-\beta E_{0}}$ is dominated by the vacuum energy
and the instantons determine how much the vacuum-energy depends on the
boundary condition. For $\theta=0$ all $Z_{n}$ add and $E_{0}$ is minimized,
as we expect for a periodic boundary, whereas for $\theta=\pi$ $E_{0}$ takes
its maximal value because the vacuum wave function is forced to have a node at
$\varphi=\pi$. At higher temperature the behavior is more complex. Even-parity
(at $\varphi=0$) states have the lowest energy for $\theta=0$, but odd-parity
states have the highest energy for $\theta=0$ and the lowest energy for
$\theta=\pi$. At temperatures well above $T=\omega_{0}$, the free energy is
determined by a mixture of these states and the effects of $\theta$ tend to
cancel out. In the Euclidean path integral, this happens because the
characteristic $t$-width (Euclidean time width) over which $\varphi$ varies in
an instanton is $t\sim 1/\omega_{0}$. When $\beta\omega<1$, the $\varphi$
value must change very rapidly in order to complete an instanton, such that
$S_{\mathrm{inst}}\geq\int_{0}^{\beta}\mathrm{d}t\>m(\partial_{t}\varphi)^{2}/2\geq
2\pi^{2}mT$. The instanton action rises linearly at large $T$, rendering
$\theta$ irrelevant. This point is illustrated in Figure 1.

Figure 1:  Instantons as a function of temperature. At low temperatures,
$\varphi$ transitions from $0$ to $2\pi$ along a smooth characteristic curve.
With increasing $T$, as the range of $t$-periodicity is narrowed, the
instanton must “hurry” between $\varphi=0$ and $\varphi=2\pi$, increasing the
associated action.

Now consider the same system, but distinguishing between $\varphi=0$ and
$\varphi=2\pi$ as physically distinct configurations. In this case, we can
consider a setup where the initial conditions are close to thermal, but with
$\varphi$ always close to the minimum at 0, and ask the completely different
question:

> Sphalerons: How fast does $\varphi$ diffuse over the range of available
> $\varphi=2\pi n$ minima over very long periods of Minkowski time?

We can define this rigorously in terms of the long-time diffusion of
$\varphi$’s value, distinguishing between $2\pi$-copies, as
$\Gamma=\lim_{t_{\mathrm{m}}\to\infty}\frac{1}{(2\pi)^{2}t_{\mathrm{m}}}\langle(\varphi(t_{\mathrm{m}})-\varphi(0))^{2}\rangle$.
This is controlled by transitions over $\varphi=\pi$ and its periodic copies
because the potential has a maximum there and this represents the principal
barrier to $\varphi$ evolution.333The definition of $\Gamma$ described in
footnote 2 is the rate of forward transitions, $\varphi<\pi$ transitioning to
$\varphi>\pi$. The analogous rate is the relevant one for the bubble
nucleation problem and is therefore the one considered by Affleck
Affleck:1980ac . Again, the rate we consider is for transitions in both
directions and is precisely twice as large. This factor emerges when computing
$\langle|\mathrm{d}\varphi/\mathrm{d}t_{\mathrm{m}}|\rangle$, which is the
mean rate of crossing in either direction. To compute the mean rate of forward
crossings only, one should average this mean field velocity only over the
forward-moving cases.

Intuitively, the sphaleron rate is the product of three factors (see figure
2):

  1. 1.

How likely is $\varphi$ to take a value near the top of the barrier? That is,
what is $\mathrm{d}P(\varphi)/\mathrm{d}\varphi$ at $\varphi=\pi$, where
$\mathrm{d}P(\varphi)/\mathrm{d}\varphi$ is the probability density to find
$\varphi$ in some narrow range about a specific value. In the toy model we
expect $\mathrm{d}P/\mathrm{d}\varphi\sim\exp(-2V_{0}/T)$.

  2. 2.

How fast does $\varphi$ change values when it is at the top of the barrier?
For our problem this is
$\langle|\mathrm{d}\varphi/\mathrm{d}t_{\mathrm{m}}|\rangle$. Note that a
rescaling of the variable $\varphi$ changes both 1 and 2, but the change
cancels such that the product of the two terms is rescaling invariant. In the
toy model we expect $m\dot{\varphi}^{2}\sim T$, so
$\langle|\mathrm{d}\varphi/\mathrm{d}t_{\mathrm{m}}|\rangle\sim\sqrt{T/m}$.

  3. 3.

How correlated are successive crossings of $\varphi=\pi$? This is a pure
number, usually close to 1, which we will call $\mathbf{d}$.

The total transition rate is then

$\Gamma=\left.\frac{\mathrm{d}P(\varphi)}{\mathrm{d}\varphi}\right|_{\varphi=\pi}\times\left\langle\left|\frac{\mathrm{d}\varphi}{\mathrm{d}t_{\mathrm{m}}}\right|\right\rangle\times\mathbf{d}\,.$
(8)

Figure 2:  Cartoon of thermal activation in the pendulum toy model. Left: when
one basin starts filled and one starts empty, the rate to move across is the
flux across the top, which is the probability to be near the top times the
rate of crossing. Right: for typical damping, an evolution bounces in one
basin, gains enough energy to cross, and then bounces in the next basin.

Figure 3:  Cartoon illustrating the dynamical prefactor $\mathbf{d}$. For
extremely strong damping (left), the trajectory can cross the barrier several
times in succession, leading to only one (or zero) true crossings. For
underdamping (right), it can cross several barriers in succession.

The first two points are intuitively clear, but are explained in proper detail
in the seminal paper of Affleck. The last point appears to have been first
noticed by Arnold and McLerran Arnold:1987mh , and is illustrated in Figure 3.
By very weak damping, the system can cross several barriers in succession.
Crossings are positively correlated and the number of individual crossings
underestimates the square of the change in $\varphi/2\pi$. In fact, in the
absence of damping as considered by Arnold and McLerran Arnold:1987zg ,
$\mathbf{d}$ is infinite and the mean-squared $\varphi$ value rises
quadratically, rather than linearly, with time. But for very strong
overdamping, several crossings may occur in a single “macroscopic” crossing
event – as appears to occur for the electroweak sphaleron rate Arnold:1996dy ;
Bodeker:1999gx . In this case, the true rate is smaller than the count of
individual crossings. In a broad intermediate range of damping strengths, the
$\dot{\varphi}$ direction does not change as it crosses the top of the barrier
(say, the region where $V_{\mathrm{max}}-V<T$), but there is enough damping by
the time it reaches the next barrier that it cannot cross. In this case we
expect $\mathbf{d}\simeq 1$.

We should emphasize that the technique we develop here cannot determine
$\mathbf{d}$. Instead, we will concentrate on determining the product of the
first two factors, and will assume that $\mathbf{d}$ is close to 1. This is a
weakness of our approach, and it is not clear to us how it can be overcome.

In a more general context, the exact location of and distance from the barrier
may be harder to quantify. Therefore, a successful thermal activation
calculation typically requires a few things:

  1. 1.

Identify the “separatrix” which divides configurations into those closer to
one minimum and those closer to another.

  2. 2.

Identify some measure of distance $\varphi$ from the separatrix.

  3. 3.

Use a Euclidean path integral to find the probability to be within some narrow
tolerance $|\varphi|<\Delta$ of the separatrix.

  4. 4.

Determine the mean velocity of the $\varphi$ variable when it is at or very
near the separatrix.

We will show how to combine the last 3 steps into a single step, which does
not require determining the value of $\varphi$.

###  2.2 Crossing speed and Euclidean fluctuations

How do we estimate $\langle|\mathrm{d}\varphi/\mathrm{d}t|\rangle$ (or its
equivalent in more complex problems) using only Euclidean data? The answer
lies in the size of the fluctuations in $\varphi(t)$ as a function of the
Euclidean time $t$. The amplitude of these fluctuations is proportional to the
rate at which $\varphi$ changes in real time. After all, the real-time
evolution can be determined by analytical continuation from the Euclidean-time
dependence, as calculated in detail by Arnold and McLerran Arnold:1987zg .
Consider again our toy example, with Lagrange function given in Eq. (6). The
Hamiltonian is

$H=\frac{m}{2}\dot{\varphi}^{2}+V_{0}(1-\cos(\varphi))$ (9)

and classically, the thermal distribution of $\dot{\varphi}$ values is given
by equipartition,

$\displaystyle P(\dot{\varphi})$
$\displaystyle\propto\exp(-m\dot{\varphi}^{2}/2T)$
$\displaystyle\left\langle\left|\frac{\mathrm{d}\varphi}{\mathrm{d}t_{\mathrm{m}}}\right|\right\rangle$
$\displaystyle=\frac{\int_{0}^{\infty}\mathrm{d}\dot{\varphi}\;\dot{\varphi}e^{-m\dot{\varphi}^{2}/2T}}{\int_{0}^{\infty}\mathrm{d}\dot{\varphi}\;e^{-m\dot{\varphi}^{2}/2T}}=$
(10) $\displaystyle=\sqrt{\frac{2T}{\pi m}}\,.$

For comparison, in Euclidean space, the dominant $t$-dependence in
$\varphi(t)$ arises from fluctuations with Matsubara frequency $\omega=2\pi
T$:

$\varphi(t)\simeq\varphi_{0}+c\,\cos(2\pi tT)+s\,\sin(2\pi tT)$ (11)

with $\varphi_{0}$ the average of $\varphi(t)$ over $t$ values and $c,s$
coefficients indicating two independent $t$-dependent fluctuations with this
frequency. The action associated with these fluctuations is approximately

$S(c,s)=\int_{0}^{\beta}\\!\\!\\!\mathrm{d}t\;\frac{m}{2}(\partial_{t}\varphi)^{2}=\int_{0}^{\beta}\\!\\!\mathrm{d}t\;2\pi^{2}mT^{2}\left(c^{2}\cos^{2}(2\pi
tT)+s^{2}\sin^{2}(2\pi tT)\right)=\pi^{2}mT(c^{2}+s^{2})\,.$ (12)

We can characterize $\varphi$’s $t$ dependence by looking at the RMS
difference between $\varphi(t=\beta/2)$ and $\varphi(t=0)$, which equals $2c$:

$\langle(\varphi(\beta/2)-\varphi(0))^{2}\rangle=\langle(2c)^{2}\rangle=\frac{\int_{-\infty}^{\infty}\mathrm{d}c(2c)^{2}e^{-\pi^{2}mTc^{2}}}{\int_{-\infty}^{\infty}\mathrm{d}ce^{-\pi^{2}mTc^{2}}}=\frac{2}{\pi^{2}mT}\;\Rightarrow\;\sqrt{\langle(2c)^{2}\rangle}=\sqrt{\frac{2}{\pi^{2}mT}}$
(13)

which is smaller than
$\langle|\mathrm{d}\varphi/\mathrm{d}t_{\mathrm{m}}|\rangle$ by a factor of
$1/T\sqrt{\pi}$. So the side-to-side fluctuations in $\varphi$ are smaller
than $|\mathrm{d}\varphi/\mathrm{d}t|$ by a factor of $1/T$, which one could
guess on dimensional grounds, and an order-1 numerical factor.

Is this relation special? We claim it is quite general – the same Hamiltonian
controls real-time evolution and Euclidean-time fluctuations, since the
Hamiltonian $H$ equals the Euclidean-time Lagrangian
$\mathcal{L}_{\mathrm{E}}$.

To determine the numerical factor precisely, we have to decide precisely what
question we want to ask about the Euclidean $\varphi$ fluctuations. The point
$\varphi=\pi$ separates values closer to the minimum at 0 and the values
closer to the minimum at $2\pi$, so we follow standard usage and call it a
separatrix. In the following, it will be useful to ask how often
$\varphi(t=0)$ and $\varphi(t=\beta/2)$ lie on opposite sides of this
separatrix. In particular, we expect on the above reasoning that the first two
terms in Eq. (8) can be recast as:

$\left.\frac{\mathrm{d}P(\varphi)}{\mathrm{d}\varphi}\right|_{\varphi=\pi}\times\left\langle\left|\frac{\mathrm{d}\varphi}{\mathrm{d}t_{\mathrm{m}}}\right|\right\rangle=\frac{\mbox{order-1
coeff}}{T}\times
2\,\mathrm{Prob}\Big{(}\varphi(t=0)<\pi<\varphi(t=\beta/2)\Big{)}\,.$ (14)

Here $2\,\mathrm{Prob}(...)$ is the chance that $\varphi$ is on opposite sides
of $\pi$ at the Euclidean times $t=0$ and $t=\beta/2$; the factor of $2$
accounts for $\varphi(t=0)>\pi>\varphi(t=\beta/2)$.

To determine the order-1 coefficient, we must repeat the quick calculation of
$(2c)_{\mathrm{RMS}}$ above with more care. We do so in Appendix A, accounting
for all Matsubara modes, finding that the real-time rate is $2T$ times the
probability for $\varphi(t=0)$ and $\varphi(t=\beta/2)$ to be on opposite
sides of the separatrix.444Here is a quick way to see why. Compare Eq. (10) to
Eq. (13). The latter was computed with only the lowest Matsubara frequency,
but all odd frequencies contribute, multiplying it by a factor of
$\sqrt{\sum_{n=1,3,\ldots}n^{-2}=\pi^{2}/8}$. Also, Eq. (13) is based on an
RMS average, but we should compute the mean absolute value as used in Eq.
(10), which is smaller by $\sqrt{2/\pi}$ for Gaussian distributions. The two
effects correct Eq. (13) to $\sqrt{2/(\pi^{2}mT)\times\pi^{2}/8\times
2/\pi}=(1/2)\sqrt{2/\pi mT}$ which must be multiplied by $2T$ to get Eq. (10).
The appendix also shows how to include the effects of gradient flow, which we
will need later.

The advantage of this formulation is that it is only necessary to determine
which side of the “separatrix” $\varphi=\pi$ the field $\varphi(t)$ is on at
the two (Euclidean) time values $t=0$ and $t=\beta/2$. When we try to apply
these ideas to the sphaleron rate of QCD, we don’t know how to measure the
distance from the separatrix, but it is perfectly feasible to determine which
side of the separatrix one is on.

One might worry that, in higher-dimensional problems like quantum field
theory, the separatrix is generally a codimension-1 surface, rather than a
single point. Does this approach still work? We argue that it does. At any
point on the separatrix, one can find a local orthogonal direction and ask
about motion, and fluctuations, in that direction. At each such point, the
speed a real-time configuration moves through the separatrix is determined by
the same kinetic Hamiltonian as the Euclidean fluctuations orthogonal to the
separatrix. Therefore, their relative amplitudes are proportional. The real-
time Affleck-type calculation involves an integral over the separatrix surface
of a local probability density times the local
$\langle|\mathrm{d}\varphi/\mathrm{d}t|\rangle$, while our calculation
involves the same surface average of $\langle\sqrt{\Delta c^{2}}\rangle$,
which is proportional with fixed proportionality constant $1/2T$. Therefore,
to the extent that configurations close to the separatrix really control the
transition rate, the two calculations will be equivalent when the factor
$1/2T$ is included.

##  3 Application to QCD

Let us see how to apply these ideas to the sphaleron rate in QCD. First, we
will consider one more time an evaluation in the toy model:

###  3.1 $\Gamma$ in the toy model

Before turning to QCD, suppose we were trying to determine $\Gamma$ in the toy
model, and we had a nonperturbative method for studying the finite-temperature
Euclidean path integral. However, suppose we somehow could not access the
exact value of $\varphi$, but only which side of $\pi$ it lies on. We could
perform our study as follows:

  1. 1.

Generate a sample of typical Euclidean configurations:

  2. 2.

For each, determine whether $\varphi$ is on the opposite side of $\varphi=\pi$
at $t=0$ as it is at $t=\beta/2$. This can be done by evolving $\varphi(t)$
under gradient-flow or relaxation dynamics and seeing whether it rolls down to
$\varphi=0$ or to $\varphi=2\pi$, for instance.

  3. 3.

Evaluate the sphaleron rate as $2T$ times the fraction of configurations which
make such a transition.

To understand the approach, consider Figure 4, which shows examples of how two
configurations would be evaluated in this procedure.

Figure 4: Cartoon of how two configurations would be evaluated using the
gradient-flow approach described in the text. A segment of a Euclidean
simulation, from $t=0$ to $t=\beta/2$, is cut out of the simulation, and the
starting/ending configurations are gradient-flowed to a zero-energy state. The
whole path (zero-energy to $t=0$ to $t=\beta/2$ to zero-energy) either has
$\varphi$ change by $\pm 2\pi$ (left) or by 0 (right); in the former case it
counts as a sphaleron.

The section from $t=0$ to $t=\beta/2$ is cut from a simulation, and one
“grafts” the gradient-flow path from the $t=0$ configuration to the vacuum and
from the $t=\beta/2$ configuration to the vacuum onto the two ends. The full
configuration starts and ends at the zero-energy solution, $\varphi=0$, modulo
$2\pi$. When the change is by $\pm 2\pi$, it is a sphaleron; otherwise it is
not.

###  3.2 QCD: Chern-Simons number

Our approach in QCD will be the same: to look for configurations which cross a
separatrix between $t=0$ and $t=\beta/2$, and therefore settle into distinct
vacua when acted on by gradient flow. There is no easily-measured value
$\varphi$ with a separatrix at $\varphi=\pi$ in QCD. However, we can use
gradient flow and topology to define a separatrix and to give a measure which
determines whether the $t=0$ slice and the $t=\beta/2$ slice are on opposite
sides.

Let’s start by defining gradient flow on a 3D configuration, in analogy to its
usual 4D definition Narayanan:2006rf ; Luscher:2010iy :

$B_{i}(t,\vec{r},{\tau_{\\!{}_{\mathrm{F}3}}}=0)=G_{i}(t,\vec{r})\,,\qquad\frac{\partial}{\partial{\tau_{\\!{}_{\mathrm{F}3}}}}B_{i}(t,\vec{r},{\tau_{\\!{}_{\mathrm{F}3}}})=-\frac{\delta\int\mathrm{d}^{3}x\,\mathrm{Tr}\,G_{jk}G_{jk}[B]}{\delta
B_{i}(\vec{r})}\,.$ (15)

The fields $B_{i}(t,\vec{r},{\tau_{\\!{}_{\mathrm{F}3}}})$, as a function of
the 3D coordinate $\vec{r}$ and the 3D flow-time coordinate
${\tau_{\\!{}_{\mathrm{F}3}}}$, can be considered a 4D gauge-field
configuration in the temporal gauge. One can then integrate the topological
charge density over this 4D configuration:

$\displaystyle Q(t)$
$\displaystyle\equiv-\int_{0}^{\infty}\mathrm{d}{\tau_{\\!{}_{\mathrm{F}3}}}\int\mathrm{d}^{3}\vec{r}\,\,q(t,{\tau_{\\!{}_{\mathrm{F}3}}},\vec{r})\,,$
$\displaystyle q(t,{\tau_{\\!{}_{\mathrm{F}3}}},\vec{r})$
$\displaystyle=\frac{\epsilon_{ijk}}{8\pi^{2}}\,\mathrm{Tr}\,G_{ij}\>\frac{\partial}{\partial{\tau_{\\!{}_{\mathrm{F}3}}}}B_{k}(t,\vec{r},{\tau_{\\!{}_{\mathrm{F}}}})\,.$
(16)

The result is equivalent to the definition of Chern-Simons number adopted in
Eq. (5) of Ref. Moore:1998swa , which predates papers discussing 4D gradient
flow Narayanan:2006rf ; Luscher:2010iy .

To understand $Q(t)$ better, consider an idealized 4D instanton solution
centered at the origin, and consider a series of time slices through the
instanton. For $t<0$ slices, gradient flow will make the 3D gauge field roll
down towards the vacuum described by the $t\to-\infty$ slice. This will yield
$Q(t)$ values which change from 0 at very negative $t$ towards $1/2$ for $t$
close to zero. But on any $t>0$ slice, the gauge field will gradient-flow down
towards the vacuum described by the $t\to+\infty$ slice. Therefore, $Q(t)$ as
defined above will have a sharp jump at $t=0$ (the middle slice of the
instanton), from a value around $+1/2$ to a value around $-1/2$. The slice
where this jump occurs is our separatrix, and the property of the separatrix
is that such a 3D slice, under gradient flow, moves to a saddlepoint between
3D vacua and sticks there.

Now return to the periodic (thermal) Euclidean path integral of time extent
$\beta$. We want to know whether the $t=0$ and $t=\beta/2$ slices are on
opposite sides of this separatrix. We could evaluate $Q(t)$ as defined above
at every time slice from $t=0$ to $t=\beta/2$ and see if it makes a sharp
jump. But this is inefficient. If we only want to know whether the $t=0$ and
$t=\beta/2$ slices are on opposite sides, we can consider the following
combination:

$\displaystyle Q_{\mathrm{S}}$
$\displaystyle=Q(t=0)+\left[\int_{0}^{\beta/2}\mathrm{d}t\int\mathrm{d}^{3}x\>q(x,t)\right]-Q(t=\beta/2)$
(17) $\displaystyle\equiv Q_{0}+Q_{\mathrm{half}}-Q_{\beta/2}\,.$

Here $Q_{\mathrm{S}}$ stands for the “sphaleron topology,” not to be confused
with $\int_{0}^{\beta}\mathrm{d}t\int\mathrm{d}^{3}x\>q(x,t)=Q_{\mathrm{I}}$
the instanton topology of Eq. (1). $Q_{\mathrm{S}}$ is “almost” the difference
between $Q(t)$, defined above, at the two slices of interest. If we cross the
separatrix, we expect $Q_{0}$ and $Q_{\beta/2}$ to be approximately $\pm 1/2$
and $\mp 1/2$ so the difference is $\pm 1$. But by adding the $q$-integral
over the intervening segment $Q_{\mathrm{half}}$, we make the measurement
actually topological. Namely, $Q_{\mathrm{S}}$ is the integral of $q(x,t)$
over a domain with no space boundaries (our box has periodic boundary
conditions), which begins on a vacuum configuration at
$t=0,{\tau_{\\!{}_{\mathrm{F}3}}}=\infty$ and ends on a vacuum configuration
at $t=\beta/2,{\tau_{\\!{}_{\mathrm{F}3}}}=\infty$, and has no field
discontinuities in between. Therefore the integral of $q$, $Q_{\mathrm{S}}$,
is an integer. When this integer is nonzero, we define the configuration to be
a sphaleron. This idea is illustrated in Figure 5.

Figure 5: How we define and identify sphalerons in QCD. The topological
density is integrated between $t=0$ and $t=\beta/2$, and supplemented by its
integration as these two boundaries are gradient-flowed to the vacuum. When
the configuration crosses a separatrix, the result is $\pm 1$.

Lastly, note that the sphaleron rate is technically defined as the number of
topology transitions per unit time and volume. Therefore, in defining the
sphaleron rate, we must find the mean-squared value of $Q_{\mathrm{S}}$ from
Eq. (17), multiply by $2T$, and divide by the space volume which was studied –
that is, $\langle Q_{\mathrm{S}}^{2}\rangle$ is expected to be extensive in
the lattice volume, provided the lattice is large enough to see the large-
volume behavior.

##  4 Lattice calculation of the sphaleron rate

As discussed in the last section, we will evaluate whether a given lattice
gauge-field configuration is a sphaleron by evaluating $Q_{\mathrm{S}}$
defined in Eq. (17), determining

$\frac{\Gamma_{\mathrm{Eucl}}}{T^{3}}=\frac{\expectationvalue{Q_{\mathrm{S}}^{2}}}{\left(N_{s}/N_{t}\right)^{3}}\,,$
(18)

with $N_{s}$ and $N_{t}$ the number of lattice points in the space and time
directions respectively, and then rescaling by $2T$ to convert this into the
real-time rate as described in the appendix, specifically Eq. (39):

$\Gamma_{\mathrm{sphal},s}=2T\Gamma_{\mathrm{Eucl}}\,.$ (19)

Figure 6: Sketch of the calculation performed on the lattice to determine the
sphaleron rate. One dimension of the lattice is left not drawn for
convenience.

There is one additional complication which we must overcome to make the
procedure work. Integrating $q(\vec{r},t)$ on a lattice gauge configuration
typically returns a highly noisy result, because the lattice definition of
$q(\vec{r},t)$ is contaminated with high-dimension operators which are not
topological.555We use the clover definition of $q(\vec{r},t)$. Therefore it is
essential to apply some lattice gradient flow (Wilson flow Luscher:2010iy ) on
the 4D configuration before the evaluation. This 4D gradient flow also reduces
the size of the fluctuations in the temporal direction which we discussed in
subsection 2.2 and in appendix A. This will reduce how many configurations
cross the separatrix and are counted as sphalerons. Appendix A.2 contains a
calculation of the size of this effect, which we can use to correct for the
effects of gradient flow.

Just to summarize the procedure, the steps that we need to take are the
following:

  1. 1.

Generate a valid lattice QCD configuration through a Hybrid Monte-Carlo
algorithm.

  2. 2.

Apply a certain amount ${\tau_{\\!{}_{\mathrm{F}}}}$ of 4D gradient flow in
order to reduce the UV noise.

  3. 3.

Calculate the topology enclosed in half the lattice.

  4. 4.

Extract the 3D slices at $t=0$ and $t=\beta/2$, and compute $Q(0)$ and
$Q(\beta/2)$ according to 16, using a 3D version of gradient flow all the way
to the vacuum.

  5. 5.

Calculate $Q_{\mathrm{S}}$ with these quantities, see Eq. (17). Repeat for
many configurations to extract $\expectationvalue{Q_{\mathrm{S}}^{2}}$.

The procedure for determining $Q_{\mathrm{S}}$ is illustrated in figure 6.

In this paper, we have used openQCD openQCD  for generating the configurations
and applying the 4D gradient flow, and then implemented the calculation of 16
separately. We update with the HMC algorithm with trajectories of length 1.5
(in lattice units), measuring every 60 updates. The autocorrelation between
configurations is then acceptable for the $N_{s}$ values we use. For the
statistical analysis we have used the $\Gamma$-method algorithm as provided in
DEPALMA2018  as a Python module.

When integrating the 3D slices along the flow direction, we have also used
blocking-techniques (as were described in Moore:1998swa ) after the energy had
been lowered below a threshold (chosen such that the blocking did not change
the final answer). It is also important to note that while the integration of
the 4D gradient flow discretization needs to be done as precisely as possible
(and for that, we have used the 3rd order Runge-Kutta algorithm already
implemented in openQCD), the integration for the 3D slices can be done with a
much simpler Euler algorithm with bigger step-size, as the most important part
is that we end up in vacuum (where all links are close to unity), and not the
precise configuration we pass through in the middle.

For the scale setting, we have used the fitting parameters found in
Francis:2015lha ; Burnier:2017bod  which allows us to find the
$\beta_{\text{latt}}$ necessary to equilibrate our lattice to our desired
temperature, given the number of sites in the time direction $N_{t}$. For
reference, we have included in table 1 the value of the lattice coupling for
different cases we have used in this work.

|  | $N_{t}$
---|---|---
|  | 6 | 8 | 10 | 12 | 14 | 16
$T/T_{\mathrm{c}}$ | 1.3 | 6.04979 | 6.23700 | 6.39495 | 6.53054 | 6.64889 | 6.75371
2 | 6.33718 | 6.54976 | 6.72273 | 6.86803 | 6.99311 | 7.10283
3 | 6.64031 | 6.86803 | 7.04966 | 7.20049 | 7.32939 | 7.44187
4 | 6.86803 | 7.10283 | 7.28847 | 7.44187 | 7.57254 | 7.68632
5 | 7.04966 | 7.28847 | 7.47640 | 7.63126 | 7.76294 | 7.87746
7 | 7.32939 | 7.57254 | 7.76294 | 7.91939 | 8.05216 | 8.16748
10 | 7.63126 | 7.87746 | 8.06960 | 8.22716 | 8.36070 | 8.47657
20 | 8.22716 | 8.47657 | 8.67054 | 8.82926 | 8.96359 | 9.08003
Table 1:  Value of the lattice coupling $\beta_{\text{latt}}$ for different
values of $T/T_{\mathrm{c}}$ (where $T_{\mathrm{c}}$ is the critical
temperature of the pure-glue Yang-Mills theory, given by
$T_{\mathrm{c}}=287.4(70)$ MeV) and number of sites in the temporal direction
$N_{t}$. Figure 7: Numerical data obtained from 2560 configurations on a
lattice with $10\times 32^{3}$ sites, equilibrated at a temperature
$T=5\,T_{\mathrm{c}}$, where the measurement has been performed at a gradient
flow depth of ${\tau_{\\!{}_{\mathrm{F}}}}=1.2a^{2}$. The top-left plot is
$Q_{\mathrm{S}}$ of each configuration, and is approximately an integer. The
top-right, bottom-left, and bottom-right plots show the $Q_{\mathrm{half}}$,
$Q_{0}$, and $Q_{\beta/2}$, respectively.

Before trying to understand the continuum limit of this method, and to be sure
we are on the right track, we first test that the procedure really returns
integer values of $Q_{\mathrm{S}}$. We do this on a single lattice size and
spacing and a single chosen value of 4D gradient flow depth and show the
results in 7. We see that, while the individual components
$Q_{0},\,Q_{\mathrm{half}},\,Q_{\beta/2}$ are definitely not topological and
take a range of values, the combination $Q_{\mathrm{S}}$ is indeed topological
up to modest fluctuations caused by $\mathcal{O}(a^{2})$ high-dimension
operator contamination in our definition of topological density. This effect
gets smaller as we apply more 4D gradient flow (see below). In practice we
must apply enough 4D gradient flow to cleanly separate the different integer
values, which can be checked by making a histogram of $Q_{\mathrm{S}}$ values
and seeing that they separate into clear peaks. Then one assigns integer
values by projecting to the closest integer. We illustrate how this works in
figure 8, where we see the histogram of $Q_{\mathrm{S}}$ at different gradient
flow depths, and precisely check that the peaks get sharper as we go down the
flow path.

The quantitative criterion we use to determine whether enough 4D gradient flow
has been used is the following. First, we consider all configurations with
$|Q_{\mathrm{S}}|<0.5$ and we evaluate the standard deviation
$\sigma^{2}_{\text{data}}=\langle
Q_{\mathrm{S}}^{2}\rangle_{|Q_{\mathrm{S}}|<0.5}$. Assuming that the true
distribution is Gaussian, we ask how much of the distribution is expected to
lie outside the range $|Q_{\mathrm{S}}|<0.5$, and we require that this be less
than 2%:

$\int_{1/2}^{\infty}\frac{2\mathrm{d}{z}}{\sqrt{2\pi\sigma^{2}_{\text{data}}}}e^{-\frac{z^{2}}{2\sigma^{2}_{\text{data}}}}=\text{erfc}\left(\frac{1}{2\sqrt{2}\sigma_{\text{data}}}\right)<0.02\quad\to\quad\sigma_{\text{data}}<0.215\,.$
(20)

Figure 8: Histograms of $Q_{\mathrm{S}}$ for a lattice with $10\times 32^{3}$
sites, at temperature $T=5\,T_{\mathrm{c}}$, and different gradient flow
depths. The amount of flow corresponds to
${\tau_{\\!{}_{\mathrm{F}}}}=[0.2,0.4,0.6]\,a^{2}$, top to bottom. We see that
as we go down the flow path, the peaks are more defined. While for
${\tau_{\\!{}_{\mathrm{F}}}}/a^{2}=0.2$ we cannot see the topological nature
of the jumps, that becomes clear for ${\tau_{\\!{}_{\mathrm{F}}}}/a^{2}=0.4$.
As a remark, the criterion Eq. (20) fails for the top case, but is satisfied
for the other two.

##  5 Volume, lattice spacing, and gradient-flow depth dependence

Now that we have formulated the theoretical framework of the calculation and
shown that it works numerically, we need to analyze how it behaves on the
lattice. The three main topics we should study are:

  1. 1.

how to correct the rate for the effects of 4D gradient flow,

  2. 2.

how it changes when we vary the spatial volume of our lattices, and

  3. 3.

how to take the continuum limit, that is, how it changes when we go to finer
lattices (by varying the number of sites in the time direction, and varying
$\beta_{\text{latt}}$ such that the temperature is constant).

###  5.1 Gradient flow dependence

As explained in the previous section, some 4D gradient flow is needed to
ensure that $Q_{\mathrm{S}}$ is really topological up to fluctuations which
are small enough to fix by projecting to the nearest integer. However,
applying this gradient flow can move individual time-slices across the
separatrix and reduces how many configurations are determined to be
sphalerons. We have to correct for this effect. The analytical dependence in
the case of a single degree of freedom has already been explored in section
A.2, and here we will check that the same relationship holds for the case of
QCD on the lattice. By taking the quotient between equations Eq. (37) and Eq.
(44), we can find that

$\frac{\Gamma_{\mathrm{sphal},s}^{\tau_{\\!{}_{\mathrm{F}}}}}{\Gamma_{\mathrm{sphal},s}}=\frac{2}{\pi}\sqrt{2\sum_{\begin{subarray}{c}n=1\\\
n\text{
odd}\end{subarray}}^{\infty}\frac{e^{-{\overline{\tau}_{\\!{}_{\mathrm{F}}}}n^{2}}}{n^{2}}}$
(21)

where we have defined ${\overline{\tau}_{\\!{}_{\mathrm{F}}}}\equiv
8\pi^{2}{\tau_{\\!{}_{\mathrm{F}}}}/\left(a^{2}N_{t}^{2}\right)$.

Figure 9: Left: measurement of the sphaleron rate at different gradient-flow
depths, where the points in red (representing different values of
$\Gamma_{\mathrm{sphal},s}^{\tau_{\\!{}_{\mathrm{F}}}}$) have been used to
find the best fit for $\Gamma_{\mathrm{sphal},s}$ (see equation Eq. (21)). The
points in blue do not satisfy the criterion explained in a previous section,
equation Eq. (20). Right: same plot but where each value has been corrected
according to Eq. (21), and the horizontal line is the best fit described
before. We see that the correction successfully accounts for the effects of 4D
gradient flow.

We can then explore the accuracy of this correction by examining the same
ensemble of configurations using different amounts of gradient flow. The
result can be seen in figure 9, where the points in the left represent the
sphaleron rate $\Gamma_{\mathrm{sphal},s}^{\tau_{\\!{}_{\mathrm{F}}}}$, while
on the right they have been corrected according to Eq. (21). We see that, for
gradient flow depths above a threshold of about
${\tau_{\\!{}_{\mathrm{F}}}}=0.5a^{2}$, the correction accurately describes
the effect of gradient flow.

We expect the correction to break down when ${\tau_{\\!{}_{\mathrm{F}}}}$
comes of order $\beta^{2}$, that is, $a^{2}N_{t}^{2}$, when the gradient flow
starts erasing all information stored in the gauge fields. We have not
explored this regime in detail because we don’t need to – as long as there is
a range of ${\tau_{\\!{}_{\mathrm{F}}}}$ values where the correction works
well, then we can use this to establish the sphaleron rate in a
${\tau_{\\!{}_{\mathrm{F}}}}$-independent way. This should be possible
whenever $N_{t}$ is sufficiently large, which is the same as a requirement
that the lattice spacing be sufficiently small.

###  5.2 Volume dependence

The sphaleron rate has been computed at very high temperatures (or weak gauge
couplings) using classical-field, real-time methods, and it was always found
that the rate per unit 4-volume is suppressed in small boxes but grows to an
infinite-volume asymptotic value above some box size Ambjorn:1995xm ;
Moore:1999fs ; Moore:2010jd . Intuitively, this means that sphalerons have a
typical intrinsic size, and the box must be larger than this size for a
sphaleron to “fit.” At weak coupling, the size is expected to be
parametrically of order $1/g^{2}T$ Kuzmin:1985mm ; Arnold:1987mh . In the 4D
Euclidean calculation we expect something similar to happen: in lattices with
$N_{t}$ sites in the temporal direction and $N_{s}$ sites in the spatial
direction, for small $N_{s}$ we will not observe any jumps. As we make the
spatial volume bigger (keeping $\beta_{\text{latt}}$ and $N_{t}$ fixed), we
will start observing more and more jumps. Eventually, we would expect $\langle
Q_{\mathrm{S}}^{2}\rangle/V$ to approach a constant.

Figure 10: Sphaleron rate as a function of the aspect ratio $N_{s}/N_{t}$, at
fixed $N_{t}$ and temperature.

We can see the result for different values of $N_{s}$ in figure 10, where we
have used a lattice with $N_{t}=10$ and equilibrated it at
$T=10\,T_{\mathrm{c}}$. As expected, we see the increase in the sphaleron rate
as we increase the volume until it reaches the peak. But from that point
forward, we see a decrease in the rate. We did not expect this and it is not
analogous to what occurs in the real-time case. To show this more clearly, we
have revisited the volume dependence of the classical real-time $SU(2)$ rate
using the same code as in Ref. Moore:1999fs . Fixing the lattice spacing to
$a=0.25g^{2}T$ and varying the box size, figure 11 shows that the classical
real-time rate saturates to a flat asymptotic value to within a few percent,
without the clear post-peak decrease we observe in the Euclidean method.

Figure 11: Volume dependence of the $SU(2)$ sphaleron rate when the rate is
determined using classical-field, real-time techniques as in Ref. Moore:1999fs
. As the volume increases, the rate rises and saturates to a very clear
plateau.

We do not yet understand why our sphaleron-rate volume dependence exhibits a
peak and a decrease towards very large volumes. Note that the configurations
at larger volumes than the peak frequently exhibit $|Q_{\mathrm{S}}|>1$, which
means that there are multiple sphalerons in a box or that our single-
saddlepoint picture has broken down. Since it is unclear if our approach is
really appropriate in this regime, we will use the peak in the volume-
dependence curve as our sphaleron-rate estimate in the remainder of this work.
Clearly it would be valuable to revisit this decision and to find a clearer
understanding of the large-volume behavior and the meaning of this peak and
falloff. We will assume that our results suffer an $\mathcal{O}(30\%)$
systematic uncertainty because of this assumption.

###  5.3 Continuum limit

Now the only thing we need to check is how $\Gamma_{\mathrm{sphal},s}$
approaches the continuum limit $a\to 0$ at fixed temperature and aspect ratio.
We do this by performing the calculation at a series of lattice spacings (or
$N_{t}$ values) and seeing how they scale with $a^{2}$. As has been explained
in section 4, we use the scale setting as described in Francis:2015lha ;
Burnier:2017bod  in order to relate the lattice spacing $a$ to the lattice
coupling $\beta$. Then we can fix the temperature while increasing $N_{t}$ and
calculate the sphaleron rate at the peak in each case666We also checked that
the aspect ratio $N_{s}/N_{t}$ where the rate reaches its peak is independent
of lattice spacing., using the previously established method to correct for 4D
gradient-flow effects.

Figure 12: Sphaleron rate calculated at the peak of each volume-dependence
curve, extrapolated to zero gradient-flow-depth, at different values of
lattice spacing $a$, for $T$ = 10 $T_{\mathrm{c}}$. The zero-spacing
extrapolation yields a value of $\Gamma_{\mathrm{sphal},s}/T^{4}=0.0035\pm
0.0002$.

The result can be seen in figure 12. As frequently happens, lattice spacing
effects are large at $N_{t}=6$, but for larger $N_{t}$ we find a good
convergence towards a continuum limit.

Rather than repeat this continuum limit procedure at every temperature, we
will fix to $N_{t}=10$ at all other temperatures, which we see will introduce
at most $\sim 20\%$ errors in our results.

##  6 Low and high temperatures

Before we proceed to show the final results for the sphaleron rate at
different temperatures, we also need to comment on the limits of our method.

###  6.1 Low temperatures and instantons

Consider the toy model one more time. In the high-temperature regime, the
$\varphi$ field typically reaches $\varphi=\pi$ as a result of a Euclidean
history which loiters near $\varphi=\pi$ (a sphaleron), rather than a history
which transitions from $\varphi=0$ to $\varphi=2\pi$ (an instanton). But as
the temperature is lowered, the instantons become more important. In this
case, the basic picture – that the real-time rate is controlled by classical
thermal fluctuations up to the barrier and the rate can be computed from the
frequency of such transitions over the barrier – comes into question.
Therefore we should not trust our approach if instantons are as common as
sphalerons.

One feature of our approach, as one sees in Figure 4 for instance, is that a
sphaleron will give the same answer whether we include the Euclidean path
between $t=0$ and between $t=\beta/2$ as we would find if we include the
Euclidean path between $t=\beta/2$ and $t=\beta$. We see this more explicitly
in Eq. (17). If we had chosen $t=\beta/2$ and $t=\beta$ as the lower and upper
boundaries of our region, we would define

$\displaystyle Q_{\mathrm{S}}^{\prime}$
$\displaystyle=Q(t=\beta/2)+\left[\int_{\beta/2}^{\beta}\mathrm{d}t\int\mathrm{d}^{3}x\>q(x,t)\right]-Q(t=\beta)$
$\displaystyle=Q_{\beta/2}+Q_{\mathrm{second-half}}-Q_{0}$
$\displaystyle=-(Q_{0}+Q_{\mathrm{half}}-Q_{\beta/2})+\left[\int_{0}^{\beta}\mathrm{d}t\int\mathrm{d}^{3}x\>q(x,t)\right]$
$\displaystyle=-Q_{\mathrm{S}}+Q_{\mathrm{I}}\,.$ (22)

In passing from the first to second line we use periodicity, which means that
the configuration at $t=\beta$ is the same as at $t=0$ so $Q_{\beta}=Q_{0}$.
The integral on the third line is the total topology in the box,
$Q_{\mathrm{I}}$. At high temperatures $Q_{\mathrm{I}}$ is almost always zero,
and therefore $Q_{\mathrm{S}}^{\prime}=-Q_{\mathrm{S}}$; the topology we find
in the upper half of the box is minus the topology we find in the lower half.
But for configurations with an actual instanton, $Q_{\mathrm{I}}=\pm
1=Q_{\mathrm{S}}+Q_{\mathrm{S}}^{\prime}$. Therefore, true sphalerons
($Q_{\mathrm{S}}=\pm 1$ but $Q_{\mathrm{I}}=0$) will count as a sphaleron
whether we use the top or the bottom half of our box, whereas instantons will
either count as a sphaleron if we use the top half, or if we use the bottom
half.

As we just emphasized, the whole approach of counting sphalerons and equating
them to classical transitions becomes suspect when a large fraction of what we
find are actually instantons. Therefore it is useful to investigate the ratio
of the instanton density to the sphaleron rate. When this ratio becomes of
order $1/2$, we have left the regime where classical activation makes sense.
When the ratio is small, sphalerons dominate instantons and our approach
should make sense.

Figure 13: Ratio between the amount of instantons ($Q_{\mathrm{I}}\neq 0$) and
sphalerons ($Q_{\mathrm{S}}\neq 0$) as a function of temperature (in
$T_{\mathrm{c}}$ units). Only the sphaleron rate has been extrapolated to zero
gradient-flow-depth, as described in the previous section.

We present this comparison in Figure 13. For values of temperatures above 1.75
$T_{\mathrm{c}}$ the appearance of instantons is exponentially suppressed and
finding one is very rare. For values between 1.3 and 1.75 $T_{\mathrm{c}}$,
while we can find some instantons, they are rare when compared to the amount
of sphalerons (sphalerons are at least 10 times more common than instantons).
For values lower than 1.3 $T_{\mathrm{c}}$, the instantons are very common and
therefore we cannot trust our calculation for the sphaleron rate. Therefore,
our method will only be reliable above 1.3 $T_{\mathrm{c}}$, and we will
restrict ourselves to this domain.

###  6.2 High temperatures and classical fields

At very high temperatures the gauge coupling becomes so weak that $g^{2}T\ll
T$. In this case, sphalerons are physically large objects which should be well
described in terms of classical fields, both thermodynamically and (within
limits) dynamically. In this regime, both our methods and previous, classical-
field real-time methods should be applicable. Therefore we will attempt to
make contact with the 3D real-time results of Ref. Moore:2010jd .

Real-time results are calculated in terms of the effective 3D gauge coupling
$g_{3}^{2}$, which equals $g^{2}T$ of the 4D theory at tree level. At the loop
level, it can be related to the 4D $\overline{\text{MS}}$ coupling via a
perturbative matching calculation, which has been carried out to 2 loops by
Laine and Schroeder in Ref. Laine:2005ai , see Eq. (2.34) to (2.37) of that
reference.

The $\overline{\text{MS}}$ coupling can in turn can be related to the lattice
coupling $\beta_{\text{latt}}$ by measuring the gradient-flowed expectation
value of the squared field strength
$E\equiv\,\mathrm{Tr}\,G_{\mu\nu}G^{\mu\nu}$ and using its NNLO relation to
the $\overline{\text{MS}}$ coupling at
$\mu=1/\sqrt{8{\tau_{\\!{}_{\mathrm{F}}}}}$ Harlander:2016vzb :

$t^{2}\expectationvalue{E(t)}=\frac{3\alpha_{s}}{4\pi}\left(1+k_{1}\alpha_{s}+k_{2}\alpha_{s}^{2}\right)$
(23)

where $k_{1}\approx 1.098$ and $k_{2}\approx-0.982$ for SU(3) theory without
fermions. We use the Wilson action, Zeuthen flow Ramos:2015baa , and clover
definition of the field strength, checking that the result is stable over a
range of flow depths after the 2-loop running of
$g^{2}(\mu_{\overline{\text{MS}}})$ is taken into account.

Equipped with this method, we study a lattice with $N_{t}=8$ and
$\beta_{\text{latt}}=15$, finding that it corresponds to
$\frac{g_{3}^{2}}{4\pi}=0.0406$. Since the coupling is small, the box size
must be correspondingly large before the sphaleron rate “turns on.” For this
choice of $N_{t}$, we find that we need $N_{s}=96$ to reach the peak in the
volume dependence of the sphaleron rate. Generating many independent
configurations of this size is challenging, so our statistical errors are
rather large; we find that
$\Gamma_{\mathrm{sphal},s}^{\text{new}}/(\alpha_{s}^{4}T^{4})=19\pm 7$. Based
on figure 12, we expect the real answer to be $20\%$ to $25\%$ lower once the
continuum limit is performed. Going back to Moore:2010jd , we find that for
$\alpha_{s}\text{(0 flavor)}=0.039$, the value is
$\Gamma_{\mathrm{sphal},s}^{\text{old}}/(\alpha_{s}^{4}T^{4})=11.5\pm 0.3$.
Therefore, our method is compatible with the previous 3D calculations, given
the downwards continuum corrections, the rather large statistical errors, and
the theoretical errors present in the calculation in Moore:2010jd .

##  7 Results, Discussion, and Conclusions

Figure 14: Results for the sphaleron rate for different temperatures. The left plot presents the dimensionless ratio found by scaling with $T^{4}$, while on the right we have reintroduced dimensions, using $T_{\mathrm{c}}\approx 280$ MeV for SU(3) theory. Errors are statistical only; all data are at $N_{t}=10$ without a continuum extrapolation. $T/T_{\mathrm{c}}$ | $\Gamma_{\mathrm{sphal},s}/T^{4}$ | $\Gamma_{\mathrm{sphal},s}[\mathrm{GeV}^{4}]$ | Num. configs.
---|---|---|---
1.3 | $0.061\pm 0.002$ | $0.00119\pm 0.00003$ | 5520
2.0 | $0.029\pm 0.002$ | $0.0031\pm 0.0002$ | 1970
3.0 | $0.016\pm 0.001$ | $0.0089\pm 0.0006$ | 1280
4.0 | $0.0115\pm 0.0006$ | $0.020\pm 0.001$ | 2260
5.0 | $0.0080\pm 0.0003$ | $0.034\pm 0.001$ | 7680
7.0 | $0.0069\pm 0.0007$ | $0.11\pm 0.01$ | 1280
10.0 | $0.0043\pm 0.0002$ | $0.29\pm 0.01$ | 4560
20.0 | $0.0027\pm 0.0003$ | $2.9\pm 0.3$ | 1280
Table 2:  Results presented in figure 14 in tabular form. Errors are
statistical only.

The results obtained with our method for the sphaleron rate in the range of
temperatures $[1.3-20]\,T_{\mathrm{c}}$ are shown in figure 14 and summarized
in table 2. We emphasize that these results are at $N_{t}=10$; we have not
taken a continuum limit.

The main application of these rates is to use them in Eq. (4) to see how
quickly or slowly a chiral imbalance relaxes via thermal sphaleron processes.
Using $\Gamma_{\mathrm{sphal},s}/T^{3}\simeq 0.04T$ over the range
$[1.3T_{\mathrm{pc}},2T_{\mathrm{pc}}]$ and treating the up and down flavors
as light ($n_{f}=2$), we estimate that the exponential decay lifetime for
axial number is of order $12/T$. This is slow enough to allow axial number to
play a role in early dynamics during heavy ion collisions, but marginally fast
enough for axial number to relax before the system hadronizes. Note that this
estimate involves using pure-glue results in full QCD, and it therefore still
has rather large theory errors.

We have also calculated the sphaleron rate at the electroweak scale, where the
coupling is approximately $\alpha_{s}\sim 0.1$. For that we have used the same
techniques as in section 6.2 to find the 4D coupling from gradient flow. For
$\beta_{\text{latt}}=9$, we have found that it corresponds to
$\alpha_{s}=0.108\pm 0.002$ (where the uncertainty shown is not statistical
but comes from having to pick a flow depth such that equation 23 is valid but
lattice errors are small). In this case, we obtain

$\Gamma_{\mathrm{sphal},s}(\alpha_{s}\sim 0.108)=14\pm
1\,\alpha_{s}^{4}T^{4}.$ (24)

We emphasize again that the error is statistical only; we have not taken a
continuum limit and that the result is for pure glue, so there are still much
larger systematic effects to be included.

The main breakthrough of this paper is to show how real-time processes
mediated by semiclassical saddlepoints, like thermal bubble nucleation or the
sphaleron rate, can be calculated using exclusively Euclidean and fully
nonperturbative tools. We have applied this technique to the sphaleron rate in
$SU(3)$ gauge theory partly because of its phenomenological interest and
partly because it is particularly clear in this case how to define the
separatrix and how to determine whether the $t=0$ and $t=\beta/2$ time slices
of the Euclidean path integral lie on opposite sides of the separatrix. The
same approach should be applicable to other saddle-controlled tunneling
processes such as bubble nucleation, though identifying the saddlepoint in
this case may be more challenging. Our approach also has the weakness that it
does not give us a way to determine a dynamical prefactor $\mathbf{d}$ which
appears in the tunneling rate. Note that this same defect is present in
analytical Euclidean approaches such as Affleck’s method as well. We currently
don’t see a way to overcome this limitation.

It should be straightforward to extend the current approach to full QCD with
2+1 light flavors of fermions. The only change needed is to include the
fermions when performing the HMC lattice updates, so that the configurations
are drawn from the full, rather than pure-glue, lattice ensemble. We have
already begun work on this project.

There are a few points where it would be useful to deepen our understanding of
the technique and our results. The unexpected scaling with lattice volume
needs to be better understood. And it would also be interesting to better
understand lower temperatures where instantons start to play a role. To what
extent can we incorporate instantons into a determination of axial number
relaxation? And to what extent does axial quark number make sense in a regime
where the system is moving towards a description in terms of hadronic degrees
of freedom? We consider these to be interesting open issues.

## Acknowledgments

The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the CRC-TR 211 ’Strong-interaction matter
under extreme conditions’– project number 315477589 – TRR 211. Calculations
were performed on the high-performance computer Lichtenberg at the NHR Centers
NHR4CES at TU Darmstadt.

##  Appendix A Size of Euclidean fluctuations

We have argued that the sphaleron rate can be determined by finding which
Euclidean configurations “span” the separatrix, in the sense that the spatial
slice at $t=0$ and the slice at $t=\beta/2$ are on different sides of the
separatrix.

Here we present the details of the calculation comparing the likelihood for a
Euclidean configuration to satisfy this condition to the real-time rate for
transitions as determined by the Affleck method. We also incorporate the
effect of 4D gradient flow applied to the Euclidean configuration before the
measurements are carried out.

The main assumption behind the method is that the sphaleron, critical bubble,
or other object responsible for tunneling is a robust high-action object. Even
though the transition rate may not be dominated by configurations close to the
saddle point, it will be dominated by configurations in the vicinity of the
(codimension-1) separatrix surface. At every point on the separatrix there is
a direction orthogonal to the separatrix, and we assume that fluctuations
along this direction, when comparing 3D slices at different $t$ values, are
small relative to the overall size of the sphaleron. In this case the degree
of freedom $\phi(t)$ associated with this direction can be treated as a single
variable with a kinetic term whose strength is approximately
$\phi$-independent and a potential which can be approximated as quadratic. (We
do not need this kinetic term or this potential to be the same everywhere on
the separatrix – we only need to know that there is little variation
orthogonal to the separatrix, over the field-space distances explored by the
higher Matsubara modes of this degree of freedom.) If these assumptions are
far from correct, the entire saddlepoint approach we use will not be
applicable and we have no theoretical tools to establish the rate from
Euclidean methods (and probably real-time methods will also not be applicable
since they also rely on a scale separation between the size and evolution time
scale of sphalerons and the intrinsic thermal size and time scales).

###  A.1 Calculation without gradient flow

We are interested in a configuration close to some point on a separatrix
surface, and in the degree of freedom orthogonal to the separatrix. We will
write this degree of freedom as $\phi(t)$ with $\phi=0$ representing the
separatrix, and we take it to have a canonically normalized kinetic term
(recall that a rescaling of $\phi(t)$ does not affect our answers). For the
toy model, $\phi(t)=\sqrt{m}(\varphi(t)-\pi)$. The probability that the
Euclidean fluctuations carry $\phi$ across $\phi=0$ between $t=0$ and
$t=\beta/2$ are given by:

$\Gamma_{\mathrm{Eucl}}\equiv\int_{\phi(t=0)=\phi(t=\beta)}\mathcal{D}\phi\;e^{-S_{E}}\;\left[\vphantom{\Big{|}}\theta(\phi(t{=}0))\theta(-\phi(t{=}\beta/2))+\theta(-\phi(t{=}0))\theta(\phi(t{=}\beta/2))\right].$
(25)

The two $\theta(..)\theta(..)$ factors represent cases where $\phi$ goes from
positive to negative and where it goes from negative to positive,
respectively. Here the Euclidean action $S_{E}$ is

$S_{E}=\int_{0}^{\beta}\mathrm{d}{t}\left(\frac{1}{2}\left(\partial_{t}\phi\right)^{2}-\frac{m^{2}}{2}\phi^{2}\right),$
(26)

with $m^{2}$ representing a possible potential, favoring $\phi$ to leave the
separatrix (the curvature in the unstable direction of the saddlepoint). We
insert the following Matsubara decomposition, which satisfies the periodicity
condition but is otherwise general:

$\phi(t)=\phi_{0}+\sum_{n=1}^{\infty}\left(c_{n}\cos(2\pi
nt/\beta)+s_{n}\sin(2\pi nt/\beta)\right).$ (27)

With this substitution, the measure of the path integral becomes simply

$\mathcal{D}\phi=\mathcal{N}\mathrm{d}{\phi_{0}}\prod_{k=1}^{\infty}\mathrm{d}{c_{k}}\mathrm{d}{s_{k}}$
(28)

where $\mathcal{N}$ is a constant, and we can explicitly evaluate the action

$\displaystyle S_{E}$
$\displaystyle=-\frac{\beta\,m^{2}}{2}\phi_{0}^{2}+\sum_{n=1}^{\infty}(c_{n}^{2}+s_{n}^{2})A(n,m)\,,$
(29) $\displaystyle A(n,m)$
$\displaystyle\equiv\frac{\beta}{4}\left[\left(\frac{2\pi
n}{\beta}\right)^{2}-m^{2}\right].$

Introducing $c_{\mathrm{even}}$ and $c_{\mathrm{odd}}$ as the sum over the
even and odd $c_{k}$ coefficients through

$1=\int\mathrm{d}c_{\mathrm{even}}\,\mathrm{d}c_{\mathrm{odd}}\;\delta\left(c_{\mathrm{even}}-\sum_{k=2,4,\ldots}c_{k}\right)\delta\left(c_{\mathrm{odd}}-\sum_{k=1,3,\ldots}c_{k}\right)$
(30)

we can rewrite the Heaviside theta-functions as:

$\begin{split}\theta(\phi(t=0))&=\theta\left(\phi_{0}+\sum_{k}c_{k}\right)=\theta(\phi_{0}+c_{\mathrm{even}}+c_{\mathrm{odd}})\,,\\\
\theta(-\phi(t=\beta/2))&=\theta\left(-\phi_{0}-\sum_{k}(-1)^{k}c_{k}\right)=\theta(-\phi_{0}+c_{\mathrm{even}}-c_{\mathrm{odd}})\,.\end{split}$
(31)

so that the path integral we have to perform is given by

$\displaystyle\Gamma_{\mathrm{Eucl}}$
$\displaystyle=2\int\mathrm{d}{\phi_{0}}\int\mathrm{d}c_{\mathrm{even}}\mathrm{d}c_{\mathrm{odd}}\prod_{k=1}^{\infty}\mathrm{d}{c_{k}}\mathrm{d}{s_{k}}\,\exp(\frac{m^{2}}{2}\phi_{0}^{2}\beta-\sum_{n=1}^{\infty}(c_{n}^{2}+s_{n}^{2})A(n,m))$
(32)
$\displaystyle{}\times\delta\left(c_{\mathrm{even}}{-}\\!\\!\sum_{k=2,4,\ldots}\\!\\!c_{k}\right)\delta\left(c_{\mathrm{odd}}{-}\\!\\!\sum_{k=1,3,\ldots}\\!\\!c_{k}\right)\theta(\phi_{0}{+}c_{\mathrm{even}}{+}c_{\mathrm{odd}})\theta(-\phi_{0}{+}c_{\mathrm{even}}{-}c_{\mathrm{odd}}).$

The factor of 2 accounts for the other $\theta$-function combination, which
gives precisely the same answer as the one shown.

Exponentiating the delta functions

$\delta\left(c_{\mathrm{even}}-\sum_{k}c_{k}\right)=\int\frac{\mathrm{d}{\lambda}}{2\pi}\;\exp\left(i\lambda
c_{\mathrm{even}}-i\sum_{k}\lambda c_{k}\right)$ (33)

makes all integrals either Gaussians or modified Gaussians, and it is
straightforward to evaluate them.

This calculation suffers from an unknown overall scaling, but recall that we
are interested in the ratio of this result to the standard result, which is
$\langle|\mathrm{d}\phi/\mathrm{d}t_{\mathrm{m}}|\rangle$ times the
probability density for $\phi_{0}=0$. Equipartition tells us that
$\langle|\mathrm{d}\phi/\mathrm{d}t_{\mathrm{m}}|\rangle=\sqrt{\frac{2T}{\pi}}$,
see Eq. (10). To evaluate the probability density, we compute the same path
integral, but with the two Heaviside functions replaced by a single delta
function $\delta(\phi_{0})$ – essentially, this determines the normalization
factor $\mathcal{N}$:

$\mathcal{N}=\int_{\phi(t=0)=\phi(t=\beta)}\mathcal{D}\phi\exp(-\int_{0}^{\beta}\mathrm{d}{t}\frac{1}{2}\left(\partial_{t}\phi\right)^{2})\delta(\phi_{0})\,,$
(34)

which means that

$\mathcal{N}=\prod_{n=1}^{\infty}\frac{A(n,0)}{\pi}.$ (35)

With this, by defining

$\begin{dcases}S_{\text{odd (even)}}&\equiv\sum_{n\,\text{odd
(even)}}\frac{1}{4A(n,m)}\\\
C&\equiv\frac{1}{4S_{\text{even}}}-\frac{m^{2}\beta}{2}\\\
b&\equiv\frac{m^{2}\beta}{2}\left(1+\frac{m^{2}\beta}{8C}\right)\end{dcases}$
(36)

the final result for the integral Eq. (25) is

$\Gamma_{\mathrm{Eucl}}=e^{-\frac{m^{2}\beta}{2}}\,\frac{m\beta}{2\sin(m\beta/2)}\sqrt{\frac{1}{\pi
S_{\text{even}}bC}}\operatorname{arctanh}\left(2\sqrt{S_{\text{odd}}b}\right)\,.$
(37)

When the saddlepoint is a strongly semiclassical configuration, the action
should not change too quickly as we move away from it, meaning that $m/2\pi
T\ll 1$. If we expand the result in this limit, it simplifies:

$\Gamma_{\mathrm{Eucl}}=4e^{-m^{2}\beta/2}\,\sqrt{\frac{S_{\text{odd}}}{\pi}}+\mathcal{O}(\sqrt{\beta}(m\beta)^{2})=\sqrt{\frac{\beta}{2\pi}}+\mathcal{O}(\sqrt{\beta}(m\beta)^{2}).$
(38)

This is the rate divided by the probability density to be on the separatrix.
To fully normalize it by the real-time rate, we should divide by
$\sqrt{2T/\pi}$ which we found above, to find:

$\displaystyle\frac{\Gamma_{\mathrm{Eucl}}}{\mbox{real-time
rate}}=\sqrt{\frac{1}{2\pi
T}}\left(\sqrt{\frac{2T}{\pi}}\right)^{-1}=\frac{1}{2T}\,,$ or
$\displaystyle\boxed{\mbox{Real-time rate}=2T\;\Gamma_{\mathrm{Eucl}}\,.\;}$
(39)

This provides the desired relation between the Affleck real-time rate and the
Euclidean sphaleron density which we have defined.

###  A.2 Calculation with flow

In practice it is necessary to modify the $\phi(t)$ field inside the Euclidean
path integral through the application of gradient flow, before performing any
evaluations on it. How does this affect the probability that $\phi(t=0)$ and
$\phi(t=\beta/2)$ are on opposite sides of the separatrix? We will answer this
now. Gradient flow changes the measurables, not the path integral or path
integral variables. In our case, it damps the $t$-dependent fluctuations
through application of the heat equation. Define the depth of heat-equation
flow to be ${\tau_{\\!{}_{\mathrm{F}}}}$. Then we set the
${\tau_{\\!{}_{\mathrm{F}}}}=0$ boundary conditions to be
$\phi({\tau_{\\!{}_{\mathrm{F}}}}=0,t)=\phi(t)$, and the gradient flow itself
is described by

$\partial_{\tau_{\\!{}_{\mathrm{F}}}}\phi({\tau_{\\!{}_{\mathrm{F}}}},t)-\partial_{t}^{2}\phi({\tau_{\\!{}_{\mathrm{F}}}},t)=0.$
(40)

Using the same Fourier-series expansion for $\phi(t)$ as before, one easily
finds that

$\phi({\tau_{\\!{}_{\mathrm{F}}}},t)=\phi_{0}+\sum_{n=1}^{\infty}\left(c_{n}\cos(2\pi
nt/\beta)+s_{n}\sin(2\pi
nt/\beta)\right)\exp(-{\tau_{\\!{}_{\mathrm{F}}}}\left(2\pi
n/\beta\right)^{2}).$ (41)

Now we can proceed as we did in the unflowed case, but taking the flow into
account on our theta functions, while the action remains the same as in Eq.
(29). The Heaviside functions become

$\begin{split}&\theta(\phi({\tau_{\\!{}_{\mathrm{F}}}},0))=\theta\left(\phi_{0}+\sum_{k}c_{k}\exp(-{\tau_{\\!{}_{\mathrm{F}}}}\left(2\pi
n/\beta\right)^{2})\right)\\\
&\theta(-\phi({\tau_{\\!{}_{\mathrm{F}}}},\beta/2))=\theta\left(-\phi_{0}-\sum_{k}\left(-1\right)^{k}c_{k}\exp(-{\tau_{\\!{}_{\mathrm{F}}}}\left(2\pi
n/\beta\right)^{2})\right).\end{split}$ (42)

Proceeding as we have done before, and defining in this case

$S^{\tau_{\\!{}_{\mathrm{F}}}}_{\text{odd (even)}}\equiv\sum_{n\,\text{odd
(even)}}\frac{\exp\left(-2{\tau_{\\!{}_{\mathrm{F}}}}\left(2\pi
n/\beta\right)^{2}\right)}{4A(n,m)}$ (43)

we find

$\begin{split}\Gamma_{\mathrm{Eucl}}^{\tau_{\\!{}_{\mathrm{F}}}}&=e^{-\frac{m^{2}}{2}\beta}\frac{m\beta}{2\sin(m\beta/2)}\sqrt{\frac{1}{\pi
S^{\tau_{\\!{}_{\mathrm{F}}}}_{\text{even}}bC}}\operatorname{arctanh}\left(2\sqrt{S^{\tau_{\\!{}_{\mathrm{F}}}}_{\text{odd}}b}\right)\\\
&\Rightarrow\;\sqrt{\frac{\beta}{2\pi}}\sqrt{\frac{8}{\pi^{2}}\sum_{n=1,3,\ldots}\frac{e^{-8\pi^{2}{\tau_{\\!{}_{\mathrm{F}}}}T^{2}n^{2}}}{n^{2}}}\\\
&\simeq\sqrt{\frac{\beta}{2\pi}}\times\sqrt{1-8\sqrt{\frac{2{\tau_{\\!{}_{\mathrm{F}}}}}{\pi\beta^{2}}}}.\end{split}$
(44)

As before, the quantity after $\Rightarrow$ is the result when taking $m\to
0$. The last line is the result of a small ${\tau_{\\!{}_{\mathrm{F}}}}$
expansion, with errors which are exponentially small provided that
${\tau_{\\!{}_{\mathrm{F}}}}\ll 1/(8\pi^{2}T^{2})$. In practice we will use
the unexpanded expression shown in the middle line. The final expression is
useful because it shows that the correction goes to 1 as
${\tau_{\\!{}_{\mathrm{F}}}}\to 0$ but only as a square root, that the
behavior is monotonic, and that the domain over which the gradient flow can be
viewed as “small” is roughly ${\tau_{\\!{}_{\mathrm{F}}}}\ll
1/(8\pi^{2}T^{2})$.

The best way to use this result is:

$\boxed{\mbox{Real-time
rate}=\frac{2T\;\Gamma_{\mathrm{Eucl}}^{\tau_{\\!{}_{\mathrm{F}}}}}{\sqrt{1-8\sqrt{\frac{2{\tau_{\\!{}_{\mathrm{F}}}}T^{2}}{\pi}}}}\,.\;}$
(45)

## References

  * [1] A. A. Belavin, Alexander M. Polyakov, A. S. Schwartz, and Yu. S. Tyupkin.  Pseudoparticle Solutions of the Yang-Mills Equations.  Phys. Lett. B, 59:85–87, 1975.
  * [2] Gerard ’t Hooft.  Symmetry Breaking Through Bell-Jackiw Anomalies.  Phys. Rev. Lett., 37:8–11, 1976.
  * [3] Gerard ’t Hooft.  Computation of the Quantum Effects Due to a Four-Dimensional Pseudoparticle.  Phys. Rev. D, 14:3432–3450, 1976.  [Erratum: Phys.Rev.D 18, 2199 (1978)].
  * [4] Stephen L. Adler.  Axial vector vertex in spinor electrodynamics.  Phys. Rev., 177:2426–2438, 1969.
  * [5] J. S. Bell and R. Jackiw.  A PCAC puzzle: $\pi^{0}\to\gamma\gamma$ in the $\sigma$ model.  Nuovo Cim. A, 60:47–61, 1969.
  * [6] Frans R. Klinkhamer and N. S. Manton.  A Saddle Point Solution in the Weinberg-Salam Theory.  Phys. Rev. D, 30:2212, 1984.
  * [7] V. A. Kuzmin, V. A. Rubakov, and M. E. Shaposhnikov.  On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe.  Phys. Lett. B, 155:36, 1985.
  * [8] V. A. Rubakov and M. E. Shaposhnikov.  Electroweak baryon number nonconservation in the early universe and in high-energy collisions.  Usp. Fiz. Nauk, 166:493–537, 1996.
  * [9] James M. Cline.  Baryogenesis.  In Les Houches Summer School - Session 86: Particle Physics and Cosmology: The Fabric of Spacetime, 9 2006.
  * [10] Larry D. McLerran, Emil Mottola, and Mikhail E. Shaposhnikov.  Sphalerons and Axion Dynamics in High Temperature QCD.  Phys. Rev. D, 43:2027–2035, 1991.
  * [11] G. F. Giudice and Mikhail E. Shaposhnikov.  Strong sphalerons and electroweak baryogenesis.  Phys. Lett. B, 326:118–124, 1994.
  * [12] Kenji Fukushima, Dmitri E. Kharzeev, and Harmen J. Warringa.  The Chiral Magnetic Effect.  Phys. Rev. D, 78:074033, 2008.
  * [13] Dam T. Son and Piotr Surowka.  Hydrodynamics with Triangle Anomalies.  Phys. Rev. Lett., 103:191601, 2009.
  * [14] M. A. Stephanov and Y. Yin.  Chiral Kinetic Theory.  Phys. Rev. Lett., 109:162001, 2012.
  * [15] A. Bazavov, H. T. Ding, P. Hegde, F. Karsch, C. Miao, Swagato Mukherjee, P. Petreczky, C. Schmidt, and A. Velytsky.  Quark number susceptibilities at high temperatures.  Phys. Rev. D, 88(9):094021, 2013.
  * [16] Szabocls Borsanyi, Zoltan Fodor, Christian Hoelbling, Sandor D. Katz, Stefan Krieg, and Kalman K. Szabo.  Full result for the QCD equation of state with 2+1 flavors.  Phys. Lett. B, 730:99–104, 2014.
  * [17] A. Bazavov et al.  Equation of state in ( 2+1 )-flavor QCD.  Phys. Rev. D, 90:094503, 2014.
  * [18] Dmitri Yu. Grigoriev and V. A. Rubakov.  Soliton Pair Creation at Finite Temperatures. Numerical Study in (1+1)-dimensions.  Nucl. Phys. B, 299:67–78, 1988.
  * [19] Dietrich Bodeker, Larry D. McLerran, and Andrei V. Smilga.  Really computing nonperturbative real time correlation functions.  Phys. Rev. D, 52:4675–4690, 1995.
  * [20] Peter Brockway Arnold, Dam Son, and Laurence G. Yaffe.  The Hot baryon violation rate is O (alpha-w**5 T**4).  Phys. Rev. D, 55:6264–6273, 1997.
  * [21] D. Bodeker, Guy D. Moore, and K. Rummukainen.  Chern-Simons number diffusion and hard thermal loops on the lattice.  Phys. Rev. D, 61:056003, 2000.
  * [22] Jan Ambjorn, T. Askgaard, H. Porter, and M. E. Shaposhnikov.  Sphaleron transitions and baryon asymmetry: A Numerical real time analysis.  Nucl. Phys. B, 353:346–378, 1991.
  * [23] Jan Ambjorn and A. Krasnitz.  Improved determination of the classical sphaleron transition rate.  Nucl. Phys. B, 506:387–403, 1997.
  * [24] Guy D. Moore, Chao-ran Hu, and Berndt Muller.  Chern-Simons number diffusion with hard thermal loops.  Phys. Rev. D, 58:045001, 1998.
  * [25] Peter Brockway Arnold.  Hot B violation, the lattice, and hard thermal loops.  Phys. Rev. D, 55:7781–7796, 1997.
  * [26] Dietrich Bodeker.  On the effective dynamics of soft nonAbelian gauge fields at finite temperature.  Phys. Lett. B, 426:351–360, 1998.
  * [27] Guy D. Moore.  The Sphaleron rate: Bodeker’s leading log.  Nucl. Phys. B, 568:367–404, 2000.
  * [28] Peter Brockway Arnold and Laurence G. Yaffe.  High temperature color conductivity at next-to-leading log order.  Phys. Rev. D, 62:125014, 2000.
  * [29] Guy D. Moore.  Computing the strong sphaleron rate.  Phys. Lett. B, 412:359–370, 1997.
  * [30] Guy D. Moore and Marcus Tassler.  The Sphaleron Rate in SU(N) Gauge Theory.  JHEP, 02:105, 2011.
  * [31] Luis Altenkort, Alexander M. Eller, Olaf Kaczmarek, Lukas Mazur, Guy D. Moore, and Hai-Tao Shu.  Sphaleron rate from Euclidean lattice correlators: An exploration.  Phys. Rev. D, 103(11):114513, 2021.
  * [32] Luis Altenkort, Alexander M. Eller, O. Kaczmarek, Lukas Mazur, Guy D. Moore, and H. T. Shu.  Lattice QCD noise reduction for bosonic correlators through blocking.  Phys. Rev. D, 105:094505, 2022.
  * [33] J. S. Langer.  Statistical theory of the decay of metastable states.  Annals Phys., 54:258–275, 1969.
  * [34] Ian Affleck.  Quantum Statistical Metastability.  Phys. Rev. Lett., 46:388, 1981.
  * [35] Andrei D. Linde.  Fate of the False Vacuum at Finite Temperature: Theory and Applications.  Phys. Lett. B, 100:37–40, 1981.
  * [36] Andrei D. Linde.  Decay of the False Vacuum at Finite Temperature.  Nucl. Phys. B, 216:421, 1983.  [Erratum: Nucl.Phys.B 223, 544 (1983)].
  * [37] Sidney R. Coleman.  The Fate of the False Vacuum. 1. Semiclassical Theory.  Phys. Rev. D, 15:2929–2936, 1977.  [Erratum: Phys.Rev.D 16, 1248 (1977)].
  * [38] Curtis G. Callan, Jr. and Sidney R. Coleman.  The Fate of the False Vacuum. 2. First Quantum Corrections.  Phys. Rev. D, 16:1762–1768, 1977.
  * [39] Peter Brockway Arnold and Larry D. McLerran.  The Sphaleron Strikes Back.  Phys. Rev. D, 37:1020, 1988.
  * [40] Peter Brockway Arnold and Larry D. McLerran.  Sphalerons, Small Fluctuations and Baryon Number Violation in Electroweak Theory.  Phys. Rev. D, 36:581, 1987.
  * [41] R. Narayanan and H. Neuberger.  Infinite N phase transitions in continuum Wilson loop operators.  JHEP, 03:064, 2006.
  * [42] Martin Lüscher.  Properties and uses of the Wilson flow in lattice QCD.  JHEP, 08:071, 2010.  [Erratum: JHEP 03, 092 (2014)].
  * [43] Guy D. Moore.  Measuring the broken phase sphaleron rate nonperturbatively.  Phys. Rev. D, 59:014503, 1999.
  * [44] http://luscher.web.cern.ch/luscher/openQCD/index.html.
  * [45] Barbara De Palma, Marco Erba, Luca Mantovani, and Nicola Mosco.  A python program for the implementation of the Γ-method for monte carlo simulations.  Computer Physics Communications, 2018.
  * [46] A. Francis, O. Kaczmarek, M. Laine, T. Neuhaus, and H. Ohno.  Critical point and scale setting in SU(3) plasma: An update.  Phys. Rev. D, 91(9):096002, 2015.
  * [47] Y. Burnier, H. T. Ding, O. Kaczmarek, A. L. Kruse, M. Laine, H. Ohno, and H. Sandmeyer.  Thermal quarkonium physics in the pseudoscalar channel.  JHEP, 11:206, 2017.
  * [48] Jan Ambjorn and A. Krasnitz.  The Classical sphaleron transition rate exists and is equal to 1.1 (alpha(w) T)**4.  Phys. Lett. B, 362:97–104, 1995.
  * [49] Guy D. Moore and Kari Rummukainen.  Classical sphaleron rate on fine lattices.  Phys. Rev. D, 61:105008, 2000.
  * [50] M. Laine and Y. Schroder.  Two-loop QCD gauge coupling at high temperatures.  JHEP, 03:067, 2005.
  * [51] Robert V. Harlander and Tobias Neumann.  The perturbative QCD gradient flow to three loops.  JHEP, 06:161, 2016.
  * [52] A. Ramos and S. Sint.  Symanzik improvement of the gradient flow in lattice gauge theories.  Eur. Phys. J. C, 76(1):15, 2016.