# Hydrodynamic fluctuations and topological susceptibility in chiral
magnetohydrodynamics

Arpit Das<EMAIL_ADDRESS>School of Mathematics, University of Edinburgh,
Edinburgh EH9 3FD, U.K.
Higgs Centre for Theoretical Physics, University of Edinburgh, Edinburgh EH8
9YL, U.K.
Nabil Iqbal<EMAIL_ADDRESS>Centre for Particle Theory, Department
of Mathematical Sciences, Durham University, South Road, Durham DH1 3LE, UK
Napat Poovuttikul<EMAIL_ADDRESS>Department of Physics, Faculty of
Science Chulalongkorn University, Bangkok 10330, Thailand

###### Abstract

Chiral magnetohydrodynamics is devoted to understanding the late-time and
long-distance behavior of a system with an Adler-Bell-Jackiw anomaly at finite
temperatures. The non-conservation of the axial charge is determined by the
topological density $\vec{E}\cdot\vec{B}$; in a classical hydrodynamic
description this decay rate can be suppressed by tuning the background
magnetic field to zero. However it is in principle possible for thermal
fluctuations of $\vec{E}\cdot\vec{B}$ to result in a non-conservation of the
charge even at vanishing $B$-field; this would invalidate the classical
hydrodynamic effective theory. We investigate this by computing the real-time
susceptibility of the topological density at one-loop level in
magnetohydrodynamic fluctuations, relating its low-frequency limit to the
decay rate of the axial charge. We find that the frequency-dependence of this
susceptibility is sufficiently soft as to leave the axial decay rate
unaffected, validating the classical hydrodynamic description. We show that
the susceptibility contains non-analytic frequency-dependence which is
universally determined by hydrodynamic data. We comment briefly on possible
connections to the recent formulation of the ABJ anomaly in terms of non-
invertible symmetry.

###### Contents

  1. I Background
  2. II One-loop hydrodynamic fluctuations
    1. II.1 Kinematics of plasma correlations
      1. II.1.1 Tensor structure of correlators
    2. II.2 Computation of correlator
    3. II.3 Correlation of topological density
  3. III Discussion
  4. A Memory matrix review
  5. B Finite temperature conventions and Matsubara sums
    1. B.1 Performing Matsubara sums
  6. C Details of the 1-loop integration
    1. C.1 $\Omega>0$
    2. C.2 $\Omega<0$

##  I Background

We begin by briefly stating the problem. A chiral plasma belongs to the same
universality class (in the context of global symmetry) as that of massless
Dirac fermions coupled to dynamical QED at finite temperatures Joyce:1997uy ;
Rogachevskii:2017uyc ; Boyarsky:2015faa ; DelZanna:2018dyb ; Hattori:2017usa ;
Yamamoto:2016xtu . The global symmetry structure of the chiral plasma can be
summarized as follows Das:2022auy ,

$\displaystyle\partial_{\mu}J^{\mu\nu}=0,\qquad\qquad\qquad
J^{\mu\nu}\equiv{{\frac{1}{2}}}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma},$ (1)
$\displaystyle\partial_{\mu}j^{\mu}_{A}=k\,\epsilon_{\mu\nu\rho\sigma}J^{\mu\nu}J^{\rho\sigma}$
(2)

where the first equation states that magnetic field lines are conserved (as
$J^{\mu\nu}$ measures the magnetic flux) and the second equation is the usual
Adler-Bell-Jackiw anomaly expressed in terms of the conserved 2-form current
$J^{\mu\nu}$. Here $j^{\mu}_{A}$ is the non-conserved axial $1$-form current
and $k$ is the anomaly coefficient, which is $k\equiv\frac{1}{16\pi^{2}}$ in
the precise case of a single Dirac fermion. We can define the electric and
magnetic fields, in the usual way as in electrodynamics, in terms of the
components of the $2$-form current as follows: $J^{0i}=B_{i}$. Thus,
$E_{l}=\frac{1}{2}\epsilon_{ijl}J^{ij}$. Let us denote the topological density
by the operator

$Q(x)=\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\ .$ (3)

The non-conservation equation for the axial charge then reads

$\partial_{\mu}j^{\mu}_{A}=-kQ(x)$ (4)

The fact that the right-hand side of this expression involves a dynamical
operator $Q(x)$ rather than a fixed external source (as in the case of a ’t
Hooft anomaly) distinguishes chiral MHD from the well-understood problem of
the hydrodynamics of systems with ’t Hooft anomalies Son:2009tf ;
Neiman:2010zi  (see Landsteiner:2016led  for a review). In particular, this
non-conservation equation leads one to expect that in a thermal state the
axial charge should decay exponentially in time as $n_{A}\sim
e^{-\Gamma_{A}t}$.

In a formal limit where the anomaly coefficient $k$ is taken to be small, the
decay rate for the axial charge is given by the following formula:

$\Gamma_{A}=\lim_{\Omega\to
0}\frac{k^{2}}{\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}\Omega}\textrm{Im}\,G^{R}_{QQ}(\Omega,\vec{p}=0)$
(5)

where $\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}$ is the axial charge
susceptibility and $G^{R}_{QQ}(\Omega,\vec{p})$ is the retarded correlation
function of the topological density111We are grateful to L. Delacretaz for
suggesting this route for calculation Luca:rep22 .. This can be obtained from
the memory matrix formalism; see Appendix A for a brief review and references.
This expression is perturbative in $k$ but makes no assumptions on the
dynamics of the degrees of freedom entering the topological charge density.
Let us now explore some implications of this formula.

First, we note that in elementary language, $Q(x)=8\vec{E}\cdot\vec{B}$. Let
us consider a state with a background magnetic field $\vec{B}\neq 0$ pointing
along the $z$ axis, and study only the fluctuations of the electric field
about the equilibrium (i.e. assuming that $\vec{B}$ does not fluctuate). We
then find the following expression for the axial charge decay rate, in terms
of the retarded correlation function of the electric field operator $E$.

$\Gamma_{A}=\frac{64k^{2}B^{2}}{\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\lim_{\Omega\to
0}\frac{1}{\Omega}\textrm{Im}\,G^{R}_{E_{z},E_{z}}(\Omega)$ (6)

This formula was derived from an effective theory for chiral MHD in
Das:2022fho . In particular, we note that the resistivity $\rho$ of the plasma
is defined in terms of a Kubo formula for the electric field, as explained in
Grozdanov:2016tdf . Thus this expression states that the axial charge
relaxation rate is:

$\Gamma_{A}=\frac{64k^{2}B^{2}\rho}{\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\
.$ (7)

This expression can be understood from elementary arguments involving a
quasiparticle description Figueroa:2017hun , and has also been verified in a
holographic model Das:2022auy . It has also been subject to numerical
investigation in classical simulations Figueroa:2017hun ; Figueroa:2019jsi ,
where in particular the recent work Das:2023nwl  displays precise agreement
with the effective field theory formula (7).

Importantly, this expression states that as the magnetic field $B$ is taken to
zero, the lifetime of the axial charge is arbitrarily long. Indeed it is this
parametric separation of scales that implicitly lies behind the extensive
literature treating the axial charge density using hydrodynamics Joyce:1997uy
; Rogachevskii:2017uyc ; Boyarsky:2015faa ; DelZanna:2018dyb ; Hattori:2017usa
; Yamamoto:2016xtu . We note the recent works Das:2022fho ; Landry:2022nog
treat this problem from an effective theory viewpoint.

Upon reflection, however, this is a somewhat strong statement – in particular,
the right hand side of (5) does not obviously appear to vanish at zero $B$.
One may imagine that the thermal fluctuations of the topological density
$\vec{E}\cdot\vec{B}$ would create a nonzero decay rate which would become the
dominant decay channel when the applied magnetic field is sufficiently small.
This would lead to a crossover between the classical result (5) at moderately
strong $B$ and some nonzero fluctuation-driven effect at small $B$. If the
decay rate does not vanish at zero $B$, it would have significant
consequences: it would mean that at sufficiently long time scales the axial
charge simply does not exist as a hydrodynamic degree of freedom. In
particular, it would suggest that classical discussions of chiral
magnetohydrodyanmics – including the effective field theories described in
Das:2022fho ; Landry:2022nog  – are unstable towards the inclusion of
fluctuations.

One might be tempted to argue that the right-hand side of (5) is likely to
vanish as follows: the infrared limit of the correlation function is
presumably related to the integral of the Euclidean correlation function of
$Q(\tau,\vec{x})$ over all Euclidean spacetime. But we know that $Q$ is a
total derivative:

$Q=\partial_{\mu}K^{\mu}\qquad
K^{\mu}\equiv\epsilon^{\mu\nu\rho\sigma}A_{\nu}F_{\rho\sigma}$ (8)

and thus the integral will receive contributions only from field
configurations with nontrivial topological structure. However in an Abelian
gauge theory with vanishing background $B$ field there are no $U(1)$
instantons, and thus the integral is zero. The argument phrased above is
somewhat heuristic and it seems to us that it depends sensitively on boundary
conditions at infinity. We were unable to formulate a completely satisfactory
version of this argument, and indeed the art of extrapolating real-time
dynamical physics from Euclidean non-perturbative data is quite subtle (see
e.g. Arnold:1987zg ).

In this work we thus directly compute the leading contribution to the decay
rate $\Gamma_{A}$ from thermal fluctuations by evaluating the Kubo formula (5)
in the state with zero background magnetic field from magnetohydrodynamics. In
particular, we will evaluate the contribution to the retarded correlation
function $Q(x)=8\vec{E}\cdot\vec{B}$ arising from a one-loop calculation where
the propagating degrees of freedom are diffusive MHD waves, the leading low-
frequency degrees of freedom in the MHD plasma. Importantly, we will
demonstrate explicitly that these fluctuations do result in a nontrivial real-
time correlation function $G^{R}_{QQ}(\Omega)$ for the topological density –
which we calculate as a function of frequency $\Omega$ – but that this
function vanishes quickly enough at small frequency that it does not result in
a non-vanishing fluctuation-driven decay rate. This is consistent with the
heuristic argument above, and (to this order) is consistent with safe use of
the effective hydrodynamic description.

The work performed in this paper has previously appeared in the PhD thesis of
Arpit Das (see Chapter 7 of Das:2023vls ).

##  II One-loop hydrodynamic fluctuations

In this section, we compute the finite-frequency real-time topological
susceptibility arising from magnetohydrodynamic fluctuations. In particular,
we are interested in computing the retarded correlation function

$G^{R}_{QQ}(\Omega)=-i\int dtd^{3}xe^{i\Omega t}\mathop{\rm Tr}\left(e^{-\beta
H}[Q(\vec{x},t),Q(0)]\right){\theta}(t)\ .$ (9)

where the operator $Q=8\vec{E}\cdot\vec{B}$. We will now write the retarded
correlation function above in terms of correlation function of $\vec{E}$ and
$\vec{B}$; those can then be evaluated using an appropriate model for the
dynamics.

A convenient way to proceed is to express the retarded correlation function in
terms of the Euclidean finite-temperature correlation function
$G^{E}_{QQ}(i\Omega_{l})$, which is defined on a discrete set of Euclidean
Matsubara frequncies $i\Omega_{l}=\frac{2\pi i\mathbb{Z}}{\beta}$. The
retarded correlation function at real $\Omega$ is related to the Euclidean one
by the usual formula

$G^{R}_{QQ}(\Omega)=G^{E}_{QQ}(i\Omega_{l}=\Omega+i\epsilon)$ (10)

The Euclidean correlator can be explicitly written as:

$G^{E}_{QQ}(i\Omega_{l})=\int_{0}^{\beta}d\tau\int
d^{3}xe^{-i\Omega_{l}\tau}\langle Q(\tau,\vec{x})Q(0)\rangle$ (11)

To proceed, we use $Q(x)=8\vec{E}\cdot\vec{B}$. We also assume all
correlations in the fluid are Gaussian. This is a reasonable starting point,
as generally in hydrodynamics it is expected that the current densities –
which in this case are precisely the components of $\vec{B}$ and $\vec{E}$ are
themselves weakly coupled at long distances so that a classical treatment is
valid222In low-dimensional hydrodynamics there are known examples where non-
linearities are relevant in the IR (see e.g. Kovtun:2012rj  for a review) but
to our knowledge this is not expected to be the case for (3+1)d MHD..

This assumption means that we can factorize the correlators in Euclidean space
as follows:

$\displaystyle\left\langle\left(\vec{E}\cdot\vec{B}\right)(x)\left(\vec{E}\cdot\vec{B}\right)(0)\right\rangle=\delta^{pq}\delta^{rs}\left[\langle
E_{p}(x)E_{r}(0)\rangle\,\langle B_{q}(x)B_{s}(0)\rangle+\langle
E_{p}(x)B_{s}(0)\rangle\,\langle B_{q}(x)E_{r}(0)\rangle\right],$ (12)

Note that in this expression we have assumed that there is no background
topological density, i.e. $\langle\vec{E}\cdot\vec{B}\rangle=0$; this follows
from CP invariance of the thermal state.

This can be conveniently interpreted as a Feynman diagram bubble evaluated in
Euclidean spacetime, where the propagators of the bubble are the two-point
correlators $\langle E_{p}(x)E_{r}(0)\rangle$ etc., see Fig. 1. Indeed,
expressed in this form the problem has a great deal of formal similarity to
the classic problem of evaluating (e.g.) a one-loop conductivity in terms of
the propagators of the microscopic charged degrees of freedom – in both cases
we are interested in determining the correlation function of an operator (i.e.
$Q$) which is a bilinear in terms of fields with quadratic – and known –
correlations (i.e. $\vec{E}$ and $\vec{B}$).

Figure 1: A bubble diagram representing the frequency sum in (13). Here
$G(\omega)$ schematically represents either $G_{E_{i}E_{j}}(\omega)$ or
$G_{E_{i}B_{j}}(\omega)$.

Following the standard approach for that problem333See e.g. mahan2000many  or
a review in a holographic context in Faulkner:2013bna ., it is convenient to
write the expression in (Euclidean) frequency space, where we find that it
becomes

$G^{E}_{QQ}(i\Omega_{l})=64T\sum_{i\omega_{m}}\int\frac{d^{3}p}{(2\pi)^{3}}\left(G^{E}_{E_{i}E_{j}}(i\Omega_{l}+i\omega_{m})G^{E}_{B_{i}B_{j}}(i\omega_{m})+G^{E}_{E_{i}B_{j}}(i\Omega_{l}+i\omega_{m})G^{E}_{B_{i}E_{j}}(i\omega_{m})\right)$
(13)

The above expression includes the same sum over the index structure that is
present in (12).

We now note that we do not have access to the Euclidean correlation functions
of the electric and magnetic fields $G^{E}_{E_{i}B_{j}}$ etc. in any suitable
form. However from magnetohydrodynamics we do have access to the Lorentzian
spectral densities for these correlations at small real frequencies444It is in
general not trivial to analytically continue approximate expressions from real
to Euclidean frequencies.. It is thus convenient to rewrite this expression in
terms of these spectral densities, which can be done using standard finite-
temperature techniques (reviewed in Appendix B) to obtain:

$\displaystyle G^{E}_{QQ}(i\Omega_{l})$
$\displaystyle=-64\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d\omega_{1}}{2\pi}\frac{d\omega_{2}}{2\pi}\left[\frac{f(\omega_{1})-f(\omega_{2})}{\omega_{1}-\omega_{2}-i\Omega_{l}}\right](\rho_{E_{i}E_{j}}(\omega_{1})\rho_{B_{i}B_{j}}(\omega_{2})+\rho_{E_{i}B_{j}}(\omega_{1})\rho_{B_{i}E_{j}}(\omega_{2}))$
(14)

where $\rho_{EE}$ and $\rho_{EB}$ are the spectral densities associated with
the retarded correlation functions of $\vec{E}$ and $\vec{B}$, i.e.

$\rho_{E_{i}E_{j}}(\omega,p)=-\frac{1}{\pi}\textrm{Im}\,G^{R}_{E_{i}E_{j}}(\omega,p),$
(15)

and similarly. Here $f(x)$ is the Bose distribution function:

$f(\omega)=\frac{1}{e^{\beta\omega}-1}$ (16)

We have reduced the problem to evaluating integrals over these spectral
densities. We now discuss these correlation functions.

###  II.1 Kinematics of plasma correlations

In the remainder of this section, we discuss the tensor structure of the
correlators $\left<BB\right>$, $\left<EE\right>$ and $\left<EB\right>$,
expressing them in terms of scalar functions of momenta and frequencies and
describing what is known about them from magnetohydrodynamics.

####  II.1.1 Tensor structure of correlators

In this work we are interested in fluctuations about the plasma at finite
temperature with zero background magnetic field $\vec{B}_{0}=0$. We will
restrict attention to a parity-invariant theory. Here we record the
constraints on the correlation functions arising from parity invariance and
magnetic flux conservation.

The most general possible tensor decomposition of the retarded correlation
functions is,

$\begin{split}&\left<E_{i}E_{j}\right>=A(\omega,|\vec{p}|)\delta_{ij}+X(\omega,|\vec{p}|)\frac{p_{i}p_{j}}{p^{2}}+U(\omega,|\vec{p}|)\epsilon_{ijk}\frac{p^{k}}{|\vec{p}|},\\\
&\left<B_{i}B_{j}\right>=C(\omega,|\vec{p}|)\delta_{ij}+Y(\omega,|\vec{p}|)\frac{p_{i}p_{j}}{p^{2}}+V(\omega,|\vec{p}|)\epsilon_{ijk}\frac{p^{k}}{|\vec{p}|},\\\
&\left<E_{i}B_{j}\right>=M(\omega,|\vec{p}|)\delta_{ij}+N(\omega,|\vec{p}|)\frac{p_{i}p_{j}}{p^{2}}+K(\omega,|\vec{p}|)\epsilon_{ijk}\frac{p^{k}}{|\vec{p}|},\end{split}$
(17)

where we remind ourselves that $p^{2}$ above denote the square of the norm of
the 3-vector: $\vec{p}$. At times we will use the short-hand notation for the
above scalar functions, $Z^{\omega}$ to denote $Z(\omega,|\vec{p}|)$.

Noting that $\vec{E}$ is a vector and $\vec{B}$ a pseudo-vector under parity,
the scalar coefficient functions $U^{\omega}$, $V^{\omega}$, $M^{\omega}$ and
$N^{\omega}$ are all odd under parity and thus vanish in a parity-invariant
state. Thus, $A^{\omega},X^{\omega}$ and $C^{\omega},Y^{\omega}$ will be
expressed in terms of the diagonal components of the $\langle
E_{i}E_{j}\rangle$ and $\langle B_{i}B_{j}\rangle$ correlators, respectively.
On the contrary, $K^{\omega}$ will be expressed in terms of the off-diagonal
components of the $\langle E_{i}B_{j}\rangle$.

To impose the constraints from conservation of $J^{\mu\nu}$, let us assume a
plane wave basis of the form $e^{i(\vec{p}\cdot\vec{x}-\omega t)}$. Without
loss of generality, in this section we take the spatial momentum to be aligned
along the $z$-direction, that is, $p_{z}\neq 0$ implying, $p^{2}=p_{z}^{2}$
and $|\vec{p}|=p_{z}$. Now let us look at the various components of Eq. (1).
Individually, the temporal and $z$ components give: $J^{0z}=0$. The remaining
spatial components, that is for $i\in\\{x,y\\}$, we get:

$J^{zi}=\frac{\omega}{p_{z}}J^{0i},$ (18)

Using Eq. (18) and the definition of $\vec{E},\vec{B}$ in terms of the
components of the $2$-form current, we can find expressions for the scalar
functions, given in Eq. (17), in terms of the relevant two-point functions of
the form $\langle J^{\mu\nu}J^{\rho\sigma}\rangle$.

Let us first look at the $\langle B_{i}B_{j}\rangle$ correlator. Since
$J^{0z}=0$, we have $\langle B_{z}B_{z}\rangle=0$. This implies,
$C^{\omega}+Y^{\omega}=0$. From the diagonal $x,y$-components, we find
$C^{\omega}=\langle B_{i}B_{i}\rangle$, where $i\in\\{x,y\\}$. Following a
similar analysis for the $\langle E_{i}E_{j}\rangle$ correlator, we get from
its diagonal $z$-component: $A^{\omega}+X^{\omega}=\langle E_{z}E_{z}\rangle$.
Its remaining diagonal components give: $A^{\omega}=\langle
E_{i}E_{i}\rangle$. Finally, we can use Eq. (18) to relate $A^{\omega}$ and
$C^{\omega}$ as follows,

$\displaystyle C^{\omega}=\frac{p^{2}}{\omega^{2}}A^{\omega}.$ (19)

The value of the scalar function $K^{\omega}$ is yet to be determined. This
can be done by looking at the cross correlator: $\left<E_{i}B_{j}\right>$,
which we do next. First of all, note that for $i\neq j$ we have,

$\displaystyle\left<E_{i}B_{j}\right>=K\epsilon_{ijk}\frac{p^{k}}{|\vec{p}|}=\left<B_{j}E_{i}\right>,\qquad\left<B_{i}E_{j}\right>=-K\epsilon_{ijk}\frac{p^{k}}{|\vec{p}|}=-\left<E_{i}B_{j}\right>.$

Now since any correlator with $B_{z}$ in it vanishes, as $J^{0z}=0$, we have
from above that any correlator with $E_{z}$ in it should also vanish. This is
because these correlators just differ from each other by a minus sign. So, for
non-vanishing cross correlators of the form: $\left<E_{i}B_{j}\right>$, we
must have $i,j\neq z$. Choosing $i,j=x,y$, we get,

$\displaystyle\left<E_{x}B_{y}\right>$
$\displaystyle=K^{\omega}=\left<J^{yz}J^{0y}\right>=-\frac{\omega}{p_{z}}\left<J^{0y}J^{0y}\right>=-\frac{\omega}{p_{z}}\left<B_{y}B_{y}\right>=-\frac{\omega}{p_{z}}C^{\omega}=-\frac{p_{z}}{\omega}A^{\omega},$
(20)

where the third equality results from Eq. (18) and the last equality comes
from Eq. (19).

Thus, to conclude, we have the following expressions for the correlators of
$\vec{E}$ and $\vec{B}$:

$\begin{split}&\left<E_{i}E_{j}\right>=A(\omega,|\vec{p}|)\delta_{ij}+X(\omega,|\vec{p}|)\frac{p_{i}p_{j}}{p^{2}},\\\
&\left<B_{i}B_{j}\right>=\frac{p^{2}}{\omega^{2}}A(\omega,|\vec{p}|)\left\\{\delta_{ij}-\frac{p_{i}p_{j}}{p^{2}}\right\\},\\\
&\left<E_{i}B_{j}\right>=-A(\omega,|\vec{p}|)\epsilon_{ijk}\frac{p^{k}}{\omega}.\end{split}$
(21)

So far our considerations have been purely kinematical. We now turn to the
specific dynamics of the finite-temperature plasma; now the low-frequency
limits of these functions can be obtained from magnetohydrodynamics.

A general formulation of magnetohydrodynamics in terms of higher-form symmetry
was given in Grozdanov:2016tdf . In particular, the transverse channel
$A^{\omega}$ contains the physics of diffusion of magnetic field lines, and
the precise correlator needed was recorded in Hofman:2017vwr .

$\displaystyle
G^{R}_{J^{zx},J^{zx}}(\omega,p_{z})_{\text{MHD}}=A(\omega,|\vec{p}|)_{\text{MHD}}=\frac{-i\omega^{2}\rho}{\omega+iDp^{2}},$
(22)

Here $\rho$ is the resistivity of the plasma, and $D$ the diffusion constant
for magnetic field lines. It can be expressed in terms of the resistivity and
magnetic permeability555The magnetic permeability can formally be thought of
as the 1-form charge susceptibility, i.e. the thermodynamic quantity that
measures the amount of magnetic field created by an applied field. $\Xi$ of
the plasma as

$D=\frac{\rho}{\Xi}\ .$ (23)

The longitudinal channel $X$ appears only in the electric field channel and
controls the physics of Debye screening. It is not expected to have any
universal hydrodynamic structure, and is presumably analytic in frequency and
momenta at low frequencies. We will see explicitly that it does not contribute
at this order to the correlations of the topological density
$J_{\mu\nu}\tilde{J}^{\mu\nu}$.

###  II.2 Computation of correlator

Now we are set to compute the factors in the momentum space version of Eq.
(12). The first term in this factor which consists of same pairing correlators
is given as,

$\displaystyle\delta^{ij}\delta^{kl}\left[\langle
E_{i}E_{k}\rangle_{\omega_{1}}\,\langle B_{j}B_{l}\rangle_{\omega_{2}}\right]$
$\displaystyle=\frac{p^{2}}{\omega_{2}^{2}}\delta^{ij}\delta^{kl}\left[\left\\{A^{\omega_{1}}\delta_{ik}+X^{\omega_{1}}\frac{p_{i}p_{k}}{p^{2}}\right\\}\left\\{A^{\omega_{2}}\left(\delta_{jl}-\frac{p_{j}p_{l}}{p^{2}}\right)\right\\}\right]$
$\displaystyle=\frac{2p^{2}}{\omega_{2}^{2}}A^{\omega_{1}}A^{\omega_{2}},$
(24)

A similar computation for the second term, which consists of cross-
correlators, gives,

$\displaystyle\delta^{ij}\delta^{kl}\left[\left<E_{i}B_{l}\right>_{\omega_{1}}\left<B_{j}E_{k}\right>_{\omega_{2}}\right]$
$\displaystyle=\delta^{ij}\delta^{kl}\left\\{-A^{\omega_{1}}\epsilon_{ilm}\frac{p^{m}}{\omega_{1}}\right\\}\left\\{-A^{\omega_{2}}\epsilon_{kjn}\frac{p^{n}}{\omega_{2}}\right\\}$
$\displaystyle=A^{\omega_{1}}A^{\omega_{2}}\frac{p_{m}p^{n}}{\omega_{1}\omega_{2}}\epsilon^{jkm}\epsilon_{kjn}=-\frac{2p^{2}}{\omega_{1}\omega_{2}}A^{\omega_{1}}A^{\omega_{2}}.$
(25)

Now we can compute the product of spectral densities as given in Eq. (14),
using the above eqautions and Eq. (15).

$\rho_{E_{i}E_{j}}(\omega_{1})\rho_{B_{i}B_{j}}(\omega_{2})=\frac{1}{\pi^{2}}\frac{2p^{2}}{\omega_{2}^{2}}a^{\omega_{1}}a^{\omega_{2}},\qquad\rho_{E_{i}B_{j}}(\omega_{1})\rho_{B_{i}E_{j}}(\omega_{2})=-\frac{1}{\pi^{2}}\frac{2p^{2}}{\omega_{1}\omega_{2}}a^{\omega_{1}}a^{\omega_{2}}.$
(26)

where $a^{\omega}=\text{Im}\,A^{\omega}$.

###  II.3 Correlation of topological density

We now compute the correlator of the topological density (14). More precisely,
we will explicitly compute the following frequency-dependent quantity

$\Gamma_{A}(\Omega)=\frac{k^{2}}{\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}\Omega}\textrm{Im}\,G^{R}_{QQ}(\Omega,\vec{p}=0)$
(27)

Evaluated at $\Omega=0$ this determines the decay rate of the axial charge, as
explained in (5). However in this section we will compute its full frequency
dependence.

From Eq. (14) we have,

$\displaystyle\text{Im}\,G^{R}_{QQ}(\Omega)$
$\displaystyle=-64\int\limits_{-\infty}^{\infty}\frac{d^{3}p}{(2\pi)^{3}}\int\limits_{-\infty}^{\infty}\frac{d\omega_{1}}{2\pi}\frac{d\omega_{2}}{2\pi}\left[\frac{f(\omega_{1})-f(\omega_{2})}{\pi}\right]\delta(\omega_{1}-\omega_{2}-\Omega)\frac{2p^{2}}{\omega_{1}\omega_{2}^{2}}(\omega_{1}-\omega_{2})a^{\omega_{1}}a^{\omega_{2}},$
(28)

To obtain the above we used $\Omega_{l}=-i\Omega+\varepsilon$ to go to real
frequencies in Eq. (14) and used the identity,

$\displaystyle\text{Im}\left(\frac{1}{\omega_{1}-\omega_{2}-\Omega-i\varepsilon}\right)=\pi\delta(\omega_{1}-\omega_{2}-\Omega).$

Next we evaluate the $\omega_{1}$ integral and replace $\omega_{2}$ by
$\omega$ for notational simplicity. Then we have (see Eq. (5)):

$\displaystyle\Gamma_{A}(\Omega)$
$\displaystyle=-\frac{64k^{2}}{\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\frac{1}{\Omega}\left\\{\int\limits_{0}^{\infty}\frac{4\pi
p^{2}dp}{(2\pi)^{3}}\int\limits_{-\infty}^{\infty}\frac{d\omega}{4\pi^{3}}\left[f(\omega+\Omega)-f(\omega)\right]\frac{2p^{2}}{(\omega+\Omega)\omega^{2}}\Omega\,a^{\omega+\Omega}a^{\omega}\right\\},$
(29)

where we have changed to polar coordinates in momentum space using
$d^{3}p=4\pi p^{2}dp$.

Given an explicit expression for $a^{\omega}$, the above expression is in
principle exact (again, assuming only Gaussian correlations in the plasma). To
obtain an explicit answer, we now compute the contribution arising from
hydrodynamic fluctuations alone, i.e. we set $a^{\omega}$ equal to its results
from the MHD correlator from Eq. (22). In a given UV complete theory, this is
not the full answer, as $a^{\omega}$ will not agree with the MHD result at
high frequencies; however we expect that it should capture the dominant
infrared contribution. Plugging in the MHD value (22) for $a^{\omega}$, we
find,

$\displaystyle\Gamma_{A}(\Omega)=$
$\displaystyle\,-\frac{64k^{2}}{\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\frac{1}{\Omega}\left\\{\int\limits_{0}^{\infty}\frac{4\pi
p^{2}dp}{(2\pi)^{3}}\int\limits_{-\infty}^{\infty}\frac{d\omega}{4\pi^{3}}\left[f(\omega+\Omega)-f(\omega)\right]\frac{2p^{2}}{(\omega+\Omega)\omega^{2}}\Omega\right.$
$\displaystyle\times\left.\left(\frac{\omega^{2}\rho}{\omega^{2}+D^{2}p^{4}}\right)\left(\frac{(\omega+\Omega)^{2}\rho}{(\omega+\Omega)^{2}+D^{2}p^{4}}\right)\right\\}$
(30)

. Simplifying this result we find:

$\displaystyle\Gamma_{A}(\Omega)=\,-\frac{16k^{2}\rho^{2}}{\pi^{5}\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\int\limits_{-\infty}^{\infty}d\omega\left(\left[f(\omega+\Omega)-f(\omega)\right]\omega(\omega+\Omega)^{2}\int\limits_{0}^{\infty}dp\frac{p^{4}}{[(\omega+\Omega)^{2}+D^{2}p^{4}][\omega^{2}+D^{2}p^{4}]}\right).$
(31)

This integral is even in $\Omega$, as can be seen from the change of variables
$\omega\to-\omega$ and the identity $f(x)+f(-x)=-1$. This is expected from (5)
and the fact that the imaginary part of the retarded correlator of the bosonic
operator $Q$ is an odd function of frequency.

Next, evaluating the $p$-integral we find

$\displaystyle\Gamma_{A}(\Omega)=-\frac{16k^{2}\rho^{2}}{2\sqrt{2}\pi^{4}D^{\frac{5}{2}}\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\int\limits_{-\infty}^{\infty}d\omega\,\left[f(\omega+\Omega)-f(\omega)\right]\frac{\omega(\omega+\Omega)^{2}}{\Omega(2\omega+\Omega)}\left(\sqrt{|\omega+\Omega|}-\sqrt{|\omega|}\right).$
(32)

The remaining integral over $\omega$ is not trivial and displays interesting
structure in $\Omega$. For concreteness we first study the case $\Omega>0$.
Let us denote the integrand in (32) by $L(\omega,\Omega)$. We begin by noting
that the integrand is now non-analytic as a function of $\omega$ and
$\omega+\Omega$, arising from singularities in the integrand of (31) at small
$p$ when either of these frequencies vanish. We must thus separate the
$\omega$ integral into the following three ranges

$\omega\in(-\infty,-\Omega)\cup(-\Omega,0)\cup(0,\infty)\ .$ (33)

First let us perform the integration for $\omega\in(0,\infty)$, i.e. the
integral of interest is:

$\displaystyle
I^{+}(\Omega)=\int\limits_{0}^{\infty}d\omega\,\left[f(\omega+\Omega)-f(\omega)\right]\frac{\omega(\omega+\Omega)^{2}}{\Omega(2\omega+\Omega)}\left(\sqrt{\omega+\Omega}-\sqrt{\omega}\right).$
(34)

We wish to extract the dependence on $\Omega$ in the limit that
$\Omega\ll\beta$. This integral is convergent and can readily be done
numerically; however obtaining an analytical handle on the small $\Omega$
limit of this integral is subtle. To see this note that expansion of the
integrand in powers of $\Omega$ leads to an expression which is analytic in
$\Omega$. Naively, proceeding this way one is led to believe that $\Gamma_{A}$
is also analytic in $\Omega$. However, this is not the case; if we attempt to
proceed naively, the integral over each term in the Taylor expansion in
$\Omega$ fails to converge near $\omega=0$, indicating that the integral
itself not analytic as a function of $\Omega$ though the integrand is.

We thus need to carefully extract this non-analytic dependence of the decay
rate on $\Omega$. To do this we use the fact that there is a hierarchy of
scales $\Omega\ll\beta^{-1}$ to introduce a cut-off $\Lambda$ such that
$\Omega\ll\Lambda\ll\beta^{-1}$. With this, we can separate the integral in
Eq. (34) into an IR part $\omega\in(0,\Lambda)$ and a UV part
$\omega\in(\Lambda,\infty)$. The non-analytic dependence on $\Omega$ will come
from the IR part. At the end of the calculation we will show that nothing
depends on $\Lambda$.

To do the IR part, we expand the integrand about $\beta\to 0$, and work to all
orders in $\Omega$. We then integrate the resulting expansion for the range:
$\omega\in(0,\Lambda)$. We find the following result, which we have presented
in a series expansion in $\Omega$.

$\displaystyle I^{+}(\Omega)_{\text{IR}}$
$\displaystyle=\left(-\frac{\sqrt{\Lambda}}{2\beta}+\frac{\beta\Lambda^{5/2}}{120}-\frac{\beta^{3}\Lambda^{9/2}}{4320}+\mathcal{O}\left(\Lambda^{13/2}\right)\right)\Omega+\left(\frac{5}{6\beta}-\frac{\pi}{4\sqrt{2}\beta}-\frac{\sinh^{-1}(1)}{2\sqrt{2}\beta}\right)\Omega^{3/2}$
$\displaystyle+\left(\frac{1}{8\beta\sqrt{\Lambda}}+\frac{5\beta\Lambda^{3/2}}{288}-\frac{3\beta^{3}\Lambda^{7/2}}{4480}+\mathcal{O}\left(\Lambda^{11/2}\right)\right)\Omega^{2}$
$\displaystyle+\left(-\frac{1}{24\beta\Lambda^{3/2}}-\frac{19\beta^{3}\Lambda^{5/2}}{28800}+\frac{\beta^{5}\Lambda^{9/2}}{24192}+\mathcal{O}\left(\Lambda^{13/2}\right)\right)\Omega^{3}$
$\displaystyle+\left(\frac{139\beta}{10080}-\frac{\beta\pi}{192\sqrt{2}}-\frac{\beta\sinh^{-1}(1)}{96\sqrt{2}}\right)\Omega^{7/2}$
$\displaystyle+\left(\frac{11}{{640\beta\Lambda^{5/2}}}+\frac{5\beta}{{1536\sqrt{\Lambda}}}-\frac{13\beta^{3}\Lambda^{3/2}}{55296}+\mathcal{O}\left(\Lambda^{7/2}\right)\right)\Omega^{4}+\mathcal{O}(\Omega^{5}),$
(35)

Note the presence of a non-analytic series of terms starting at
${{\mathcal{O}}}(\Omega^{\frac{3}{2}})$. An interesting thing to note is that,
in the above equation, the analytic terms in ${{\mathcal{O}}}$ depend upon the
cut-off $\Lambda$, while the non-analytic pieces are independent of $\Lambda$,
and are exact expressions which are functions of $\beta$. This suggests that
the analytic pieces will receive contributions from the UV part while the non-
analytic pieces are exact. Also, in the above equation there seems to be IR
divergences in the analytic pieces about $\Lambda=0$. As we will see below
these will be cancelled from the UV part of the integral; of course the final
answer cannot depend on $\Lambda$.

We now turn to the UV part of the integral $\omega\in(\Lambda,\infty)$. Here
we can simply expand the integrand in powers of $\Omega$; the integral over
each term will converge, with any putative IR divergences cutoff by $\Lambda$.
So, the UV limit of integration is simpler; expanding the integrand we find:

$\displaystyle L^{+}(\omega,\Omega)\xrightarrow{\Omega\to
0}-\frac{e^{\beta\omega}\beta\omega^{3/2}}{4(1-e^{\beta\omega})^{2}}\Omega+\frac{\beta\sqrt{\omega}}{64}\left(2\beta\omega\coth\left(\frac{\beta\omega}{2}\right)-5\right)\text{csch}\left(\frac{\beta\omega}{2}\right)^{2}\Omega^{2}$
$\displaystyle-\frac{\beta^{2}\sqrt{\omega}}{768}\text{csch}\left(\frac{\beta\omega}{2}\right)^{4}\left(4\beta\omega\left(2+\cosh\left(\beta\omega\right)\right)-15\sinh\left(\beta\omega\right)\right)\Omega^{3}$
$\displaystyle+\frac{\beta\left(4\beta^{3}\omega^{3}\left(11\cosh\left(\frac{\beta\omega}{2}\right)+\cosh\left(\frac{3\beta\omega}{2}\right)\right)-5\left(3+16\beta^{2}\omega^{2}+\left(-3+8\beta^{2}\omega^{2}\right)\cosh\left(\beta\omega\right)\right)\sinh\left(\frac{\beta\omega}{2}\right)\right)}{192\,e^{\frac{-5\beta\omega}{2}}\left(1-e^{\beta\omega}\right)^{5}\omega^{3/2}}\Omega^{4}$
$\displaystyle+\mathcal{O}(\Omega^{5}).$ (36)

The linear in $\Omega$ piece in the above expansion is free of any IR
divergences and can be immediately integrated over the full range,
$\omega\to(0,\infty)$. The term in $\Omega^{2}$ is slightly more subtle;
integrating it over $\omega\in(\Lambda,\infty)$ we find that the leading
dependence on $\Lambda$ is $-\frac{1}{8\beta\Lambda^{1/2}}$, which indeed
precisely cancels the $\Lambda$-dependent divergent term in the quadratic
piece in Eq. (35). We have explicitly verified that a similar cancellation
takes place for the terms up to $\mathcal{O}(\Omega^{4})$, and on general
grounds it must happen to all orders in the $\Omega$ expansion. Next we can
numerically integrate the UV part, term by term in the range
$\omega\in(\Lambda,\infty)$ and finally take $\Lambda\to 0$.

Note that, the non-analytic pieces are controlled only by the IR integral
while the analytic pieces received contribution both from the IR and the UV
parts of the integral. Hence, in this sense, the non-analytic pieces are
universal.

We may treat the remaining two pieces of the integral in (33) in the same way;
we leave the details of these integrals to the Appendix C. Now gathering
everything together we find the $\Omega$-dependence of the integral for
$\Omega>0$ to be:

$\displaystyle\Gamma_{A}(\Omega)=\frac{16k^{2}\rho^{2}}{2\sqrt{2}\pi^{4}D^{\frac{5}{2}}\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\left\\{\frac{\pi}{2\sqrt{2}\beta}\Omega^{3/2}-\frac{0.3236}{\sqrt{\beta}}\Omega^{2}+\frac{\pi\beta}{96\sqrt{2}}\Omega^{7/2}-0.00518\,\beta^{3/2}\Omega^{4}+\mathcal{O}(\Omega^{6})\right\\},$
(37)

where the numerical coefficients of the analytic pieces are obtained by
numerical integration and hence are approximate values. On the contrary, the
numerical coefficients of the non-analytic pieces are exact values which do
not receive UV corrections, and their dependence on the IR data $\rho$ and $D$
is expected to be universal.

We now recall that our computation above assumed $\Omega>0$. Note that, due to
the non-analytic dependence on $\Omega$, strictly speaking we do not know the
behavior of the integral for $\Omega<0$, though the answer is constrained by
known transformation properties of the spectral density. We explicitly compute
the $1$-loop integral for $\Omega<0$ in Appendix C.2 and show that indeed
$\Gamma_{A}(-\Omega)=\Gamma_{A}(\Omega)$ as required. So, for
$\Omega\in\mathbb{R}$, the decay rate is given as,

$\displaystyle\Gamma_{A}(\Omega)=\frac{16k^{2}\rho^{2}}{2\sqrt{2}\pi^{4}D^{\frac{5}{2}}\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\left\\{\frac{\pi}{2\sqrt{2}\beta}|\Omega|^{3/2}-\frac{0.3236}{\sqrt{\beta}}\Omega^{2}+\frac{\pi\beta}{96\sqrt{2}}|\Omega|^{7/2}-0.00518\,\beta^{3/2}\Omega^{4}+\mathcal{O}(\Omega^{6})\right\\},$
(38)

which is now manifestly even.

Finally, we may note that Eq. (38) implies that $\Gamma_{A}(\Omega\to 0)=0$.
We see that the $1$-loop contribution to the decay rate itself vanishes in the
vanishing magnetic field limit. We discuss the implications of this result
further in the conclusion.

##  III Discussion

Above we presented an explicit calculation of the real-time topological
susceptibility – i.e. the retarded correlation function $G^{R}_{QQ}$ of the
operator $Q=8\vec{E}\cdot\vec{B}$ – arising from hydrodynamic fluctuations
about an equilibrium with vanishing magnetic field $\vec{B}=0$ in a
magnetohydrodynamic plasma.

In the presence of a finite $B$ field a classical calculation leads to
$\Gamma_{A}\sim\frac{B^{2}\rho}{\chi_{A}}$ at zero frequencies, as shown in
(7). The goal of this calculation was to determine the leading fluctuation-
induced contribution to the decay rate at vanishing $B$ field. If this is non-
zero, then strictly speaking in the infrared limit axial charge should not be
considered a hydrodynamic variable.

Importantly, however, we found after a calculation that the resulting
correlation function vanishes at low frequencies, as shown explicitly in (37).
This has the immediate implication that in a chiral plasma, the decay rate of
axial charge (as computed to one-loop order in hydrodynamics) remains zero if
background magnetic field is zero, i.e. the classical analysis here is
trustworthy. In our analysis we have also computed the first few terms in a
small frequency expansion in $\Omega$; it is interesting to note that the
presence of gapless diffusive modes results in a non-analytic dependence on
$\Omega$, though this dependence begins at
${{\mathcal{O}}}(\Omega^{\frac{3}{2}})$ and so is soft enough not to
contribute to the decay rate itself.

It is interesting to compare this to corresponding results for a non-Abelian
plasma, where the quantity analogous to (5) is the Chern-Simons diffusion
rate, i.e. the low-frequency limit of the correlation function of the non-
Abelian topological density $\mathop{\rm
Tr}(F^{a}_{\mu\nu}\tilde{F}^{a\mu\nu})$. This quantity has been extensively
studied both at weak-coupling Bodeker:1998hm ; Arnold:1996dy ; Huet:1996sh
and from holography Son:2002sd , and is certainly not zero.

A universal way to understand the difference between the Abelian and non-
Abelian case is the following: in the Abelian case studied here the
topological density in question can be understood as a bilinear in a conserved
2-form current
$J^{\mu\nu}={{\frac{1}{2}}}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$, i.e.
the anomaly equation (2) reads:

$\partial_{\mu}j^{\mu}_{A}=k\,\epsilon_{\mu\nu\rho\sigma}J^{\mu\nu}J^{\rho\sigma}$
(39)

This 2-form current is associated with the continuous $U(1)$ 1-form symmetry
that protects magnetic flux conservation in electrodynamics Gaiotto:2014kfa .
The presence of this continuous 1-form symmetry gives a great deal of extra
structure to this problem. This structure is not present for non-Abelian gauge
theory, where at most we have a discrete 1-form symmetry. At a calculational
level, it was the presence of this continuous symmetry current (and its
subsequent realization in thermal equilibrium) which allowed us to obtain non-
trivial constraints on the infrared physics from magnetohydrodynamics.

More generally, it has recently been shown that a precise characterization of
the anomaly (39) is possible in terms of non-invertible symmetries
Choi:2022jqy ; Cordova:2022ieu  (see Gomes:2023ahz ; McGreevy:2022oyu ;
Brennan:2023mmt ; Shao:2023gho ; Luo:2023ive ; Bhardwaj:2023kri ; Schafer-
Nameki:2023jdn  for reviews on this subject): i.e. there still exist
topological operators that count axial charge, but these operators no longer
obey a standard group composition law, and there is no longer a simple
conserved current. Our understanding of the dynamical consequences of such
non-invertible symmetries is still in its infancy. However it seems that one
way to understand the calculation above is that the finite temperature
dynamics of a charge density protected by a non-invertible symmetry is
somewhat constrained – for example, it will not relax to nothingness unless an
external magnetic field is applied. This result is philosophically consistent
with GarciaEtxebarria:2022jky , which showed that at zero temperature a form
of Goldstone’s theorem applies to such non-invertible symmetries in that there
is a protected gapless mode when the symmetry is spontaneously broken.

This suggests that the calculation above could be reorganized to make the role
played by the non-invertible symmetry more manifest. One possible way to do so
would be to use at one-loop level the effective theories constructed in
Das:2022fho ; Landry:2022nog , which realize the symmetries more directly. It
would be very interesting to obtain a robust argument for the vanishing of
this relaxation rate to all loop order in hydrodynamics. Another direction for
future work is to compare our results for $G^{R}_{QQ}(\omega)$ to real-time
lattice computations such as those in Figueroa:2017hun ; Figueroa:2019jsi ,
where one might hope that the non-analytic dependence on $\Omega$ in (37) –
which are in principle fully determined by hydrodynamic data – could be
verified from the lattice. We hope to return to this in the future.

Acknowledgments: We are extremely grateful to Luca Delacrétaz for suggesting
this calculation Luca:rep22 . We acknowledge helpful discussions with Adrien
Florio on related issues. We express our gratitude to Aristomenis Donos, Sašo
Grozdanov, Giorgio Frangi and Benjamin Withers for many fruitful discussions.
AD is supported by the STFC Consolidated Grant ST/T000600/1 – “Particle Theory
at the Higgs Centre”. NI is supported in part by STFC grant number
ST/T000708/1. The work of NP is supported by the grant for development of new
faculty (Ratchadapiseksomphot fund) and Sci-Super IX_66_004 from Chulalongkorn
University.

##  Appendix A Memory matrix review

For a system at finite temperature $T$ that undergoes a generic time
evolution, it is common to assume that most of the operators ${\cal O}$ will
decay away at late times as encoded by the retarded Green’s function
$G^{R}_{{\cal O}{\cal O}}(t)\sim\exp(-\Gamma_{\cal O}t)$. For a theory with a
hydrodynamic description, it is implicitly assumed that there are only a few
‘long-lived’ operators where $\Gamma_{\cal O}^{-1}$ is much larger than a
typical time scale (set by temperature and other macrocopic quantities),
namely the conserved charge densities. Often, it is assumed that all other,
non-conserved, operators have already decayed away by the late times at which
the hydrodynamic descriptions is applicable. However, it is sometimes possible
to have a situation where the lifetimes $\Gamma^{-1}_{\cal O}$ of non-
conserved operators are parametrically large enough such that they can
interfere with the hydrodynamic modes.

When there are only a few of these long-lived operators, the memory matrix
formalism forster2018hydrodynamic  (see also Hartnoll:2012rj ; Lucas:2015pxa ;
Grozdanov:2018fic  for a more recent discussion) is a powerful tools to
understand the correlation functions of such systems. In particular, it allows
one to extract the decay rate, or inverse-lifetime, of long-lived operators in
terms of simple 2-point correlation functions:

$\Gamma_{\cal O}=\frac{1}{\chi_{{\cal O}{\cal O}}}\lim_{\omega\to
0}\frac{1}{\omega}\text{Im}G^{R}_{\dot{\cal O}\dot{\cal O}}(\omega,\vec{p}=0)$
(40)

with $\dot{\cal O}\equiv\partial_{t}{\cal O}$ and $\chi_{{\cal O}{\cal O}}$ is
the susceptibility of the operator ${\cal O}$.

For completeness, we review the derivation of this formula. The starting point
is to realize that one can formally define an inner product between operators

$C_{AB}(t-t^{\prime})=(A(t)|B(t^{\prime})=(A|e^{-i(t-t^{\prime})L}|B)\equiv
T\int^{1/T}_{0}d\lambda\langle A(t)B(t^{\prime}+i\lambda)\rangle$ (41)

with $\langle...\rangle$ is the average over either quantum or thermal
fluctuation and $L=[H,\circ]$ is the Liouville operator. The relations between
$C_{AB}$, and its Laplace transformed $\tilde{C}_{AB}(\omega)$, can be shown
related to the standard correlation functions via

$C_{AB}(0)=(A|B)\equiv
T\chi_{AB}\,,\qquad\tilde{C}_{AB}(\omega)=(A|\frac{i}{\omega-L}|B)=\frac{T}{i\omega}\left(G^{R}_{AB}(\omega)-G^{R}_{AB}(0)\right)$
(42)

This translation to the inner product allows one to employ standard linear
algebra techniques to show that, for a long-lived operator ${\cal O}$

$\tilde{C}_{{\cal O}{\cal O}}(\omega)=\sum_{AB}\chi_{{\cal
O}A}\left(\frac{1}{-i\omega\chi+(N+M)}\right)_{AB}\chi_{B{\cal O}}$ (43)

where the sum is done over sets of long-lived operators $A,B$ and where the
matrix $N_{AB}=(A|\dot{B})$ vanishes in a theory with time-reversal symmetry.
The matrix $M_{AB}$ is the memory matrix defined via

$M_{AB}(\omega)=\frac{i}{T}(\dot{A}|\mathfrak{q}\frac{1}{\omega-\mathfrak{q}L\mathfrak{q}}|\dot{B})$
(44)

where $\mathfrak{q}$ is the projector that satisfies $\mathfrak{q}|A)=0$ when
$A$ is the long-lived operator and $\mathfrak{q}|A)=|A)$ when it is not. By
relating the Laplace-transformed $\tilde{C}_{AB}(\omega)$ in (43) to (42), one
can see that the decay rate encoded by a pole $G_{{\cal O}{\cal
O}}(\omega)\sim(i\omega-\Gamma_{\cal O})^{-1}$ can be obtained by
diagonalising the matrix $(\chi^{-1}M)_{AB}$.

This formalism can be readily applied to a system with ABJ anomaly, provided
that we assume the hydrodynamic description so that the energy $E$, momentum
$\vec{P}$ and the magnetic field $\vec{B}$ are the only conserved quantities
and that the axial charge density $n_{A}$ is long-lived. At zero magnetic
field and zero axial charge density, one can show that the only non-vanishing
‘overlap’ of these operators are

$(n_{A}|n_{A})=T\chi_{A}\,,\qquad(B^{i}|B^{j})=T\chi_{B}\delta^{ij}$ (45)

where $\chi_{A},\chi_{B}$ are the axial charge and magnetic susceptibility
respectively. Thus, we easily obtain the correlation of axial charge density

$\tilde{C}_{n_{A}n_{A}}(\omega)=\frac{\chi_{A}}{-i\omega+M_{n_{A}n_{A}}/\chi_{A}}$
(46)

It is clear that decay rate is controlled by the memory matrix
$M_{n_{A}n_{A}}(\omega\to 0)$ which can be written as

$M_{n_{A}n_{A}}(\omega\to 0)=\frac{1}{T}\lim_{\omega\to
0}(\dot{n}_{A}|\frac{i}{\omega-L}|\dot{n}_{A})=\frac{1}{T}\lim_{\omega\to
0}\tilde{C}_{\dot{n}_{A}\dot{n}_{A}}(\omega)=\lim_{\omega\to
0}\frac{1}{\omega}\text{Im}G^{R}_{\dot{n}_{A}\dot{n}_{A}}(\omega)$ (47)

where we use the fact that $\dot{n}_{A}=kQ=8k(\vec{E}\cdot\vec{B})$ is not a
long-lived operator (due to the r.h.s. containing the non-conserved electric
field $\vec{E}$) and thus $\mathfrak{q}|\dot{n}_{A})=|\dot{n}_{A})$. Upon
converting (46) back to the retarded correlation function via (42) and extract
the decay rate $\Gamma_{A}$ from the pole
$G^{R}_{n_{A}n_{A}}(\omega)\sim(\omega-\Gamma_{A})^{-1}$, we are then able to
express the decay rate in terms of the correlation function of $Q$ as in (5).

##  Appendix B Finite temperature conventions and Matsubara sums

Here we collect some identities that are useful for performing the frequency
integrals in the main text. All of these results are standard, and further
background can be found e.g. in Bellac:2011kqa .

Consider a quantum field theory with a bosonic Hermitian operator
${{\mathcal{O}}}(t,\vec{x})$. We study the theory in the thermal state with
temperature $\beta^{-1}$. There are various basic two-point functions for
${{\mathcal{O}}}$, including the Euclidean correlation function,

$G^{E}(\tau,\vec{x})\equiv\langle{{\mathcal{O}}}(\tau,\vec{x}){{\mathcal{O}}}(0)\rangle,$
(48)

and the retarded real-time correlation function

$G^{R}(t,\vec{x})=-i\theta(t)\mathop{\rm Tr}\left(e^{-\beta
H}[{{\mathcal{O}}}(t,\vec{x}),{{\mathcal{O}}}(0)]\right)\ .$ (49)

The Euclidean correlation function in frequency space can be written in terms
of the spectral density $\rho(\Omega)$:

$G^{E}(i\omega_{n})=\int\frac{d\Omega}{2\pi}\frac{\rho(\Omega)}{i\omega_{n}+\Omega}$
(50)

We may also obtain the retarded correlator from the Euclidean one by
evaluating the latter at a real frequency:

$G^{R}(\Omega)=G^{E}(i\omega_{l}=\Omega+i\epsilon)$ (51)

Inserting (50) into (51) and using the identity

$\textrm{Im}\,\left(\frac{1}{x-i\epsilon}\right)=\pi\delta(x)$ (52)

we conclude that the imaginary part of $G^{R}(\omega)$ directly measures the
spectral density.

$\textrm{Im}\,G^{R}(\omega)=-\pi\rho(\omega)$ (53)

It is shown in Bellac:2011kqa  that for a bosonic operator $\rho(\omega)$ is
an odd function of $\omega$, and furthermore is positive for positive
$\omega$, i.e. $\omega\rho(\omega)>0$.

###  B.1 Performing Matsubara sums

We will need to perform a loop sum over Euclidean frequencies. Here we review
a standard trick to express such sums in terms of the corresponding spectral
densities, following the discussion in Faulkner:2013bna . Consider summing
over a set of discrete Matsubara frequencies $i\Omega_{m}=\frac{2\pi
m}{\beta}$, $m\in\mathbb{Z}$. We can express this in terms of a contour
integral over a contour $C$ in the complex $\omega$ plane, i.e.

$T\sum_{i\omega_{m}}\to\frac{1}{2\pi
i}\int_{C}d\omega{{\frac{1}{2}}}\coth\left(\frac{\beta\omega}{2}\right)$ (54)

Here the hyperbolic function in the integrand has poles at each of the
discrete Matsubara frequencies along the imaginary axis, and the contour $C$
is a series of disjoint circles which encircles each of these poles.

To see an application of this, consider evaluating the following sum, where
$i\Omega_{l}$ is a Matsubara frequency and $\omega_{1,2}$ are two real
frequencies:

$S(i\Omega_{l},\omega_{1},\omega_{2})=T\sum_{i\omega_{m}}\frac{1}{i(\omega_{m}+\Omega_{l})-\omega_{1}}\frac{1}{i\omega_{m}-\omega_{2}}$
(55)

From above we see that this can be written as the following contour integral:

$S(i\Omega_{l},\omega_{1},\omega_{2})=\frac{1}{2\pi
i}\int_{C}d\omega{{\frac{1}{2}}}\coth\left(\frac{\beta\omega}{2}\right)\frac{1}{\omega+i\Omega_{l}-\omega_{1}}\frac{1}{\omega-\omega_{2}}$
(56)

Now consider deforming the contour $C$ into two parallel lines, one running
down the imaginary $\omega$ axis at infinitesimal real positive $\omega$ and
the other running up the imaginary $\omega$ axis at infinitesimal real
negative $\omega$. We can now attempt to deform these lines away to infinity.
At large $|\omega|$ the integrand behaves as $|\omega|^{-2}$. The contribution
to the integrand at infinity can be neglected, and the full integral arises
from the contribution at the non-Matsubara poles of the integrand, which
appear only at $\omega=\omega_{1}-i\Omega_{l}$ and $\omega=\omega_{2}$.
Performing the integral by residues we find

$S(i\Omega_{l},\omega_{1},\omega_{2})=-{{\frac{1}{2}}}\left(\coth\left(\frac{\beta(\omega_{1}-i\Omega_{l})}{2}\right)-\coth\left(\frac{\beta\omega_{2}}{2}\right)\right)\frac{1}{\omega_{1}-i\Omega_{l}-\omega_{2}}$
(57)

which after some algebra can be seen to be equal to

$S(i\Omega_{l},\omega_{1},\omega_{2})=-\frac{f(\omega_{1})-f(\omega_{2})}{\omega_{1}-i\Omega_{l}-\omega_{2}}$
(58)

where we have used the fact that $e^{i\beta\Omega_{l}}=1$ on a Matsubara
frequency. Here $f(\omega)$ is the Bose distribution function:

$f(\omega)=\frac{1}{e^{\beta\omega}-1}$ (59)

Let us now use this form to perform a frequency sum. In the bulk of the text
we will find ourselves needing to calculate sums of the form

$F(i\Omega_{l})=T\sum_{i\omega_{m}}G_{1}^{E}(i\Omega_{l}+i\omega_{m})G_{2}^{E}(i\omega_{m})$
(60)

where here $G^{E}_{1,2}(i\Omega)$ are two (possibly different) Euclidean
propagators. It is very convenient to express this in terms of the spectral
densities $\rho_{1,2}(\omega)$ associated with these propagators. To do this,
we first use (50) and then perform the sum over $i\omega_{m}$ using (58) to
find

$F(i\Omega_{l})=-\int\frac{d\omega_{1}}{2\pi}\frac{d\omega_{2}}{2\pi}\left(f(\omega_{1})-f(\omega_{2})\right)\frac{\rho_{1}(\omega_{1})\rho_{2}(\omega_{2})}{\omega_{1}-i\Omega_{l}-\omega_{2}}\
.$ (61)

This expression is used to obtain (14) in the main text.

##  Appendix C Details of the 1-loop integration

###  C.1 $\Omega>0$

Here we give details of the remaining parts of the $1$-loop integration with
$\Omega>0$.

As described around (33), the frequency integral in (32) must be split up into
three parts

$\omega\in(-\infty,-\Omega)\cup(-\Omega,0)\cup(0,\infty)\ .$ (62)

The last integral was performed in detail in the bulk of the text. In this
Appendix we perform the other two using the same methods. We begin with
$\omega\in(-\Omega,0)$. To do this integral let us perform the following
change of variable: $\omega\to-\omega$ and then perform the integration over
the positive range: $\omega\in(0,\Omega)$. The integral of interest is:

$\displaystyle
I^{-}_{1}(\Omega)=\int\limits_{0}^{\Omega}d\omega\,\left[f(\Omega-\omega)-f(-\omega)\right]\frac{\omega(\Omega-\omega)^{2}}{\Omega(2\omega-\Omega)}\left(\sqrt{\Omega-\omega}-\sqrt{\omega}\right),$
(63)

To do the above integral, we expand the integrand about $\beta\to 0$ and then
integrate term by term. We get,

$\displaystyle I^{-}_{1}(\Omega)$
$\displaystyle=\frac{(-2+\sqrt{2}\sinh^{-1}(1))}{2\beta}\Omega^{3/2}+\frac{\beta(-26+15\sqrt{2}\sinh^{-1}(1))}{1440}\Omega^{7/2}$
$\displaystyle+\frac{\beta^{3}(214-105\sqrt{2}\sinh^{-1}(1))}{2419200}\Omega^{11/2}.$
(64)

As expected the above integral is non-analytic in $\Omega$ and these pieces,
as discussed before, do not receive any UV corrections.

Next let us move on to performing the integration over the range:
$\omega\in(-\infty,-\Omega)$. As before let us do the variable change:
$\omega\to-\omega$ and then the integration has to be performed over the
range: $\omega\in(\Omega,\infty)$. The integral of interest is:

$\displaystyle
I^{-}_{2}(\Omega)=\int\limits_{\Omega}^{\infty}d\omega\,\left[f(\Omega-\omega)-f(-\omega)\right]\frac{\omega(\Omega-\omega)^{2}}{\Omega(2\omega-\Omega)}\left(\sqrt{\omega-\Omega}-\sqrt{\omega}\right),$
(65)

We can perform the above integration using the methods employed to do the
integration in Eq. (34). We obtain the following results.

$\displaystyle I^{-}_{2}(\Omega)_{\text{IR}}$
$\displaystyle=\left(\frac{\sqrt{\Lambda}}{2\beta}-\frac{\beta\Lambda^{5/2}}{120}+\frac{\beta^{3}\Lambda^{9/2}}{4320}+\mathcal{O}\left(\Lambda^{13/2}\right)\right)\Omega+\left(\frac{1}{6\beta}-\frac{\pi}{4\sqrt{2}\beta}-\frac{\sinh^{-1}(1)}{2\sqrt{2}\beta}\right)\Omega^{3/2}$
$\displaystyle+\left(\frac{1}{8\beta\sqrt{\Lambda}}+\frac{5\beta\Lambda^{3/2}}{288}-\frac{3\beta^{3}\Lambda^{7/2}}{4480}+\mathcal{O}\left(\Lambda^{11/2}\right)\right)\Omega^{2}$
$\displaystyle+\left(\frac{1}{24\beta\Lambda^{3/2}}+\frac{19\beta^{3}\Lambda^{5/2}}{28800}-\frac{\beta^{5}\Lambda^{9/2}}{24192}+\mathcal{O}\left(\Lambda^{13/2}\right)\right)\Omega^{3}$
$\displaystyle+\left(\frac{43\beta}{10080}-\frac{\beta\pi}{192\sqrt{2}}-\frac{\beta\sinh^{-1}(1)}{96\sqrt{2}}\right)\Omega^{7/2}$
$\displaystyle+\left(\frac{11}{{640\beta\Lambda^{5/2}}}+\frac{5\beta}{{1536\sqrt{\Lambda}}}-\frac{13\beta^{3}\Lambda^{3/2}}{55296}+\mathcal{O}\left(\Lambda^{7/2}\right)\right)\Omega^{4}+\mathcal{O}(\Omega^{5}).$
(66)

Comparing Eq. (66) with Eq. (35) we see that, for the analytic pieces: the
terms with odd powers of $\Omega$ are of opposite signs and the terms with
even powers of $\Omega$ have the same sign.

###  C.2 $\Omega<0$

Here we work out the $1$-loop integration with $\Omega<0$. For simplicity, let
us define $t=-\Omega$ with $t>0$. From Eq. (32) we get,

$\displaystyle\Gamma_{A}(-t)=-\frac{16k^{2}\rho^{2}}{2\sqrt{2}\pi^{4}D^{\frac{5}{2}}\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\int\limits_{-\infty}^{\infty}d\omega\,\left[f(\omega-t)-f(\omega)\right]\frac{\omega(\omega-t)^{2}}{(-t)(2\omega-t)}\left(\sqrt{|\omega-t|}-\sqrt{|\omega|}\right).$
(67)

From the structure of the square-root above we see that the integral should be
integrated over the following intervals,

$\displaystyle\omega\in(-\infty,0)\cup(0,t)\cup(t,\infty)$

Contrast this to the intervals in the $\Omega>0$ case. Both are mirror images
of each other. Now let us do the integral over the range,
$\omega\in(-\infty,0)$. This integral, after a change of variable:
$\omega\to-\omega$ becomes,

$\begin{split}\Gamma_{A}(-t)&=\frac{16k^{2}\rho^{2}}{2\sqrt{2}\pi^{4}D^{\frac{5}{2}}\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\int\limits_{0}^{\infty}d\omega\,\left[f(-\omega-t)-f(-\omega)\right]\frac{\omega(\omega+t)^{2}}{t(2\omega+t)}\left(\sqrt{\omega+t}-\sqrt{\omega}\right)\\\
&=-\frac{16k^{2}\rho^{2}}{2\sqrt{2}\pi^{4}D^{\frac{5}{2}}\raisebox{0.0pt}[2.79857pt][2.79857pt]{$\chi$}_{A}}\int\limits_{0}^{\infty}d\omega\,\left[f(\omega+t)-f(\omega)\right]\frac{\omega(\omega+t)^{2}}{t(2\omega+t)}\left(\sqrt{\omega+t}-\sqrt{\omega}\right)\end{split}$
(68)

where to get to the second equality we used the identity: $f(-x)=-1-f(x)$.
Note that the above integral is the same as the integral in Eq. (34).

Similarly, using the above identity, one can show that $\Gamma_{A}(-t)$ for
$\omega\in(0,t)$ matches with the integral in Eq. (63) and $\Gamma_{A}(-t)$
for $\omega\in(t,\infty)$ matches with the integral in Eq. (65). Thus, we find
$\Gamma_{A}(\Omega)=\Gamma_{A}(-\Omega)$, as expected.

## References

  * (1) M. Joyce and M. E. Shaposhnikov, “Primordial magnetic fields, right-handed electrons, and the Abelian anomaly,” Phys. Rev. Lett. 79 (1997) 1193–1196, arXiv:astro-ph/9703005.
  * (2) I. Rogachevskii, O. Ruchayskiy, A. Boyarsky, J. Fröhlich, N. Kleeorin, A. Brandenburg, and J. Schober, “Laminar and turbulent dynamos in chiral magnetohydrodynamics-I: Theory,” Astrophys. J. 846 no. 2, (2017) 153, arXiv:1705.00378 [physics.plasm-ph].
  * (3) A. Boyarsky, J. Frohlich, and O. Ruchayskiy, “Magnetohydrodynamics of Chiral Relativistic Fluids,” Phys. Rev. D 92 (2015) 043004, arXiv:1504.04854 [hep-ph].
  * (4) L. Del Zanna and N. Bucciatini, “Covariant and 3 + 1 equations for dynamo-chiral general relativistic magnetohydrodynamics,” Mon. Not. Roy. Astron. Soc. 479 no. 1, (2018) 657–666, arXiv:1806.07114 [astro-ph.HE].
  * (5) K. Hattori, Y. Hirono, H.-U. Yee, and Y. Yin, “MagnetoHydrodynamics with chiral anomaly: phases of collective excitations and instabilities,” Phys. Rev. D 100 no. 6, (2019) 065023, arXiv:1711.08450 [hep-th].
  * (6) N. Yamamoto, “Scaling laws in chiral hydrodynamic turbulence,” Phys. Rev. D 93 no. 12, (2016) 125016, arXiv:1603.08864 [hep-th].
  * (7) A. Das, R. Gregory, and N. Iqbal, “Higher-form symmetries, anomalous magnetohydrodynamics, and holography,” SciPost Phys. 14 (2023) 163, arXiv:2205.03619 [hep-th].
  * (8) D. T. Son and P. Surowka, “Hydrodynamics with Triangle Anomalies,” Phys. Rev. Lett. 103 (2009) 191601, arXiv:0906.5044 [hep-th].
  * (9) Y. Neiman and Y. Oz, “Relativistic Hydrodynamics with General Anomalous Charges,” JHEP 03 (2011) 023, arXiv:1011.5107 [hep-th].
  * (10) K. Landsteiner, “Notes on Anomaly Induced Transport,” Acta Phys. Polon. B47 (2016) 2617, arXiv:1610.04413 [hep-th].
  * (11) L. Delacrétaz, “Report on arXiv:2205.03619v2, delivered 2022-09-12, doi:10.21468/SciPost.Report.5679,”.
  * (12) A. Das, N. Iqbal, and N. Poovuttikul, “Towards an effective action for chiral magnetohydrodynamics,” arXiv:2212.09787 [hep-th].
  * (13) S. Grozdanov, D. M. Hofman, and N. Iqbal, “Generalized global symmetries and dissipative magnetohydrodynamics,” Phys. Rev. D95 no. 9, (2017) 096003, arXiv:1610.07392 [hep-th].
  * (14) D. G. Figueroa and M. Shaposhnikov, “Anomalous non-conservation of fermion/chiral number in Abelian gauge theories at finite temperature,” JHEP 04 no. 04, (2018) 026, arXiv:1707.09967 [hep-ph]. [Erratum: JHEP 07, 217 (2020)].
  * (15) D. G. Figueroa, A. Florio, and M. Shaposhnikov, “Chiral charge dynamics in Abelian gauge theories at finite temperature,” JHEP 10 (2019) 142, arXiv:1904.11892 [hep-th].
  * (16) A. Das, A. Florio, N. Iqbal, and N. Poovuttikul, “Higher-form symmetry and chiral transport in real-time lattice $U(1)$ gauge theory,” arXiv:2309.14438 [hep-th].
  * (17) M. J. Landry and H. Liu, “A systematic formulation of chiral anomalous magnetohydrodynamics,” arXiv:2212.09757 [hep-ph].
  * (18) P. B. Arnold and L. D. McLerran, “The Sphaleron Strikes Back,” Phys. Rev. D 37 (1988) 1020.
  * (19) A. Das, Generalised symmetries, anomalous magnetohydrodynamics and holography.  PhD thesis, Durham U., Durham University, 8, 2023.
  * (20) P. Kovtun, “Lectures on hydrodynamic fluctuations in relativistic theories,” J. Phys. A45 (2012) 473001, arXiv:1205.5040 [hep-th].
  * (21) G. D. Mahan, Many-particle physics.  Springer Science & Business Media, 2000.
  * (22) T. Faulkner, N. Iqbal, H. Liu, J. McGreevy, and D. Vegh, “Charge transport by holographic Fermi surfaces,” Phys. Rev. D 88 (2013) 045016, arXiv:1306.6396 [hep-th].
  * (23) D. M. Hofman and N. Iqbal, “Generalized global symmetries and holography,” SciPost Phys. 4 (2018) 005, arXiv:1707.08577 [hep-th].
  * (24) D. Bodeker, “On the effective dynamics of soft nonAbelian gauge fields at finite temperature,” Phys. Lett. B 426 (1998) 351–360, arXiv:hep-ph/9801430.
  * (25) P. B. Arnold, D. Son, and L. G. Yaffe, “The Hot baryon violation rate is O (alpha-w**5 T**4),” Phys. Rev. D 55 (1997) 6264–6273, arXiv:hep-ph/9609481.
  * (26) P. Huet and D. T. Son, “Long range physics in a hot nonAbelian plasma,” Phys. Lett. B 393 (1997) 94–102, arXiv:hep-ph/9610259.
  * (27) D. T. Son and A. O. Starinets, “Minkowski space correlators in AdS / CFT correspondence: Recipe and applications,” JHEP 09 (2002) 042, arXiv:hep-th/0205051.
  * (28) D. Gaiotto, A. Kapustin, N. Seiberg, and B. Willett, “Generalized Global Symmetries,” JHEP 02 (2015) 172, arXiv:1412.5148 [hep-th].
  * (29) Y. Choi, H. T. Lam, and S.-H. Shao, “Non-invertible Global Symmetries in the Standard Model,” arXiv:2205.05086 [hep-th].
  * (30) C. Cordova and K. Ohmori, “Non-Invertible Chiral Symmetry and Exponential Hierarchies,” arXiv:2205.06243 [hep-th].
  * (31) P. R. S. Gomes, “An Introduction to Higher-Form Symmetries,” arXiv:2303.01817 [hep-th].
  * (32) J. McGreevy, “Generalized Symmetries in Condensed Matter,” arXiv:2204.03045 [cond-mat.str-el].
  * (33) T. D. Brennan and S. Hong, “Introduction to Generalized Global Symmetries in QFT and Particle Physics,” arXiv:2306.00912 [hep-ph].
  * (34) S.-H. Shao, “What’s Done Cannot Be Undone: TASI Lectures on Non-Invertible Symmetry,” arXiv:2308.00747 [hep-th].
  * (35) R. Luo, Q.-R. Wang, and Y.-N. Wang, “Lecture Notes on Generalized Symmetries and Applications,”  7, 2023.  arXiv:2307.09215 [hep-th].
  * (36) L. Bhardwaj, L. E. Bottini, L. Fraser-Taliente, L. Gladden, D. S. W. Gould, A. Platschorre, and H. Tillim, “Lectures on Generalized Symmetries,” arXiv:2307.07547 [hep-th].
  * (37) S. Schafer-Nameki, “ICTP Lectures on (Non-)Invertible Generalized Symmetries,” arXiv:2305.18296 [hep-th].
  * (38) I. García-Etxebarria and N. Iqbal, “A Goldstone theorem for continuous non-invertible symmetries,” arXiv:2211.09570 [hep-th].
  * (39) D. Forster, Hydrodynamic Fluctuations, Broken Symmetry, And Correlation Functions.  CRC Press, 2018.  https://books.google.co.th/books?id=sl0PEAAAQBAJ.
  * (40) S. A. Hartnoll and D. M. Hofman, “Locally Critical Resistivities from Umklapp Scattering,” Phys. Rev. Lett. 108 (2012) 241601, arXiv:1201.3917 [hep-th].
  * (41) A. Lucas and S. Sachdev, “Memory matrix theory of magnetotransport in strange metals,” Phys. Rev. B 91 no. 19, (2015) 195122, arXiv:1502.04704 [cond-mat.str-el].
  * (42) S. Grozdanov, A. Lucas, and N. Poovuttikul, “Holography and hydrodynamics with weakly broken symmetries,” arXiv:1810.10016 [hep-th].
  * (43) M. L. Bellac, Thermal Field Theory.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 3, 2011.