Revisiting relativistic magnetohydrodynamics from quantum electrodynamics

We provide a statistical mechanical derivation of relativistic magnetohydrodynamics on the basis of the $(3+1)$-dimensional quantum electrodynamics; the system endowed with the magnetic one-form symmetry. The conservation laws and the constitutive relations are presented in a manifestly covariant way with respect to the general coordinate transformation. The method of the local Gibbs ensemble (or nonequilibrium statistical operator) combined with the path-integral formula for the thermodynamic functional enables us to obtain an exact form of the constitutive relations. Applying the derivative expansion to the exact formula, we derive the first-order constitutive relations for the relativistic magnetohydrodynamics. The result for the QED plasma preserving the parity and charge-conjugation symmetries is equipped with two electrical resistivities and five (three bulk and two shear) viscosities. We also show that those transport coefficients satisfy the Onsager's reciprocal relation and a set of inequalities, indicating the semi-positivity of the entropy production rate consistent with the local second law of thermodynamics.


Introduction
Establishment of quantum electrodynamics (QED) is a milestone at the dawn of diverse developments in quantum field theories. The success of the covariant perturbative renormalization a la Tomonaga, Schwinger, Feynmann, and Dyson [1][2][3][4][5][6][7][8] (together with later developments in the renormalization-group method [9][10][11][12][13][14]) has deepened our fundamental theoretical backgrounds in modern physics. The further success of the perturbation theory for the electron (lepton) anomalous magnetic moment [15] has reached the state-of-art technology as the most accurate theoretical calculation presented by humankind [16][17][18][19]. Besides, QED has played a role of the prototype of quantum gauge theories: generalization to non-Abelian theories has led us to quantum chromodynamics (QCD) and electroweak theory-other basic building blocks of the standard model-that contain richer nonperturbative dynamics induced by, e.g., the asymptotic freedom and the Higgs mechanism.
In applications of QED to many-body systems, real-time dynamics of the plasma has been one of the central issues and still occupies one theoretical hot place. Such plasmas are composed of the electron and proton (and/or positron) in astrophysics [20,21] and high-intensity laser physics [22,23] and of the electron and hole in table-top materials in condensed matter physics. In the weak-coupling regime, real-time dynamics of plasmas may be described with a quasi-particle approximation, leading to the well-established kinetic theory, or the Boltzmann equation [24]. However, this is not always the case because plasmas in, e.g., a certain kind of materials has an effectively large electromagnetic coupling, which invalidates the weak-coupling approximation. A general formalism for describing the strong-coupling plasma has not been established thus far (See, e.g., Refs. [25][26][27] and references therein for reviews on strong-coupling phenomena beyond the quasi-particle approximations and condensed matter applications of the holographic principle). Therefore, it is temping to develop the hydrodynamic framework that works even in the strong-coupling regime when we focus on the low-energy (long-time/long-length) dynamics [28]. The presence of a dynamical electromagnetic field modifies the usual hydrodynamic equation to that of magnetohydrodynamics as the low-energy effective theory.
The conventional formulation of magnetohydrodynamics is simply to combine two famous equations, i.e., the Navier-Stokes equation for a fluid and the Maxwell equation for an electromagnetic field (See, e.g., Refs. [29][30][31][32] for representative textbooks). In the relativistic formulation, they are given as [33,34] ∂ µ T µν matt = F νµ j µ , (1.1a) where T µν matt and F µν are the energy-momentum tensor of the matter (fluid) and the electromagnetic field strength tensor. The right-hand sides of those equations are the Jouleheat/Lorentz-force term and the electrical current j ν which provide sources of the energy/momentum and the electric field, respectively. They are "non-conservation equations" and contain a gapped mode that dissipates in time, i.e., the electric field damped out by the Debye screening effect. Therefore, the hydrodynamic variables (or, conserved charge densities) followed from global symmetries of the system, are not appropriately identified in such a conventional formulation. An electric field should be, instead, induced by the dynamics of hydrodynamic variables, and thus, obey a constitutive relation. Another aspect of this issue is that energy-momentum tensors of the matter T µν matt and the electromagnetism are introduced separately 1 . However, what the translational symmetry of the system tells us is the conservation of the total energy and momentum composed of a mixture of matter and electromagnetic fields, instead of the matter component alone. Those fundamental issues have been overlooked in the conventional formulation 2 .
Recent progress in reformulation of magnetohydrodynamics has clarified those theoretical issues [35][36][37][38][39]. The vital point is to recognize the Bianchi identity as the conservation 1 Namely, eliminating the electric current j µ in Eqs. (1.1a) and (1.1b) and using the Bianchi identity, one can render them in the conservation equation: ∂µ(T µν matt + T µν Maxwell ) = 0, where the Maxwell stress tensor is given by T µν Maxwell = F µα F ν α − η µν F αβ F αβ /4 with the Minkowski metric ηµν = diag(−1, +1, +1, +1). Clearly, one finds a separation of the electromagnetic part from the matter part T µν matt . 2 Nevertheless, the presence of those theoretical issues do not necessarily means that the conventional approach fails to describe magnetohydrodynamic behaviors. Eliminating the electric field, one can indeed extract magnetohydrodynamic equations as long as the separation between matter and electromagnetic parts of the energy-momentum tensor works as a good approximation. law that describes time-evolution of a magnetic flux density. Moreover, motivated by the theoretical derivation of nonlinear relativistic dissipative magnetohydrodynamics in the exact form that systematically generates the derivative expansions of constitutive relations on an order-by-order basis. This development builds a solid foundation for the preceding phenomenological derivations in the literature.
The organization of the paper is as follows: In Sec. 2, we briefly review the conservation laws for the QED plasma from its global symmetries equipped with background gauge fields. In Sec. 3, we present a basic strategy of the local Gibbs ensemble method to derive the magnetohydrodynamic equation. Section 4 is devoted to the exact results both for the local equilibrium and off-equilibrium parts. In Sec. 5, we perform a derivative expansion and obtain the first-order constitutive relations supplemented with the Green-Kubo formulas for transport coefficients. The Onsager's reciprocal relation and the inequalities are also explicitly shown there (The comparison of those results to the previous works is given in Appendix A.). In Sec. 6, we summarize the results and discuss prospects.

Symmetries of QED and conservation laws
We consider (3 + 1)-dimensional QED plasmas such as an electron-positron plasma. Its microscopic dynamics is governed by the following Lagrangian: where ψ denotes a charged fermionic field with a mass m and electric charge q (q is negative for a negatively charged fermion). For simplicity, we assume that ψ represents a single component Dirac field (with the gamma matrices γ µ ), but generalization to, e.g., multicomponent cases with or without charged scalar fields is straightforward. We also introduced the U(1) gauge field A µ with the field strengh tensor F µν ≡ ∂ µ A ν − ∂ ν A µ . The QED coupling constant e is included in the definition of charge q, and η µν is the Minkowski metric with the mostly-plus convention, η µν = diag(−1, +1, +1, +1). For brevity, we will use a collective notation for dynamical variables (fermion and photon fields) as ϕ = {ψ,ψ, A µ }. The goal of this paper is to derive the hydrodynamic equations on the basis of the QED Lagrangian (2.1). Since a set of hydrodynamic equations consists of conservation laws, the first step is to specify the corresponding continuous global symmetry of QED. One can immediately identify the Poincaré invariance, i.e., the spacetime translation symmetry and Lorentz symmetry, as such symmetries resulting in the conservation law for the energymomentum and total angular momentum. As shown shortly, only the energy-momentum conservation serves as a dynamical equation, whereas the angular-momentum conservation turns out to be a constraint equation.
Besides, the QED Lagrangian (2.1) enjoys the so-called magnetic one-form symmetry (often called the generalized global symmetry). While our original dynamical variables ϕ are inert under the magnetic one-form symmetry, it acts on a dual gauge field A µ defined by ∂ µ A ν − ∂ ν A µ ≡ 1 2 µνρσ F µν with a totally anti-symmetric tensor µνρσ ( 0123 = 1). The magnetic one-form symmetry generates a shift of the dual gauge field A µ → A µ + θ µ , and the corresponding conservation law is found to be ∂ µ F µν = 0. (2.2) This is nothing but the Bianchi identity. This equation means that one can regard the magnetic flux as a conserved quantity associated with the magnetic one-form symmetry. Notice that we distinguish the magnetic flux from a magnetic field. The latter will be introduced as a conjugate variable in the subsequent section. In contrast, the QED Lagrangian (2.1) is not invariant under the electric one-form symmetry, which generates a shift of the original gauge field A µ . The explicit breaking of the electric one-form symmetry is owing to the presence of a charged matter, so that the Maxwell equation with the electric current (1.1b) describes non-conservation of the electric flux 5 . One thus only needs to take into account dynamics of the energy-momentum density and magnetic flux density in magnetohydrodynamics.
Here is another remark: Despite its conserving property, the electric charge density does not serve as a hydrodynamic (or gapless) variable in the presence of dynamical electromagnetic fields. Thus, the U(1) gauge symmetry does not yield a gapless hydrodynamic mode due to its non-local nature. In other words, a finite electric charge density induces a gapped mode, i.e., an electric flux density. This clearly causes a conflict with the absence of an electric flux in equilibrium states without spontaneous symmetry breaking such as the charge density wave state. The local equilibrium is realized only after they are damped out 6 , so that a net electric charge density does not appear in magnetohydrodynamics. Both an electric current (including an electric charge density) and electric flux are thus induced by dynamics of hydrodynamic variables as already mentioned below Eq. (1.1).
As widely known, global symmetries of the system can be tied to conservation laws in a systematic and useful way by introducing the corresponding background fields 7 . For this purpose, we introduce the background fields associated with the Poincaré and magnetic one-form symmetries, which results in the following gauged action: (2.3) We summarize generalities of the gauged action (2.3) as follows. To gauge the spacetime symmetries, we introduced a vierbein field e a µ and its inverse e µ a . Throughout this paper, we use the Greek letters (µ, ν, . . . = 0, 1, 2, 3) for the curved spacetime indices, and the Latin letters (a, b, . . . = 0, 1, 2, 3) for the local Lorentz indices. The vierbein specifies the local relation between the curved metric g µν (x) and the Minkowski metric η ab = diag(−1, +1, +1, +1) as g µν (x) = e a µ (x)e b ν (x)η ab , so that its determinant gives a correct measure with g ≡ det g µν < 0. 5 The free Maxwell theory in the absence of a charged matter also enjoys the electric one-form symmetry, and the electric flux is conserved. 6 Intuitively, one may not put like-sign charges in a finite distance, which costs a finite Coulomb energy.
Therefore, the QED plasma in an infinite volume system cannot reach a thermal (or even steady) state with a finite electric charge density coupled to the dynamical U(1) gauge field. 7 Besides, it will serve as a solid basis to discuss transport phenomena taking place in the local thermal equilibrium as discussed in Sec. 4.1.
We also defined the covariant derivative of the Dirac field as where we introduced a spin connection ω ab µ upon a requirement of the covariance under the local Lorentz transformation. Recalling the transformation property of the Dirac field, one can identify its representation matrix Σ ab as Σ ab = i[γ a , γ b ]/4. Within the vierbein postulate, D µ e a ν = ∇ µ e a ν + ω a µ b e b ν = 0, and the torsinonless assumption, one can express the spin connection by the vierbein as where we introduced the Ricci rotation coefficient C µνρ ≡ e c µ (∂ ν e ρc − ∂ ρ e νc ). We also introduced a background two-form field b µν minimally coupled to the magnetic current F µν defined by Here, we used a normalized totally anti-symmetric tensor ε µνρσ = µνρσ / √ −g in the curved spacetime. Note that ε µνρσ (instead of µνρσ ) behaves as a tensor under the general coordinate transformation. In the following, we will also use a collective notation to express a set of background fields as j ≡ {e a µ , b µν }. Now, based on the gauged action (2.3), let us derive the covariant conservation laws in the presence of background fields. The crucial point here is that the gauged action (2.3) remains invariant (δ χ S QED = 0) under the local transformation acting on both the dynamical and background fields: Here, we introduced a set of infinitesimal local transformation parameters χ ≡ {ξ µ , α a b , θ µ } for a general coordinate (ξ µ ), local Lorentz (α ab = −α ba ), and magnetic one-form gauge transformation (θ µ ), respectively. On the other hand, (re)introducing a set of gauge currents in terms of functional derivatives of the action (2.3) as we obtain another expression for the variation of the action (for notational simplicity, we hereafter express the magnetic current as J µν ≡ F µν ): To reach the second line, we performed an integration by parts and used the equation of motion δS QED /δϕ = 0. Since the action is invariant (δ χ S QED = 0) under the transformation with an arbitrary χ = {ξ µ , α ab , θ ν }, we find conservation laws of the energy-momentum tensor and magnetic flux together with the following constraint on the energy-momentum tensor: The source term J αβ H ναβ appearing on the right-hand side of the energy-momentum equation (2.10) is an analogue of the Lorentz-force term induced by the background electromagnetic field. We have thus identified the conservation laws (2.10) associated with global symmetries of the QED plasma. Those two equations govern the low-energy dynamics of the conserved charge densities, that is, the energy-momentum density and magnetic flux density. In the following sections, we will derive hydrodynamic equations based on operator versions of the conservation laws (2.10), which finally constitute magnetohydrodynamics.

Local Gibbs ensemble method for magnetohydrodynamics
In this section, we present our setup and oveview of the local Gibbs ensemble method in deriving relativistic magnetohydrodynamics. After introducing a crucial assumption on the initial density operator, we explain the difficulty of the problem and how the use of the local Gibbs ensemble method overcomes that.

Local Gibbs ensemble
Based on the preparation in the last section, we here describe the basic framework of the local Gibbs ensemble method [63][64][65] applied to our formulation of magnetohydrodynamics. Our starting point is operator versions of the conservation laws in Eq. (2.10): where H ναβ remains a classical external field defined in Eq. (2.10). While these equations provide operator identities, which serves as a basis of hydrodynamics, we still need to specify a state (or density operator) to derive hydrodynamic equations. Therefore, we, here, make a crucial assumption that the system is in a local thermal equilibrium at the initial time t 0 . One can describe a local thermal equilibrium state with the local Gibbs distribution where we defined the entropy operatorŜ[t; λ t ] aŝ Namely, we assume that the initial density operator take the formρ 0 =ρ LG [t 0 ; λ t 0 ]. In the following, we express an expectation value of an arbitrary operatorÔ over the local Gibbs The local Gibbs distribution, a generalization of the Gibbs distribution in a global equilibrium, fixes average values of the conserved charge densities with the help of conjugate Lagrange multipliers, which we collectively express as λ t ≡ {β µ (t, x), H ν (t, x)}. Since the QED plasma supports the energy-momentum and magnetic flux as conserved charges, we introduced corresponding parameters β µ and H ν , which are identified as the fluid fourvelocity u µ and magnetic field H µ multiplied by the local inverse temperature as β ν = βu ν with u ν u ν = −1 and H µ = βH µ . In this paper, we call the conjugate parameter H ν a "reduced magnetic field" after the presence of the inverse temperature prefactor. Since the zeroth component of the magnetic flux density identically vanishes (recallĴ 00 = 0), the reduced magnetic field has only the spatial components; in other words, H 0 = 0. Apart from a first-order phase transition, these Lagrange multipliers have one-to-one correspondences to the average conserved charge densities. Note that the first argument t in Eqs. (3.2)-(3.3) corresponds to the operator argument, whereas the suffix t for the parameter λ specifies its configuration at time t. We also used the covariant description by introducing the spacelike hypersurface Σ t with its normal vector dΣ tµ proportional to a 3-dimensional volume element on the hypersurface (See, e.g., Refs. [63,64] and explanation around Eq. (4.2) in detail).
To normalize the local Gibbs distribution as Trρ LG [t; λ t ] = 1, we introduced This functional, defined in the local thermal equilibrium state, serves as a generalization of the thermodynamic potential for the global equilibrium state and is called the Massieu-Planck functional. As in a global thermal equilibrium case, taking a variation with respect to the parameter λ, one can obtain the local Gibbs average of the conserved charge densities: .
Thus, the Massieu-Planck functional keeps full information on the local Gibbs average of conserved charge densities. We shall again introduce a collective notation for the average charge densities as LG t } 8 . Besides, one finds that the local Gibbs average of the entropy functional operator, S[c t ] ≡ Tr ρ LG [t; λ t ]Ŝ[t; λ t ] , can be identified as the Legendre transformation from the parameter λ t to the average charge densities c t : Accordingly, we obtain a direct relation between the value of the parameter λ t and the average charge density-T 0 , (3.8) whose global equilibrium limit matches with the definition of thermodynamic parameters following from the thermodynamic entropy. This is the reason why the functional operator S[t; λ t ] is referred to as the entropy operator.

Setup of the problem
Let us now set up the problem. We first define an expectation value of a Heisenberg operator O(t) at an arbitrary time t: where we employed the Heisenberg picture, and thus, the density operatorρ 0 does not evolve in time. The above assumption on the initial local Gibbs state will be explicit if one writes Ô (t) = Ô (t) LG t 0 . Then, taking the expectation values of the operator identities (3.1), we obtain conservation laws for the averaged quantities at an arbitrary time t (≥ t 0 ): The goal of this paper is to show that these equations finally result in a set of magnetohydrodynamic equations. To be more specific, we will express the spatial components of currents as functionals of the expectation values of conserved charge densities, so that the averaged conservation laws (3.1) form a closed set of equations. Those relations are called the constitutive relations: where the right-hand side represents functional dependences of the average currents on the average charge densities. Thus, one of our main tasks is to derive the constitutive relations from the field-theoretical framework. We will explain that the derivation of constitutive relations is far from trivial even if we put a crucial assumption on the local thermal equilibrium for the initial state. Before presenting challenging points of the problem, we shall check the numbers of the dynamical equations and degrees of freedom. For the averaged conservation laws (3.10) to form a closed set of hydrodynamic equations, the number of independent equations should match that of the conserved charges { T 0 ν , Ĵ 0ν }, namely seven (recall thatĴ µν is anti-symmetric so thatĴ 00 = 0). One, however, finds that the temporal component of the second equation, ∇ µ Ĵ µ0 = 0, is a nondynamical constraint equation, implying an apparent mismatch between the numbers of dynamical equations and independent variables. However, it is crucial to notice that we have a trivial identity ∇ ν (∇ µ Ĵ µν ) ∝ R µν Ĵ µν = 0 (we used the fact that Ricci tensor R µν is symmetric: R µν = R νµ ). Therefore, both the numbers of independent equations and of conserved charges are six, and indeed match each other. In the following, we do not explicitly solve the constraint equation.

Time evolution as renormalized/optimized perturbation theory
We here explain challenging points of the problem and sketch our strategy to derive the constitutive relations. At the initial time t 0 , the expectation values are simply specified by the initial local Gibbs distributionρ LG [t 0 ; λ t 0 ] and are provided as functionals of the thermodynamic parameters λ t 0 at that time. Recalling the one-to-one correspondences between the values of parameters λ and those of charge densities T 0 ν and J 0ν , one can deduce the presence of constitutive relations (3.11) though its explicit form needs further investigation. At time t (> t 0 ) of our interest in the future, the situation becomes more complicated. If we express the average currents T µν (t, x) and Ĵ µν (t, x) by applying the derivative expantion of the parameter at the initial time t 0 , we soon encounter a meaningless divergent result due to a lack of the backreaction to thermodynamic variables (currents endlessly flow according to the initial parameter profile). Then, one would expect that the average currents T µν (t, x) and Ĵ µν (t, x) are expressed by a new parameter set λ t = {β µ (t, x), H µ (t, x)} at time t instead of the initial one λ 0 . However, the density operator in the Heisenberg picture does not evolve in time and, furthermore, we have not yet defined such a new parameter set.
To resolve the above difficulty, we turn our attention to what we practically do in solving hydrodynamic equations. Suppose that one knows the configuration of charge densities at time t and, moreover, is equipped with constitutive relations in terms of the inverse temperature, fluid velocity, and magnetic field. Then, one can solve the hydrodynamic equations within an infinitesimal time step dt to obtain a new configuration of the average charge densities, T 0 ν (t + dt, x) and Ĵ 0ν (t + dt, x) . To keep running hydrodynamic equations with the same algorithm that relies on the constitutive relations in terms of the parameter λ, we then need to feedback information on T 0 ν (t + dt, x) and Ĵ 0ν (t + dt, x) to update the corresponding parameters λ t+dt (or define λ t+dt from the average charge densities). One standard way is to use the local thermodynamic relation, which enable us to relates the charge densities to their conjugate parameters. From our perspective, this assumption can be interpreted as the matching condition for the full average charge density to a "hypothetical" local Gibbs average 9 : The right-hand sides of these equations contain a new parameter set λ t+dt , which one can determine with the left-hand sides which have been obtained by solving the hydrodynamic equation. Note that the number of the matching condition coincides with that of parameters λ(t). This matching condition completes our algorithm to solve the averaged conservation laws relying on the constitutive relations expressed by the local thermodynamic parameter λ. In short, we generally employ the matching condition These conditions enable us to define the local thermodynamic parameters λ t at a later time t (> t 0 ). Based on the new parameter set λ t , we are now ready to explain how we express the average currents T µν (t, x) and Ĵ µν (t, x) in terms of λ t . As we mentioned, the initial density operator does not evolve in time in the Heisenberg picture, and it is not clear how the new parameter set λ t manifests itself in expectation values of the conserved charge densities. However, a simple trick of the renormalized/optimized perturbation [71][72][73] resolves the problem: we formally rewrite the initial density operator aŝ Note that the operatorΣ[t, t 0 ; λ] gives the entropy difference between time t and t 0 , so that we will call it the entropy production operator. This arrangement is always allowed since the right-hand side is just identical to the original definition ofρ LG [t 0 ; λ t 0 ]. Keeping the noncommutative properties of the operators in mind, we factorize the exponential factor aŝ where T τ denotes the time-ordering operator for τ . We defined a shorthand notation O τ ≡ e τŜ[t;λ]Ô e −τŜ[t;λ] for an arbitrary operatorÔ. Thanks to this rearranged factorization, we obtain an identity, which expresses the expectation value of an arbitrary operatorÔ(t) at time t by the local Gibbs average at the same time: This observation motivates us to evaluate Ô (t) by expanding the exponential operator U [t, t 0 ; λ], which leads to an expansion on top of the "hypothetical" local equilibrium state at each time slice. Namely, one could expand the "evolution operator"Û [t, t 0 ; λ] with respect toΣ[t, t 0 ; λ], in which the leading term reduces to the thermal expectation value Ô (t) LG t in line with our expectation. This prescription gives the renormalized/optimized perturbation theory for the hydrodynamic time evolution.
In order for this renormalized/optimized peturbation to work, the entropy production operatorΣ[t, t 0 ; λ] needs to be a small quantity so that it works as a controlled expansion parameter. One can confirm this is indeed the case in the sense of the derivative expansion.
By rewriting the entropy production operator aŝ 16) one finds that the right-hand side consists only of the parameter (and background gauge field) derivatives. Here, we introduced LG t for an arbitrary Heisenberg operatorÔ(t). We also used the Stokes' theorem ∂ t dΣ tµ f µ = d 3 x √ −g∇ µ f µ and the covariant conservation laws (3.1), together with the following expression for the time derivative of the Massieu-Planck functional Ψ[λ t ]: The subtraction of ∂ t Ψ[λ t ] leads to displacements δT µν and δĴ µν in Eq. (3.16). The derivative quantities on the right-hand side of Eq. (3.16) are expected to take small values when hydrodynamics works in describing the low-energy behavior (i.e., low-frequency and long-wavelength limits) of conserved charge densities. In such a situation, one may systematically organize a derivative expansion of the constitutive equations via the expansion of Our basic strategy has been implemented with the above rearrangements. Based on those crucial observations with introduction ofÛ [t, t 0 ; λ], we now decompose the problem into two separate parts: the local equilibrium average and deviation from it. This decomposition is accomplished by arranging the average currents as where the displacements are given by Those terms contain at least oneΣ[t, t 0 ; λ], and thus, lead to derivative corrections to the constitutive relations. In the subsequent sections, we will separately investigate the first and second terms of Eq. (3.18) in detail to derive the constitutive relations for relativistic magnetohydrodynamics.
Before closing this section, we provide a comment on the matching condition used to update the parameter set λ in the time increment. Using the above notations, one can rewrite the matching conditions as δT 0 ν (t, x) = 0 and δĴ0 ν (t, x) = 0, indicating the absence of the derivative corrections to conserved charges. However, this does not mean that the system is in a local equilibrium state at every time slice t since the density operator is always given by the initial local Gibbs distribution. Putting it differently, the matching condition only means that we have matched the conserved charge densities at time t to those in the hypothetical local equilibrium forms by adjusting the parameter set λ t (note that the full averages of other operators do not agree with the newly introduced local Gibbs average of those operators). Thus, it is legitimate to say that introduction of a new parameter set λ t with the matching condition (3.12) should be understood as one particular approximation scheme such as the one employed in the renormalized/optimized perturbation theory [71][72][73]. From this viewpoint, one can identify the matching condition (3.12) as a counterpart of the fastest apparent convergence condition used in the optimized perturbation theory [71]. It is also worth mentioning that our matching condition (3.12) with the decomposition β µ = βu µ (u µ u µ = −1) leads to one natural definition of a "fluid-velocity" u µ , which is known as that defined in the so-called beta frame [62]. This fluid velocity is, in general, different from that of the Landau-Lifshitz frame (or, of courase, that of the Eckart frame). This "frame choice" is usually implicitly assumed in the phenomenological derivation of relativistic hydrodynamics, but is not at all a unique choice (see, e.g., Refs. [74,75]).

Constitutive relations: Exact results
In this section, we derive an exact formula for the constitutive equations (3.11) for both the nondissipative and dissipative components. Before starting a discussion, we briefly summarize our notations for the coordinate system (see Refs. [63,64,76] for more systematic explanations). While we have introduced a covariant notation to describe constant time hypersurface Σ t in the previous section, other parts of the coordinate system have been left implicit. In order to make our formulation to be manifestly covariant under diffeomorphism, we put the QED plasma in a curved spacetime equipped with the background vierbein (metric) as before. However, we emphasize that even in a flat spacetime there is a benefit to maintain the general covariance by using the curvlinear coordinate system-as was the case in the covariant formulation (super-many-time theory) of QED [1][2][3][4][5][6][7][8]-so that we can choose a suitable coordinate system to describe the local equilibrium state with a nonuniform flow velocity.
Letxμ be our choice of the coordinate system and assume existence of an inverse mapping x µ (x) to another coordinate system x µ = (t, x) in general. Using a scalar time functiont(x), we introduce a foliation speficying equal-time (or constant time) hypersurfaces Σt, on each of which the time function takes a constant value:t(x) = const. Expressing its component asxμ(x) ≡ t (x),x(x) , one findsx defines the coordinate on the spatial hypersurface Σt at each time slice parametrized by a value of the time coordinatet. Then, we define two timelike vectors which specify the time direction of the local coordinate system and the normal direction to the hypersurface Σt. Explicitly, they are given by The time vector t µ and the normal vector n µ are proportional to each other if we use the Cartesian coordinate system. The normal vector has a unit norm n µ n µ = −1, and the normalization function N (x) is called the lapse function. On the other hand, the time vector is not necessarily normalized, and its norm represents how we measure the scale along local time directions. We also require that the normal vector satisfy a condition n∧dn = 0 so that the hypersurface covers all the space without intersecting another hypersurface thanks to the Frobenius theorem. By using the normal vector, we can write the infinitesimal surface element vector as On the rightmost side, we introduced an (invariant) spatial volume element on the hypersurface by using the determinant of the induced metric: γ = det γ µν with γ µν ≡ g µν + n µ n ν .
Decomposing the time vector with the normal vector n µ and its orthogonal component, where N µ is called a shift vector. Using this (3 + 1)-decomposition, we can express the (inverse) metric gμν (gμν) in the coordinate systemx µ as where we introduced Nī = γījNj and the inverse γīj of the induced metric satisfying γīkγkj = δj i . Note that the determinant of spacetime metric g µν is given by product of the Lapse N and spatial determinant γ as √ −g = N √ γ. Since we will basically use only the coordinate system ofx µ in the rest of this paper, we will omit all the overbars to express our coordinate systemx µ , which is thus simply denoted as x µ .

Nondissipative part from the Massieu-Planck functional
In this subsection, we examine the nondissipative (or local equilibrium) part given by the first term in Eq. (3.18). First, we show that the Massieu-Planck functional (3.5) works as a generating functional of the conserved currents-not only for the charge density as given in Eq. (3.6). This means that all the transport phenomena taking place in the local equilibrium state could be derived from one functional Ψ[λ t ]. Then, we provide a path-integral representation of the Massieu-Planck functional, and show that the local equilibrium configuration can be generally interpreted as an emergent curved spacetime geometry/two-form gauge connection in the imaginary-time (or Matsubara) formalism [77,78]. Notion of thermally induced curved spacetime/gauge connection, in particular its symmetry properties, allows us to give a systematic construction of the Massieu-Planck functional. Since this formalism can be constructed in parallel to previous works [63,64,66], we will highlight differences from the previous works that are induced by the one-form symmetry newly introduced in this work.
which leads to a compact expression of the Massieu-Planck functional as Ψ[t; λ, j] = log Tr e −K[t;λ,j] . Since the background field j = {e a µ , b µν } plays a crucial role in the following discussion, we explicitly show it in the functional argument.
On the basis of this definition, let us consider a variation of the Massieu-Planck functional under an infinitesimal transformation acting on both the local Gibbs parameter λ and background field j. To be explicit, we consider a combination of the infinitesimal general coorinate transformation and 1-form gauge transformation acting on them as with a specific choice of the transformation parameter ξ µ = β µ and θ µ = (H µ − β ν b νµ ) ( denotes an infinitesimal constant). Note that the background two-form gauge field has a contribution from the one-form gauge transformation, and that the time function t(x) specifying the equal-time hypersurface also transforms under the general coordinate transformation. Direct computation leads to the following explicit forms: where we used the vierbein postulate ∇ µ e a ν +ω a µ b e b ν = 0 for the third equation, and defined β ≡ −β µ n µ for the last equation.
From the definition of the Massieu-Planck functional (Ψ[t; λ, j] = log Tr e −K[t;λ,j] ), one immediately finds a simple variational relation LG t , where δ para λ means that the transformation only acts on the paramters inK[t; λ, j]. The vital point here is that the functional operatorK[t; λ, j] remains invariant under the simultaneous transformation acting both on the parameter and operator; namely, δ tot λK ≡ δ para λK + δ op λK = 0. Since the operator variation is generated by taking a self commutation relation, it identically vanishes: δ op λK = [iK,K] = 0. Therefore, the above variational relation leads to the following identity for the Massieu-Planck functional: Let us then specify an explicit form of δ λ Ψ[t; λ, j] and derive the variational formula. Recalling that Ψ[t; λ, j] depends on the time function t(x), which also transforms as a scalar, we evaluate the variation of the Massieu-Planck functional with respect to the time function as follows: where we used the conservation laws at the operator level (3.1). Then, using the transformation rule (4.7) together with the operator version of the constraint (2.11), we obtained the rightmost side of Eq. (4.9) (see Refs. [64,66] for more details). Taking into account that contribution, we obtain δ λ Ψ as follows: To obtain the third line, recalling Eq. (3.6), we performed an integration by parts as where we neglected the surface term and used the zeroth-component of the conservation law forĴ 0µ . As a consequence, the identity δ λ Ψ = 0 results in (4.12) Here, the variations δ λ e a µ and δ λ b µν are still left arbitrary since any parameter configuration λ should satisfy δ λ Ψ = 0 irrespective of the hydrodynamic evolution of β µ and H µ in the transformation parameters. Therefore, we eventually find the following variational formulas: where we defined a shorthand notation for anti-symmetrization as It is worth emphasizing that this gives an exact formula without relying on the derivative expansion. Using those exact formulas, one can extract all the transport phenomena taking place in the local thermal equilibrium state from the single functional Ψ[λ t ].
The derivative expansion of each local equilibrium constitutive relation is reduced from that of the Massieu-Planck functional on an order-by-order basis. In short, construction of the constitutive relation for the nondissipative parts amounts to computing the Massieu-Planck functional Ψ[λ t ].
Hydrostatic gauge. While we have derived the general variational formulas without specifying the coordinate system, there is a useful gauge (or choice of the coordinate system), which we call hydrodstatic gauge [64]. To get motivated to introduce the hydrostatic gauge, we turn our attention to the fact that each fluid element has its local rest frame. On the other hand, we have acquired the freedom to choose the local time-direction of our coordinate system thanks to the covariant formulation. Therefore, we can make a suitable choice of the local coordinate system, in that all the fluid elements in the system look as if static.
Here, we introduce the hydrostatic gauge and simplify the variational formulas (4.13). The hydrostatic gauge is defined by the following gauge conditions: where we introduced an arbitrary reference inverse temperature β 0 which has the same mass dimension as that of the inverse temperature. Here, the subscripts stand for the quantities given in the hydrostatic gauge. The first condition for the time vector t µ indicates that the local time direction is taken along the local fluid flow and the scale of the time axis is specified by the local inverse temperature [63]. Besides, we put the second condition to extend the hydrostatic gauge for magnetohydrodynamics, which means that we interpret the reduced magnetic field H µ as a temporal component of the background two-form gauge field 10 .
Thanks to the hydrostatic condition, we obtain H µ (x) − β ν b νµ hs = 0 so that one-form gauge transformation in Eq. (4.7) becomes trivial as θ ν hs = (H µ −β ν b νµ ) hs = 0. Thus, we can simplify a little bit complicated transformation δ λ in Eq. (4.7) to just the lie derivative along the time direction £ β = β 0 £ t hs . Furthemore, the variation with respect to H µ is absorbed into that with respect to b µν so that we only need to take account of the variation with respect to background fields. Therefore, the identity δ λ Ψ = 0 in the hydrostatic gauge takes the following simple form: where we used β = β 0 N hs . Therefore, we obtain simpler variational formulas in the hydrostatic gauge This expression means that the Massieu-Planck functional in the hydrostatic gauge serves as an analogue of the action functional, which implies a possibility of the action principle description of hydordynamics (See, e.g., Ref. [44] for recent discussions).

Path-integral formalism for the Massieu-Planck functional
According to the variational formulas shown in the previous section, we only need to compute the Massieu-Planck functional to evaluate the local equilibrium averages of current operators. In the derivation of the variational formulas, we only used the consequences resulting from symmetries, specifically the spacetime and magnetic one-form symmetries, which are independent of microscopic details of the system. We now specify the QED action (2.3) as the underlying microscopic theory 11 , and derive the path-integral representation of the Massiue-Planck functional (3.5). We take the the operator (Hamiltonian) formalism on the basis of the (3 + 1) decomposition [76] as the starting point. The relevant dynamical degrees of freedom are the sum of the Dirac field and U(1) gauge field separately considered in Ref. [64], so that we can apply the path-integral formulas shown therein. Putting an emphasis on the specific points arising in the presence of the magnetic one-form symmetry, we here sketch our derivation (see Ref. [64] for the detailed calculation).
Derivation of path-integral formula. We shall elaborate the QED action (2.3) equipped with the background field. Since the fermion sector is almost same as the one in Ref. [64], we here concentrate on the gauge field sector, which contains the two-form background field. The (3 + 1) decomposition of the photon sector leads to the following Lagrangian: from which one can define the canonical momentum of the gauge field A i as Thus, the presence of the background two-form field generates a c-number shift for the conjugate momentum, which does not affect the canonical commutation relation. It is useful to define a shifted momentuḿ which satisfies the same commutation relation as Π i does.
On the other hand, we need to pay attention to the Gauss's law since the Maxwell equation contains the contribution from the source term: 11 While we focus on QED in (3 + 1) dimension, we expect that the emergence of the same thermal background and the geometrical interpretation of the local thermal equilibrium state holds in any system endowed with one-form symmetry.
where we introduced an electric current j µ by a variation of the fermionic part of the QED action S mat QED with respect to the gauge field. Expressing the contribution from the background field as j µ bkd (x) ≡ 1 2 ∇ ν (ε µνρσ b ρσ ), one can regard the insertion of the background two-form field as putting the background electric current (see also Ref. [35] for more discussions). The Gauss's law from the temporal component of the Maxwell equation reads which contains a c-number shift if we useΠ i . Note that the temporal component of the canonical momentum vanishes as usual, Π 0 = 0, due to the gauge redundancy, so that there is no contribution from the background field. As a next step, let us write down the conserved charge densities based on their definitions (2.8). Recalling that the coupling term induced by the two-form external field is independent of the vierbein [see the second line of Eq. (4.17)], we find that the background two-form field does not contribute to the energy-momentum tensor [recall also the energy-momentum tensor (2.8) is defined via the variation of the action with respect to the vierbein]. Therefore, the energy-momentum tensor has exactly the same form as the conventional case. In terms ofΠ i defined in Eq. (4.19), the explicit forms of the photon contribution to the conserved charge densities read Note that (J 0i ) photon coinsides with the total magnetic flux densities J 0i since the charged matter sector does not contribute. We have specified all the ingredients induced by the two-form background field, which enables us to write down the path-integral formula for the Massieu-Planck functional. With the help of the previous result for the path-integral formula (see the detailed account in Ref. [64]), we obtain the following expression in the axial gauge (A 3 = 0): where we introduced the phase-space Lagrangian L H as where we used a decomposed form of the fluid vector β µ = βu µ /β 0 and the reduced magnetic field H µ = βH µ with a normalized four-velocity u µ (u µ u µ = −1) and the ratio of the inverse temperature to the constant reference e σ(x) = β(x)/β 0 . A new vierbeinẽ a µ and its inversẽ e µ a as well as a new covariant derivative D µ will be explained shortly. We also defined the field strength tensor including the imaginary-time direction as F 0i = i(∂ τ A i − ∂ i A 0 ). The use of the axial gauge A 3 = 0 leads to the insertion of δ(A 3 ) det(∇ 3 ), which is replaced by δ(F ) det(∂F/∂α) in a general gauge satisifying F = 0 with the gauge parameter α.
Performing the functional integration with respect to the bilinear form of the modified conjugate momentumΠ i , we obtain the path-integral formula for the Massieu-Planck functional: with the imaginary-time action This Lagrangian density has the same form as the original one (2.3). This indicates that information on the local thermal equilibrium is thoroughly encoded in the thermally induced background field with the aforementioned vierbein and two-form gauge fieldb µν defined bỹ from which one finds the orthogonality ofẽ i a to the fluid velocity as u aẽ i a = 0. In the same way as the original spacetime metric, the emergent thermal metric and its inverse satisfỹ g µν = η abẽ a µẽ b ν andg µν = η abẽ µ aẽ ν b , respectively. Their explicit forms are obtained by replacing the Lapse function N and shift vector N i in the (3 + 1) parametrization (4.4) with the following thermal ones: Here, we used the previously defined dilaton-like parametrization e σ(x) = β(x)/β 0 . From this parametrization, one can easily extract the measure factor √ −g =Ñ √ γ = e σ u 0 √ −g. We also introduced a partial derivative in the thermal spacetime as∂ µ = (i∂ τ , ∂ i ) and a covariant derivative for the Dirac field by replacing the spin connection in Eqs. (2.4)-(2.5) with that of thermal spacetime composed of the thermal vierbein (4.26). We note that the constraint (2.11) coming from the local Lorentz symmetry is crucial to obtain the correct expression for the imaginary-time spin connection (see Ref. [64] for the details).
Symmetries of the Massieu-Planck functional. Based on the path-integral formula, we demonstrate symmetries of the Massieu-Planck functional. The above result indicates that one needs to perform the path integral in the presence of the emergent background to investigate the transport phenomena taking place in the locally equilibrated QED plasma. Notice that one can summarize information on the background field in the form of the emergent line element ds 2 and two-form gauge connectionb in the Kaluza-Klein parametrizations: where we introduced (dt, dx) ≡ (−idτ, dx) and identified the parameters

(4.32)
This Kaluza-Klein gauge symmetry results from the independence of induced background fieldsj ≡ {ẽ a µ ,b µν } in the imaginary-time coordinate (note, however, that they are still real-time dependent quantities taking values defined on each equal-time hypersurface of the foliation while its dependence is not explicit in our notations). Investigating the transformation rule induced by the Kaluza-Klein gauge transformation (4.32), one finds that the lower-time and upper-spatial indices are inert while the upper-time and lower-spatial indices transform [41,63,64]. For example, the Kaluza-Klein gauge transformation acts on arbitrary vectors A µ andB µ in thermal spacetime as Noting a i =g 0i /g 00 , one sees the Kaluza-Klein gauge field a i indeed obeys this transformation rule. Using the Kaluza-Klein gauge field, one can construct the Kaluza-Klein gauge invariant lower-spatial (upper-temporal) component A i (B 0 ) as One sees that γ ij andb ij in Eq. (4.31) are indeed constructed according to this procedure. As will be shown shortly, the Kaluza-Klein gauge invariance restricts possible forms of the Massieu-Planck functional.
In addition to the Kauluza-Klein gauge invariance, the Massieu-Planck functional also enjoys spatial differmorphism and two-form gauge symmetry applied to the original two-form gauge field b µν as We note that only the spatial component of the two-form gauge field,b ij , transforms under the original gauge transformation, andb 0i is inert under that (recall thatb 0i does not contains b 0i ). The emergent line element is the same as the previously studied one, and a new ingredient for the QED plasma is the presence of the two-form background connectioñ b induced by the reduced magnetic field H i . Furthemore, we also note that there is another restrction coming from discrete symmetry. For instance, the insertion of charge conjugation C and time reversal T leads to the following identities: where Θ λ C = ±1 (Θ λ T = ±1) denotes an eigenvalue of the charge density operatorsĉ under charge conjugation (time reversal): Cĉ(t, x) C −1 = Θ λ Cĉ (t, x) and Tĉ(t, x) T −1 = Θ λ Tĉ (t, x). Here, we note that the time reversal transformation is defined as the reflection with respect to the given hypersurface Σ t under consideration, so that the time-argument of the operator c(t, x) is unchanged. We identify C and T eigenvalues as In summary, we have obtained the path-integral formula for the Massieu-Planck functional for the locally equilibrated QED plasma in Eq. (4.24). We can fully capture the effects resulting from inhomogeneous thermodynamic parameters such as the local temperature by the use of the emergent backgrounds with the line element (4.29) and two-form gauge connection (4.30). The background data tells us the symmetries of the Massieu-Planck functional; the Kaluza-Klein gauge symmetry (4.32), spatial diffeormorphism (4.35), and spatial one-form gauge symmetry (4.36). These symmetries together with the discretesymmetry argument (4.37) constrain possible forms of the Massieu-Planck functional. For instance, one can use e σ and the magnitude ofb 0i as basic building blocks, which give the local inverse temperature and the magnitude of the reduced magnetic field since they are Kaluza-Klein and one-form gauge invarint quantites. On the other hand, possible forms of the dependences on the Kaluza-Klein gauge field a i and the spatial component of the two-form gauge fieldb ij are restricted by the symmetries (4.32)-(4.36). Consequently, dependences on a i andb ij are allowed only in the form of their field strengths, which are accompanied by, at least, one spatial derivative. We can now find an exhaustive list of the Kaluza-Klein and one-form gauge invariant quantities up to the first-order derivative: Since the Massieu-Planck functional enjoys spatial diffeomorphism too, we can construct spatial scalar quantities from these building blocks. In this way, at a given order of the derivative expansion, one can identify the finite number of the invariants compatible with all symmetries. We will explicitly use the symmetry constraints clarified here when we organize a derivative expansion of the Massieu-Planck functional in the subsequent section.

Dissipative part
Let us now look back on the real-time (not imaginary-time) evolution discussed in Sec. 3.3. As we found there, the dissipative part δT µ ν and δĴ µν are associated with the entropy production operatorΣ[t, t 0 ; λ] given in Eq. (3.16). In this section, we focus on this operator and rewrite it in a formal but more suitable form for the derivation of dissipative corrections. According to our matching condition (3.12), dissipative corrections are purely spatial, i.e., Therefore, the dissipative effects occur in such a way that the conserved charges distributed on the hypersurface Σ t diffuse in the purely spatial direction perpendicular to n µ . It is thus useful to introduce a projection matrix, which decomposes the curved spacetime index into the time and spatial components as where v µ ≡ t µ /N is the normalized time vector satisfying v µ n µ = −1. Clearly, this matrix has desired orthogonal properties P µ ν n µ = 0 and P µ ν v ν = 0. When one acts the projection matrix on an arbitrary tensor index, a purely spatial index is returned in our coordinate system. Now, by the use of the projection matrix, the derivative operator is decomposed as Due to the matching condition (4.40), only the spatial derivative ∇ ⊥µ will be shown to appear in the leading-order dissipative corrections. This decomposition allows us to rewrite the entropy production operator (3.16) aŝ where we defined a short-hand notation The crucial point here is that the first two terms in Eq. (4.44) contain the time derivatives of parameters ∇ t λ ≡ {∇ t β µ , ∇ t H µ }, although it is not obvious at this stage whether those time-derivative terms cause any problem. It will turn out that further manipulation is necessary to eliminate them so that one can obtain the well-defined Green-Kubo formulas for the transport coefficients. In fact, without eliminating those time-derivative terms, they would propagate through the following formulation and finally manifest themselves in the constitutive relations. Then, the current operator contains the contributions from the (linear) hydrodynamic modes, and the gapless hydrodynamic modes contribute to current-current correlators in the resulting Green-Kubo formulas. Consequently, when taking the hydrodynamic limit of correlators, one will suffer from the divergence, and cannot obtain the correct transport coefficients. We will, therefore, appropriately eliminate the time derivatives by the use of the equations of motion (3.10).
Following the procedure developed in Ref. [63], we formally solve the equation of motion (3.10) and express the time derivatives ∇ t λ in terms of the spatial derivatives ∇ ⊥µ λ. Using the collective (vector-like) notation-the conserved charge densityĉ a = {T 0 µ ,Ĵ 0µ }, conserved currentĴ µ a = {T µ ν ,Ĵ µν } satisfying ∇ µĴ µ a =Ĉ a , and the parameter derivatives ∇ µ Λ a = {∇ µ β ν , ∇ µ H µ }-we obtain the following formal expression (see Refs. [64,65] for the details): We introduced inverse susceptibilities as a second variation of the entropy functional: , (4.47) and the Kubo-Mori-Bogoliubov (KMB) inner product:  Here, we defined a projected operator This projection operatorP gives the local Gibbs version of Mori's projection operator [79]. We also defined thermodynamically conjugate operators Equation (4.49) gives the exact expression for the entropy production operator. As a result, the combination of Eqs. (3.19) and (4.49) gives the exact form of constitutive relations induced by the deviation from the local thermal equilibrium. Note that Eq. (4.49) contains δT µ ν and δĴ µν , so that we need to solve Eqs. (3.19) and (4.49) in a self-consistent manner. However, as discussed in the next section, the terms in the second line of Eq. (4.49) are found to be the second-order derivative corrections with the help of an appropriate power counting scheme. This observation simplifies the procedure to obtain a self-consistent constitutive relation within the first-order derivative expansion.

Constitutive equations: Derivative expansion
In this section, applying a derivative expansion to the exact results developed in the previous section, we derive the constitutive relations for relativistic magnetohydrodynamic up to the first derivative order.
We first present our power counting scheme. In this paper, we employ the simplest power counting scheme, where all the thermodynamic parameters and background fields are counted as zeroth-order quantities 12 : As a result, one finds that the first line of the entropy production operator (4.49) contains O(∂ 1 ) contribution, while the second line provides only the higher-order terms (recall that δT µ ν and δĴ µν contain, at least, one derivative and H ναβ = O(∂ 1 ) in our counting). Thus, to derive the constitutive relations up to the first derivative, it is sufficient to use the simplified form of the entropy production operator:

Leading order: Ideal magnetohydrodynamics
In the leading-order derivative expansion, we only need to evaluate the expectation value over the local Gibbs distribution [recall that the entropy production (5.24) is accompanied by at least one derivative]. As shown in Eq. where we used the fact that H 0 = 0 to express γ ij H i H j in the covariant manner. Here, we also introduced an upper-component induced metric and normal vector in thermal spacetime as γ µν ≡g µν +ñ µñν =ẽ µ aẽ ν b (η ab + n a n b ), n µ ≡g µνñ µ =ẽ µ a n a , (5.4) whose components can be directly read off from Eq. (4.28). Note that n a in the second expression is not the thermal vector but the original one with the local Lorentz index. The normalized property ofñ µ in the thermal spacetime, orñ µñνg µν = −1, immediately leads toñ µγ µν = 0. Using these zeroth-order invariants, we identify the leading-order expression of the Massieu-Planck functional: with β = −g µν β µ β ν and H = γ µν H µ H ν . We note again that the metric contracting their indices is different from each other: g µν for β, andγ µν for H. Also, note that the measure factor √ −g is necessary to make the volume element dτ d 3 Recalling that u aẽ i a = 0 as given just after Eq. (4.27) together with H 0 = 0, one finds that the former is orthogonal to the fluid vector: u a H a = 0.
Combining the result in Eq. (5.5) with the variational formula shown in the previous section, we can derive the leading-order constitutive relations that consist ideal magnetohydrodynamics 13 . For that purpose, let us then consider the variation with the fixed hypersurface. Useful variational formulas of the various quantities are summarized as 13 Here, we define our ideal magnetohydrodynamics with the constitutive relations at the non-derivative (or zeroth) order. The conventional definition of the "ideal" magnetohydrodynamics goes in a different way (see, e.g., Refs. [29][30][31][32]) and relies on an illegitimate notation of an "infinite" electrical conductivity that is a dimensionfull quantity and is, moreover, not defined a priori in the formulation of hydrodynamics.
where, in addition to β µ = βu µ , we defined the following decompositions of H a and H µ : H a = Hh a , H µ = Hb µ with b µẽ a µ h a = 1.
For the computation of δp, we used the variation of the thermal vierbein given by µ e a ν δβ ν + δe a µ (5.9) and the orthogonal properties: u a h a = 0 and H µñ µ = 0. Equation (5.9) follows from the fact that the fluid-vector β µ generally gives the time vectort µ in thermal spacetime as β µ = β 0t µ (when we use the hydrostatic gauge, it will also match the time vector along the real-time direction). With the help of the above formulas, one can quickly evaluate the variation of the leading-order Massieu-Planck functional (5.5) as Based on the variational formula (4.13), we obtain the constitutive relations of ideal magnetohydrodynamics While the first two terms of the energy-momentum tensor (5.11a) are common to the usual relativistic perfect fluid without one-form symmetry, the last term is induced by magnetic one-form symmetry inherent in QED. Although it is not mandatory, one can also extract the conserved charge densities from Eqs. (3.6) and (5.10), which allows us to express ∂p/∂β and ∂p/∂ H in terms of conserved charge densities. The constitutive relation (5.11b) describes the local electromagnetic field in the QED plasma flowing with a finite velocity, sinceĴ 0i andĴ ij are the magnetic flux density and electric field, respectively. One then finds only the magnetic flux sits in the fluid rest frame at the ideal order, i.e., as long as the system is in the local thermal equilibrium state. The new term in Eq. (5.11a) describes a pressure anisotropy induced by the magnetic flux density in the QED plasma. To see this, it is helpful to rewrite the energy-momentum tensor (5.11a) into the following decomposed form: Here, we defined energy density e ≡ u µ u a T µ a (t, x) LG  with magnitudes of the magnetic flux density B ≡ β∂p/∂ H and effective magnetic field H = β H. The difference between the parallel and perpendicular components is evident in this expression: ∆p ≡ p ⊥ −p = HB, so that the parallel pressure is effectively reduced (see Fig. 1 for HB > 0). Note again that the matter and magnetic-field contributions are never separated, and p itself contains the contribution of the magnetic flux such as the Maxwell stress and magnetization. Remarkably, T µ a (t, x) LG t u a = −eu µ automatically follows from the orthogonal property u a h a = 0. Thus, we find that our fluid velocity u µ agrees with that in the Landau frame in ideal magnetohydrodynamics.
On the other hand, the energy-momentum tensor in the conventional approach is given by (see Eqs. (14) and (15) of Ref. [33]) where b µ is a normalized spatial vector along the magnetic field B µ conv = B conv b µ . This expression suggests the following clear identification between our result (5.12) and the conventional one: H ↔ H conv , B ↔ B conv . According to the tensor structures, the energy density e conv and perpendicular pressure p conv are also identified with ours. However, since the definition of their variables goes in a different way, we shall briefly summarize the points (see Ref. [33] for detailed accounts). First of all, most of variables in the conventional approach are defined in the additive forms: e conv = ε + 1 2 B 2 conv and p conv = p matt + 1 2 B 2 conv − M conv B conv . The first and second terms come from the separate contributions of the matter and Maxwell-stress parts. The perpendicular pressure as well as ε contain a contribution of the medium response, that is, the magnetization potential. Indeed, a part of their energy density ε is defined with the inclusion of the magnetization potential, i.e, ε = ε matt − M conv B conv , where M conv denotes the magnitude of a magnetization vector as mentioned below Eq. (19) of Ref. [33]. The cooperative contributions of the Maxwell stress and magnetization potential lead to the resultant anisotropic pressure ∆p conv = p conv − (p matt − 1 2 B 2 conv ) = H conv B conv with the (in-medium) magnetic field H conv = B conv − M conv . All these identifications clarify the breakdown for our energy den-sity and pressure under the assumption of a separation among the matter, Maxwell-stress and magnetization contributions.

Next-to-leading order: Dissipative magnetohydrodynamics
Derivation of constitutive relations. Let us proceed to the derivative expansion in the next-to-leading order. One can examine possible derivative corrections from the derivative expansions of the Massieu-Planck functional and the "evolution operator"Û [t, t 0 ; λ] in Eq. (3.19). However, the charge-neutral QED plasma under consideration is free from the first-order derivative corrections to the Massieu-Planck functional owing to the discretesymmetry constraint (4.37). Therefore, the local equilibrium parts of constitutive relations are the same as those of ideal magnetohydrodynamics in Eq. (5.11). In the following, we will identify the dissipative contributions to constitutive relations supplemented with the relevant Green-Kubo formula for transport coefficients.
As shown in Eq. (3.16), the "evolution operator"Û [t, t 0 ; λ] containsΣ τ [t, t; λ] in that all the terms have at least one derivative [see Eq. (5. 24)]. Therefore, we may writê Here, we maintain the term up to the first derivative in the current working accuracy, which allows us to use Eq. (5.24) for the entropy production operator. Inserting the above first-order expansion into Eq. (3.19), we have where we used the KMB inner product defined in Eq. (4.48).
It is now useful to decompose the projected current operatorsδT µν andδĴ µν by the use of available tensors at hand. The transversality of those dissipative terms to n µ greatly restricts possible tensor structures, so that one should utilize the projection matrix P µ ν defined in (4.41). Besides, we wish to introduce the projection matrix that allows for the decomposition with respect to the direction of the magnetic field. For that purpose, we introduce the lower-component vector h µ ≡ẽ a µ h a from h a , which satisfies b µ h µ = 1 thanks to Eq. (5.8). As was mentioned when we introduced the local Gibbs distribution, the zeroth component of the reduced magnetic field H 0 is set to be zero since the conserved chargê J 0ν does not have the temporal componentĴ 00 . Then, noting h µ = H µ / H according to its definition, we find that h µ is a spatial vector: Furthermore, using the expression of the original normal vector in terms of the thermal one, n µ = e σ u 0ñ µ , one also finds the upper-index vector b µ satisfies b µ n µ = 0. Thus, the vectors b µ and h µ serve as basic building blocks to prepare the desired projection matrix. Based on this observation, we define another projection matrix which satisfies the useful properties: This projection matrix enables us to extract a tensor index pependicular to n µ , v µ , h µ , and b µ . In short, we find two relevant tensors We define Ξ µν as an "inverse matrix" such that Ξ µρ Ξ ρν = ∆ µ ν . By using those tensors, we can decomposeδT µν andδĴ µν as where we used A (µν) = (A µν + A νµ )/2 and A [µν] = (A µν − A νµ )/2 to express symmetric and anti-symmetric components. Here, the projected operators are given bỹ Note that all the projected vectors and tensors are transverse; that is,δπ µ n µ =δπ µ h µ = 0, δÊ µ n µ =δÊ µ h µ = 0,δτ µν n ν =δτ µν h ν = 0, andδD µν n ν = δD µν h ν = 0 are satisfied. The rank-two tensors are symmetricδτ µν =δτ νµ and anti-symmetricδD µν = −δD νµ , respectively. Likewise, the flow gradient can be decomposed as where we defined θ ≡ b σ b σ ∇ ⊥ρ β σ and θ ⊥ ≡ Ξ ρσ ∇ ⊥ρ β σ . These scalars correspond to the expansion/compression flow in parallel and perpendicular to the magnetic field, respectively. Substituting the decomposed forms (5.21) and (5.23) into Eq. (5.24), we rewrite the entropy production operator aŝ

(5.24)
We then use this result in Eq. (5.16) to find Here, we assumed that the integration kernel (or the inner product for the projected current operators) behaves moderately in spacetime. In other words, we assume that (δĴ µ a (t, x),δĴ ν b (t , x )) t shows a sufficiently rapid decay with respect to spacetime coordinates (e.g. an exponential decay), so that we can apply the Markov approximation. This assumption, expected from the scale separation between hydrodynamic variables (charge densities) and non-hydrodynamic variables (currents), is crucial to obtain the local form of the constitutive relation. It is worth emphasizing that the use of the projection operator (4.50) is essentially important here since it eliminates the linear hydrodynamic mode from the integration kernel. In other words, without the projection operator, we have a contribution from the linear hydrodynamic mode, which causes an ill-mannered behavior of the integration kernel, and thus, prevents us from applying the Markov approximation employed in Eq. (5.25). The subtraction of the linear hydrodynamic mode motivates us to expect the well-mannered behavior described above, but it is fair to say that Eq. (5.25) is still an assumption that we cannot completely verify in this paper 14 .
With the help of Eqs. (5.16), (5.24), and (5.25), we finally specify the first-order dissipative corrections to constitutive relations as where we introduced the symmetric and traceless projection of a tensor as Here, we also defined a set of transport coefficients in Eq. (5.26) in the form of the spacetime 14 In fact, the low-dimensional hydrodynamics suffers from the nonlinear hydrodynamic fluctuation, which forbids us to derive the local form of the constitutive relation. In this case, we do not have the well-defined transport coefficient in the thermodynamic limit but have the scale-dependent one in the same manner as the running coupling constant in quantum field theory.
integral of the KMB inner product: (5.28) These equations give the Green-Kubo formula [67][68][69] for the transport coefficients of relativistic dissipative magnetohydrodynamics. From their roles in the constitutive relations (5.26), we identify ζ and η as the bulk and shear viscosities, respectively. Recalling that 0ijkĴ jk represents the electric field while 0ijk ∇ j H k 0ijk β ρ H ρjk ∝ j i bkd does the background electric current, we can identify ρ with the resistivity. Since the magnetic flux lines break a spatial rotational symmetry, we have two distinct coefficients for the bulk and shear viscosities and the resistivity, i.e., the parallel components ζ , η , ρ and perpendicular ones ζ ⊥ , η ⊥ , ρ ⊥ (see Figs. 2 and 3). Besides, we have two additional cross bulk viscosities ζ × , ζ × , which represent deviation of the pressure in parallel (perpendicular) to the magnetic flux in response to the compression/expansion in the other direction, i.e., the perpendicular (parallel) direction (see Fig. 4). We will show that they should take the same values within the first-order derivative expansion just below.
Onsager's reciplocal relation. Since ζ × and ζ × correspond to the reciprocal processes to each other, one can show that Onsager's reciprocal relation constrains their values [70]. To figure out Onsager's relation for transport coefficients within the first-order derivative expansion, one can reduceK[t; λ t ] in the Green-Kubo formula (5.28) to a global equilibrium one given byK where we introduced a total HamiltonianĤ and total magnetic fluxΦ i along   . 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 (compare with Fig. 3). Those two cases are reciprocal processes to each other.
We also reduce a spacetime background to the flat spacetime with the Minkowski metric. Then, the Green-Kubo formulas for the cross bulk viscosities read dτ Tr e −Keq e τKeqδp (t, x)e −τKeqδp ⊥ (0, 0) , dτ Tr e −Keq e τKeqδp ⊥ (t, x)e −τKeqδp (0, 0) . where we also used the Hermitian properties of the operators. We then reparametrize the spatial coordinate with the opposite sign, which leads to Therefore, we have shown Onsager's relation ζ × (β, H i ) = ζ × (β, H i ) instead of the naïve one with the flipped magnetic field. As a consequence, we have five (three bulk and two shear) viscosities and two resistivities for the QED plasma preserving the parity and chargeconjugation symmetries. This result from the nonequilibrium statistical operator method provides a basis for that from the phenomenological formulation [35].
Inequalities for the transport coefficients. The transport coefficients in Eq. (5.28) satisfy a set of inequalities as a consequence of the positive-definite property of the KMB inner product in the global thermal equilibrium: (Â,Â) eq ≥ 0 (the equality holds forÂ = 0). Since all the parallel/perpendicular components of transport coefficients are given by the inner products between the same operators, we immediately find the following constraints: On the other hand, the cross viscosity ζ × (= ζ × ) is not necessarily a positive-definite quantity since it is given by the inner product between the different operators. However, we can find a constraint on the cross viscosity ζ × by paying attention to the positive-definiteness of the inner product for the following linear combination with an arbitrary real-valued parameter y: Let us first consider the case with ζ ⊥ > 0, putting aside the case ζ ⊥ = 0. This inequality is satisfied for an arbitrary y if and only if the cross viscosity ζ × satisfies an inequality When ζ ⊥ = 0, the inequality (5.36) is consistent with the inequality for the other bulk viscosity ζ ≥ 0 only when ζ × = 0. This case can be combined with the inequality (5.37), which implies ζ × = 0 when ζ ⊥ = 0 (as long as ζ × is a real-valued quantity). The inequality (5.37) implies that ζ × is nonzero only when both of ζ ,⊥ are finite. As wrap up, we have shown the following set of inequalities If one wishes to find the minimum set of independent inequalities, one may remove either ζ ≥ 0 or ζ ⊥ ≥ 0 from the above list, because one of them, e.g., ζ ≥ 0 together with ζ ζ ⊥ − ζ 2 × ≥ 0 immediately implies the other inequality ζ ⊥ ≥ 0, and vice versa. These inequalities are consistent with those found from the phenomenological analysis in Ref. [35], in which the authors demand the local second law of thermodynamics; that is, the semi-positivity of the local entropy production rate ∇ µ ŝ µ (t, x) ≥ 0 with the entropy current ŝ µ (t, x) . We shall verify this assumption from our statistical mechanical viewpoint. Since our entropy production operator Σ[t; t 0 ; λ] gives a spacetime integral of the local entropy production rate, we can identify ∇ µ ŝ µ (t, x) with the integrand of Eq. (5.24) (see Ref. [63] for a definition of the entropy-current operator). Inserting the first-order constitutive relations (5.26) into the integrand of Σ[t; t 0 ; λ] in Eq. (5.24), we obtain the local entropy production rate as Note that the right-hand-side has the quadratic form, of which the (matrix) coefficient is given by the set of transport coefficients. Therefore, as a corollary of the set of inequalities (5.38), we can easily verify the semi-positive definiteness of the local entropy production rate-i.e., the local second law of thermodynamics-within the first-order derivative expansion: While this local positivity constraint is usually postulated as the starting point of the phenomenological formulation of hydrodynamics [28], we have verified it from our statistical mechanical formulation. For reader's convenience, we further compare our viscous coefficients and inequalities to those introduced in Refs. [33,34]. As discussed in Appendix A, we find a redundancy in the set of inequalities shown in Ref. [34] while we found consistent results in the other parts.

Summary and Discussions
In this paper, we have investigated a field-theoretical formulation of relativistic magnetohydrodynamics with the local Gibbs ensemble method (also known as the nonequilibrium statistical operator method). Starting from the QED Lagrangian in the (3 + 1) dimensional spacetime, we identified the symmetries of the system (including magnetic one-form symmetry) and the local Gibbs density operators relevant for describing the locally equilibrated QED plasma with a dynamical magnetic field. We have provided the exact results for both the local equilibrium and off-equilibrium parts of the constitutive relations, and then performed the derivative expansions for both of them up to the first order. Besides, we have clarified the Green-Kubo formulas, Onsager's reciprocal relation for the cross bulk viscosities ζ × , and the inequalities for transport coefficients without relying on the phenomenological assumptions. Here is a short summary of our assumption used in deriving hydrodynamic equations in this paper: (1) Special initial condition: The initial density operatorρ 0 is assumed to be the local Gibbs distribution,ρ 0 =ρ LG [t 0 ; λ t 0 ] defined in Eqs.  (3) Applicability of derivative expansion 2 (dissipatve part): Configurations of thermodynamic parameters are slowly varying functions of spacetime coordinates, and allow us to expand the "evolution operator"Û [t, t 0 ; λ] with respect to the entropy productionΣ[t, t 0 ; λ] as given in Eq. (5.15).
(4) Scale separtion and Markov approximation: The two-point current correlator (δĴ µ a (t, x),δĴ ν b (t , x )) t is assumed to show a sufficiently rapid decay with respect to spacetime coordinates, so that we can apply the Markov approximation (5.25).
While we have focused on the QED plasma, the same set of hydrodynamic equations may emerge as the universal low-energy effective theory of systems endowed with the same symmetries, i.e., the Poincaré invariance and one-form symmetry.
As concluding remarks, we present some interesting outlooks in order.
Extension to n-group symmetry. We first note that our framework is applicable to the recently proposed n-group symmetry [80]. As a simple example, let us consider the Abelian two-group symmetry equipped with the Poincaré invariance. The conservation laws are generalized as Here, we have an additional zero-form symmetry current J µ 5 , which, together with J µν , forms the Abelian two-group symmetry with a certain numberκ A . We also introduced a background field A 5 µ and its field strength F 5 µν ≡ ∂ µ A 5 ν − ∂ ν A 5 µ coupled with the current J µ 5 . Despite its anomaly-like appearance, the source term in the last equation is actually a consequence of the two-group symmetry (see Ref. [80] in detail). In this setup, one can introduce the conjugate parameter ν 5 to the new charge density J 0 5 , and the entropy functional takes a generalized form S[t; λ t ] = − dΣ tµ T µ ν (t, x)β ν (t, x) +Ĵ µν (t, x)H ν (t, x) +Ĵ µ 5 (t, x)ν 5 (t, x) + Ψ[λ t ]. (6.2) Repeating the same procedure presented in this paper, one can, for instance, find the firstorder constitutive relation δĴ µν for the Abelian two-group symmetry by the following replacement in Eq. (5.26b): This indicates that the background field of zero-form symmetry induces the two-form current. It is notable that the induced two-form current is proportional to the chemical potential ν 5 , showing a property similar to the chiral magnetic effect [81][82][83][84][85] (see Refs. [36,41,[86][87][88][89][90][91] for hydrodynamic derivations). The complete analysis, including nondissipative transport phenomena, may deserve further studies.
Spin transport. There are also several interesting prospects more relevant to familiar physical systems. The first direction is to clarify magnetohydrodynamics with a nonrelativistic component (like the proton in the astrophysical QED plasma) and to investigate roles of their spin degrees of freedom. In the relativistic case, one can interpret Eq. (2.11) as the total angular momentum conservation in the flat spacetime limit. Due to the mutual conversion between the spin and orbital parts of the angular momentum, the spin polarization itself is not a conserved quantity, and does not serve as a strict hydrodynamic variable (see discussions in, e.g., Ref. [92]). Nevertheless, the nonrelativistic (or large-mass) limit leads to an emergent internal SU(2) spin symmetry (like the heavy-quark symmetry [93]), and the spin density for the heavy fermion serves as an emergent hydrodynamic variable. While the interplay between the magnetic field and spin takes place via the Zeeman coupling as the leading correction, which leads to the approximate internal SU(2) symmetry, it is interesting to investigate the coupled dynamics of the nonrelativistic spin density and dynamical magnetic field as in Ref. [94]. Such a hydrodynamic framework will be applied to various systems such as the astrophysical QED plasma composed of (heavy) proton and (light) electron.
Evaluation of the transport coefficients. It is also an important direction to investigate the physical properties of the QED plasmas with the help of field-theoretical techniques. While we have clarified the universal forms of the constitutive relations, they contain several physical parameters; the equation of state p(β, H) and transport coefficients. As we have presented all the field-theoretical formulas necessary for determining those parameters, one can systematically evaluate them by using, e.g., the finite-temperature perturbation theory [95][96][97][98][99][100] (see also Refs. [51,52,[101][102][103][104] for the phenomenological treatments of the collisional effects) and/or the strong-coupling methods (see Refs. [105][106][107][108] for holographic calculations and Ref. [109] for a recent lattice Monte Carlo simulation).
In evaluation of the Green-Kubo formulas with the finite-temperature perturbation theory [95][96][97][98][99][100], the magnetic field is often treated as a background field B bkd and is coupled to the charged matter (which gives rise to the Landau quantization in the strong magnetic field). This is in sharp contrast to the current formulation since we fix the value of the dynamical magnetic flux density by using the (reduced) magnetic field, which is not coupled to the charged matter but to the magnetic flux. As a result, effects of the magnetic field on the charged matter are encoded in a different way as the conventional formulation. In many of the above studies, the resummation for the coupling between the charged fermions and magnetic field plays crucial roles. Therefore, it is interesting to push forward with the current "dual" formulation for the evaluation of physical parameters, which could open an avenue for a new resummation scheme. We leave those issues as future works. The coefficients with the superscripts "HK" and "HSR" are introduced in Refs. [34] and [33], respectively, while those without superscripts are introduced in this paper. We note that the viscous coefficients in the current paper agree with those in Eqs. (3.9), (3.10), (3.13), and (3.14) of Ref. [35]. According to the above correspondences (A.1), it turns out that the quite involved inequalities in Eq. (B.2) of HK [34] can be rewritten with our conventions in drastically simple forms: Since one may remove one inequality ζ ≥ 0 from our list of inequalities (5.38) as mentioned there, an essential difference between the set of inequalities in Eqs. (5.38) and (A.2) is only the existence of the last inequality in Eq. (A.2), i.e., ζ + ζ ⊥ + 2ζ × ≥ 0. This inequality can be, however, deduced from the other inequalities as follows. According to the inequality ζ ζ ⊥ ≥ ζ 2 × , one can immediately show that (ζ +ζ ⊥ ) 2 −(2ζ × ) 2 = (ζ −ζ ⊥ ) 2 +4(ζ ζ ⊥ −ζ 2 × ) ≥ 0. Since we have ζ ,⊥ ≥ 0, we can get rid of the square to find the inequality ζ + ζ ⊥ + 2ζ × ≥ 0, regardless of the sign of ζ × . This proof suggests that the fourth inequality in Eq. (A.2), and thus in Eq. (B.2) of HK [34], is a redundant inequality. Therefore, the minimal set of inequalities is given by Eq. (5.38) without either ζ ≥ 0 or ζ ⊥ ≥ 0. While all the transport coefficients in Ref. [33] are assumed to be semi-positive definite, this requirement is too strong; that is, not a necessary condition but sufficient one.
Comparison to the viscosities defined in Li and Yee (LY) [98] may be also of reader's interest. In this reference, the authors investigated hydrodynamics in an external nondynamical magnetic field and identified the hydrodynamic modes and relevant viscosities. Since the perpendicular components of momentum charges are not conserved quantities due to the Lorentz force [cf. the discussion below Eq. (1.1)], the number of hydrodynamic modes and of available tensor structures are reduced in the setup of Ref. [98]. Although we should stress the difference between the dynamical and nondynamical magnetic fields, the correspondences to the three viscous coefficients in Eq. (A5) of Ref. [98] are found to be η LY = η , ζ LY = ζ , ζ LY = ζ × , (A. 3) where the superscripts LY denote the viscosities introduced in Ref. [98]. It is clear that all of them quantify the responses to flow perturbations along the external magnetic field. The cross viscosity ζ × = ζ LY , however, contributes only to the transverse pressure, i.e., the transverse components of the energy-momentum tensor that are not conserved currents in the setup of Ref. [98]. Therefore, the cross viscosity may not be allowed in the strict hydrodynamic limit. Nevertheless, there could be a potential relevance of the cross viscosity if one considers the transverse momentum as a quasi-hydrodynamic mode. An entropycurrent analysis with a decent order-counting scheme for this specific setup is necessary for clarifying this point. From our point of view, if the cross viscosity contributes to the local entropy production at the first order, the cross bulk viscosity ζ × = ζ LY should vanish according to the inequality (5.37) with ζ ⊥ = 0 for the semi-positive definite entropy production.