# New developments in relativistic magnetohydrodynamics

Koichi Hattori<EMAIL_ADDRESS>Zhejiang Institute of Modern
Physics, Department of Physics, Zhejiang University, Hangzhou, 310027, China
Research Center for Nuclear Physics, Osaka University, Osaka 567-0047 Japan.
Masaru Hongo<EMAIL_ADDRESS>Department of Physics, Niigata
University, Niigata 950-2181, Japan  RIKEN iTHEMS, RIKEN, Wako 351-0198, Japan
Xu-Guang Huang<EMAIL_ADDRESS>Physics Department and Center for
Field Theory and Particle Physics, Fudan University, Shanghai 200438, China
Key Laboratory of Nuclear Physics and Ion-beam Application (MOE), Fudan
University, Shanghai 200433, China  Shanghai Research Center for Theoretical
Nuclear Physics, National Science Foundation of China and Fudan University,
Shanghai 200438, China

###### Abstract

Relativistic magnetohydrodynamics (RMHD) provides an extremely useful
description of the low-energy long-wavelength phenomena in a variety of
physical systems from quark-gluon plasma in heavy-ion collisions to matters in
supernova, compact stars, and early universe. We review the recent theoretical
progresses of RMHD, such as a formulation of RMHD from the perspective of
magnetic flux conservation using the entropy-current analysis, the
nonequilibrium statistical operator approach applied to quantum
electrodynamics, and the relativistic kinetic theory. We discuss how the
transport coefficients in RMHD are computed in kinetic theory and perturbative
quantum field theories. We also explore the collective modes and instabilities
in RMHD with a special emphasis on the role of chirality in a parity-odd
plasma. We also give some future prospects of RMHD, including the interaction
with spin hydrodynamics and the new kinetic framework with magnetic flux
conservation.

###### Contents

  1. I Introduction
  2. II Macroscopic approach: the entropy-current analysis
    1. II.1 Primer to the entropy-current analysis
    2. II.2 Relativistic MHD from the magnetic flux conservation
      1. II.2.1 Entropy production rate with magnetic flux
      2. II.2.2 Zeroth-order in derivatives: Nondissipative RMHD
      3. II.2.3 First-order derivative corrections: Dissipative RMHD
  3. III Nonequilibrium statistical operator method for relativistic MHD
    1. III.1 Optimized perturbation with local Gibbs distribution
    2. III.2 Evaluating the Local Gibbs averages
    3. III.3 Evaluating the dissipative corrections
  4. IV Interlude: Connection to the conventional MHD
    1. IV.1 Anisotropic pressure
    2. IV.2 First-order constitutive relations including Hall transports
  5. V Towards relativistic MHD from kinetic equations
    1. V.1 Chapman-Enskog method
    2. V.2 Grad’s method of moments
  6. VI Perturbative evaluation of transport coefficients
    1. VI.1 Transport coefficients in strong magnetic fields
    2. VI.2 Roles of the chiral symmetry
    3. VI.3 Contributions of higher Landau levels
  7. VII Chiral magnetohydrodynamics
    1. VII.1 Entropy-current analysis with chiral anomaly
    2. VII.2 Linear waves in relativistic MHD
    3. VII.3 Helical instabilities in the chiral MHD
  8. VIII Summary and future prospects
  9. A Electric charges in MHD
  10. B Matching and frame conditions in relativistic MHD
  11. C Derivation of Eqs. (73a)-(73c)
  12. D Viscous coefficients and inequalities

##  I Introduction

Relativistic hydrodynamics provides an incredibly successful description of
the macroscopic dynamics of interacting many-body systems in relativistic
arena since its birth Eckart (1940); Kluitenberg  _et al._ (1953a, b);
Kluitenberg and De Groot (1954a); Landau and Lifshitz (1959). Its
applicability ranges from very small systems such as the quark gluon plasma
(QGP) created in high-energy heavy-ion collisions to very large systems like
the expanding universe and the explosive supernovas Yagi  _et al._ (2008);
Romatschke and Romatschke (2019). From a theoretical point of view,
relativistic hydrodynamics is a typical example of effective field theories
which is valid at low-energy and long-wavelength limit. The dynamical
variables in relativistic hydrodynamics are the coarse-grained conserved
quantities stemmed from the underlying symmetries. As such, there exist a
variety of ways to systematically construct relativistic hydrodynamics when an
appropriate derivative expansion is employed.

Strong magnetic fields exist and play critical roles in a number of systems
where the relativistic hydrodynamics can be applied. Examples range from high-
energy heavy-ion collisions to supernovas, neutron stars, and the early
universe. In heavy-ion collisions, the colliding nuclei induce very strong
transient magnetic fields exerting on the produced quark-gluon plasma (QGP)
Skokov  _et al._ (2009); Voronyuk  _et al._ (2011); Bzdak and Skokov (2012);
Deng and Huang (2012); Bloczynski  _et al._ (2013); Tuchin (2013a); Yan and
Huang (2021); Wang  _et al._ (2022). The peak strength of such magnetic fields
can reach $10^{19}$ Gauss in Au + Au collisions at Relativistic Heavy Ion
Collider (RHIC) and $10^{20}$ Gauss in Pb + Pb collisions at Large Hadron
Collider (LHC). Given that the QGP is also an electromagnetic (EM) plasma,
such magnetic fields can strongly influence the dynamics of QGP and induce
novel transport phenomena in QGP such as the chiral magnetic effect which
provides a valuable machinery to access the quantum chromodynamics (QCD)
topological sector in an experimentally feasible way Kharzeev  _et al._
(2008); Fukushima  _et al._ (2008). (For reviews of strong magnetic-field
effects in heavy-ion collisions, see Refs. Miransky and Shovkovy (2015); Huang
(2016); Hattori and Huang (2017); Kharzeev  _et al._ (2016); Kharzeev and Liao
(2021).)

In compact stellar objects such as the neutron stars, the surface magnetic
fields can reach the order of $10^{12}–10^{13}$ Gauss Shapiro and Teukolsky
(1983), with a subclass of neutron stars (called magnetars) having surface
magnetic fields of the order of $10^{14}–10^{15}$ Gauss Duncan and Thompson
(1992) (Even stronger transient magnetic fields may be created in binary
neutron star mergers Price and Rosswog (2006); Kiuchi  _et al._ (2015)). The
internal magnetic fields could be by orders of magnitude stronger than the
surface magnetic fields. Although the origin of magnetar magnetic fields is
not fully understood, one widely accepted theory is that they result from the
magnetohydrodynamic dynamo processes in the interior fluids of some newly born
neutron stars Thompson and Duncan (1993). Observationally, soft gamma ray
repeaters (SGRs) and anomalous X-ray pulsars (AXPs) are commonly identified as
magnetars. The strong magnetic fields are crucial for understanding the
phenomenology of neutron stars, particularly magnetars, such as radio
emission, cooling properties, equations of state, shape deformation and
gravitational-wave emission, merging processes and post-merger evolution of
binary neutron stars, among other things. (See Refs. Lai (2001); Harding and
Lai (2006); Mereghetti  _et al._ (2015); Turolla  _et al._ (2015); Kaspi and
Beloborodov (2017); Enoto  _et al._ (2019); Ciolfi (2020) for reviews.)

Although a conclusive evidence for primordial magnetic fields in the early
universe is still lacking, astrophysical studies of the large-scale
intergalactic magnetic fields in the current universe strongly support their
existence. The origin of the primordial magnetic fields has been pursued over
the last three decades with possible scenarios stemming from the big-bang era
Turner and Widrow (1988); Carroll  _et al._ (1990); Garretson  _et al._ (1992)
and cosmic phase transitions like the QCD phase transition and the electroweak
phase transition Hogan (1983); Quashnock  _et al._ (1989); Vachaspati (1991);
Cheng and Olinto (1994); Baym  _et al._ (1996); Son (1999). Besides, the
strong primordial magnetic fields may have played an important role in
understanding the cosmological structure formation, the thermal spectrum and
polarization anisotropies of the cosmic microwave background, the big-bang
nucleosynthesis and baryogenesis, and so on. (For reviews of the primordial
magnetic fields, see Refs. Grasso and Rubinstein (2001); Kandus  _et al._
(2011); Subramanian (2016).)

The relativistic magnetohydrodynamics (relativistic MHD or RMHD) is often used
as a standard tool to analyze the physical processes in the systems mentioned
above. It provides a macroscopic framework to self-consistently describe the
evolution of the matter coupled with either dynamical electromagnetic (EM)
fields or in external EM fields. There is a long history of astrophysical and
cosmological applications of RMHD (typically in the presence of gravity), and
there are already a number of excellent reviews and textbooks (see, for
example, books Lichnerowicz (1994); Anile (2005); Rezzolla and Zanotti (2013);
Goedbloed  _et al._ (2019); Kato and Fukue (2020)). Over the last decade, due
to the realization of strong magnetic fields in heavy ion collisions, various
aspects and applications of RMHD in the context of heavy ion collisions have
been extensively investigated, including the study of evolution of magnetic
fields in the QGP Tuchin (2013b); Li  _et al._ (2016), the computation of RMHD
transport coefficients in QCD matter using perturbative field theory Hattori
_et al._ (2017a); Hattori and Satow (2016); Hattori  _et al._ (2017b, c); Li
and Yee (2018), kinetic theories Denicol  _et al._ (2018, 2019); Dey  _et al._
(2021a); Chen  _et al._ (2020); Dash  _et al._ (2020); Singh  _et al._ (2020);
Panda  _et al._ (2021a); Ghosh and Haque (2022), and holographic models
Critelli  _et al._ (2014); Finazzo  _et al._ (2016); Li  _et al._ (2018);
Fukushima and Okutsu (2022), the anisotropic evolution of the QGP coupled with
magnetic fields Tuchin (2012); Gursoy  _et al._ (2014); Roy  _et al._ (2015);
Pu and Yang (2016); Pu  _et al._ (2016); Inghirami  _et al._ (2016); Gürsoy
_et al._ (2018); She  _et al._ (2019); Inghirami  _et al._ (2020); Kord  _et
al._ (2022); Emamian  _et al._ (2020), the simulation of chiral magnetic
effect and other anomalous transport phenomena Hongo  _et al._ (2017); Yee and
Yin (2014); Hirono  _et al._ (2014); Yin and Liao (2016); Huang  _et al._
(2016); Jiang  _et al._ (2018); Shi  _et al._ (2018); Guo  _et al._ (2017);
Shi  _et al._ (2020); Siddique  _et al._ (2019), among other things. In
particular, the presence of the axial charge induced by a topological nature
of the system motivated people to extend RMHD to the chiral
magnetohydrodynamics (chiral MHD) Hattori  _et al._ (2019a), which gives a
potential theoretical tool to investigate the dynamical origin of the
primordial large-scale magnetic field Tashiro  _et al._ (2012); Boyarsky  _et
al._ (2015); Rogachevskii  _et al._ (2017); Schober  _et al._ (2018);
Brandenburg  _et al._ (2017); Boyarsky  _et al._ (2021) (see also Refs. Joyce
and Shaposhnikov (1997); Field and Carroll (2000); Giovannini (2004); Semikoz
and Sokoloff (2005); Laine (2005); Boyarsky  _et al._ (2012) for earlier
approaches to this problem based on RMHD).

Instead of exploring into the details of RMHD applications in many subfields
of physics, the goal of this article is to present a two-fold overview of the
theoretical aspects of the special relativistic MHD: One based on the recent
formulation of RMHD motivated by the generalized symmetry viewpoint, and the
other with the conventional approach in which the matter components and the EM
fields are separated discussed. We will present our view to relate these two
formulations with each other.

As a branch of relativistic hydrodynamics, the RMHD can also be considered as
an effective theory of low-energy long-wavelength modes of the system. Such
modes are usually the conserved modes like the energy and momentum modes. When
the EM fields are dynamical, it has been suggested recently that the Bianchi
identity for EM fields can be regarded as a conservation law (associated with
a one-form magnetic $U(1)$ symmetry) Grozdanov  _et al._ (2017); Hattori  _et
al._ (2019a); Glorioso and Son (2018); Armas and Jain (2019, 2020); Hongo and
Hattori (2021), allowing a formulation of RMHD based solely on symmetry
argument and derivative expansions of conserved quantities (in this case,
called hydrodynamic variables) 111See also Refs. Grozdanov and Poovuttikul
(2018); Armas  _et al._ (2018); Gralla and Iqbal (2019); Delacrétaz  _et al._
(2020); Iqbal and Poovuttikul (2020); Landry (2021) for several applications
of higher-form (and higher-group) symmetry to hydrodynamics.. We will give
detailed discussion about the construction of RMHD in this manner in Sec. II
using a phenomenological method based on entropy-current analysis and in Sec.
III using a nonequilibrium statistical operator method. A crucial observation
of such a construction is that the magnetic field, like temperature and fluid
velocity, persists at thermal equilibrium and can thus be assigned as a
leading-order variable in derivative expansion. This is known as a strong
magnetic field, and is the scenario that we are concerned with 222We will not
discuss the scenario with weak magnetic fields at sub-leading orders in
derivative expansion because the results are less novel, though it is also
very useful in practical applications.. As a result, the constitutive
relations exhibit anisotropies at both the ideal and dissipative levels, which
is a distinguished feature of RMHD (also non-relativistic MHD, in fact)
Braginskii (1965); Lifschitz and Pitaevskii (1981); Huang  _et al._ (2010,
2011); Hernandez and Kovtun (2017).

We will then examine how this new formulation is connected to the conventional
approach in Sec. IV, and proceed with the evaluations of transport
coefficients in kinetic theory and in perturbative quantum field theories
(QED, as an example) in Sec. V and Sec. VI, respectively. More novel features
can appear in RMHD when the system allows parity violation, leading to the
chiral MHD. This will be explored in Sec. VII. We note that there is an
significant difference between the treatments of EM fields in Secs. II, III,
VII and that in Secs. V and VI. In Secs. II, III, and VII, we treat the EM
fields as dynamical variables following a recent formulation of RMHD. On the
other hand, the EM fields in Secs. V and VI are considered as backgrounds, and
formulations are closer to the conventional approach with the EM fields
treated separately from the matter components. Considering this difference, we
compare the result of Secs. II-III to the conventional one to clarify these
two formulations in the intermediate section IV. Finally, we discuss the
future prospects of RMHD in Sec. VIII.

Here is a summary of our notations. We use the natural units,
$\hbar=c=k_{B}=1$. For the electromagnetism, we use the Heaviside-Lorentz
convention $\epsilon_{0}=\mu_{0}=1$ and $\alpha=e^{2}/(4\pi)$=1/137. Other
notations are:

  * •

Minkowski metric: $\eta_{\mu\nu}=\eta^{\mu\nu}={\rm diag}(1,-1,-1,-1)$. A
curved-spacetime metric is denoted by $g_{\mu\nu}$.

  * •

Levi-Civita tensor in Minkowski spacetime: $\varepsilon^{\mu\nu\rho\sigma}$
with $\varepsilon^{0123}=-\varepsilon_{0123}=1$.

  * •

Fluid velocity four vector: $u^{\mu}=\gamma\,(1,\bm{v})^{T}$ with
$\gamma=1/\sqrt{1-\bm{v}^{2}}$ and $u^{2}=1$.

  * •

Direction of the magnetic field: $b^{\mu}=B^{\mu}/B$ with
$B=\sqrt{-B_{\mu}B^{\mu}}$. Note that $b^{2}=-1$ and $b\cdot u=0$.

  * •

Projector transverse to $u^{\mu}$:
$\Delta^{\mu\nu}=\eta^{\mu\nu}-u^{\mu}u^{\nu}$.

  * •

Projector transverse to both $u^{\mu}$ and $b^{\mu}$:
$\Xi^{\mu\nu}=\Delta^{\mu\nu}+b^{\mu}b^{\nu}$.

  * •

Cross projector:
$b_{\star}^{\mu\nu}=-b_{\star}^{\nu\mu}=\varepsilon^{\mu\nu\rho\sigma}u_{\rho}b_{\sigma}$.
Note that $u_{\mu}b_{\star}^{\mu\nu}=0=b_{\mu}b_{\star}^{\mu\nu}$,
$b_{\star}^{\mu\lambda}b^{\nu}_{\star\lambda}=\Xi^{\mu\nu}$.

  * •

Co-moving derivative (or material derivative or proper-time derivative) of
$A$: $DA=\dot{A}=dA/d\tau=u^{\mu}\partial_{\mu}A$.

  * •

Spatial gradient of $A$: $\nabla_{\mu}A=\Delta_{\mu\nu}\partial^{\nu}A$.

  * •

Symmetrization, anti-symmetrization, and traceless symmetrization of a rank-
two tensor $A^{\mu\nu}$: $A^{(\mu\nu)}=(A^{\mu\nu}+A^{\nu\mu})/2$,
$A^{[\mu\nu]}=(A^{\mu\nu}-A^{\nu\mu})/2$, and
$A^{\langle\mu\nu\rangle}=(1/2)\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}(A^{\alpha\beta}+A^{\alpha\beta})-(1/3)\Delta^{\mu\nu}\Delta_{\alpha\beta}A^{\alpha\beta}$.

  * •

Decomposition of velocity gradient:
$\partial_{\mu}u_{\nu}=u_{\mu}Du_{\nu}+\omega_{\mu\nu}+w_{\mu\nu}$ where
$\omega_{\mu\nu}=\nabla_{[\mu}u_{\nu]}$ is the vorticity tensor and
$w_{\mu\nu}=\nabla_{(\mu}u_{\nu)}=\nabla_{\langle\mu}u_{\nu\rangle}+\frac{1}{3}\Delta_{\mu\nu}\theta$
with $\nabla_{\langle\mu}u_{\nu\rangle}$ the shear tensor and
$\theta=\Delta_{\alpha\beta}w^{\alpha\beta}$ the expansion rate of the fluid.

##  II Macroscopic approach: the entropy-current analysis

In this section, we review a phenomenological way to derive the RMHD equations
based on the second law of local thermodynamics. In Sec. II.1, we first give a
brief overview of how one can build relativistic hydrodynamics without
coupling to the EM fields before going into the discussion of RMHD. This can
help us understand some basic ingredients that are essential to the
construction of RMHD. In Sec. II.2, we generalize our discussion to RMHD based
on the symmetry associated with the Bianchi identity.

As mentioned in Sec. I, hydrodynamics is regarded as a low-energy effective
theory that describes the dynamics of gapless modes. In other words,
hydrodynamics describes the macroscopic behavior of conserved charges that do
not dissipate away. Therefore, it is crucial to identify symmetries of the
system and associated conservation laws that serve as the equations of motion
(EOMs) in hydrodynamics. This observation allows us to derive hydrodynamics
from the viewpoint of the symmetry and thermodynamics of irreversible
processes Landau and Lifshitz (1959); De Groot and Mazur (1962) without going
into details of the system. This macroscopic (phenomenological) approach is
often called the entropy-current analysis since the constitutive relations
(see definition below Eq. (3)) are derived by requiring the semi-positivity of
the entropy-current divergence, or the second law of local thermodynamics.

###  II.1 Primer to the entropy-current analysis

For the sake of simplicity, let us consider a relativistic system that only
enjoys spacetime translational symmetries, and thus, respects the energy-
momentum conservation law:

$\displaystyle\partial_{\mu}T^{\mu\nu}=0,$ (1)

where $T^{\mu\nu}$ denotes the energy-momentum tensor.Throughout this paper,
we always consider relativistic systems respecting the Lorentz symmetry, so
that one can write down the hydrodynamic equations in a Lorentz covariant
manner. One may think that the Lorentz symmetry itself leads to the angular
momentum conservation law, which should give independent hydrodynamic
equations as well. However, in the strict hydrodynamic limit (i.e., the long
wavelength and low frequency limit), the angular momentum conservation law
reduces to the constraint equation that forces the energy-momentum tensor to
be symmetric under the exchange of its two Lorentz indices. We thus use the
symmetric energy-momentum tensor $T^{\mu\nu}=T^{\nu\mu}$ in the following
discussion 333 In the transient time scale, called the spin hydrodynamic
regime in Refs. Hongo  _et al._ (2021, 2022); Cao  _et al._ (2022), the spin
density of microscopic constituents may show its intrinsic dynamics, and the
energy-momentum tensor has the anti-symmetric components. This is the main
topic of relativistic spin hydrodynamics which we will briefly discuss in Sec.
VIII.. We also note that the discrete charge-conjugation (C), time reversal
(T), and parity (P) symmetries also impose strong constraints on the transport
phenomena. As will be seen in Sec. VII, breaking one or two of them will allow
the appearance of new transport terms in the constitutive relations. In this
section, however, we assume that C, P, T symmetries are not violated.

Dynamical variables of hydrodynamics are the conserved energy-momentum
densities, or the corresponding thermodynamic conjugate variables given by the
fluid four-velocity $u^{\mu}$ — normalized as $\eta_{\mu\nu}u^{\mu}u^{\nu}=1$
with the Minkowski metric $\eta_{\mu\nu}=\mathrm{diag}(+1,-1,-1,-1)$ — and the
local inverse temperature $\beta=1/T$. In this paper, we employ the Landau-
Lifshitz frame Landau and Lifshitz (1959) to define the fluid velocity by

$T^{\mu}_{~{}\nu}u^{\nu}=\epsilon u^{\nu},$ (2)

with the energy density $\epsilon$. The local inverse temperature is related
to the energy density $\epsilon$ through the usual thermodynamic relation. We
will later give a concrete definition of the inverse temperature in Eq. (8)
when we start to consider the entropy density.

Then, as a long-wavelength low-frequency effective theory, one can formulate
hydrodynamics based on an expansion of $T^{\mu\nu}$ with respect to the small
derivatives of $u^{\mu}$ and $\beta$ as

$\displaystyle T^{\mu\nu}=T^{\mu\nu}_{(0)}+T^{\mu\nu}_{(1)}+O(\partial^{2}),$
(3)

where the subscripts $(0)$ and $(1)$ denote numbers of derivatives. This
equation, which expresses $T^{\mu\nu}$ in terms of the dynamical variable
$u^{\mu}$ and $\beta$ (or $\epsilon$) is called the constitutive relation. We
will construct the constitutive relation on an order-by-order basis.

In the leading-order derivative expansion, one can uniquely decompose the
zeroth-order term $T^{\mu\nu}_{(0)}$ as

$\displaystyle T^{\mu\nu}_{(0)}=\epsilon u^{\mu}u^{\nu}-p\Delta^{\mu\nu},$ (4)

with the projection tensor $\Delta_{\mu\nu}=\eta_{\mu\nu}-u_{\mu}u_{\nu}$
satisfying $u^{\mu}\Delta_{\mu\nu}=0$. This is because no other rank-two
symmetric tensors can be constructed from algebraic combinations of available
building blocks $u^{\mu},\eta^{\mu\nu}$, and Levi-Civita tensor
$\varepsilon^{\mu\nu\rho\sigma}$ (normalized to be
$\varepsilon^{0123}=-\varepsilon_{0123}=1$). The scalar coefficient functions
$\epsilon$ and $p$ are regarded as the energy density and pressure in the
fluid rest frame. This interpretation becomes manifest, $T_{(0)}^{\mu\nu}={\rm
diag}(\epsilon,p,p,p)$, when we take $u^{\mu}=(1,\bm{0})$. We note that we
have already imposed the matching condition for $\epsilon$,
$T^{\mu\nu}u_{\mu}u_{\nu}=T_{(0)}^{\mu\nu}u_{\mu}u_{\nu}$, in Eq. (4), which
enforces that $T^{\mu\nu}_{(1)}u_{\nu}=0$ following Eqs. (2). On the other
hand, no such matching condition is required for $p$, and there could be
$O(\partial^{1})$ correction to the physical pressure in the fluid.

At this stage, Eq. (4) is just a parameterization of $T^{\mu\nu}_{(0)}$ using
the unknown function $p$. In order to organize Eq. (1) in a solvable leading-
order hydrodynamic equation with $T^{\mu\nu}_{(0)}$, we have to relate the
pressure $p$ to $\epsilon$. If our system is close to the local thermal
equilibrium, we can expect that such a relation is provided by the equation of
state (EOS), $p=p(\epsilon)$, so that $p$ is identified as the thermodynamic
pressure. As we will see shortly, $p(\epsilon)$ has to satisfy a certain
thermodynamic relation in order to respect the second law of thermodynamics.
Then, substituting the leading-order constitutive relation (4) into the
conservation law (1), one obtains the relativistic Euler equations:

$\displaystyle(\epsilon+p)Du^{\mu}-\nabla^{\mu}p=0\,,$ (5a) $\displaystyle
D\epsilon+(\epsilon+p)\theta=0\,,$ (5b)

where we defined the material (or co-moving time) derivative
$D=u^{\mu}\partial_{\mu}$, the spatial gradient
$\nabla_{\mu}=\Delta_{\mu\nu}\partial^{\nu}$, and the expansion rate of the
fluid $\theta=\partial_{\mu}u^{\mu}$.

In the same manner, we can continue to perform the tensor decomposition of the
first-order term $T^{\mu\nu}_{(1)}$. Noting that $T^{\mu\nu}_{(1)}$ satisfies
$T^{\mu\nu}_{(1)}u_{\nu}=0$ enforced by the Landau-Lifshitz frame condition
and the matching condition for $\epsilon$, we use $\Delta^{\mu\nu}$ to
decompose the rank-two symmetric tensor $T^{\mu\nu}_{(1)}$ into the trace part
$\Pi$ and symmetric traceless part $\pi^{\mu\nu}$ (or
$\pi^{\mu\nu}\Delta_{\mu\nu}=0$) as

$\displaystyle T^{\mu\nu}_{(1)}=-\Delta^{\mu\nu}\Pi+\pi^{\mu\nu}.$ (6)

This equation gives the decomposition of the rank-two symmetric tensor into
irreducible representations of the $\mathrm{SO}(3)$ rotation, and thus,
$\Delta^{\mu\nu}\Pi$ and $\pi^{\mu\nu}$ do not mix under the $\mathrm{SO}(3)$
rotation. Here, we can regard $\Pi$ as the derivative correction to the
pressure because it appears at the same place as the thermodynamic pressure
$p$ does. As we will see, $\Pi$ describes a viscous correction of the pressure
proportional to the expansion rate $\theta$, while $\pi^{\mu\nu}$ gives the
shear viscous contribution describing a friction-like process induced by the
velocity gradient. As being $O(\partial)$-order quantities, $\Pi$, and
$\pi^{\mu\nu}$ should be structured linearly in the gradients of $u^{\mu}$ and
$\beta$ (or equivalently, $\epsilon$).

The entropy current analysis gives an elegant way to work out such linear
structures based only on the thermodynamic laws. First, we write down the
entropy current in a similar manner as that for $T^{\mu\nu}$,

$\displaystyle s^{\mu}=su^{\mu}+s_{(1)}^{\mu}+O(\partial^{2}),$ (7)

where $s=s(\epsilon)$ is the thermodynamic entropy density in the fluid rest
frame and $s_{(1)}^{\mu}$ is a derivative correction transverse to $u^{\mu}$.
Note that the thermodynamic entropy density is a function of the energy
density $\epsilon$, whose derivative with respect to $\epsilon$ gives a
definition of the local inverse temperature $\beta(x)=1/T(x)$ as

$\beta(x):=\left.\frac{\partial
s(\epsilon)}{\partial\epsilon}\right|_{\epsilon=\epsilon(x)}~{}\Leftrightarrow~{}T(x)\mathrm{d}s(x)=\mathrm{d}\epsilon(x).$
(8)

Then, using the chain rule for the material derivative of the entropy density,
and the contracted equation of motion $u_{\nu}\partial_{\mu}T^{\mu\nu}=0$ with
Eqs. (4) and (6), we find the divergence of the entropy current to be

$\displaystyle\partial_{\mu}s^{\mu}=\beta(Ts-\epsilon-p)\theta-\beta\Pi\theta+\beta\partial_{\langle\mu}u_{\nu\rangle}\pi^{\mu\nu}+\partial_{\mu}s_{(1)}^{\mu}+O(\partial^{3}),$
(9)

where we introduced the symmetric traceless projection of the velocity
gradient $\partial_{\langle\mu}u_{\nu\rangle}$ (called shear tensor) as

$\displaystyle\partial^{\langle\mu}u^{\nu\rangle}=\nabla^{\langle\mu}u^{\nu\rangle}=\frac{1}{2}\left(\nabla^{\mu}u^{\nu}+\nabla^{\nu}u^{\mu}\right)-\frac{1}{3}\Delta^{\mu\nu}\theta.$
(10)

Since the divergence of the entropy current gives a local entropy production
rate, the second law of local thermodynamics requires
$\partial_{\mu}s^{\mu}\geq 0$ for any configuration of the hydrodynamic
variables $\beta$ and $u^{\mu}$. This is achieved by identifying the first-
order entropy current as $s_{(1)}^{\mu}=0$ and the following conditions:

$\displaystyle Ts$ $\displaystyle=$ $\displaystyle\epsilon+p,$ (11a)
$\displaystyle\Pi$ $\displaystyle=$ $\displaystyle-\zeta\,\theta,$ (11b)
$\displaystyle\pi^{\mu\nu}$ $\displaystyle=$ $\displaystyle
2\eta\,\partial^{\langle\mu}u^{\nu\rangle}.$ (11c)

Here, Eq. (11a) gives the thermodynamic relation, which restricts
$p(\epsilon)$ appearing in the leading-order constitutive relation to be the
thermodynamic pressure. Equations (11b)-(11c) complete the first-order
constitutive relations with the phenomenological parameters $\zeta,\eta$
identified as the transport coefficients called bulk viscosity and shear
viscosity, respectively. The second law of local thermodynamics
$\partial_{\mu}s^{\mu}\geq 0$ requires the semi-positivity of those transport
coefficients.

With the identification of $\Pi$ and $\pi^{\mu\nu}$ given in Eqs. (11b)-(11c),
we obtain the hydrodynamic equations up to $O(\partial^{2})$ by substituting
the derived constitutive relations into the conservation law (1). By
projecting the obtained equation into the spatial and temporal directions, we
have

$\displaystyle(\epsilon+p-\zeta\theta)Du^{\mu}-\nabla^{\mu}(p-\zeta\theta)+2\Delta^{\mu}_{\nu}\partial_{\rho}(\eta\partial^{\langle\nu}u^{\rho\rangle})=0\,,$
(12a) $\displaystyle
D\epsilon+(\epsilon+p-\zeta\theta)\theta-2\eta\partial_{\langle\mu}u_{\nu\rangle}\partial^{\langle\mu}u^{\nu\rangle}=0\,.$
(12b)

These are the relativistic Navier-Stokes equations.

One can continue the above procedure to higher orders in derivative expansion.
The resultant constitutive relations and the hydrodynamic equations become
complicated, and we do not report these higher order results here. The readers
who are interested in such results can find excellent discussions in, e.g.,
Refs. Israel and Stewart (1979a); Baier  _et al._ (2008); Romatschke and
Romatschke (2019).

###  II.2 Relativistic MHD from the magnetic flux conservation

We now proceed to the discussion of RMHD. Here, what is different from the
usual relativistic hydrodynamics is the existence of the $\mathrm{U}(1)$ gauge
field, or dynamical electric and magnetic fields. We thus first need to
examine whether or not those fields deserve to be qualified as hydrodynamic
variables that persist in the low-energy long-wavelength limit without being
damped out in a fluid.

Firstly, one finds that the electric field is not a hydrodynamic variable. To
see this, notice that electric fields are subject to the Debye screening,
which is a static screening effect in the long spacetime limit Kapusta and
Gale (2011); Bellac (2011). Namely, an electric-charge density is
redistributed so that the Coulomb field offsets the electric field in an
equilibrium state. Therefore, an electric field is a gapped excitation due to
the Debye screening mass, and does not deserve a hydrodynamic variable. One
also finds the same conclusion from the Maxwell equation
$\partial_{\mu}F^{\mu\nu}=J^{\nu}$, where
$F_{\mu\nu}:=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\nu}$ is the field
strength tensor and $J^{\mu}$ is the electric current. The Maxwell equation
indicates that the electric flux is not conserved due to the presence of the
electric current. With a given constitutive relation for $J^{\mu}$ describing
the Ohmic current, as shown in Appendix A, the electric field behaves like a
decay mode in neutral plasma.

In contrast, there is no static screening effect on a magnetic field, meaning
that magnetic field enjoys its intrinsic dynamics in the macroscopic scale
even in a medium. This qualitative difference from the electric fields stems
from the absence of a magnetic monopole. Namely, there is no “magnetic charge”
distribution that can screen magnetic fields. This is a consequence of the
Bianchi identity $\partial_{\mu}\widetilde{F}^{\mu\nu}=0$, where we defined
the Hodge dual of $F_{\mu\nu}$ by
$\widetilde{F}^{\mu\nu}=\frac{1}{2}\varepsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$.

The novel formulation of RMHD reviewed in Secs. II and III of this paper is
based on this simple observation that the magnetic flux is a conserved
quantity Schubring (2015); Grozdanov  _et al._ (2017); Hattori  _et al._
(2019a); Glorioso and Son (2018); Armas and Jain (2019, 2020); Hongo and
Hattori (2021) due to the Bianchi identity. This observation allows us to
straightforwardly derive RMHD according to the philosophy of hydrodynamics
that motivates us to keep only the conserved quantities resulting from the
associate symmetries 444The corresponding symmetry associated with the
conserved magnetic flux is called a one-form magnetic $\mathrm{U}(1)$ symmetry
Gaiotto  _et al._ (2015). This gives a generalization of the usual global
symmetry because the conserved magnetic-flux acts on a one-dimensional
extended object, or a ’t Hooft loop, rather than a conventional local operator
having a point-like charge. . In this section, we formulate RMHD by plugging
the magnetic-flux conservation law into the entropy-current analysis. In Sec.
III, we identify the symmetry for the magnetic-flux conservation inherent in
QED and provide the derivation of RMHD on the basis of the nonequilibrium
statistical method.

Before diving into the entropy-current analysis, it is useful to pay our
attention to a basic difference with the conventional formulation of MHD
Landau  _et al._ (1984); Kluitenberg and De Groot (1954b, c); Lichnerowicz
(1967); Israel (1978); Gedalin (1991). In the conventional formulation of MHD,
the starting point is to couple the Navier-Stokes (or Euler) equation with the
Maxwell equation and Bianchi identity. This means that we also keep the
electric field as a dynamical variable, which shows a non-hydrodynamic
relaxational behavior according to our identification. Thus, the conventional
MHD describes the transient dynamics including the electric field. On the
other hand, the new formulation only keeps conserved quantities as the
dynamical variable, and gives a low-energy effective theory in a strict
hydrodynamic limit.

####  II.2.1 Entropy production rate with magnetic flux

Let us perform the entropy-current analysis to derive the RMHD equation. This
is carried out as a direct extension of the simple analysis in Sec. II.1 by
additionally considering the magnetic-flux conservation law.

First of all, the energy-momentum and magnetic-flux conservation laws are
given as

$\displaystyle\partial_{\mu}T^{\mu\nu}=0\,,\quad\partial_{\mu}\tilde{F}^{\mu\nu}=0\,.$
(13)

Here, we note that $T^{\mu\nu}$ is the total energy-momentum tensor of the
system that includes not only the matter part but also the contribution from
the EM fields. The second equation, or the Bianchi identity, indicates the
conservation of magnetic flux, which corresponds to the Gauss’s law for the
magnetic flux in the absence of a magnetic monopole. Following the philosophy
of hydrodynamics discussed above, we do not include another Maxwell equation
into the set of hydrodynamic equations because it describes the time evolution
of the gapped electric fields.

We then introduce the dynamical variables of RMHD. Associated with the
conserved energy-momentum density, we again employ the Landau-Lifshitz frame
to define the energy density $\epsilon$ and the normalized fluid four-velocity
$u^{\mu}$ by Eq. (2). Besides, we define a magnetic flux density by
$B^{\mu}=\tilde{F}^{\mu\nu}u_{\nu}$ in a covariant manner. This four vector
reduces to $B^{\mu}=(0,{\bm{B}})$ in the local rest frame with
$u^{\mu}=(1,0,0,0)$. Thus, the number of components of it is not four but
three, which is manifest by noting that $B^{\mu}$ is transverse to the fluid
velocity, i.e., $u_{\mu}B^{\mu}=0$.

Let us next introduce the conjugate variables based on the first law of local
thermodynamics. Relying on the thermodynamic entropy density
$s(\epsilon,B^{\mu})$ in the rest frame, which is a function of the energy
density $\epsilon$ and magnetic flux density $B^{\mu}$, we define the local
inverse temperature $\beta$ and (in-medium) magnetic field 555We note that in
Sec. II, III, and VII we call $H_{\mu}$ the magnetic field and $B^{\mu}$ the
magnetic flux density, but in other sections, we call $B^{\mu}$ the magnetic
field without confusion. $H_{\mu}$ as

$\displaystyle\beta:=\frac{\partial
s(\epsilon,B^{\mu})}{\partial\epsilon},\quad\beta H_{\mu}:=-\frac{\partial
s(\epsilon,B^{\mu})}{\partial
B^{\mu}}~{}\Leftrightarrow~{}T\mathrm{d}s=\mathrm{d}\epsilon-
H_{\mu}\mathrm{d}B^{\mu}\,.$ (14)

With the help of these variables, we will derive the constitutive relations
for $T^{\mu\nu}$ and $\widetilde{F}^{\mu\nu}$.

To perform the entropy-current analysis with Eqs. (13) and (14), it is useful
to decompose $T^{\mu\nu}$ and $\widetilde{F}^{\mu\nu}$ into possible tensor
structures. In the present setup, we have the magnetic flux vector $B^{\mu}$
in addition to the flow vector $u^{\mu}$ as available zeroth-order vectors. It
is then convenient to introduce a normalized vector $b^{\mu}=B^{\mu}/B$ with
$B=\sqrt{-B_{\mu}B^{\mu}}$ such that $\eta_{\mu\nu}b^{\mu}b^{\nu}=-1$ and
$u_{\mu}b^{\mu}=0$. With this vector, we also introduce a projector
$\Xi^{\mu\nu}=\Delta^{\mu\nu}+b^{\mu}b^{\nu}$ transverse to both $u^{\mu}$ and
$b^{\mu}$ ($\Xi^{\mu\nu}u_{\mu}=0=\Xi^{\mu\nu}b_{\mu}$). By the use of those
tensors, we parameterize the constitutive relations as

$\displaystyle T^{\mu\nu}$ $\displaystyle=$ $\displaystyle\epsilon
u^{\mu}u^{\nu}-p_{\perp}\Xi^{\mu\nu}+p_{\parallel}b^{\mu}b^{\nu}+T^{\mu\nu}_{(1)}\,,$
(15a) $\displaystyle\tilde{F}^{\mu\nu}$ $\displaystyle=$ $\displaystyle
B^{\mu}u^{\nu}-B^{\nu}u^{\mu}+\tilde{F}^{\mu\nu}_{(1)}\,.$ (15b)

Here, we assume that $p_{\perp},\,p_{\parallel}$ are zeroth-order in
derivatives. The first-order corrections are collectively denoted as
$T^{\mu\nu}_{(1)}$ and $\tilde{F}_{(1)}^{\mu\nu}$. Considering a charge-
neutral and parity-even plasma, we assume the charge-conjugation symmetry and
require the energy-momentum tensor $T^{\mu\nu}$ not to have a charge-
conjugation odd term proportional to $u^{(\mu}b^{\nu)}$ at the zeroth-order in
derivatives. Likewise, $\tilde{F}^{\mu\nu}$ could not have a term
$\varepsilon^{\mu\nu\alpha\beta}u_{\alpha}b_{\beta}$ at the zeroth-order
because the resulting scalar coefficient for such a term is parity odd.
Plugging those expressions into the equations of motion (13) and contracting
those with $\beta u_{\nu}$ and $H_{\nu}$, we obtain

$\displaystyle 0$
$\displaystyle=u_{\nu}\partial_{\mu}T^{\mu\nu}=D\epsilon+(\epsilon+p_{\perp})\theta+(p_{\perp}-p_{\parallel})b_{\nu}b^{\mu}\partial_{\mu}u^{\nu}+u_{\nu}\partial_{\mu}T_{(1)}^{\mu\nu}\,,$
(16a) $\displaystyle 0$
$\displaystyle=H_{\nu}\partial_{\mu}\tilde{F}^{\mu\nu}=-H_{\mu}DB^{\mu}-B^{\mu}H_{\mu}\theta+B^{\mu}H_{\nu}\partial_{\mu}u^{\nu}+H_{\nu}\partial_{\mu}\tilde{F}_{(1)}^{\mu\nu}\,.$
(16b)

Here, we note that our definitions of $\epsilon$, $u^{\mu}$, and $B^{\mu}$
lead to the matching conditions
$u_{\nu}T_{(1)}^{\mu\nu}=0=u_{\nu}\tilde{F}_{(1)}^{\mu\nu}$ (see Appendix B).

Now, let us compute the divergence of the entropy current
$s^{\mu}=su^{\mu}+s^{\mu}_{(1)}$. As in the previous section, we use the chain
rule and the equations of motion (16) to eliminate $D\epsilon$ and
$-H_{\mu}DB^{\mu}$. Then, we obtain

$\displaystyle\partial_{\mu}s^{\mu}$ $\displaystyle=$ $\displaystyle
s\theta+Ds+\partial_{\mu}s_{(1)}^{\mu}$ (17) $\displaystyle=$ $\displaystyle
s\theta+\beta D\epsilon-\beta H_{\mu}DB^{\mu}+\partial_{\mu}s_{(1)}^{\mu}$
$\displaystyle=$ $\displaystyle\beta(Ts-\epsilon-
p_{\perp}+B^{\mu}H_{\mu})\theta-\beta
b^{\mu}[(p_{\perp}-p_{\parallel})b^{\nu}+H^{\nu}]\partial_{\mu}u_{\nu}$
$\displaystyle+T_{(1)}^{\mu\nu}\partial_{\mu}(\beta
u_{\nu})+\tilde{F}_{(1)}^{\mu\nu}\partial_{\mu}(\beta
H_{\nu})+\partial_{\mu}(s_{(1)}^{\mu}-\beta u_{\nu}T_{(1)}^{\mu\nu}-\beta
H_{\nu}\tilde{F}_{(1)}^{\mu\nu})\,.$

Applying the second law of local thermodynamics $\partial_{\mu}s^{\mu}\geq 0$
to Eq. (17), we will derive the RMHD equation in the leading and next-to-
leading orders in derivatives.

####  II.2.2 Zeroth-order in derivatives: Nondissipative RMHD

At the leading order in derivatives, we find a set of constraints by requiring
the absence of entropy production. For this requirement to be satisfied in any
hydrodynamic configuration, all the leading-order terms should vanish
independently because they are proportional to the irreducible tensor
decomposition of $\partial_{\mu}u_{\nu}$ and a non-derivative form of
$u_{\nu}$ in Eq. (17). We thus find the following four constraints

$\displaystyle Ts+B^{\mu}H_{\mu}=\epsilon+p_{\perp}\,,$ (18a)
$\displaystyle(p_{\perp}-p_{\parallel})b^{\nu}+BH^{\nu}=0\,.$ (18b)

The first constraint serves as an extension of the thermodynamic relation with
the magnetic flux. The second constraint indicates that there is a pressure
anisotropy $p_{\perp}-p_{\parallel}=B^{\mu}H_{\mu}(>0)$ induced by the finite
magnetic flux (cf. Fig. 1). Besides, noting that the magnetic field $H^{\mu}$
is parallel to $B^{\mu}$ and recalling that the (in-medium) magnetic field is
parameterized by $H^{\mu}=-\mu_{m}^{-1}B^{\mu}$ with the magnetic permeability
$\mu_{m}$, one identifies that the pressure difference gives the magnetic
permeability, $\mu_{m}^{-1}=(p_{\perp}-p_{\parallel})/|B|^{2}$.

Figure 1: Anisotropic pressure induced by the magnetic field and the
magnetization.

We conclude that the leading-order analysis leads to the following
constitutive relations

$\displaystyle T^{\mu\nu}_{(0)}$ $\displaystyle=$ $\displaystyle\epsilon
u^{\mu}u^{\nu}-p_{\perp}\Xi^{\mu\nu}+(p_{\perp}-B^{\alpha}H_{\alpha})b^{\mu}b^{\nu}\,,$
(19a) $\displaystyle\tilde{F}^{\mu\nu}_{(0)}$ $\displaystyle=$ $\displaystyle
B^{\mu}u^{\nu}-B^{\nu}u^{\mu}\,.$ (19b)

Note that the derived constitutive relation for $\tilde{F}^{\mu\nu}_{(0)}$
indicates that the electric field is absent at the leading order because it is
given by
$E^{\mu}_{(0)}:=-\frac{1}{2}\epsilon^{\mu\nu\alpha\beta}u_{\nu}\tilde{F}_{{(0)}\alpha\beta}=0$.
This result is consistent with our starting point treating the electric field
as a non-hydrodynamic gapped variable. We will find that an electric field
(and corresponding electric current) appears as a first-order derivative
correction. This will clearly show that they are slaved to the true
hydrodynamic variables in the strict hydrodynamic limit as shown below in the
form of the first-order constitutive relation.

####  II.2.3 First-order derivative corrections: Dissipative RMHD

We proceed to the derivation of the first-order corrections
$T_{(1)}^{\mu\nu},\,\tilde{F}_{(1)}^{\mu\nu}$. Being equipped with the
identified constraints (18), we can simplify the entropy production rate (17)
as

$\displaystyle\partial_{\mu}s^{\mu}$ $\displaystyle=$ $\displaystyle\beta
T_{(1)}^{\mu\nu}\partial_{\mu}u_{\nu}+\tilde{F}_{(1)}^{\mu\nu}\partial_{\mu}(\beta
H_{\nu})+\partial_{\mu}\big{(}s_{(1)}^{\mu}-\beta
u_{\nu}T_{(1)}^{\mu\nu}-\beta H_{\nu}\tilde{F}_{(1)}^{\mu\nu}\big{)}\,,$ (20)

where we used $T_{(1)}^{\mu\nu}u_{\nu}=0$. Noting that $u_{\nu}$ and $\beta
H_{\nu}$ are independent variables, one can ensure the semi-positive entropy
production, or the local second law $\partial_{\mu}s^{\mu}\geq 0$, in any
hydrodynamic configuration by requiring that

$\displaystyle{\mathcal{R}}_{u}$ $\displaystyle:=\beta
T_{(1)}^{\mu\nu}\partial_{\mu}u_{\nu}\geq 0\,,$ (21a)
$\displaystyle{\mathcal{R}}_{E}$
$\displaystyle:=\tilde{F}_{(1)}^{\mu\nu}\partial_{\mu}(\beta H_{\nu})\geq
0\,,$ (21b) $\displaystyle s_{(1)}^{\mu}$ $\displaystyle=\beta
u_{\nu}T_{(1)}^{\mu\nu}+\beta H_{\nu}\tilde{F}_{(1)}^{\mu\nu}\,.$ (21c)

As we will see, the first two conditions are satisfied by requiring each term
to be a semi-positive bilinear form, which determines the first-order
constitutive relations. With the resulting constitutive relations, we can also
find the derivative corrections to the entropy current from Eq. (21c). Below,
we will separately analyze the constitutive relations for $T_{(1)}^{\mu\nu}$
and $\tilde{F}_{(1)}^{\mu\nu}$.

##### Electric field and resistivities:

Let us first derive the constitutive relation for $\tilde{F}_{(1)}^{\mu\nu}$,
focusing on Eq. (21b). To make ${\mathcal{R}}_{E}$ a semi-positive bilinear,
we can express $\tilde{F}_{(1)}^{\mu\nu}$ as

$\tilde{F}_{(1)}^{\mu\nu}=T\rho^{\mu\nu\rho\sigma}\partial_{[\rho}(\beta
H_{\sigma]}).$ (22)

Here, we introduced a rank-four tensor $\rho^{\mu\nu\rho\sigma}$, which will
be identified as a resistivity tensor. As we will see, the semi-positivity of
${\mathcal{R}}_{E}$ and the anti-symmetric property in the Lorentz indices
restrict possible structures of the resistivity tensor
$\rho^{\mu\nu\rho\sigma}$.

To identify the tensor structure of $\rho^{\mu\nu\rho\sigma}$, we first recall
that $\tilde{F}_{(1)}^{\mu\nu}$ is transverse to $u^{\nu}$,
$\tilde{F}_{(1)}^{\mu\nu}u_{\nu}=0$. Thus, we can use $b^{\mu}$ and
$\Xi^{\mu\nu}$ to perform the tensor decomposition of
$\rho^{\mu\nu\rho\sigma}$. Moreover, we note that it is anti-symmetric with
respect to the exchange of its Lorentz indices $\mu\leftrightarrow\nu$ and
$\rho\leftrightarrow\sigma$. Taking into account these properties, we can
write down the most general form for $\rho^{\mu\nu\rho\sigma}$ in a neutral
plasma as

$\rho^{\mu\nu\rho\sigma}=2\rho_{\perp}\left(b^{\mu}\Xi^{\nu[\rho}b^{\sigma]}-b^{\nu}\Xi^{\mu[\rho}b^{\sigma]}\right)+2\rho_{\parallel}\Xi^{\mu[\rho}\Xi^{\sigma]\nu}.$
(23)

The two coefficients $\rho_{\perp,\parallel}$ will be identified with two
components of the electric resistivity (see Sec. IV.2). Those tensor
structures project out the gradient of the magnetic field in parallel and
perpendicular to $b^{\mu}$ as opposed to what the subscripts of
$\rho_{\perp,\parallel}$ denote. This is because the electric field is defined
with $\widetilde{F}^{\mu\nu}_{(1)}$ and the antisymmetric tensor that swaps
the directions (see Eq. (25) below). Note that $\rho^{\mu\nu\rho\sigma}$
cannot have a tensor structure including
$b_{\star}^{\mu\nu}:=\varepsilon^{\mu\nu\alpha\beta}u_{\alpha}b_{\beta}$
because any such term would violate the charge-conjugation symmetry. To
respect the second law of local thermodynamics, we require semi-positivity of
the two resistivities, $\rho_{\parallel}\geq 0$ and $\rho_{\perp}\geq 0$. In
fact, substituting Eqs. (22) and (23) into Eq. (21b), we find

${\mathcal{R}}_{E}=8T\rho_{\perp}\left(b^{[\nu}\Xi^{\mu]\rho}b^{\sigma}\partial_{[\rho}\big{(}\beta
H_{\sigma]}\big{)}\right)^{2}+2T\rho_{\parallel}\left(\Xi^{\mu[\rho}\Xi^{\sigma]\nu}\partial_{[\rho}\big{(}\beta
H_{\sigma]}\big{)}\right)^{2},$ (24)

which is positive semi-definite. Equations (22) and (23) complete the
constitutive relation for $\tilde{F}_{(1)}^{\mu\nu}$. The first-order
corrections to $\tilde{F}_{(1)}^{\mu\nu}$ gives rise to an induced electric
field

$\displaystyle
E_{(1)}^{\mu}:=-\frac{1}{2}\varepsilon^{\mu\nu\alpha\beta}u_{\nu}\tilde{F}_{{(1)}\alpha\beta}=(\rho_{\parallel}b^{\mu}b^{\nu}-\rho_{\perp}\Xi^{\mu\nu})T\varepsilon_{\nu\lambda\alpha\beta}u^{\lambda}\partial^{\alpha}(\beta
H^{\beta})\,,$ (25)

where we used an identity
$2b_{\star}^{\mu[\alpha}b^{\beta]}=\Xi^{\mu}_{\rho}\epsilon^{\rho\nu\alpha\beta}u_{\nu}$
obtained from the Schouten identity.

##### Stress tensor and viscosities:

One can derive the first-order corrections to the energy-momentum tensor
$T^{\mu\nu}_{(1)}$ in the same manner. To ensure the semi-positivity of (21a),
the left-hand side should also be a positive semi-definite bilinear.
Introducing the rank-four viscous tensor $\eta^{\mu\nu\rho\sigma}$, we now
express $T_{(1)}^{\mu\nu}$ as

$T_{(1)}^{\mu\nu}=\eta^{\mu\nu\rho\sigma}\partial_{\rho}u_{\sigma},$ (26)

and identify the tensor structure of $\eta^{\mu\nu\rho\sigma}$.

We first recall that we employ the Landau-Lifshitz frame (see Appendix B for a
detailed discussion on the frame choice), in which $T_{(1)}^{\mu\nu}$ is
transverse to $u^{\nu}$ as $T_{(1)}^{\mu\nu}u^{\nu}=0$. Thus, we can only use
$b^{\mu}$ and $\Xi^{\mu\nu}$ as a possible vector and tensor to decompose the
viscous tensor $\eta^{\mu\nu\rho\sigma}$. Note that $b_{\star}^{\mu\nu}$
cannot be used to decompose $\eta^{\mu\nu\rho\sigma}$ because of charge-
conjugation symmetry. Moreover, recalling that $T^{\mu\nu}$ is the symmetric
energy-momentum tensor, one finds that the viscous tensor should be symmetric
with respect to the exchanges between its Lorentz indices
$\mu\leftrightarrow\nu$ and $\rho\leftrightarrow\sigma$. These properties
allow us to perform the tensor decomposition as

$\displaystyle\eta^{\mu\nu\rho\sigma}=$
$\displaystyle\zeta_{\parallel}b^{\mu}b^{\nu}b^{\rho}b^{\sigma}+\zeta_{\perp}\Xi^{\mu\nu}\Xi^{\rho\sigma}-\zeta_{\times}(b^{\mu}b^{\nu}\Xi^{\rho\sigma}+\Xi^{\mu\nu}b^{\rho}b^{\sigma})$
$\displaystyle-4\eta_{\parallel}b^{(\mu}\Xi^{\nu)(\rho}b^{\sigma)}+\eta_{\perp}\left(\Xi^{\mu\alpha}\Xi^{\nu\beta}+\Xi^{\nu\alpha}\Xi^{\mu\beta}-\Xi^{\mu\nu}\Xi^{\alpha\beta}\right),$
(27)

where we introduced five viscosities — three bulk viscosities
$\zeta_{\parallel},\zeta_{\perp},\zeta_{\times}$ and two shear viscosities
$\eta_{\parallel},\eta_{\perp}$. As we will specify, these viscosities must
satisfy a semi-positivity constraint to ensure the second law of local
thermodynamics.

To get a physical intuition of the dissipative processes and find the semi-
positivity constraint attached to each viscosity, it is useful to decompose
the velocity gradient as

$\displaystyle\partial_{(\mu}u_{\nu)}=-\theta_{\parallel}b_{\mu}b_{\nu}+\frac{1}{2}\theta_{\perp}\Xi_{\mu\nu}-2b_{\alpha}b_{(\mu}\Xi_{\nu)\beta}\partial^{(\alpha}u^{\beta)}+\partial_{\\{\mu}u_{\nu\\}},$
(28)

where we defined $\theta_{\parallel}:=(-b^{\mu}b^{\nu})\partial_{\mu}u_{\nu}$,
$\theta_{\perp}:=\Xi^{\mu\nu}\partial_{\mu}u_{\nu}$, and

$\partial_{\\{\mu}u_{\nu\\}}:=\frac{1}{2}\left(\Xi_{\alpha\mu}\Xi_{\nu\beta}+\Xi_{\alpha\nu}\Xi_{\mu\beta}-\Xi_{\mu\nu}\Xi_{\alpha\beta}\right)\partial^{\alpha}u^{\beta}.$
(29)

Using this decomposition, we find the first-order derivative corrections to
the constitutive relation (26) to be

$\displaystyle T_{(1)}^{\mu\nu}=$
$\displaystyle-(\zeta_{\parallel}\theta_{\parallel}+\zeta_{\times}\theta_{\perp})b^{\mu}b^{\nu}+(\zeta_{\perp}\theta_{\perp}+\zeta_{\times}\theta_{\parallel})\Xi^{\mu\nu}-4\eta_{\parallel}b^{(\mu}\Xi^{\nu)(\rho}b^{\sigma)}\partial_{\rho}u_{\sigma}+2\eta_{\perp}\partial^{\\{\mu}u^{\nu\\}}.$
(30)

Equation (30) allows us to get a physical intuition on the dissipative process
attached to each viscosity. Without the magnetic flux, an
expansion/compression rate is given as $\theta=\nabla_{\mu}u^{\mu}$. In the
presence of the magnetic flux,the expansion/compression in the parallel and
perpendicular directions should be distinguished from each other, and thus, we
have $\theta_{\parallel}$ and $\theta_{\perp}$. Two viscosities
$\zeta_{\parallel}$ and $\zeta_{\perp}$ describes a resistance to such two
expansion/compression, respectively. Similarly, the flow gradient for the
shear deformation is projected into the parallel and perpendicular directions,
which leads to the two friction-like processes described by $\eta_{\parallel}$
and $\eta_{\perp}$ (cf. Figs. 2 and 3). Based on these observations, we
identify two coefficients $\zeta_{\parallel}$ and $\zeta_{\perp}$ with bulk
viscosities and another two coefficients $\eta_{\parallel}$ and $\eta_{\perp}$
with shear viscosities. Besides, there is an off-diagonal cross response
proportional to $\zeta_{\times}$. We also identify $\zeta_{\times}$ as one of
the bulk viscosities since the associated term describes the cross response of
the anisotropic pressures to the two expansion/compression rates (cf. Fig. 4).
These processes are reciprocal to one another, and the associated transport
coefficients should be the same according to the Onsager’s reciprocal relation
Onsager (1931); Hooyman  _et al._ (1954); De Groot and Mazur (1962); Grozdanov
_et al._ (2017). Putting $\zeta_{\parallel}=\zeta_{\perp}=\zeta_{\times}$ and
$\eta_{\parallel}=\eta_{\perp}$, one can confirm that the anisotropic viscous
corrections in Eq. (30) reduce to the isotropic form (6) by the use of
$\theta=\theta_{\parallel}+\theta_{\perp}$ and the identity (28).

Figure 2: Longitudinal and transverse bulk viscosities that are the off-
equilibrium responses of the pressure (green arrows) in the direction of the
expansion/compression (blue arrows) of the system. Figure 3: Shear
deformations in and out of planes with respect to the magnetic field
direction. Figure 4: Cross bulk viscosities that are the response of the
pressure (green arrows) in the orthogonal direction to the
expansion/compression (blue arrows) of the system. Those two cases are
reciprocal processes to each other. Compare those two with Fig. 2.

We next investigate the positivity constraint required for the viscosities.
Substituting the decomposition (28) into Eq. (21a), we obtain

$\displaystyle{\mathcal{R}}_{u}$
$\displaystyle=\beta\eta^{\mu\nu\rho\sigma}\partial_{\mu}u_{\nu}\partial_{\rho}u_{\sigma}$
$\displaystyle=\beta\begin{pmatrix}\theta_{\parallel}&\theta_{\perp}\end{pmatrix}\begin{pmatrix}\zeta_{\parallel}&\zeta_{\times}\\\
\zeta_{\times}&\zeta_{\perp}\end{pmatrix}\begin{pmatrix}\theta_{\parallel}\\\
\theta_{\perp}\end{pmatrix}+2\beta\eta_{\parallel}\big{(}2b_{(\mu}\Xi_{\nu)\rho}b_{\sigma}\partial^{\rho}u^{\sigma}\big{)}^{2}+2\beta\eta_{\perp}\big{(}\partial_{\\{\mu}u_{\nu\\}}\big{)}^{2}\,.$
(31)

Firstly, the last two terms are positive semi-definite by requiring
$\eta_{\parallel}\geq 0$ and $\eta_{\perp}\geq 0$. On the other hand, for the
first term to respect the local second law, the matrix composed of the bulk
viscosities should be positive semi-definite. This means that the eigenvalue
of that matrix needs to be non-negative:
$\zeta_{\parallel}+\zeta_{\perp}\geq\pm\sqrt{(\zeta_{\parallel}-\zeta_{\perp})^{2}+4\zeta_{\times}^{2}}$.
Notice also that the semi-positivity should be separately ensured in a
parallel expansion/compression ($\theta_{\parallel}\not=0,\ \theta_{\perp}=0$)
and in a perpendicular one ($\theta_{\perp}\not=0,\ \theta_{\parallel}=0$) as
necessary conditions, requiring that $\zeta_{\parallel}\geq 0$ and
$\zeta_{\perp}\geq 0$. Then, we find an inequality,
$\zeta_{\parallel}\zeta_{\perp}-\zeta_{\times}^{2}\geq 0$, from the above
eigenvalues. Summarizing, we found five inequalities

$\displaystyle\eta_{\parallel}\geq 0\,,\ \ \ \eta_{\perp}\geq 0\,,\ \ \
\zeta_{\parallel}\geq 0\,,\ \ \ \zeta_{\perp}\geq 0\,,\ \ \
\zeta_{\parallel}\zeta_{\perp}-\zeta_{\times}^{2}\geq 0\,.$ (32)

Since the third and fifth (fourth and fifth) inequalities imply the fourth
(third) inequality, one can get rid of either third or fourth inequality. See
Appendix D for more discussions and comparisons to the results in the
literature. We also note that the second law of local thermodynamics does not
require the sign of $\zeta_{\times}$ to be semi-positive, but requires an
inequality among $\zeta_{\times}$ and $\zeta_{\parallel,\perp}$.

##  III Nonequilibrium statistical operator method for relativistic MHD

From the underlying quantum field theory, one can also derive RMHD equations
by generalizing the nonequilibrium statistical operator method which was
initiated by some Japanese physicists in 1950s Nakajima (1957); Mori (1958),
further developed in 1960-1970s McLennan (1960, 1988); Kawasaki and Gunton
(1973); Zubarev  _et al._ (1979, 1996, 1997), and sophisticated quite recently
Sasa (2014); Becattini  _et al._ (2015); Hayata  _et al._ (2015); Hongo (2017,
2019); Becattini  _et al._ (2019); Hongo and Hidaka (2019); Hongo and Hattori
(2021). We here review such an approach for deriving RMHD (see Ref. Hongo and
Hattori (2021) for more details).

###  III.1 Optimized perturbation with local Gibbs distribution

The vital point in the nonequilibrium statistical operator method is to
correctly identify the appropriate form of the density operator. As we already
discussed, the dynamical variables in RMHD are the energy-momentum density
$T^{0}_{~{}\mu}$ and the magnetic flux density $\widetilde{F}^{0\mu}$ in a
given coordinate system. Let us first identify these operators by considering
QED as an underlying quantum theory. The QED Lagrangian is given by

$\mathcal{L}_{\mathrm{QED}}=\mathrm{i}{\bar{\psi}}\gamma^{\mu}D_{\mu}\psi-m{\bar{\psi}}\psi-\frac{1}{4}\eta^{\mu\nu}\eta^{\alpha\beta}F_{\mu\alpha}F_{\nu\beta},$
(33)

where we introduced the Dirac field $\psi$ with electric charge $q$, its Dirac
conjugate ${\bar{\psi}}=\psi^{\dagger}\gamma^{0}$, and the $\mathrm{U}(1)$
gauge field $A_{\mu}$. We also defined the covariant derivative
$D_{\mu}\psi:=\partial_{\mu}\psi-\mathrm{i}qA_{\mu}\psi$ and the field
strength tensor $F_{\mu\nu}:=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. The
Poincaré symmetry and the Bianchi identity enables us to find that the
operators 666To keep the notations simple, we use the same symbol for operator
$\tilde{F}^{\mu\nu}$ and its expectation value,
$\langle\tilde{F}^{\mu\nu}\rangle$, the meaning should be self-explained in
the context.,

$\hat{T}^{\mu}_{~{}\nu}=\mathrm{i}{\bar{\psi}}\gamma^{\mu}D_{\nu}\psi+F^{\mu\rho}F_{\nu\rho}-\delta^{\mu}_{\nu}\mathcal{L}_{\mathrm{QED}},\quad\widetilde{F}^{\mu\nu}=\frac{1}{2}\varepsilon^{\mu\nu\rho\sigma}F_{\rho\sigma},$
(34)

satisfy the following Ward-Takahashi identities:

$\partial_{\mu}\hat{T}^{\mu}_{~{}\nu}(t,{\bm{x}})=0,\quad\partial_{\mu}\widetilde{F}^{\mu\nu}(t,{\bm{x}})=0.$
(35)

One may wonder how the Bianchi identity is related to the symmetry of QED. In
fact, it is not so obvious how the QED Lagrangian is equipped with the
corresponding symmetry — the magnetic $\mathrm{U}(1)$ one-form symmetry —
since it does not act on the local operator sitting in Eq. (33). Rather, it
acts on the line operator
$\displaystyle{\exp\left(\mathrm{i}\int_{C}\mathrm{d}x^{\mu}\widetilde{A}_{\mu}\right)}$,
called the t’ Hooft line operator composed of a dual gauge field
$\widetilde{A}_{\mu}$ defined by
$\widetilde{F}_{\mu\nu}=\partial_{\mu}\widetilde{A}_{\nu}-\partial_{\nu}\widetilde{A}_{\mu}$.
The insertion of the t’Hooft line is regarded as putting the test particle
with the magnetic charge and is analogous to the insertion of the Wilson line
Gaiotto  _et al._ (2015). Thus, considering Eq. (35) as the equations of
motion is on the canonical line to set the starting point in constructing a
symmetry-based effective theory.

We then identify $T^{\mu\nu}$ and $\widetilde{F}^{\mu\nu}$ in the previous
section with expectation values of the above quantum operators
$\hat{T}^{\mu\nu}$ and $\widetilde{F}^{\mu\nu}$. This identification motivates
us to specify the appropriate density operator as the one describing fixed
expectation values of our dynamical variables
$\langle\hat{T}^{0}_{~{}\nu}(t,{\bm{x}})\rangle$ and
$\langle\widetilde{F}^{0\mu}(t,{\bm{x}})\rangle$. The so-called local Gibbs
(LG) distribution $\hat{\rho}_{\mathrm{LG}}[\lambda_{t};t]$ realizes such an
density operator, which is parameterized by a set of the Lagrange multipliers
$\lambda_{t}:=\\{\beta^{\mu}(t,{\bm{x}}),\mathcal{H}_{\nu}(t,{\bm{x}})\\}$ as

$\hat{\rho}_{\mathrm{LG}}[\lambda_{t};t]:=\exp(-\hat{S}[\lambda_{t};t]),$ (36)

where we introduced the entropy functional operator

$\hat{S}[\lambda_{t};t]:=-\int\mathrm{d}^{3}x\left[\hat{T}^{0}_{~{}\nu}(t,{\bm{x}})\beta^{\nu}(t,{\bm{x}})+\widetilde{F}^{0\nu}(t,{\bm{x}})\mathcal{H}_{\nu}(t,{\bm{x}})\right]+\Psi[\lambda_{t}].$
(37)

Here, the first argument $\lambda_{t}$ in
$\hat{\rho}_{\mathrm{LG}}[\lambda_{t};t]$ and $\hat{S}[\lambda_{t};t]$
describes their functional dependence on the Lagrange multipliers,
$\lambda_{t}:=\\{\beta^{\mu}(t,{\bm{x}}),\mathcal{H}_{\nu}(t,{\bm{x}})\\}$,
while the second one $t$ represents the time argument for a set of the
conserved charge-density operators
$\hat{c}_{t}:=\\{\hat{T}^{0}_{~{}\nu}(t,{\bm{x}}),\widetilde{F}^{0\mu}(t,{\bm{x}})\\}$.
One finds that these Lagrange multipliers can be decomposed as
$\beta^{\mu}(t,{\bm{x}})=\beta(t,{\bm{x}})u^{\mu}(t,{\bm{x}})$ and
$\mathcal{H}_{\mu}(t,{\bm{x}})=\beta(t,{\bm{x}})H_{\mu}(t,{\bm{x}})$ with the
local inverse temperature $\beta(t,{\bm{x}})$, the fluid four-velocity
$u^{\mu}(t,{\bm{x}})$, and the magnetic field $H_{\mu}(t,{\bm{x}})$. In the
following, we express the average of a quantum operator $\hat{\mathcal{O}}$
over the LG distribution as

$\langle\hat{\mathcal{O}}\rangle^{\text{LG}}_{\lambda_{t}}:={\rm
Tr}(\hat{\rho}_{\mathrm{LG}}[\lambda_{t};t]\hat{\mathcal{O}}).$ (38)

In Eq. (37), we also defined the local thermodynamic functional
$\Psi[\lambda_{t}]$, called the Massieu-Planck functional, as a normalization
factor of the LG distribution

$\Psi[\lambda_{t}]:=\log{\rm
Tr}\exp\left(\int\mathrm{d}^{3}x\left[\hat{T}^{0}_{~{}\mu}(t,{\bm{x}})\beta^{\mu}(t,{\bm{x}})+\widetilde{F}^{0\mu}(t,{\bm{x}})\mathcal{H}_{\mu}(t,{\bm{x}})\right]\right).$
(39)

This functional is used to extract the average charge densities

$\langle\hat{T}^{0}_{~{}\mu}(t,{\bm{x}})\rangle^{\text{LG}}=\frac{\delta\Psi[\lambda_{t}]}{\delta\beta^{\mu}(t,{\bm{x}})},\quad\langle\widetilde{F}^{0\mu}(t,{\bm{x}})\rangle^{\text{LG}}=\frac{\delta\Psi[\lambda_{t}]}{\delta\mathcal{H}_{\mu}(t,{\bm{x}})}.$
(40)

Besides, we define the entropy functional $S[c_{t}]$ by taking the average of
Eq. (37) over $\hat{\rho}_{\mathrm{LG}}[\lambda_{t};t]$ as

$S[c_{t}]:=-\int\mathrm{d}^{3}x\left[\langle\hat{T}^{0}_{~{}\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}\beta^{\nu}(t,{\bm{x}})+\langle\widetilde{F}^{0\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}\mathcal{H}_{\nu}(t,{\bm{x}})\right]+\Psi[\lambda_{t}].$
(41)

In other words, the entropy functional $S[c_{t}]$ is defined by the Legendre
transform of the Massieu-Planck functional $\Psi[\lambda_{t}]$, and thus, its
argument is the averaged charge densities
$c_{t}=\\{\langle\hat{T}^{0}_{~{}\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}},\langle\widetilde{F}^{0\mu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}\\}$.
One then finds that the Lagrange multipliers, or local thermodynamic variables
conjugate to the conserved charge densities, as

$\beta^{\mu}(t,{\bm{x}})=-\frac{\delta S[c_{t}]}{\delta
T^{0}_{~{}\mu}(t,{\bm{x}})},\quad\mathcal{H}_{\mu}(t,{\bm{x}})=-\frac{\delta
S[c_{t}]}{\delta\widetilde{F}^{0\mu}(t,{\bm{x}})},$ (42)

which are consistent with Eq. (14) in the previous section.

To describe the dissipative transport phenomena with the nonequilibrium
statistical operator method, we require a crucial assumption that the density
operator $\hat{\rho}_{0}$ at the initial time $t_{0}$ be given by the LG
distribution $\hat{\rho}_{0}=\hat{\rho}_{\mathrm{LG}}[\lambda_{t_{0}};t_{0}]$
Zubarev  _et al._ (1979, 1996, 1997); Sasa (2014); Hayata  _et al._ (2015);
Hongo and Hattori (2021). Using the Heisenberg picture, we define the
expectation value of an arbitrary Heisenberg operator $\hat{\mathcal{O}}(t)$
at time $t(\geq t_{0})$ as

$\langle\hat{\mathcal{O}}(t)\rangle:={\rm
Tr}[\hat{\rho}_{0}\hat{\mathcal{O}}(t)]=\langle\hat{\mathcal{O}}(t)\rangle^{\text{LG}}_{\lambda_{t_{0}}},$
(43)

where we expressed the rightmost side using the above assumption and Eq. (38).
Then, according to our identification, the averaged Ward-Takahashi identities,

$\partial_{\mu}\langle\hat{T}^{\mu}_{~{}\nu}(t,{\bm{x}})\rangle=0,\quad\partial_{\mu}\langle\widetilde{F}^{\mu\nu}(t,{\bm{x}})\rangle=0,$
(44)

should provide the RMHD equations after an appropriate derivative expansion is
employed.

From this microscopic point of view, we have already fixed the definition of
expectation values $\langle\hat{T}^{\mu}_{~{}\nu}(t,{\bm{x}})\rangle$ and
$\langle\widetilde{F}^{\mu\nu}(t,{\bm{x}})\rangle$ by Eq. (43). Thus, the
remaining problem is to derive the constitutive relations for Eq. (43) based
on our density operator
$\hat{\rho}_{0}=\hat{\rho}_{\mathrm{LG}}[\lambda_{t_{0}};t_{0}]$. However,
recalling the result in the previous section, one realizes that this is a
tough problem; We expect that the resulting constitutive relations should be
expressed by the conjugate variables $\lambda_{t}$ at time $t$, whereas our
density operator
$\hat{\rho}_{0}=\hat{\rho}_{\mathrm{LG}}[\lambda_{t_{0}};t_{0}]$ only contains
those $\lambda_{t_{0}}$ at the initial time $t_{0}$. In fact, the expectation
value at time $t\,(>t_{0})$ in the present setup is always defined by taking
average over the initial density operator
$\hat{\rho}_{0}=\hat{\rho}_{\mathrm{LG}}[\lambda_{t_{0}};t_{0}]$, which only
contains conjugate variables at the initial time. Thus, it sounds impossible
to express $\langle\hat{T}^{\mu}_{~{}\nu}(t,{\bm{x}})\rangle$ and
$\langle\widetilde{F}^{\mu\nu}(t,{\bm{x}})\rangle$ in terms of $\lambda_{t}$
at time $t$ as we did in the previous section.

The above observation implies that the initial density operator
$\hat{\rho}_{0}=\hat{\rho}_{\mathrm{LG}}[\lambda_{t_{0}};t_{0}]$ does not give
a useful starting point to evaluate the expectation value at later time
$t(>t_{0})$. Instead, one immediately finds that the better starting point is
the local Gibbs distribution at time $t$. The question is how we can shift
into such a different distribution when we do not even have conjugate
variables other than those at initial time $t_{0}$.

We can resolve this problem by invoking the optimized (or renormalized)
perturbation theory (see, e.g., Refs. Stevenson (1981); Kleinert (2009);
Jakovác and Patkós (2015)). Suppose that we know the configuration of the
conserved charge densities at time $t$. With the help of the entropy
functional $S[c_{t}]$, we first define the conjugate variable $\lambda_{t}$ by
Eq. (42). One can show that this definition is equivalent to requiring the
following matching conditions

$\langle T^{0}_{~{}\nu}\rangle=\langle
T^{0}_{~{}\nu}\rangle^{\text{LG}}_{\lambda_{t}},\quad\langle\widetilde{F}^{0\nu}\rangle=\langle\widetilde{F}^{0\nu}\rangle^{\text{LG}}_{\lambda_{t}}.$
(45)

Using the defined conjugate variables $\lambda_{t}$, we rearrange our density
operator as

$\hat{\rho}_{0}=\mathrm{e}^{-\hat{S}[\lambda_{t};t]+\hat{\Sigma}[t,t_{0};\lambda]}=\mathrm{e}^{-\hat{S}[\lambda_{t};t]}\mathrm{T}_{\tau}\exp\left(\int_{0}^{1}\mathrm{d}\tau\mathrm{e}^{\tau\hat{S}[\lambda_{t};t]}\hat{\Sigma}[t,t_{0};\lambda]\mathrm{e}^{-\tau\hat{S}[\lambda_{t};t]}\right),$
(46)

where we just added and subtracted the entropy operator
$\hat{S}[t;\lambda_{t}]$ at time $t$ and defined the entropy production
operator $\hat{\Sigma}[t,t_{0};\lambda]$ as

$\hat{\Sigma}[t,t_{0};\lambda]:=\hat{S}[\lambda_{t};t]-\hat{S}[\lambda_{t_{0}};t_{0}].$
(47)

The rightmost side of Eq. (46) gives a useful formula for the derivative
expansion since the entropy production operator
$\hat{\Sigma}[t,t_{0};\lambda]$ will be shown to be $O(\partial^{1})$. Thus,
we can regard $\hat{\Sigma}[t,t_{0};\lambda]$ as a derivative correction, and
Eq. (46) gives a familiar perturbative expansion formula in the interacting
picture, which we learn in the elementary course of quantum mechanics. As a
result, we see that the nonequilibrium statistical operator method gives one
useful expansion scheme relying on the new optimized (or renormalized)
parameter $\lambda_{t}$.

Expanding the rightmost side of Eq. (46) at the first-order in
$\hat{\Sigma}[t,t_{0};\lambda]$, we obtain

$\begin{split}\langle\hat{T}^{\mu\nu}(t,{\bm{x}})\rangle&=\langle\hat{T}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}+\bm{(}{\hat{T}^{\mu\nu}(t,{\bm{x}}),\hat{\Sigma}[t,t_{0};\lambda]}\bm{)}_{t}+O(\hat{\Sigma}^{2}),\\\
\langle\widetilde{F}^{\mu\nu}(t,{\bm{x}})\rangle&=\langle\widetilde{F}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}+\bm{(}{\widetilde{F}^{\mu\nu}(t,{\bm{x}}),\hat{\Sigma}[t,t_{0};\lambda]}\bm{)}_{t}+O(\hat{\Sigma}^{2}),\end{split}$
(48)

where we introduced the Kubo-Mori-Bogoliubov inner product:

$\bm{(}{\hat{A},\hat{B}}\bm{)}_{t}:=\int_{0}^{1}\mathrm{d}\tau\langle\mathrm{e}^{\tau\hat{S}[\lambda_{t};t]}\hat{A}\mathrm{e}^{-\tau\hat{S}[\lambda_{t};t]}\hat{B}^{{\dagger}}\rangle^{\text{LG}}_{\lambda_{t}}.$
(49)

Equation (48) indicates that we have separated the problem into two parts: The
first one is to evaluate the expectation value of currents with the LG
distribution describing the local thermal equilibrium

$T^{\mu\nu}_{(0)}=\langle\hat{T}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}},\quad\widetilde{F}^{\mu\nu}_{(0)}=\langle\widetilde{F}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}},$
(50)

which will be shown to contain the leading-order terms in derivative. The
second one is to find the dissipative corrections

$T^{\mu\nu}_{(1)}=\bm{(}{\hat{T}^{\mu\nu}(t,{\bm{x}}),\hat{\Sigma}[t,t_{0};\lambda]}\bm{)}_{t},\quad\widetilde{F}^{\mu\nu}_{(1)}=\bm{(}{\widetilde{F}^{\mu\nu}(t,{\bm{x}}),\hat{\Sigma}[t,t_{0};\lambda]}\bm{)}_{t},$
(51)

by computing the entropy production.

###  III.2 Evaluating the Local Gibbs averages

Let us first investigate the LG-averaged currents given in Eq. (50). In
evaluating these expectation values, it is useful to put our system in the
curved spacetime described by the vierbein $e_{\mu}^{~{}a}$ and introduce a
background two-form gauge field $b_{\mu\nu}$ that couples to
$\widetilde{F}^{\mu\nu}$. In the presence of these background fields, one can
show the following variational formulas Hongo and Hattori (2021)

$\begin{split}\langle\hat{T}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}&=\frac{1}{\beta^{0}\sqrt{\gamma}}\frac{\delta\Psi[\lambda_{t}]}{\delta
e_{\mu}^{~{}a}(t,{\bm{x}})}e^{\nu}_{~{}a},\\\
\langle\widetilde{F}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}&=\frac{2}{\beta^{0}\sqrt{\gamma}}\left(\beta^{[\mu}\delta^{\nu]}_{\rho}\frac{\delta\Psi[\lambda_{t}]}{\delta\mathcal{H}_{\rho}(t,{\bm{x}})}+\frac{\delta\Psi[\lambda_{t}]}{\delta
b_{\mu\nu}(t,{\bm{x}})}\right),\end{split}$ (52)

where $e^{\nu}_{~{}a}$ is the inverse of vierbein and $\beta^{0}$ denotes the
zeroth-component of the four vector $\beta^{\mu}$. Here, we also used
$\gamma:=\det g_{ij}$ with a spatial part of the metric
$g_{ij}=e_{i}^{~{}a}e_{j}^{~{}b}\eta_{ab}$. Thus, we can derive the LG
expectation values
$\langle\hat{T}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}$ and
$\langle\widetilde{F}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}$
once we get the form of the local thermodynamic functional $\Psi[\lambda_{t}]$
under the background fields.

To specify the form of the Massieu-Planck functional $\Psi[\lambda_{t}]$, we
can rely on the path-integral formula Hongo (2017); Hongo and Hattori (2021).
Substituting Eq. (34) into the definition (39) of $\Psi[\lambda_{t}]$ and
following the usual procedure of deriving the path-integral representation for
the partition function, we obtain

$\Psi[\lambda_{t}]=\log\int\mathcal{D}\psi\mathcal{D}{\bar{\psi}}\mathcal{D}A_{\mu}\delta(F)\det(\partial
F/\partial\alpha)\exp\left(\int_{0}^{\beta_{\mathrm{ref}}}\mathrm{d}\tau\mathrm{d}^{3}x\sqrt{-\tilde{g}}\widetilde{\mathcal{L}}_{\mathrm{QED}}\right),$
(53)

where $\beta_{\mathrm{ref}}$ is an arbitrary constant reference inverse
temperature and $F$ denotes a gauge fixing condition $F=0$ with a gauge
parameter $\alpha$. Due to the inhomogeneity of the local thermodynamic
variables $\lambda_{t}$, we need to perform the path-integral not in the flat
Euclidean spacetime but in the emergent curved spacetime with the two-form
gauge field. In fact, we find the Lagrangian density in Eq. (53) to be

$\widetilde{\mathcal{L}}_{\mathrm{QED}}=\mathrm{i}{\bar{\psi}}\gamma^{a}\tilde{e}_{a}^{~{}\mu}\widetilde{D}_{\mu}\psi-m{\bar{\psi}}\psi-\frac{1}{4}\tilde{g}^{\mu\nu}\tilde{g}^{\alpha\beta}F_{\mu\alpha}F_{\nu\beta}+\frac{1}{2}\tilde{b}_{\mu\nu}\widetilde{F}^{\mu\nu},$
(54)

and the background fields — the thermal vierbein $\tilde{e}_{\mu}^{~{}a}$ and
thermal two-form gauge field $\tilde{b}_{\mu\nu}$ — are given by

$\tilde{e}_{0}^{~{}a}:=\mathrm{e}^{\sigma}u^{a},\quad\tilde{e}_{i}^{~{}a}:=e_{i}^{~{}a},\quad\tilde{b}_{0i}:=\mathrm{e}^{\sigma}H_{i},\quad\tilde{b}_{ij}:=b_{ij},$
(55)

where we defined $\mathrm{e}^{\sigma(x)}=\beta(x)/\beta_{\mathrm{ref}}$. We
also defined the thermal metric
$\tilde{g}_{\mu\nu}:=\tilde{e}_{\mu}^{~{}a}\tilde{e}_{\nu}^{~{}b}\eta_{ab}$
with $\tilde{g}:=\det\tilde{g}_{\mu\nu}$, and the covariant derivative

$\widetilde{D}_{\mu}\psi:=\tilde{\partial}_{\mu}\psi-\mathrm{i}qA_{\mu}\psi-\frac{\mathrm{i}}{2}\widetilde{\omega}_{\mu}^{~{}ab}\Sigma_{ab}\psi,$
(56)

with the representation matrix of the Dirac spinor under the Lorentz
transformation $\Sigma_{ab}:=\mathrm{i}[\gamma_{a},\gamma_{b}]/4$. The
derivative in the thermal space is given as
$\tilde{\partial}_{\mu}=(\mathrm{i}\partial_{\tau},\partial_{i})$ and the spin
connection $\widetilde{\omega}_{\mu}^{~{}ab}$ is determined by the thermal
vierbein $\tilde{e}_{\mu}^{~{}a}$ as Hongo (2017)

$\widetilde{\omega}_{\mu}^{~{}ab}=\frac{1}{2}\tilde{e}^{a\nu}\tilde{e}^{b\rho}(\widetilde{C}_{\nu\rho\mu}-\widetilde{C}_{\rho\nu\mu}-\widetilde{C}_{\mu\nu\rho})\quad{\rm
with}\quad\widetilde{C}_{\mu\nu\rho}:=\tilde{e}_{\mu}^{~{}c}\big{(}\tilde{\partial}_{\nu}\tilde{e}_{\rho
c}-\tilde{\partial}_{\rho}\tilde{e}_{\nu c}\big{)}.$ (57)

In short, the Massieu-Planck functional is described by performing the path
integral for QED in the presence of the curved and two-form backgrounds. This
result is a generalization of the background field method to the locally
thermalized system. We can read off the symmetry properties of the Massieu-
Planck functional as follows. To see this, note that the line element
$\mathrm{d}\tilde{s}^{2}$ and two-form gauge connection $\tilde{b}$,

$\displaystyle\mathrm{d}\tilde{s}^{2}$
$\displaystyle:=\tilde{g}_{\mu\nu}\mathrm{d}\tilde{x}^{\mu}\otimes\mathrm{d}\tilde{x}^{\nu}=-\mathrm{e}^{2\sigma}(\mathrm{d}\tilde{t}+a_{i}\mathrm{d}x^{i})^{2}+g^{\prime}_{ij}\mathrm{d}x^{i}\otimes\mathrm{d}x^{j},$
(58) $\displaystyle\tilde{b}$
$\displaystyle:=\frac{1}{2}\tilde{b}_{\mu\nu}\mathrm{d}\tilde{x}^{\mu}\wedge\mathrm{d}\tilde{x}^{\nu}=\tilde{b}_{0i}(\mathrm{d}\tilde{t}+a_{j}\mathrm{d}x^{j})\wedge\mathrm{d}x^{i}+\frac{1}{2}\tilde{b}_{ij}^{\prime}\mathrm{d}x^{i}\wedge\mathrm{d}x^{j},$
(59)

describe the backgrounds. Here, we defined
$(\mathrm{d}\tilde{t},\mathrm{d}\widetilde{\bm{x}})=(-\mathrm{i}\mathrm{d}\tau,\mathrm{d}{\bm{x}})$
and expressed $\mathrm{d}\tilde{s}^{2}$ and $\tilde{b}$ using the Kaluza-Klein
parameterization

$a_{i}:=-\mathrm{e}^{-\sigma}u_{i},\quad
g^{\prime}_{ij}:=g_{ij}+u_{i}u_{j},\quad\tilde{b}_{ij}^{\prime}:=b_{ij}-\tilde{b}_{0j}a_{i}-\tilde{b}_{i0}a_{j}.$
(60)

The apparently complicated Kaluza-Klein parameterization is, indeed, useful
because the Massieu-Planck functional is invariant under the Kaluza-Klein
gauge transformation

$\begin{split}\tilde{t}&\to\tilde{t}^{\prime}=\tilde{t}+\chi({\bm{x}}),\\\
{\bm{x}}&\to{\bm{x}}^{\prime}={\bm{x}},\\\ a_{i}({\bm{x}})&\to
a_{i}^{\prime}({\bm{x}})=a_{i}({\bm{x}})-\partial_{i}\chi({\bm{x}}),\end{split}$
(61)

and $\mathrm{e}^{\sigma}$, $\gamma_{ij}^{\prime}$ and $\tilde{b}_{0i}$ and
$\tilde{b}_{ij}^{\prime}$ are all invariant under the Kaluza-Klein gauge
transformation. Besides, the Massieu-Planck functional is invariant under the
spatial diffeomorphism and gauge transformation acting on $b_{ij}$,

$\displaystyle x^{i}$ $\displaystyle\to x^{i\prime}=f^{i}({\bm{x}}),$ (62)
$\displaystyle\tilde{b}_{ij}$
$\displaystyle\to\tilde{b}_{ij}+\partial_{i}\theta_{j}-\partial_{j}\theta_{i}.$
(63)

The crucial point here is that the Massieu-Planck functional has to respect
the symmetries under the transformations (61)-(63). Using this symmetry
property and relying on the derivative expansion, one can write down the most
general form of the Massieu-Planck functional in an order-by-order basis in
derivative. In the leading-order expansion, we have two zeroth-order invariant
scalars $\mathrm{e}^{\sigma(x)}=\beta(x)/\beta_{\mathrm{ref}}$ and
$\sqrt{\big{(}\tilde{b}_{0i}(x)\big{)}^{2}}=\mathrm{e}^{\sigma(x)}H(x)=\mathcal{H}(x)/\beta_{\mathrm{ref}}$,
where we decomposed the magnetic field as $H_{\mu}=-Hb_{\mu}$ and
$\mathcal{H}_{\mu}=-\mathcal{H}b_{\mu}$ with the normalized spatial vector
$b^{\mu}$. Moreover, there is no invariant scalar at $O(\partial^{1})$, and
thus, the most general form of $\Psi[\lambda]$ in the leading-order expansion
is given by

$\begin{split}\Psi[\lambda]&=\int_{0}^{\beta_{\mathrm{ref}}}\mathrm{d}\tau\mathrm{d}^{3}x\sqrt{-\tilde{g}}p(\beta,\mathcal{H})+O(\partial^{2})=\int\mathrm{d}^{3}x\sqrt{\gamma}\beta^{0}p(\beta,\mathcal{H})+O(\partial^{2}).\end{split}$
(64)

Recalling the variational formula (52), we find that the functional derivative
of Eq. (64) leads to

$\begin{split}\langle\hat{T}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}&=\epsilon
u^{\mu}u^{\nu}-p_{\perp}(\eta^{\mu\nu}-u^{\mu}u^{\nu}+b^{\mu}b^{\nu})+p_{\parallel}b^{\mu}b^{\nu},\\\
\langle\widetilde{F}^{\mu\nu}(t,{\bm{x}})\rangle^{\text{LG}}_{\lambda_{t}}&=B(b^{\mu}u^{\nu}-b^{\nu}u^{\mu}),\end{split}$
(65)

where we have taken the flat background limit and introduced a set of the
scalar functions

$\epsilon=-\frac{\partial(\beta p)}{\partial\beta},\quad p_{\perp}=p,\quad
p_{\parallel}=p-\beta^{-1}\mathcal{H}B,\quad B=\beta\frac{\partial
p}{\partial\mathcal{H}}.$ (66)

All of them can be extracted from the single function $p(\beta,\mathcal{H})$.
Equations (65)-(66) give the leading-order constitutive relations of RMHD,
which agree with Eq. (19) in the previous section (recall
$\beta^{-1}\mathcal{H}B=H^{\mu}B_{\mu}$). In contrast to the entropy-current
analysis, we now have the microscopic path-integral formula for
$p(\beta,\mathcal{H})$. One can thus, in principle, compute all coefficient
functions in Eq. (66), or the equations of state, from the underlying
microscopic theory, i.e., QED.

###  III.3 Evaluating the dissipative corrections

We next evaluate the dissipative corrections given in Eq. (51). For this
purpose, we first rewrite the entropy production operator by using the Ward-
Takahshi identity (35) and performing integration by parts. The resultant
expression reads

$\displaystyle\hat{\Sigma}[t,t_{0};\lambda]$
$\displaystyle:=\hat{S}[\lambda_{t};t]-\hat{S}[\lambda_{t_{0}};t_{0}]=\int_{t_{0}}^{t}\mathrm{d}t^{\prime}\partial_{t^{\prime}}\hat{S}[\lambda_{t^{\prime}};t^{\prime}]$
$\displaystyle=-\int_{t_{0}}^{t}\mathrm{d}t^{\prime}\partial_{t^{\prime}}\left[\int\mathrm{d}^{3}x^{\prime}\left(\beta^{\nu}(t^{\prime},{\bm{x}}^{\prime})\hat{T}^{0}_{~{}\nu}(t^{\prime},{\bm{x}}^{\prime})+\mathcal{H}_{\nu}(t^{\prime},{\bm{x}}^{\prime})\widetilde{F}^{0\nu}(t^{\prime},{\bm{x}}^{\prime})\right)+\Psi[\lambda_{t^{\prime}}]\right]$
$\displaystyle=-\int_{t_{0}}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\left[\delta\hat{T}^{\mu\nu}(t^{\prime},{\bm{x}}^{\prime})\partial_{\mu}\beta_{\nu}(t^{\prime},{\bm{x}}^{\prime})+\delta\widetilde{F}^{\mu\nu}(t^{\prime},{\bm{x}}^{\prime})\partial_{\mu}\mathcal{H}_{\nu}(t^{\prime},{\bm{x}}^{\prime})\right],$
(67)

where we defined a deviation
$\delta\hat{\mathcal{O}}(t):=\hat{\mathcal{O}}(t)-\langle\hat{\mathcal{O}}(t)\rangle^{\text{LG}}_{\lambda_{t}}$.
To obtain the third line, we also used the following identity:

$\begin{split}\partial_{t}\Psi[\lambda_{t}]&=-\Big{\langle}\partial_{t}\int\mathrm{d}^{3}x\left(\beta^{\nu}\hat{T}^{0}_{~{}\nu}+\mathcal{H}_{\nu}\widetilde{F}^{0\nu}\right)\Big{\rangle}^{\mathrm{LG}}_{\lambda_{t}}\\\
&=\int\mathrm{d}^{3}x\left(\langle\hat{T}^{\mu\nu}\rangle^{\text{LG}}_{\lambda_{t}}\partial_{\mu}\beta_{\nu}+\langle\widetilde{F}^{\mu\nu}\rangle^{\text{LG}}_{\lambda_{t}}\partial_{\mu}\mathcal{H}_{\nu}\right).\end{split}$
(68)

Note that the entropy production operator $\hat{\Sigma}[t,t_{0};\lambda]$ in
Eq. (67) contains the time derivative of the parameters
$\lambda_{t^{\prime}}$. This can be explicitly seen by using the projection
tensor, $\delta^{\mu}_{\nu}=u^{\mu}u_{\nu}+\Delta^{\mu}_{\nu}$, that
decomposes a derivative as

$\partial_{\mu}=u_{\mu}D+\nabla_{\mu}\quad{\rm with}\quad
D:=u^{\mu}\partial_{\mu},\quad\nabla_{\mu}:=\Delta^{\nu}_{\mu}\partial_{\nu}.$
(69)

Substituting this decomposition together with $\beta^{\mu}=\beta u^{\mu}$ and
$\mathcal{H}_{\nu}=-\mathcal{H}b_{\nu}$, we obtain the following expression
for the entropy production operator

$\begin{split}\hat{\Sigma}[t,t_{0};\lambda]=-\int_{t_{0}}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}&\Big{[}u_{\mu}u_{\nu}\delta\hat{T}^{\mu\nu}D\beta+u_{\mu}\delta\hat{T}^{\mu\nu}\beta
Du_{\nu}-u_{\mu}b_{\nu}\delta\widetilde{F}^{\mu\nu}D\mathcal{H}-\mathcal{H}u_{\mu}\delta\widetilde{F}^{\mu\nu}Db_{\nu}\\\
&+\beta\delta\hat{T}^{\mu\nu}\nabla_{\mu}u_{\nu}+u_{\nu}\delta\hat{T}^{\mu\nu}\nabla_{\mu}\beta+\mathcal{H}\delta\widetilde{F}^{\mu\nu}b^{\rho}u_{\nu}\nabla_{\mu}u_{\rho}+\delta\widetilde{F}^{\mu\nu}\Delta_{\nu}^{\rho}\nabla_{\mu}\mathcal{H}_{\rho}\Big{]}.\end{split}$
(70)

Here, it is important to eliminate time-derivative terms in the first line.
Otherwise, the resulting Green-Kubo formulas pick up contributions from
gapless linear hydrodynamic modes, which prevents us from obtaining convergent
time integrals for the transport coefficients. In Ref. Hongo and Hattori
(2021), we accomplish such a procedure by a formal manipulation. To complement
that formal manipulation, we will here explicitly demonstrate this procedure.

To eliminate the time-derivative terms, we solve the leading-order RMHD
equations following the leading-order constitutive relations (65). Using Eq.
(65), we find the leading-order equations of motion

$\displaystyle 0$
$\displaystyle=\partial_{\mu}T^{\mu\nu}_{(0)}=\partial_{\mu}\big{[}\epsilon
u^{\mu}u^{\nu}-p_{\perp}(\eta^{\mu\nu}-u^{\mu}u^{\nu})+B^{\mu}H^{\nu}\big{]},$
(71a) $\displaystyle 0$
$\displaystyle=\partial_{\mu}\widetilde{F}^{\mu\nu}_{(0)}=\partial_{\mu}[B^{\mu}u^{\nu}-B^{\nu}u^{\mu}],$
(71b)

where we used $p_{\parallel}=p_{\perp}-HB$ to express $p_{\parallel}$ in terms
of $B^{\mu}=Bb^{\nu}$ and $H^{\nu}=-Hb^{\nu}$. After contracting these ideal
RMHD equations with appropriate tensors such as $b_{\mu}$, we obtain the
following set of equations:

$\displaystyle D\epsilon$
$\displaystyle=-(\epsilon+p_{\perp})\theta+HB\theta_{\parallel},$ (72a)
$\displaystyle(\epsilon+p_{\perp})Du_{\nu}$
$\displaystyle=\nabla_{\nu}p_{\perp}-H_{\nu}\partial_{\mu}B^{\mu}-\Delta_{\nu\rho}B^{\mu}\partial_{\mu}H^{\rho},$
(72b) $\displaystyle DB$ $\displaystyle=-B\theta_{\perp},$ (72c)
$\displaystyle Db^{\nu}$
$\displaystyle=\Xi^{\nu\rho}b^{\mu}\partial_{\mu}u_{\rho}+B^{-1}u^{\nu}\partial_{\mu}B^{\mu},$
(72d)

where we used $\theta=\theta_{\parallel}+\theta_{\perp}$ with
$\theta_{\parallel}=-b^{\mu}b^{\nu}\partial_{\mu}u_{\nu}$ and
$\theta_{\perp}=\Xi^{\mu\nu}\partial_{\mu}u_{\nu}$. Combining these equations
with thermodynamic relations, we can further simplify $Du_{\nu}$ and also find
the time derivatives of conjugate variables $D\beta$ and $D\mathcal{H}$ as
(see Appendix C for a derivation)

$\displaystyle D\beta$ $\displaystyle=\beta\left(\frac{\partial
p_{\perp}}{\partial\epsilon}-B\frac{\partial
H}{\partial\epsilon}\right)\theta_{\parallel}+\beta\frac{\partial
p_{\perp}}{\partial\epsilon}\theta_{\perp},$ (73a) $\displaystyle
D\mathcal{H}$ $\displaystyle=-\beta\frac{\partial p_{\perp}}{\partial
B}\theta_{\perp}-\beta\left(\frac{\partial p_{\perp}}{\partial
B}-B\frac{\partial H}{\partial B}\right)\theta_{\parallel},$ (73b)
$\displaystyle Du_{\nu}$
$\displaystyle=-\beta^{-1}\nabla_{\nu}\beta-\frac{1}{\epsilon+p_{\perp}}\left[2\beta^{-1}Bb^{\mu}\nabla_{[\mu}\mathcal{H}_{\nu]}+\theta_{\parallel}HBu_{\nu}\right],$
(73c) $\displaystyle Db_{\nu}$
$\displaystyle=\Xi_{\nu}^{~{}\rho}b^{\mu}\nabla_{\mu}u_{\rho}+B^{-1}u_{\nu}\partial_{\mu}B^{\mu}.$
(73d)

These leading-order equations of motion for thermodynamic parameters enable us
to eliminate the time-derivative terms in the entropy production operator.
Substituting them into Eq. (70) and rearranging terms, we eventually obtain
the entropy production operator

$\begin{split}\hat{\Sigma}[t,t_{0};\lambda]=-\int_{t_{0}}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x&\Big{[}\tilde{\delta}\hat{p}_{\parallel}\beta\theta_{\parallel}+\tilde{\delta}\hat{p}_{\perp}\beta\theta_{\perp}+2\tilde{\delta}\hat{\pi}^{(\mu}b^{\nu)}\beta\nabla_{\mu}u_{\nu}+\tilde{\delta}\hat{\tau}^{\mu\nu}\beta\nabla_{\mu}u_{\nu}\\\
&+2\tilde{\delta}\hat{E}^{[\mu}b^{\nu]}\nabla_{\mu}\mathcal{H}_{\nu}+\tilde{\delta}\hat{D}^{\mu\nu}\nabla_{\mu}\mathcal{H}_{\nu}+O(\nabla^{2})\Big{]},\end{split}$
(74)

where we defined the projected components of the operators
$\tilde{\delta}\hat{T}^{\mu\nu}$ and
$\tilde{\delta}\widetilde{F}^{\rho\sigma}$ as

$\displaystyle\tilde{\delta}\hat{p}_{\parallel}$
$\displaystyle=-b_{\mu}b_{\nu}\delta\hat{T}^{\mu\nu}+u_{\mu}u_{\nu}\delta\hat{T}^{\mu\nu}\left(\frac{\partial
p_{\perp}}{\partial\epsilon}-B\frac{\partial
H}{\partial\epsilon}-\frac{HB}{\epsilon+p_{\perp}}+u_{\mu}b_{\nu}\delta\widetilde{F}^{\mu\nu}\frac{\partial
p_{\parallel}}{\partial B}\right),$ (75a)
$\displaystyle\tilde{\delta}\hat{p}_{\perp}$
$\displaystyle=\frac{1}{2}\delta\hat{T}^{\mu\nu}\Xi_{\mu\nu}+u_{\mu}u_{\nu}\delta\hat{T}^{\mu\nu}\frac{\partial
p_{\perp}}{\partial\epsilon}+u_{\mu}b_{\nu}\delta\widetilde{F}^{\mu\nu}\frac{\partial
p_{\perp}}{\partial B},$ (75b) $\displaystyle\tilde{\delta}\hat{\pi}^{\mu}$
$\displaystyle=\left(-b_{\sigma}\delta\hat{T}^{\rho\sigma}+H\delta\widetilde{F}^{\rho\sigma}u_{\sigma}\right)\Xi^{~{}\mu}_{\rho},$
(75c) $\displaystyle\tilde{\delta}\hat{\tau}^{\mu\nu}$
$\displaystyle=\Big{(}\Xi^{\mu}_{~{}\rho}\Xi^{\nu}_{~{}\sigma}-\frac{1}{2}\Xi^{\mu\nu}\Xi_{\rho\sigma}\Big{)}\delta\hat{T}^{\rho\sigma},$
(75d) $\displaystyle\tilde{\delta}\hat{E}^{\mu}$
$\displaystyle=b_{\rho}\delta\widetilde{F}^{\rho\sigma}\Xi_{\sigma}^{~{}\mu}+\frac{2B}{\epsilon+p_{\perp}}u_{\rho}\delta\hat{T}^{\rho\mu},$
(75e) $\displaystyle\tilde{\delta}\hat{D}^{\mu\nu}$
$\displaystyle=\Xi^{\mu}_{~{}\rho}\Xi^{\nu}_{~{}\sigma}\delta\widetilde{F}^{\rho\sigma}.$
(75f)

We then substitute the obtained entropy production operator into Eq. (51).
Assuming that the correlation of projected operators
$\tilde{\delta}\hat{T}^{\mu\nu}$ and $\tilde{\delta}\widetilde{F}^{\mu\nu}$
decays with the microscopic scales, we perform the Markovian approximation for
the integration kernel. As we emphasized, this approximation does not work if
we do not solve the ideal RMHD equations to obtain the projected operator
(75). After this procedure, we eventually obtain the dissipative corrections
to the constitutive relations

$\displaystyle T^{\mu\nu}_{(1)}=$
$\displaystyle-(\zeta_{\parallel}\theta_{\parallel}+\zeta_{\times}\theta_{\perp})b^{\mu}b^{\nu}+(\zeta_{\perp}\theta_{\perp}+\zeta_{\times}^{\prime}\theta_{\parallel})\Xi^{\mu\nu}$
$\displaystyle-2\eta_{\parallel}\left(b^{\mu}\Xi^{\nu(\rho}b^{\sigma)}+b^{\nu}\Xi^{\mu(\rho}b^{\sigma)}\right)\nabla_{\rho}u_{\sigma}+2\eta_{\perp}\nabla^{\\{\mu}u^{\nu\\}},$
(76a) $\displaystyle\widetilde{F}^{\mu\nu}_{(1)}=$ $\displaystyle
2\rho_{\perp}T\left(b^{\mu}\Xi^{\nu[\rho}b^{\sigma]}-b^{\nu}\Xi^{\mu[\rho}b^{\sigma]}\right)\nabla_{\rho}\mathcal{H}_{\sigma}+2\rho_{\parallel}T\Xi^{\mu[\rho}\Xi^{\sigma]\nu}\nabla_{\rho}\mathcal{H}_{\sigma}.$
(76b)

Those tensor structures are the same as those in Eqs. (23) and (30). We
defined a set of transport coefficients which are expressed in the form of the
spacetime integral of the Kubo-Mori-Bogoliubov inner product:

$\begin{split}\zeta_{\parallel}&=\beta(t,{\bm{x}})\int_{-\infty}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\bm{(}{\tilde{\delta}\hat{p}_{\parallel}(t,{\bm{x}}),\tilde{\delta}\hat{p}_{\parallel}(t^{\prime},{\bm{x}}^{\prime})}\bm{)}_{t},\\\
\zeta_{\perp}&=\beta(t,{\bm{x}})\int_{-\infty}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\bm{(}{\tilde{\delta}\hat{p}_{\perp}(t,{\bm{x}}),\tilde{\delta}\hat{p}_{\perp}(t^{\prime},{\bm{x}}^{\prime})}\bm{)}_{t},\\\
\zeta_{\times}&=\beta(t,{\bm{x}})\int_{-\infty}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\bm{(}{\tilde{\delta}\hat{p}_{\parallel}(t,{\bm{x}}),\tilde{\delta}\hat{p}_{\perp}(t^{\prime},{\bm{x}}^{\prime})}\bm{)}_{t},\\\
\zeta_{\times}^{\prime}&=\beta(t,{\bm{x}})\int_{-\infty}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\bm{(}{\tilde{\delta}\hat{p}_{\perp}(t,{\bm{x}}),\tilde{\delta}\hat{p}_{\parallel}(t^{\prime},{\bm{x}}^{\prime})}\bm{)}_{t},\\\
\eta_{\parallel}&=\frac{\beta(t,{\bm{x}})}{2}\int_{-\infty}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\bm{(}{\tilde{\delta}\hat{\pi}^{\mu}(t,{\bm{x}}),\tilde{\delta}\hat{\pi}^{\nu}(t^{\prime},{\bm{x}}^{\prime})}\bm{)}_{t}\Xi_{\mu\nu},\\\
\eta_{\perp}&=\frac{\beta(t,{\bm{x}})}{4}\int_{-\infty}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\bm{(}{\tilde{\delta}\hat{\tau}^{\mu\nu}(t,{\bm{x}}),\tilde{\delta}\hat{\tau}^{\rho\sigma}(t^{\prime},{\bm{x}}^{\prime})}\bm{)}_{t}\Xi_{\mu\rho}\Xi_{\nu\sigma},\\\
\rho_{\perp}&=\frac{\beta(t,{\bm{x}})}{2}\int_{-\infty}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\bm{(}{\tilde{\delta}\hat{E}^{\mu}(t,{\bm{x}}),\tilde{\delta}\hat{E}^{\nu}(t^{\prime},{\bm{x}}^{\prime})}\bm{)}_{t}\Xi_{\mu\nu},\\\
\rho_{\parallel}&=\frac{\beta(t,{\bm{x}})}{2}\int_{-\infty}^{t}\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\bm{(}{\tilde{\delta}\hat{D}^{\mu\nu}(t,{\bm{x}}),\tilde{\delta}\hat{D}^{\rho\sigma}(t^{\prime},{\bm{x}}^{\prime})}\bm{)}_{t}\Xi_{\mu\rho}\Xi_{\nu\sigma}.\end{split}$
(77)

They are the Green-Kubo formulas Green (1954); Nakano (1956); Kubo (1957) for
the seven transport coefficients — three bulk viscosities
($\zeta_{\parallel},\zeta_{\perp},\zeta_{\times}$), two shear viscosities
($\eta_{\parallel},\eta_{\perp})$, and two electric resistivities
($\rho_{\parallel},\rho_{\perp}$) in RMHD (see Sec. IV.2).

Two remarks are in order. In the previous section, we do not count
$\zeta_{\times}^{\prime}$ as an independent transport coefficient. This is
because the corresponding Green-Kubo formula in (77) respects Onsager’s
reciprocal relation Onsager (1931): $\zeta_{\times}^{\prime}=\zeta_{\times}$.
This can be shown by performing an expansion around the global equilibrium and
using the charge-conjugation and time-reversal symmetries applied to the above
Green-Kubo formula for $\zeta_{\times}^{\prime}$. We also note that the Green-
Kubo formulas (77) automatically provide a set of the simi-positivity
constraints given by Eq. (32) specified in the previous section. This stems
from the property of the Kubo-Mori-Bogoliubov inner product. Here, it is worth
emphasizing that the semi-positivity constraints (and Onsager’s reciprocal
relation) are not required but derived in the nonequilibrium statistical
operator method.

While the expressions in Eq. (77) are given in terms of the Kubo-Mori-
Bogoliubov inner product, one can derive a set of more familiar Green-Kubo
formulas in terms of the retarded Green’s functions. For this purpose, we
again expand Eq. (77) on top of the global equilibrium. By inserting the
convergence factor $\mathrm{e}^{\epsilon(t^{\prime}-t)}$ which will be
eventually turned off by taking $\epsilon\to 0$ after the whole calculation,
we generally obtain

$\begin{split}\int_{-\infty}^{t}\mathrm{d}t^{\prime}\mathrm{e}^{\epsilon(t^{\prime}-t)}\bm{(}{\hat{A}(t),\hat{B}(t^{\prime})}\bm{)}_{\mathrm{eq}}=-\frac{1}{\beta}\int\frac{\mathrm{d}\omega}{2\pi}\mathcal{P}\left(\frac{1}{\omega}\right)\frac{\partial}{\partial\omega}G_{R}^{A,B}(\omega)-\frac{\mathrm{i}}{2\beta}\lim_{\omega\to
0}\frac{\partial}{\partial\omega}G_{R}^{A,B}(\omega),\end{split}$ (78)

where $\mathcal{P}$ stands for the principal value,
$\bm{(}{\hat{A}(t),\hat{B}(t^{\prime})}\bm{)}_{\mathrm{eq}}$ is the global
equilibrium limit of the inner product
$\bm{(}{\hat{A}(t),\hat{B}(t^{\prime})}\bm{)}_{t}$, and $G_{R}^{A,B}$ is the
retarded Green’s function,
$G_{R}^{A,B}(t-t^{\prime})=-i\theta(t-t^{\prime})\langle[\hat{A}(t),\hat{B}(t^{\prime})]\rangle_{\rm
eq}$. With the help of this identity, we can replace the inner products in Eq.
(77) with the retarded Green’s functions $G_{R}^{A,B}$.

##  IV Interlude: Connection to the conventional MHD

In previous sections, we have reviewed the recent formulation of the RMHD with
only the energy-momentum conservation law and the Bianchi identity as relevant
equations of motion from the very beginning. On the other hand, the RMHD can
also be formulated in a conventional approach in which the Maxwell equation
and electric charge conservation law enter as additional dynamical equations.
In the following two sections V and VI, we will review hydrodynamics under the
strong background magnetic field that is closer to the latter formulation
though the magnetic field is non-dynamical. Thus, in this intermediate
section, we discuss the relation between the two formulations, focusing on the
anisotropic pressure and the first-order constitutive relations. The second
topic also serves as a basis for the following two sections.

###  IV.1 Anisotropic pressure

We start with the correspondence of the zeroth-order terms, especially the
anisotropic pressure, between the conventional formulation and the formulation
we discussed in Sec. II and Sec. III. In the conventional formulation, one
needs to separate the matter and EM components in the system, and connect the
two components via the Maxwell equation $\partial_{\mu}F^{\mu\nu}=J^{\nu}$ and
energy-momentum non-conservation equation $\partial_{\mu}T^{\mu\nu}_{\rm
matt}=F^{\nu\mu}J_{\mu}$, where $J^{\mu}$ is the electric current and
$T^{\mu\nu}_{\rm matt}$ is the matter energy-momentum tensor. To make a
connection with our formulation in Sec. II and Sec. III, let us decompose the
total energy density and pressure presented in Sec. II into the matter and
magnetic components as $\epsilon\sim\epsilon_{\rm
matt}+|{\bm{B}}|^{2}/2-{\bm{M}}\cdot{\bm{B}}$ and $p\sim p_{\rm
matt}+|{\bm{B}}|^{2}/2-{\bm{M}}\cdot{\bm{B}}$, respectively 777 In the above
expressions, the Lorentz scalars are given in the rest-frame expressions for
clarity. The energy density including the magnetization $\epsilon_{\rm
matt}-{\bm{M}}\cdot{\bm{B}}$ corresponds to the definition of $\epsilon$ in
Ref. Huang  _et al._ (2011) as stated there. . Then, inserting those
expressions into Eq. (19), we have

$\displaystyle T^{\mu\nu}_{(0)}$ $\displaystyle\simeq$
$\displaystyle(\epsilon_{\rm
matt}+\frac{1}{2}|{\bm{B}}|^{2}-{\bm{M}}\cdot{\bm{B}})u^{\mu}u^{\nu}$ (79)
$\displaystyle-(p_{\rm
matt}+\frac{1}{2}|{\bm{B}}|^{2}-{\bm{M}}\cdot{\bm{B}})\Xi^{\mu\nu}+(p_{\rm
matt}-\frac{1}{2}|{\bm{B}}|^{2})b^{\mu}b^{\nu}\,.$

The above expression agrees with those in Eqs. (14) and (15) of Ref. Huang
_et al._ (2011) (see also Refs. Kluitenberg and De Groot (1954b, c);
Lichnerowicz (1967); Israel (1978); Gedalin (1991) for classic works): the
terms with $|\bm{B}|^{2}$ are combined into the Maxwell energy-momentum tensor
(see below) and the rest terms give $T^{\mu\nu}_{{\rm matt}(0)}$. In this way,
the conventional “ideal MHD limit” is reproduced from the leading-order result
in the new formulation that generalizes the conventional formulation in the
following points.

  * •

The right-hand side of the energy-momentum non-conservation law,
$\partial_{\mu}T^{\mu\nu}_{\rm matt}=F^{\nu\mu}J_{\mu}$, describes the Joule
heat and Lorentz force which provide the source and/or dissipation of the
energy and momentum, respectively. Similarly, the electric current in the
Maxwell equation $\partial_{\mu}F^{\mu\nu}=J^{\nu}$ provides a source of EM
fields. This latter equation constrain the electric field as a gapped mode
excited by the source term. The new formulation does not contain such a
redundancy and can work in the strict hydrodynamic limit.

  * •

The pressures $p_{\perp}$ and $p_{\parallel}$ in Eq. (19) satisfy
$p_{\perp}-p_{\parallel}=B^{\mu}H_{\mu}(>0)$. Subtracting the magnetic
pressures $p_{B\perp}=|\bm{B}|^{2}/2$ and $p_{B\parallel}=-|\bm{B}|^{2}/2$, we
obtain the anisotropic matter pressures as $p_{{\rm matt}\parallel}=p_{\rm
matt}$ and $p_{{\rm matt}\perp}=p_{\rm matt}-{\bm{M}}\cdot{\bm{B}}$,
respectively. One thus finds that $p_{{\rm matt}\perp}-p_{{\rm
matt}\parallel}=-\bm{M}\cdot\bm{B}$ in the rest frame of the fluid. This leads
to $p_{{\rm matt}\perp}<p_{{\rm matt}\parallel}$ since the magnetic
susceptibility is usually positive.

  * •

These non-conservation equations can be combined together into the form
$\partial_{\mu}(T^{\mu\nu}_{\rm matt}+T^{\mu\nu}_{\rm Maxwell})=0$, where
$T^{\mu\nu}_{\rm Maxwell}$ is the Maxwell tensor 888 Explicitly, this means
that $F^{\nu\mu}J_{\mu}=\partial_{\mu}(F^{\mu\alpha}F^{\nu}_{\,\
\alpha}-g^{\mu\nu}F^{\alpha\beta}F_{\alpha\beta}/4)=:-\partial_{\mu}T^{\mu\nu}_{\rm
Maxwell}$. . Therefore, this equation can be reduced from the first equation
in Eq. (13) if one assumes a clear separation between the matter and
electromagnetic contributions to the energy-momentum tensor as in Eq. (79).
However, it would not be possible to separate those contributions in a
strongly coupled system where excitations are composed of mixture of matter
and electromagnetic fields. Moreover, hydrodynamic framework itself should not
care such microscopic details of the system, and the translational symmetry of
the system only tells us the conservation of the total energy-momentum. Those
facts should be respected in the formulation.

  * •

Excluding an electric field from the set of hydrodynamic variables, one does
not need to assume an “infinite electric conductivity” as in the conventional
formulation (see, e.g., Refs. Lichnerowicz (1967); Davidson (2002)). Note that
the electric conductivity is a dimensionful quantity and is, moreover, not
defined a priori in the formulation of hydrodynamics. If it implies an
infinitesimally short relaxation time, there would be also a conceptual
conflict when one tries to include (finite) dissipative effects such as a
viscosity in the derivative corrections. The formulation in Sec. II and Sec.
III is free of such a dilemma.

###  IV.2 First-order constitutive relations including Hall transports

We next discuss the correspondence between the two formulations at the first-
order in derivatives and show how the electric field and current arises in the
new formulation. We also include the charge-conjugation odd terms in
$T_{(1)}^{\mu\nu}$ and $\tilde{F}_{(1)}^{\mu\nu}$ to identify how the Hall
components could appear when a finite charge density is allowed in a finite
time scale. When charge-conjugation odd terms are allowed in the tensor
decomposition (23), we have an additional term as

$\rho^{\mu\nu\rho\sigma}=2\rho_{\perp}\left(b^{\mu}\Xi^{\nu[\rho}b^{\sigma]}-b^{\nu}\Xi^{\mu[\rho}b^{\sigma]}\right)+2\rho_{\parallel}\Xi^{\mu[\rho}\Xi^{\sigma]\nu}+2\rho_{H}\left(b^{\mu}b_{\star}^{\nu[\rho}b^{\sigma]}-b^{\nu}b_{\star}^{\mu[\rho}b^{\sigma]}\right),$
(80)

where $b^{\mu\nu}_{\star}=\varepsilon^{\mu\nu\alpha\beta}u_{\alpha}b_{\beta}$.
The charge-conjugation odd term does not create entropy in Eq. (21b). Thus,
the additional coefficient $\rho_{H}$ can be both positive and negative
values, while $\rho_{\parallel,\perp}$ should be positive semi-definite as we
have seen. The first-order correction $\tilde{F}_{(1)}^{\mu\nu}$ provides the
constitutive relation of the electric field

$\displaystyle E_{(1)}^{\mu}$ $\displaystyle=$
$\displaystyle-[\rho_{\parallel}(-b^{\mu}b^{\rho})+\rho_{\perp}\Xi^{\mu\rho}+\rho_{H}b_{\star}^{\mu\rho}]T\epsilon_{\rho\nu\alpha\beta}u^{\nu}\partial^{\alpha}(\beta
H^{\beta})\,,$ (81)

where we used identities
$b_{\star}^{\mu\lambda}b^{\nu}_{\star\lambda}=\Xi^{\mu\nu}$,
$2b_{\star}^{\mu[\alpha}b^{\beta]}=\Xi^{\mu\rho}\epsilon_{\rho\nu\alpha\beta}u^{\nu}$,
and
$2\Xi^{\mu[\alpha}b^{\beta]}=-b_{\star}^{\mu\rho}\epsilon_{\rho\nu\alpha\beta}u^{\nu}$.
This is an extension of Eq. (25) with the Hall term. On the other hand, the
in-medium Maxwell equation provides the relation between the zeroth-order
field strength tensor and the free electric current (also called conduction
current) as 999Using $H^{\mu}=-\mu_{m}^{-1}B^{\mu}$ with $\mu_{m}$ the
magnetic permeability , one can find that the free current $J_{f{(1)}}^{\mu}$
is related to the total current
$J_{(1)}^{\mu}=-(1/2)\varepsilon^{\nu\mu\alpha\beta}\partial_{\nu}(B_{\alpha}u_{\beta}-B_{\beta}u_{\alpha})$
by
$J_{f{(1)}}^{\mu}=\mu_{m}^{-1}J^{\mu}_{(1)}+\varepsilon^{\mu\nu\alpha\beta}H_{\alpha}u_{\beta}\nabla_{\nu}\ln\mu_{m}$.

$\displaystyle J^{\mu}_{f{(1)}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\varepsilon^{\nu\mu\alpha\beta}\partial_{\nu}(H_{\alpha}u_{\beta}-H_{\beta}u_{\alpha})$
(82) $\displaystyle=$
$\displaystyle-T\varepsilon^{\mu\nu\alpha\beta}u_{\nu}\nabla_{\alpha}(\beta
H_{\beta})+\varepsilon^{\mu\nu\alpha\beta}H_{\nu}\nabla_{\alpha}u_{\beta}-\varepsilon^{\mu\nu\alpha\beta}u_{\nu}H_{\alpha}(Du_{\beta}-\nabla_{\beta}\ln
T)\,.$

The first term is generated by a nonzero curl of the magnetic field. The
second term is along $u^{\mu}$ direction and thus generates a non-zero charge
density in the rest frame of the fluid,
$u_{\mu}J_{f{(1)}}^{\mu}=\varepsilon^{\mu\nu\alpha\beta}u_{\mu}H_{\nu}\nabla_{\alpha}u_{\beta}=-2\omega_{\mu}H^{\mu}$
with the vorticity vector
$\omega^{\mu}:=(1/2)\varepsilon^{\mu\nu\alpha\beta}u_{\nu}\nabla_{\alpha}u_{\beta}$.
The rest of terms are Hall-like currents that are driven by the acceleration
and the temperature gradient, which, upon using Eq. (73c), can be re-written
as

$\displaystyle-\varepsilon^{\mu\nu\alpha\beta}u_{\nu}H_{\alpha}(Du_{\beta}-\nabla_{\beta}\ln
T)$ $\displaystyle=$
$\displaystyle\frac{2T}{\epsilon+p_{\perp}}\varepsilon^{\mu\nu\alpha\beta}u_{\nu}H_{\alpha}B^{\lambda}\nabla_{[\lambda}(\beta
H_{\beta]})$ (83) $\displaystyle=$
$\displaystyle\frac{TBH}{\epsilon+p_{\perp}}\Xi^{\mu}_{\lambda}\varepsilon^{\lambda\nu\alpha\beta}u_{\nu}\nabla_{\alpha}(\beta
H_{\beta}).$

Therefore, we have

$\displaystyle-T\varepsilon^{\mu\nu\alpha\beta}u_{\nu}\nabla_{\alpha}(\beta
H_{\beta})=\frac{\epsilon+p_{\perp}}{\epsilon+p_{\perp}-BH}\Xi^{\mu}_{\lambda}J^{\lambda}_{f{(1)}}-b^{\mu}b_{\lambda}J^{\lambda}_{f{(1)}}.$
(84)

Plugging those expressions, one can rewrite the entropy production rate (21b)
as

$\displaystyle{\mathcal{R}}_{E}=-\beta
E_{(1)}^{\mu}\left[\frac{\epsilon+p_{\perp}}{\epsilon+p_{\perp}-BH}\Xi_{\mu\lambda}J^{\lambda}_{f{(1)}}-b_{\mu}b_{\lambda}J^{\lambda}_{f{(1)}}\right]+{O}(\partial^{3})\,.$
(85)

The origin of the entropy production ${\mathcal{R}}_{E}$ is identified with
the Joule heat due to the induced electric field and current.

To get a direct relation between the electric field and current, we eliminate
the magnetic fields in Eq. (81) using Eq. (82) to find 101010The structure of
Eq. (86) becomes more transparent if we write it in three-vector form in frame
$u^{\mu}=(1,\bm{0})$. The result is
$\bm{E}_{(1)}=\tilde{\rho}_{\perp}\bm{J}_{f{(1)}}+\frac{\tilde{\rho}_{\parallel}-\tilde{\rho}_{\perp}}{|\bm{B}|^{2}}(\bm{B}\cdot\bm{J}_{f{(1)}})\bm{B}+\frac{\tilde{\rho}_{H}}{|\bm{B}|}\bm{J}_{f{(1)}}\times\bm{B},$
where
$\bm{J}_{f{(1)}}=\mu_{m}^{-1}\left(\bm{J}_{(1)}-\bm{\nabla}\ln\mu_{m}\times\bm{B}\right)$
in terms of the total current $\bm{J}_{(1)}$. When order-zero effective
background charges are present, $\bm{E}_{(1)}$ may contain also terms
proportional to the gradient of the effective background charge potential.
Similar results have recently been obtained by using the method of effective
field theory Vardhan  _et al._ (2022).

$\displaystyle E^{\mu}_{(1)}$ $\displaystyle=$
$\displaystyle[-\tilde{\rho}_{\parallel}b^{\mu}b^{\nu}+\tilde{\rho}_{\perp}\Xi^{\mu\nu}+\tilde{\rho}_{H}b_{\star}^{\mu\nu}]J_{f{(1)}\nu}\,,$
(86)

where
$\tilde{\rho}_{\parallel}:=\rho_{\parallel},\tilde{\rho}_{\perp}:=\rho_{\perp}(\epsilon+p_{\perp})/(\epsilon+p_{\perp}-BH)$,
and $\tilde{\rho}_{H}:=\rho_{H}(\epsilon+p_{\perp})/(\epsilon+p_{\perp}-BH)$.
Note that, same as $\rho_{\parallel}$ and $\rho_{\perp}$,
$\tilde{\rho}_{\parallel}$ and $\tilde{\rho}_{\perp}$ are positive semi-
definite as well. It is more convenient to express the free electric current
in terms of the electric fields and we obtain

$\displaystyle J^{\mu}_{f{(1)}}$ $\displaystyle=$ $\displaystyle-2\omega\cdot
Hu^{\mu}+[-\sigma_{\parallel}b^{\mu}b^{\nu}+\sigma_{\perp}\Xi^{\mu\nu}+\sigma_{H}b_{\star}^{\mu\nu}]E_{{(1)}\nu}\,.$
(87)

The three coefficients are given by the coefficients
$\tilde{\rho}_{\perp,\parallel,H}$ as

$\displaystyle\sigma_{\parallel}=\frac{1}{\tilde{\rho}_{\parallel}},\quad\sigma_{\perp}=\frac{\tilde{\rho}_{\perp}}{\tilde{\rho}_{\perp}^{2}+\tilde{\rho}_{H}^{2}},\quad\sigma_{H}=\frac{\tilde{\rho}_{H}}{\tilde{\rho}_{\perp}^{2}+\tilde{\rho}_{H}^{2}}\,.$
(88)

The Hall term $\propto\sigma_{H}$ identically vanishes in Eq. (85) and does
not create entropy. The other coefficients $\sigma_{\parallel}$ and
$\sigma_{\perp}$ are the parallel and perpendicular components of the Ohmic
conductivity with respect to the direction of the magnetic flux. The Ohmic
conductivities should be positive semi-definite $\sigma_{\parallel,\perp}\geq
0$, while the Hall conductivity $\sigma_{H}$ can be both positive and negative
values.

Next, we discuss the Hall components in the viscous tensor. One can divide the
energy-momentum tensor into the dissipative and nondissipative components as

$\displaystyle T_{(1)}^{\mu\nu}=T_{(1){\rm
dis}}^{\mu\nu}+T_{(1)\star}^{\mu\nu}.$ (89)

In the preceding sections, we have already discussed the dissipative component
$T_{(1){\rm dis}}^{\mu\nu}$ that creates an entropy. Here, we focus on the
nondissipative component $T_{(1)\star}^{\mu\nu}$ that can be further
decomposed as

$\displaystyle
T_{(1)\star}^{\mu\nu}=\big{[}\,2\eta_{H\parallel}(-b^{\alpha}b^{(\mu})b_{\star}^{\nu)\beta}+2\eta_{H\perp}\Xi^{\alpha(\mu}b_{\star}^{\nu)\beta}\,\big{]}w_{\alpha\beta}\,,$
(90)

where $w_{\alpha\beta}=\nabla_{(a}u_{\beta)}$. This term does not create
entropy 111111 We immediately notice that
$b^{\alpha}b^{(\mu}b_{\star}^{\nu)\beta}w_{\alpha\beta}w_{\mu\nu}=2b_{\star}^{\nu\beta}(b^{\alpha}w_{\alpha\beta})(b^{\mu}w_{\mu\nu})=0$,
and then that
$\Xi^{\alpha(\mu}b_{\star}^{\nu)\beta}w_{\alpha\beta}w_{\mu\nu}=g^{\alpha(\mu}b_{\star}^{\nu)\beta}w_{\alpha\beta}w_{\mu\nu}=0$.,
and the Hall viscous coefficients $\eta_{H_{\parallel},\perp}$ can take both
positive and negative values. The presence of the antisymmetric tensor
$b_{\star}^{\mu\nu}$ implies that the direction of the stress is orthogonal to
both the flow velocity and the magnetic field.

Figure 5: Hall viscosities induced by the Lorentz force exerting on charged
fluid cells. The flow velocity and the Lorentz force are shown with green and
orange arrows, respectively. The length and thickness of the arrows denote the
magnitudes of the vectors. The direction of a magnetic field is taken along
the $z$ direction.

We briefly demonstrate the mechanisms that generate the Hall viscosities in
magnetic fields (assuming that there is a possible electric-charge density).
In Fig. 5, the direction of a magnetic field is taken along the $z$ direction,
and the flow velocity and the Lorentz force are shown with green and orange
arrows, respectively. When there is a gradient of the transverse flow along
the magnetic field $\partial_{z}u_{x}(z)>0$, the Lorentz force exerting on the
fluid volume is oriented in the $y$ direction and has a gradient along the
magnetic field (see the left panel). This configuration corresponds to the
term proportional to $\eta_{H\parallel}$. The gradient of the Lorentz force
gives rise to a shear stress in the plane orthogonal to the flow-velocity and
magnetic fields. On the other hand, when there is a flow gradient in the
transverse plane $\partial_{y}u_{x}(y)>0$, the Lorentz force is oriented in
the $y$ direction and has a gradient in the same direction (see the right
panel). Therefore, the fluid volume is stretched in this direction according
to the term proportional to $\eta_{H\perp}$. This term also induces a shear
deformation in the $x$-$y$ plane in response to an expansion/compression in
the $x$ direction, $\partial_{x}u_{x}(x)\not=0$. This effect can be understood
in a similar manner.

In the nonrelativistic theory, it has been known for some time that there are
seven viscous coefficients Hooyman  _et al._ (1954); Lifschitz and Pitaevskii
(1981) (see also Ref. De Groot and Mazur (1962)). In this paper, we have
obtained the breakdown of the seven viscous coefficients in the relativistic
extension. Including the five dissipative viscosities discussed in the
previous sections, we identified three bulk, two shear, and two Hall
composnents. The relativistic extension was also carried out earlier in Refs.
Huang  _et al._ (2011); Hernandez and Kovtun (2017); Grozdanov  _et al._
(2017); Hongo and Hattori (2021). However, those authors use different tensor
bases. For readers’ convenience, we provide an explicit comparison among those
tensor bases in Appendix D.

##  V Towards relativistic MHD from kinetic equations

In addition to hydrodynamics, kinetic theory is also often used for the
studies of many-body systems (See Refs. Lifschitz and Pitaevskii (1981); de
Groot  _et al._ (1980) for the standard textbooks of kinetic theory). It is
valid in the regime where the system permits well-defined particles (or quasi-
particles) and when the densities of these particles are low enough. More
precisely, let the typical microscopic scales be set by the interaction range
$r_{c}\sim\sqrt{\sigma}$ with $\sigma$ the scattering cross section, the
inter-particle distance $r_{d}\sim n^{-1/3}$ with $n$ the number density, and
the mean-free path $\lambda_{\rm mfp}\sim 1/(n\sigma)$. The applicability of
kinetic theory requires $\lambda_{\rm mfp}\gg r_{d}\gg r_{c}$ (dilute
condition). The scatterings of the particles drive the system to evolve
towards global thermal equilibrium, a process called kinetic thermalization.
This process is usually associated with the arising of a new, macroscopic
scale, $L$, over which the macroscopic properties of the system vary. In the
late stage of the thermalization, the scale $L$ will be clearly separated from
the microscopic scales. This gives a characteristic parameter ${\rm
Kn}=\lambda_{\rm mfp}/L$ (the Knudsen number). When ${\rm Kn}\ll 1$ or
equivalently when the derivative $\partial$ is much smaller than
$1/\lambda_{\rm mfp}$, the macroscopic properties become insensitive to the
microscopic details of the system and the hydrodynamic description is expected
to arise. Therefore, when the kinetic theory is applied to the regime where
${\rm Kn}\ll 1$, its solution is expected to be expressed by the local
hydrodynamic variables and their derivatives. This gives a systematical way to
derive the hydrodynamic constitutive relations (including expressing the
transport coefficients in terms of kinetic-theory parameters) and EOMs by
expanding the kinetic equation and the distribution functions in ${\rm Kn}$ or
$\partial$. This is the essential idea of the Chapman-Enskog method Chapman
and Cowling (1970) which we will now discuss. Another frequently utilized
approach to hydrodynamics from kinetic theory is the Grad’s method of moments
Grad (1949) which we will explore in Sec. V.2.

###  V.1 Chapman-Enskog method

Before studying RMHD using kinetic theory, in order to demonstrate the
methodology, let us consider a simpler case in which the EM fields are absent.
The starting point is the relativistic Boltzmann equation for on-shell
particles:

$\displaystyle k^{\mu}\partial_{\mu}f_{{\bm{k}}}=C[f],$ (91)

where $k^{\mu}=(k^{0}=\sqrt{{\bm{k}}^{2}+m^{2}},{\bm{k}})$ with $m$ being the
mass of the particles, $f_{\bm{k}}(x)$ is the distribution function, and all
external forces are omitted. For the sake of simplicity, we consider only one
species of particles. We assume that the collisional process conserves the
energy-momentum and also the particle numbers, so that the collision kernel
satisfies

$\displaystyle\int dKC[f]=\int dKk^{\mu}C[f]=0,$ (92)

where $dK=gd^{3}{\bm{k}}/[(2\pi)^{3}k^{0}]$ is the invariant momentum-space
measure with $g=2S+1$ the spin degenerate factor with $S$ spin of the
particles. By using Eq. (91), Eq. (92) leads to the conservation of particle-
number current $N^{\mu}$ and energy-momentum tensor $T^{\mu\nu}$,

$\displaystyle\partial_{\mu}N^{\mu}$ $\displaystyle=$ $\displaystyle 0,$ (93a)
$\displaystyle\partial_{\mu}T^{\mu\nu}$ $\displaystyle=$ $\displaystyle 0,$
(93b)

where $N^{\mu}$ and $T^{\mu\nu}$ are expressed by

$\displaystyle N^{\mu}$ $\displaystyle=$ $\displaystyle\langle
k^{\mu}\rangle=nu^{\mu}+\nu^{\mu},$ (94a) $\displaystyle T^{\mu\nu}$
$\displaystyle=$ $\displaystyle\langle k^{\mu}k^{\nu}\rangle=\epsilon
u^{\mu}u^{\nu}-\Delta^{\mu\nu}(p+\Pi)+h^{\mu}u^{\nu}+h^{\nu}u^{\mu}+\pi^{\mu\nu},$
(94b)

with $\langle\cdots\rangle=\int dK(\cdots)f_{\bm{k}}$. Here, we introduced a
timelike unit vector field $u^{\mu}$ which will be identified as the fluid
velocity. Using $u^{\mu}$, the momentum can be decomposed into two parts
(following the notations of Ref. Denicol  _et al._ (2012a)):
$k^{\mu}=E_{\bm{k}}u^{\mu}+k^{\langle\mu\rangle}$, with $E_{\bm{k}}=u\cdot k$
and $k^{\langle\mu\rangle}=\Delta^{\mu\nu}k_{\nu}$. We have also introduced
the number density $n$, the number diffusion current $\nu^{\mu}$, the energy
density $\epsilon$, the thermodynamic pressure $p$, the viscous pressure
$\Pi$, the shear viscous tensor $\pi^{\mu\nu}$, and the heat flux $h^{\mu}$,

$\displaystyle n=\langle E_{\bm{k}}\rangle,\;\nu^{\mu}=\langle
k^{\langle\mu\rangle}\rangle,\;\epsilon=\langle
E_{\bm{k}}^{2}\rangle,\;p+\Pi=-\frac{1}{3}\Delta^{\mu\nu}\langle
k_{\mu}k_{\nu}\rangle,\;\pi^{\mu\nu}=\langle
k^{\langle\mu}k^{\nu\rangle}\rangle,\;h^{\mu}=\langle
E_{\bm{k}}k^{\langle\mu\rangle}\rangle,$

where
$A^{\langle\mu\nu\rangle}=(1/2)\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}(A^{\alpha\beta}+A^{\beta\alpha})-(1/3)\Delta^{\mu\nu}\Delta_{\alpha\beta}A^{\alpha\beta}$
is the spatial traceless symmetrization of $A^{\mu\nu}$. The expression for
$p$ itself will be determined later. To specify the rest frame of the fluid,
we use the Landau-Lifshitz choice so that $T^{\mu\nu}u_{\nu}=\epsilon
u^{\mu}$. This imposes the constraint

$\displaystyle h^{\mu}=\langle E_{\bm{k}}k^{\langle\mu\rangle}\rangle=0.$ (96)

To proceed, one takes the Chapman-Enskog expansion of the distribution
function, $f_{\bm{k}}=f^{(0)}_{{\bm{k}}}+f_{{\bm{k}}}^{(1)}+\cdots$, and of
the Boltzmann equation (91) order by order in ${\rm Kn}$ or equivalently in
derivatives. At zeroth order in the derivative expansion, Eq. (91) reads

$\displaystyle 0=C[f^{(0)}].$ (97)

Its solution is called the local-equilibrium distribution
$f^{(0)}_{\bm{k}}=f_{0{\bm{k}}}$ because it saturates the local detailed
balance. We assume that $f_{0{\bm{k}}}=1/[\exp(\beta E_{\bm{k}}-\alpha)+a]$
with $a=1,0$ and $-1$ for fermions, classical Boltzmann particles, and bosons.
Here, $\beta$ and $\alpha=\beta\mu$ are the local inverse temperature and the
ratio of the chemical potential to temperature. Their values are fixed by
matching conditions Stewart (1971); de Groot  _et al._ (1980):

$\displaystyle\epsilon=\langle E_{\bm{k}}^{2}\rangle=\langle
E_{\bm{k}}^{2}\rangle_{0},\;\;n=\langle E_{\bm{k}}\rangle=\langle
E_{\bm{k}}\rangle_{0},$ (98)

where $\langle\cdots\rangle_{0}=\int dK(\cdots)f_{0{\bm{k}}}$. Substituting
$f_{0{\bm{k}}}$ into Eqs. (94a) and (94b), one obtains the zeroth-order
$N^{\mu}$ and $T^{\mu\nu}$:

$\displaystyle N_{(0)}^{\mu}$ $\displaystyle=$ $\displaystyle nu^{\mu},$ (99a)
$\displaystyle T_{(0)}^{\mu\nu}$ $\displaystyle=$ $\displaystyle\epsilon
u^{\mu}u^{\nu}-p\Delta^{\mu\nu},$ (99b)

where $p=-\frac{1}{3}\Delta^{\mu\nu}\langle k_{\mu}k_{\nu}\rangle_{0}$, and
Eqs. (93a) and (93b) become the ideal hydrodynamic equations.

At first order in derivative expansion, the collision kernel becomes a linear
integral operator acting on $f_{\bm{k}}^{(1)}$,
$C^{(1)}[f]=\hat{L}f_{\bm{k}}^{(1)}$. The Boltzmann equation (91) at $O({\rm
Kn})$ thus reads

$\displaystyle k^{\mu}\partial_{\mu}f_{0{\bm{k}}}=\hat{L}f_{\bm{k}}^{(1)}.$
(100)

Once $f_{\bm{k}}^{(1)}$ is solved out from this equation, the first-order
constitutive relations are then obtained:

$\displaystyle\nu^{\mu}=\langle
k^{\langle\mu\rangle}\rangle_{1},\;\Pi=-\frac{1}{3}\Delta^{\mu\nu}\langle
k_{\mu}k_{\nu}\rangle_{1},\;\pi^{\mu\nu}=\langle
k^{\langle\mu}k^{\nu\rangle}\rangle_{1},$ (101)

with $\langle\cdots\rangle_{1}=\int dK(\cdots)f^{(1)}_{{\bm{k}}}$. To see how
this procedure works in practice, let us consider a collision kernel in
relaxation-time approximation (RTA) 121212Note that this collision kernel does
not automatically conserve particle number and energy-momentum, Eq. (92).
However, these conservation laws are recovered once the matching conditions
(98) and Landau-Lifshitz frame constraint (96) are imposed.,

$\displaystyle C[f]=-E_{\bm{k}}\frac{f_{\bm{k}}-f_{0{\bm{k}}}}{\tau_{R}}.$
(102)

Writing $f^{(1)}_{\bm{k}}$ as
$f^{(1)}_{\bm{k}}=f_{0{\bm{k}}}\tilde{f}_{0{\bm{k}}}\phi_{\bm{k}}$ with
$\phi_{\bm{k}}$ at $O({\rm Kn})$ and $\tilde{f}_{\bm{k}}=1-af_{\bm{k}}$, one
obtains
${\hat{L}}f^{(1)}_{\bm{k}}=-\tau_{R}^{-1}E_{\bm{k}}f_{0{\bm{k}}}\tilde{f}_{0{\bm{k}}}\phi_{\bm{k}}$.
From Eq. (100), one finds

$\displaystyle\phi_{\bm{k}}=\frac{\beta\tau_{R}}{E_{\bm{k}}}\bigg{[}\left(\frac{1}{3}\Delta_{\mu\nu}k^{\mu}k^{\nu}+E^{2}_{\bm{k}}\frac{\partial
p}{\partial\epsilon}\Big{|}_{n}+E_{\bm{k}}\frac{\partial p}{\partial
n}\Big{|}_{\epsilon}\right)\theta+T\left(\frac{nE_{\bm{k}}}{\epsilon+p}-1\right)k^{\mu}\nabla_{\mu}\alpha+k^{\langle\mu}k^{\nu\rangle}\xi_{\mu\nu}\bigg{]},$

where we defined $\xi_{\mu\nu}:=w_{\langle\mu\nu\rangle}$ which is the shear
tensor 131313The symbol $\xi_{\mu\nu}$ for shear tensor is used only in this
section in order to simplify the equations. In other sections, we simply use
$w_{\langle\mu\nu\rangle}$ or $\nabla_{\langle\mu}u_{\nu\rangle}$ to denote
the shear tensor., and the terms are organized in such a way that they are
mutually orthogonal under integral $\int dK(\cdots)F(E_{\bm{k}})$ with
$F(E_{\bm{k}})$ an arbitrary converging function of $E_{\bm{k}}$. Substituting
$\phi_{\bm{k}}$ into Eq. (101), one finds

$\displaystyle\nu^{\mu}$ $\displaystyle=$
$\displaystyle\kappa\,\nabla^{\mu}\alpha,\quad{\rm
with}\quad\kappa=\tau_{R}T\left(\frac{2}{3}\frac{\epsilon-3p}{m^{2}}-\frac{n^{2}}{\epsilon+p}+\frac{m^{2}}{3}\left\langle\frac{1}{E_{\bm{k}}^{2}}\right\rangle_{0}\right),\,$
(104a) $\displaystyle\Pi$ $\displaystyle=$
$\displaystyle-\zeta\,\theta,\quad{\rm
with}\quad\zeta=\tau_{R}\left[\left(\frac{1}{3}-c_{s}^{2}\right)(\epsilon+p)-\frac{2}{9}(\epsilon-3p)-\frac{m^{4}}{9}\left\langle\frac{1}{E_{\bm{k}}^{2}}\right\rangle_{0}\right],\,$
(104b) $\displaystyle\pi^{\mu\nu}$ $\displaystyle=$ $\displaystyle
2\eta\,\xi^{\mu\nu},\quad{\rm
with}\quad\eta=\tau_{R}\left[\frac{4}{5}p+\frac{1}{15}(\epsilon-3p)-\frac{m^{4}}{15}\left\langle\frac{1}{E_{\bm{k}}^{2}}\right\rangle_{0}\right],\,$
(104c)

where $c_{s}^{2}=\partial p/\partial\epsilon|_{s/n}=\partial
p/\partial\epsilon|_{n}+(\epsilon+p)^{-1}n\,\partial p/\partial n|_{\epsilon}$
is the sound velocity squared, $\kappa,\zeta,\eta$ are conductivity of number
diffusion current, bulk viscosity, and shear viscosity. Note that when $m=0$,
the bulk viscosity vanishes because $\epsilon=3p$ and $c_{s}^{2}=1/3$ in
massless limit. Thus, at $O({\rm Kn})$, we recover the relativistic Navier-
Stokes hydrodynamics.

With the above preparation, let us now consider RMHD in the Chapman-Enskog
method. In this case, the Boltzmann equation contains the EM force term:

$\displaystyle
k^{\mu}\partial_{\mu}f_{{\bm{k}}}+qF^{\mu\nu}k_{\nu}\partial_{k^{\mu}}f_{\bm{k}}=C[f],$
(105)

where $q$ is the charge of the particles (we consider only one species of
particles and assume $q>0$). We assume that the particles are under binary
elastic collisions and the colliding processes are not interfered by the
external EM field $F^{\mu\nu}$. Hence the collision kernel preserves the
particle number (and the electric charges as well) and energy-momentum, i.e.,
Eq. (92) still holds which implies

$\displaystyle\partial_{\mu}N^{\mu}$ $\displaystyle=$
$\displaystyle\partial_{\mu}\int dKk^{\mu}f_{\bm{k}}=0,$ (106)
$\displaystyle\partial_{\mu}T^{\mu\nu}$ $\displaystyle=$
$\displaystyle\partial_{\mu}\int
dKk^{\mu}k^{\lambda}f_{\bm{k}}=qF^{\nu\lambda}N_{\lambda}=F^{\nu\lambda}J_{\lambda},$
(107)

where $J^{\mu}=qN^{\mu}$ is the charge current. We emphasize that the EM
fields appear in our kinetic approach as external fields which is very
different from Sec. II and Sec. III in which the dynamics of the EM fields
plays important role. This leads to different leading-order constitutive
relations, but they share the same structures of the dissipative constitutive
relations for the fluid (namely, the viscous and conductivity tensors).

Let us clarify a significant distinction between how the EM fields are treated
in this section and how they are treated in Sec. II and Sec. III. In Sec. II
and Sec. III, the EM fields are dynamical and the energy-momentum tensor
includes contributions from the EM fields as well, resulting in a conserved
energy-momentum tensor. The RMHD derived in Sec. II and Sec. III is therefore
a strict hydrodynamic theory. This is not the case with this section, where
the EM fields are external fields (as clearly seen in Eq. (105)), hence
$N^{\mu}$ and $T^{\mu\nu}$ obtained from the kinetic theory only include
contributions from the matter. As shown in the right-hand side of Eq. (107),
matter exchanges energy and momentum with EM fields, causing the $T^{\mu\nu}$
to be non-strictly conserved. For this reason, the pressure and energy density
include only the matter contributions.

The presence of the EM fields introduces new scales into the kinetic equation
(and also the hydrodynamic equations). So let us clarify the range of scales
that we will focus in this case. The typical microscopic scales are still set
by the interaction range $r_{c}$, inter-particle distance $r_{d}$, and the
mean-free path $\lambda_{\rm mfp}$. In terms of ${\rm Kn}$ or the gradient
$\partial$ of the conserved densities, the electric field in the rest frame of
the fluid, $E^{\mu}=F^{\mu\nu}u_{\nu}$, is considered as $O(\partial)$, while
the magnetic field in the rest frame of the fluid,
$B^{\mu}=\varepsilon^{\mu\nu\alpha\beta}u_{\nu}F_{\alpha\beta}/2$, can be an
$O(1)$ quantity, as we already emphasized in Sec. II. This amounts to the fact
that the electric field is screened by the gradient of charge distribution in
a plasma, while magnetic field is not. The magnetic field makes the motion of
charged particles curvilinear with radius $R_{L}=k_{\perp}/(qB)$ (the Larmor
radius) with $k_{\perp}$ being the momentum of the particles transverse to the
magnetic field. For hot relativistic plasma, $k_{\perp}\sim T$; so
$R_{T}=T/(qB)$ sets the magnetic cyclotron scale. We will always assume that
the thermal wavelength $\beta=1/T$ of the particles and the interaction range
$r_{c}$ are much smaller than $R_{T}$, $\beta\ll R_{T}$ and $r_{c}\ll R_{T}$.
The first inequality means that the Landau quantization effect is not
significant and we can safely use a classical treatment. The second inequality
means that in the collision process the magnetic field can be neglected so we
can use a magnetic-field independent collision operator C (but the
distribution function $f_{\bm{k}}$ can certainly depend on $B^{\mu}$). The
situation with an even stronger magnetic field, $R_{T}\ll\beta$, that can make
the Landau quantization effect significant, will be discussed in Sec. VI.

We will always use the Landau-Lifshitz frame for the fluid and use the
matching conditions for $\epsilon$ and $n$. The zeroth-order distribution
function is still chosen as $f^{(0)}_{\bm{k}}=f_{0{\bm{k}}}=1/[\exp(\beta
E_{\bm{k}}-\alpha)+a]$ with $a=1,0$ and $-1$ for fermions, classical Boltzmann
particles, and bosons. Choosing such an equilibrium distribution amounts to
assume that the magnetization pressure $-MB$ is much smaller than
thermodynamic pressure $p$ and thus is omitted. This also implies that the
magnetization current is vanishing and we thus do not distinguish the free
electric current and the total electric current in this section. Therefore,
the ideal-fluid constitutive relations are still given by Eqs. (99a) and
(99b). To obtain $f^{(1)}_{\bm{k}}$, let us again take the RTA for the
collision kernel,

$\displaystyle C[f]=-E_{\bm{k}}\frac{f_{\bm{k}}-f_{0{\bm{k}}}}{\tau_{R}}.$
(108)

First, we consider a situation in which both $E^{\mu}$ and $B^{\mu}$ in Eq.
(105) are of order $O(\partial)$, namely, when both the electric and magnetic
fields are very weak. This is not in the MHD region but let us make a case
study first. Writing $f^{(1)}_{\bm{k}}$ as
$f^{(1)}_{\bm{k}}=f_{0{\bm{k}}}\tilde{f}_{0{\bm{k}}}\phi_{\bm{k}}$ with
$\phi_{\bm{k}}$ at $O({\partial})$ and $\tilde{f}_{\bm{k}}=1-af_{\bm{k}}$, one
obtains from Eq. (105),

$\displaystyle\phi_{\bm{k}}$ $\displaystyle=$
$\displaystyle-\frac{\tau_{R}}{E_{\bm{k}}}(f_{0{\bm{k}}}\tilde{f}_{0{\bm{k}}})^{-1}\left(k^{\mu}\partial_{\mu}f_{0{\bm{k}}}+qF^{\mu\nu}k_{\nu}\partial_{k^{\mu}}f_{0{\bm{k}}}\right)$
(109) $\displaystyle=$
$\displaystyle\frac{\beta\tau_{R}}{E_{\bm{k}}}\bigg{[}T\left(\frac{nE_{\bm{k}}}{\epsilon+p}-1\right)k^{\mu}\left(\nabla_{\mu}\alpha+\beta
qE_{\mu}\right)$
$\displaystyle+\left(\frac{1}{3}\Delta_{\mu\nu}k^{\mu}k^{\nu}+E^{2}_{\bm{k}}\frac{\partial
p}{\partial\epsilon}\Big{|}_{n}+E_{\bm{k}}\frac{\partial p}{\partial
n}\Big{|}_{\epsilon}\right)\theta+k_{\perp}^{\langle\mu}k_{\perp}^{\nu\rangle}\xi_{\mu\nu}\bigg{]}.$

Compared with Eq. (V.1), the only new term is the electric-field term which is
always accompanied with $\nabla_{\mu}\alpha$ and the magnetic field drops out.
To obtain $\phi_{\bm{k}}$, we have used the EOMs at $O(1)$ order (ideal
hydrodynamic equations):

$\displaystyle Du^{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon+p}\left(\nabla^{\mu}p+nqE^{\mu}\right),$ (110)
$\displaystyle D\alpha$ $\displaystyle=$ $\displaystyle-\beta\frac{\partial
p}{\partial n}\Big{|}_{\epsilon}\theta,$ (111) $\displaystyle D\beta$
$\displaystyle=$ $\displaystyle\beta\frac{\partial
p}{\partial\epsilon}\Big{|}_{n}\theta.$ (112)

Substituting $f^{(1)}_{\bm{k}}$ into Eqs. (101), we find that the $\Pi$ and
$\pi^{\mu\nu}$ are still given by Eqs. (104b) and (104c) but $\nu^{\mu}$ is
given by $\nu^{\mu}=\kappa\left(\nabla^{\mu}\alpha+\beta qE^{\mu}\right)$ with
the number diffusion constant $\kappa$ given in Eq. (104a). The charge
diffusion current $j^{\mu}=q\langle k^{\langle\mu\rangle}\rangle_{1}$ is

$\displaystyle j^{\mu}=q\nu^{\mu}=q\kappa\left(\nabla^{\mu}\alpha+\beta
qE^{\mu}\right).$ (113)

This relation gives that the electric conductivity (the coefficient in front
of $E^{\mu}$) is determined by number diffusion conductivity, $\sigma=\beta
q^{2}\kappa$, a relation representing the Wiedemann–Franz law.

The above situation with $B^{\mu}\sim O(\partial)$ is not in the MHD region.
Let us now consider the MHD region in which $E^{\mu}$ is at $O(\partial)$ but
$B^{\mu}$ is at $O(1)$. The Boltzmann equation at $O(\partial)$ reads

$\displaystyle
k^{\mu}\partial_{\mu}f_{0{\bm{k}}}-qE^{\nu}k_{\nu}u^{\mu}\partial_{k^{\mu}}f_{0{\bm{k}}}=-qBb_{\star}^{\mu\nu}k_{\nu}\partial_{k^{\mu}}f^{(1)}_{\bm{k}}-\frac{E_{\bm{k}}}{\tau_{R}}f^{(1)}_{\bm{k}},$
(114)

where $B=\sqrt{-B^{\mu}B_{\mu}}$ is the strength of the magnetic field and
$b_{\star}^{\mu\nu}=\varepsilon^{\mu\nu\alpha\beta}u_{\alpha}b_{\beta}$ is the
cross projector. The solution for $f^{(1)}_{\bm{k}}$ is
$f^{(1)}_{\bm{k}}=f_{0{\bm{k}}}\tilde{f}_{0{\bm{k}}}\phi_{\bm{k}}$ with

$\displaystyle\phi_{\bm{k}}$ $\displaystyle=$
$\displaystyle\phi_{\parallel}+\frac{\tau_{R}}{E_{\bm{k}}}\frac{\beta}{1+\xi_{\bm{k}}^{2}}\left(V^{\mu}+\xi_{\bm{k}}b_{\star}^{\mu\nu}V_{\nu}\right)k^{\perp}_{\mu}+\frac{\tau_{R}}{E_{\bm{k}}}\frac{\beta}{1+4\xi_{\bm{k}}^{2}}\left(w_{\perp}^{\mu\nu}+2\xi_{\bm{k}}b_{\star}^{\mu\lambda}w^{\nu}_{\perp\lambda}\right)k^{\perp}_{\\{\mu}k^{\perp}_{\nu\\}},$
(115) $\displaystyle\phi_{\parallel}$ $\displaystyle=$
$\displaystyle\frac{\beta\tau_{R}}{E_{\bm{k}}}\bigg{[}\left(\frac{1}{2}\Xi_{\alpha\beta}k^{\alpha}k^{\beta}+E^{2}_{\bm{k}}\frac{\partial
p}{\partial\epsilon}\Big{|}_{n}+E_{\bm{k}}\frac{\partial p}{\partial
n}\Big{|}_{\epsilon}\right)\theta_{\perp}+\left(-k_{\parallel}^{2}+E^{2}_{\bm{k}}\frac{\partial
p}{\partial\epsilon}\Big{|}_{n}+E_{\bm{k}}\frac{\partial p}{\partial
n}\Big{|}_{\epsilon}\right)\theta_{\parallel}$ (116) $\displaystyle-
Tk_{\parallel}\left(\frac{nE_{\bm{k}}}{\epsilon+p}-1\right)b^{\mu}\left(\nabla_{\mu}\alpha+\beta
qE_{\mu}\right)\bigg{]},$ $\displaystyle V^{\mu}$ $\displaystyle=$
$\displaystyle\left[T\left(\frac{nE_{\bm{k}}}{\epsilon+p}-1\right)\Xi^{\mu\lambda}\left(\nabla_{\lambda}\alpha+\beta
qE_{\lambda}\right)-2k_{\parallel}\Xi^{\mu\lambda}w_{\lambda\rho}b^{\rho}+\frac{qB}{\epsilon+p}E_{\bm{k}}b_{\star}^{\mu\lambda}\nu_{\lambda}\right],$
(117)

where $\xi_{\bm{k}}:=qB\tau_{R}/E_{\bm{k}}$ and we have introduced

$\displaystyle\begin{split}&k_{\parallel}=b\cdot
k,\;\;\Xi^{\mu\nu}=\Delta^{\mu\nu}+b^{\mu}b^{\nu},\;\;k_{\perp}^{\mu}=\Xi^{\mu\nu}k_{\nu},\;\;k_{\perp}^{\\{\mu}k_{\perp}^{\nu\\}}=k_{\perp}^{\mu}k_{\perp}^{\nu}-\frac{1}{2}\Xi^{\mu\nu}\Xi_{\alpha\beta}k^{\alpha}k^{\beta},&\;\\\