Relativistic magnetohydrodynamics

We present the equations of relativistic hydrodynamics coupled to dynamical electromagnetic fields, including the effects of polarization, electric fields, and the derivative expansion. We enumerate the transport coefficients at leading order in derivatives, including electrical conductivities, viscosities, and thermodynamic coefficients. We find the constraints on transport coefficients due to the positivity of entropy production, and derive the corresponding Kubo formulas. For the neutral state in a magnetic field, small fluctuations include Alfven waves, magnetosonic waves, and the dissipative modes. For the state with a non-zero dynamical charge density in a magnetic field, plasma oscillations gap out all propagating modes, except for Alfven-like waves with a quadratic dispersion relation. We relate the transport coefficients in the"conventional"magnetohydrodynamics (formulated using Maxwell's equations in matter) to those in the"dual"version of magnetohydrodynamics (formulated using the conserved magnetic flux).


Introduction
In a macroscopic system, near-equilibrium phenomena can often be described by classical hydrodynamics. When the microscopic theory contains weakly coupled U(1) gauge fields, long-range correlations mediated by those fields are possible. Maxwell's equations in matter give an effective description of such correlations in terms of classical gauge fields. These equations are useful when the coupling between electromagnetic and thermal/mechanical degrees of freedom can be neglected. We would like to understand the effective description of relativistic systems in which macroscopic electromagnetic degrees of freedom are coupled to the macroscopic thermal and mechanical degrees of freedom. This amounts to coupling Maxwell's equations in matter to hydrodynamic equations. When the matter is electrically conducting and electric fields are neglected, such classical effective theory is usually called magneto-hydrodynamics (MHD). Our motivation it two-fold. From a fundamental point of view, a number of recent developments in relativistic hydrodynamics have pushed the boundaries of the "traditional" theory, as described for example in the classic textbook [1]. These include: a systematic derivative expansion in hydrodynamics [2], an equivalence between hydrodynamics and black hole dynamics [3], the manifestation of chiral anomalies in hydrodynamic equations [4], the relevance of partition functions [5,6], elucidation of the role of the entropy current [7,8], new insights into relativistic hydrodynamic turbulence [9], convergence properties of the hydrodynamic expansion [10], and a classification of hydrodynamic transport coefficients [11]. It is reasonable to expect that the above insights will also lead to an improved understanding of the "traditional" MHD. For example, there does not appear to be an agreement in the current literature on such basic question as the number of transport coefficients in MHD.
From an applied point of view, recent years have seen relativistic hydrodynamics expand from its traditional areas of astrophysical plasmas and hot subnuclear matter into the domain of condensed matter physics. Examples include transport near relativistic quantum critical points [12], in graphene [13,14] and in Weyl semi-metals [15]. For conducting matter, MHD is a natural extension of such hydrodynamic models.
In what follows, we will outline the construction of classical relativistic hydrodynamics with dynamical electromagnetic fields, starting from equilibrium thermodynamics. In order to write down the hydrodynamic equations, we will assume that the system is locally in thermal equilibrium. We will further assume that the departures from local equilibrium may be implemented through a derivative expansion such that the parameters which characterize the equilibrium (temperature, chemical potential, magnetic field, fluid velocity) vary slowly in space and time. At one-derivative order, transport coefficients such as viscosity and electrical conductivity appear in the constitutive relations. We are not aware of previous treatments that list all one-derivative terms in the constitutive relations of magnetohydrodynamics.
For parity-preserving conducting fluids in magnetic field, we find eleven transport coefficients at one-derivative order. One transport coefficient is thermodynamic, and determines the angular momentum of charged fluid induced by the magnetic field. Three transport coefficients are non-equilibrium and non-dissipative: these are the two Hall viscosities (transverse and longitudinal), and one Hall conductivity. There are also seven non-equilibrium dissipative transport coefficients: two electrical conductivities (transverse and longitudinal), two shear viscosities (transverse and longitudinal), and three bulk viscosities. The constitutive relations for the energy-momentum tensor are given in eqs. (3.1), (3.11), and for the current in eqs. (3.2), (3.12). The dissipative coefficients have to satisfy the inequalities in eq. (3.19) imposed by the positivity of entropy production, or alternatively by the positivity of the spectral function. As a simple application of the hydrodynamic equations, we study eigenmodes of small oscillations near thermal equilibrium in constant magnetic field.
We start in Section 2 with a discussion of equilibrium thermodynamics in the presence of external electromagnetic and gravitational fields. In Section 3, we will discuss hydrodynamics, again when electromagnetic and gravitational fields are external. The magnetic fields are taken as "large" and electric fields as "small" in the sense of the derivative expansion. The smallness of the electric field is due to electric screening. Our procedure will improve on existing studies by taking into account the effects of polarization (magnetic, electric, or both), electric fields, and by enumerating all transport coefficients at leading order in derivatives. In Section 4 we discuss hydrodynamics with dynamical electromagnetic fields, as an extension of hydrodynamics with fixed electromagnetic fields. As a simple example, one can study Alfvén and magnetosonic waves in a neutral state (including their damping and polarization), and waves in a dynamically charged (but overall electrically neutral) state. We compare our results with the recent "dual" formulation of MHD in Section 5, and with some of the previous studies of transport coefficients of relativistic fluids in magnetic field in the Appendix.

Thermodynamics
Let us start with equilibrium thermodynamics. For a system in equilibrium subject to an external non-dynamical gauge field A µ and an external non-dynamical metric g µν , we write the logarithm of the partition function W s = −i ln Z as and we will call F the free energy density. [Conventions: metric is mostly plus, ǫ 0123 =1/ √ −g.] For a system with short-range correlations in equilibrium and for external sources A and g which only vary on scales much longer than the correlation length, F is a local function of the external sources, and W s is extensive in the thermodynamic limit. The density F may then be written as an expansion in derivatives of the external sources [5,6]. The current J µ (defined by varying W s with respect to the gauge field) and the energy-momentum tensor T µν (defined by varying W s with respect to the metric) automatically satisfy is the generating functional of static (zero frequency) correlation functions of T µν and J µ in equilibrium. Of course, the conservation laws (2.2) are also true out of equilibrium, being a consequence of gauge-and diffeomorphism-invariance in the microscopic theory.
Being in equilibrium means that there exists a timelike Killing vector V such that the Lie derivative of the sources with respect to V vanishes. The equilibrium temperature T , velocity u α and the chemical potential µ are functions of the Killing vector and the external sources [5,6] Here β 0 is a constant setting the normalization of temperature, and Λ V is a gauge parameter which ensures that µ is gauge-invariant [16]. The electromagnetic field strength tensor where E µ ≡ F µν u ν is the electric field, and B µ ≡ 1 2 ǫ µναβ u ν F αβ is the magnetic field, satisfying u·E = u·B = 0. The decomposition (2.4) is just an identity, true for any antisymmetric F µν and any timelike unit u µ . Electric and magnetic fields are not independent, but are related by the "Bianchi identity" ǫ µναβ ∇ ν F αβ = 0, which in equilibrium becomes Here Ω µ ≡ ǫ µναβ u ν ∇ α u β is the vorticity and a µ ≡ u λ ∇ λ u µ is the acceleration. In equilibrium, the acceleration is related to temperature by ∂ λ T = −T a λ . Relations (2.5) are curved-space versions of the familiar flat-space equilibrium identities ∇·B = 0 and ∇×E = 0. In order to write down the density F in the derivative expansion, we need to specify the derivative counting of the external sources A and g. The natural derivative counting for the metric is g ∼ O(1) (assuming we are interested in transport phenomena in flat space), while the derivative counting for A depends on the physical system under consideration.
As an example, consider an insulator, such as a system made out of particles which carry electric/magnetic dipole moments, but no electric charges. In such a system, there is no conserved electric charge, and the above µ is not a relevant thermodynamic variable. If we are interested in thermodynamics of such a system subject to external electric and magnetic fields, we are free to choose B ∼ O(1) and E ∼ O(1) in the derivative expansion. The free energy density is then The leading-order term is the pressure, whose dependence on E and B encodes the electric, magnetic, and mixed susceptibilities. For the list of O(∂) contributions to F , see ref. [17]. As another example, consider a system that has electrically charged degrees of freedom (a conductor), such that µ gives a non-negligible contribution to thermodynamics. In equilibrium, ∂ λ µ = E λ − µa λ is satisfied identically, which suggests that counting µ ∼ O(1) leads to E ∼ O(∂). This is a manifestation of electric screening. The magnetic field, on the other hand, may still be counted as O(1). The counting B ∼ O(1) and E ∼ O(∂) is the relevant derivative counting for MHD. The free energy density is then where s (1) n are O(∂) gauge-and diffeomorphism-invariants, and the coefficients M n need to be determined by the microscopic theory, just like the pressure p. Following ref. [17], we list the invariants s (1) n in Table 1. The rows labeled C, P, T indicate the eigenvalue of the invariant under charge conjugation, parity, and time reversal. The last row shows the weight w of the invariant under a local rescaling of the metric: g µν →g µν = e −2ϕ g µν , and s n →s n = e wϕ s n . The invariant s (1) 3 does not transform homogeneously under the rescaling, and can not appear in a conformally invariant generating functional. Hence, we expect that in a conformal theory M 3 = 0. The coefficient M 5 is the usual magneto-electric (or electro-magnetic) susceptibility; similarly M 4 may be termed magneto-vortical susceptibility. For the rest of the paper, we will adopt the derivative counting B ∼ O(1) and E ∼ O(∂), as is appropriate for MHD.
As an example, consider a parity-invariant theory in magnetic field. The only O(∂) thermodynamic coefficient is the magneto-vortical susceptibility M Ω ≡ M 4 , which affects T µν and J µ when there is non-zero vorticity, and higher-point equilibrium correlation functions of T µν and J µ when there is no vorticity. We define static (zero frequency) correlation functions of T µν and J µ by varying the generating functional (2.1) with respect to g µν and A µ in the standard fashion. For example, in flat space at constant temperature T 0 , constant chemical potential µ 0 , and constant magnetic field B 0 in the z-direction, one finds the following static correlation functions at small momentum The first expression may be used to evaluate the magneto-vortical susceptibility M Ω in a system that is not subject to magnetic field, and is not rotating.
3 Hydrodynamics with external electromagnetic fields

Constitutive relations
Hydrodynamics is conventionally formulated as an extension of thermodynamics, in the sense that hydrodynamic variables are inherited from the thermodynamic parameters. This is a strong assumption, and we expect the hydrodynamic description only to be valid for B ≪ T 2 , otherwise new non-hydrodynamic degrees of freedom (such as those associated with Landau levels) must be taken into account. Let us start by taking E and B fields as external and nondynamical. In hydrodynamics, the thermodynamic variables T , u α , and µ are promoted to time-dependent quantities. Out of equilibrium, they no longer have a microscopic definition, but are merely auxiliary variables used to build the non-equilibrium energy-momentum tensor and the current. The expressions of T µν and J µ in terms of the auxiliary variables T , u α , and µ are called constitutive relations; they contain both thermodynamic contributions (coming from the variation of F ), and non-equilibrium contributions (such as the viscosity). It is worth noting that thermodynamic contributions and non-equilibrium contributions to the constitutive relations may appear at the same order in the derivative expansion. The constitutive relations are then used together with the conservation laws (2.2) to find the energy-momentum tensor and the current. While in thermodynamics Eqs. (2.2) are mere identities reflecting the symmetries of W s , solving Eqs. (2.2) in hydrodynamics can be a challenging endeavour leading to rich physics. We will write the energy-momentum tensor using the decomposition with respect to the timelike velocity vector u µ , where ∆ µν ≡ g µν + u µ u ν is the transverse projector, Q µ is transverse to u µ , and T µν is transverse to u µ , symmetric, and traceless. Explicitly, the coefficients are E ≡ u µ u ν T µν , Similarly, we will write the current as where the charge density is N ≡ −u µ J µ , and the spatial current is J µ ≡ ∆ µλ J λ . Using the equilibrium free energy (2.7), one can isolate O(1) and O(∂) contributions to the energy-momentum tensor and the current: where ǫ = −p + T (∂p/∂T ) + µ(∂p/∂µ), Π = p − 2 3 α BB B 2 , n = ∂p/∂µ, and the magnetic susceptibility

Field redefinitions
Out of equilibrium, the variables T , u α , and µ may be redefined. Such a redefinition is often referred to as a choice of "frame", see e.g. ref. [18] for a discussion. Consider changing the hydrodynamic variables to T ′ = T + δT , u ′α = u α + δu α , µ ′ = µ + δµ, where δT , δu α , and δµ are O(∂). The same energy-momentum tensor and the current may be expressed either in terms of T , u α , µ, or in terms of T ′ , u ′α , µ ′ (note that B 2 = B ′2 + O(∂ 2 )). Physical transport coefficients must be derived from O(∂) quantities which are invariant under such changes of hydrodynamic variables. A direct evaluation shows that the following combinations are invariant under "frame" transformations: Here B µν ≡ ∆ µν − B µ B ν /B 2 is the projector onto a plane orthogonal to both u µ and B µ , all thermodynamic derivatives are evaluated at fixed B 2 , and B ≡ √ B 2 . When the magnetic susceptibility α BB is T -and µ-independent, the stress f µν T is frame-invariant. As an example, one can choose δT and δµ such that E ′ = ǫ(T ′ , µ ′ , B ′2 ), N ′ = n(T ′ , µ ′ , B ′2 ), and further choose δu α such that Q ′ α = 0. This corresponds to the Landau-Lifshitz frame [1]. The components of energy-momentum tensor and the current take the following form in the Landau-Lifshitz frame: where the frame invariants are given by eq. (3.3). In the Landau-Lifshitz frame, a non-zero value of the pseudoscalar frame-invariant ℓ indicates a current flowing along the magnetic field. In a constant external magnetic field such currents arise as consequences of chiral anomalies [4]; in an inhomogeneous external field, an electric current flowing along the magnetic field can arise without chiral anomalies, owing to a non-zero magnetic susceptibility.

Thermodynamic frame
The energy-momentum tensor and the current derived from the static generating functional W s correspond to a different frame, termed in [6] the thermodynamic frame. Taking the variation of the free energy (2.7), one finds the following equilibrium O(∂) contributions in the thermodynamic frame: where the bar signifies equilibrium contributions, and the coefficients ǫ n , π n , φ n , γ n , δ n , θ n are all O(1) functions of the five thermodynamic coefficients M n (T, µ, B 2 ) and of the magnetic susceptibility α BB = 2∂p/∂B 2 . The explicit expressions are given in Appendix A. The onederivative scalars s (1) n are given in  is the Poynting vector. Bottom: Non-zero symmetric transverse traceless O(∂) tensors that appear in the equilibrium stress T µν . For any two transverse vectors X µ and Y µ , the angular brackets stand for In the vector invariant, we have defined v (1)µ 5 ≡ s (1) 2 B µ . The subscript "non-eq" denotes non-equilibrium contributions which by definition vanish in equilibrium. The functions Φ n (T, µ, B 2 ), Λ n (T, µ, B 2 ), Γ n (T, µ, B 2 ), Θ n (T, µ, B 2 ) are non-dissipative thermodynamic transport coefficients. Explicitly, We see that the constitutive relations for energy-momentum tensor and the current contain twenty-one thermodynamic transport coefficients Φ n , Λ 2 , Γ n , Θ n . These twenty-one coefficients are not independent, but can all be expressed in terms of only five parameters M n of the equilibrium generating functional.
Let us now write down the constitutive relations in the thermodynamic frame that is a natural generalization of the Landau-Lifshitz frame. We will define the thermodynamic frame (primed variables) by redefinitions of T , µ, and u α that give In other words, in this thermodynamic frame the coefficients E, N , and Q α in the decompositions (3.1), (3.2) take their equilibrium values, derived from the equilibrium generating functional W s . The other coefficients take the following form in the thermodynamic frame:

Non-equilibrium contributions
With the equilibrium contributions out of the way, the next task is to find the non-equilibrium terms in the constitutive relations (3.6). This amounts to finding one-derivative scalars, vectors (orthogonal both to B µ and to u µ ), and transverse traceless symmetric tensors that vanish in equilibrium. Note that non-equilibrium contributions (those that vanish in equilibrium) are not the same as dissipative contributions (those that contribute to hydrodynamic entropy production). Every dissipative contribution is non-equilibrium, but not every nonequilibrium contribution is dissipative. The six independent non-equilibrium one-derivative scalars are given in Table 3. The scalar u λ ∂ λ B 2 is not independent as a consequence of the electromagnetic Bianchi identity, and can be expressed as a combination of ∇·u and B µ B ν ∇ µ u ν . Three scalar equations of motion ∇ µ J µ = 0, u ν ∇ µ T µν + E µ J µ = 0, and B ν ∇ µ T µν + (E·B)(u·J) = 0 taken at zeroth order provide three relations among the scalars. We choose to eliminate s (1) 1 non-eq. , s (1) 2 non-eq. , and s (1) 6 non-eq. and write the scalar and pseudo-scalar constitutive relations as f non-eq. = c 1 s (1) 3 non-eq. + c 2 s (1) 4 non-eq. + c 3 s (1) 5 non-eq. , ℓ non-eq. = c 4 s (1) 3 non-eq. + c 5 s (1) 4 non-eq. + c 6 s (1) 5 non-eq. , with some undetermined transport coefficients c n .
The independent non-equilibrium transverse one-derivative vectors are given in Table 3, where the shear tensor is We use the vector equation of motion (2.2a) projected with B µν at zeroth order to eliminate one of the vectors,   Table 3. Non-equilibrium scalars and transverse non-equilibrium vectors at O(∂), written in terms of b µ ≡ B µ /B. In addition to the vectors listed in the table, there are corresponding transverse nonequilibrium vectorsṽ (1)µ non-eq. ≡ ǫ µνρσ u ν b ρ v (1) non-eq. σ . The table also shows the parity of non-equilibrium scalars and vectors. Under time-reversal, the scalars s (1) n non-eq. are T-odd, the vectors v (1)µ n non-eq. are T-even, and the vectorsṽ (1)µ n non-eq. are T-odd.
and write the vector constitutive relation as The tilded vectors are defined asṽ There is a number of symmetric transverse traceless non-equilibrium one-derivative tensors besides the shear tensor σ µν . One such tensor is Other tensors can be formed by B µ B ν s (1) n non-eq. , or by symmetrizing B µ with a transverse non-equilibrium vector. Again, we eliminate three scalars and one vector by the zeroth order equations of motion and write the tensor constitutive relation in terms of b µ ≡ B µ /B as with some undetermined transport coefficients c n . Thus there are five equilibrium functions M n (T, µ, B 2 ), and nineteen non-equilibrium functions c n (T, µ, B 2 ) that determine onederivative contributions to the energy-momentum tensor and the current in strong magnetic field. If the microscopic system is parity-invariant, all thermodynamic coefficients M n vanish except for M 4 . In addition, the dynamical coefficients c 3 , c 4 , c 5 , c 8 , c 10 , c 14 , c 15 , c 17 must vanish by parity invariance. Thus a conducting parity-invariant system in magnetic field has one thermodynamic coefficient M 4 , three "electrical conductivities" c 6 , c 7 , and c 9 , and eight "viscosities" c 1 , c 2 , c 11 , c 12 , c 13 , c 16 , c 18 , and c 19 . We will see later that the Onsager relations impose a relation between c 2 , c 12 , and c 13 , plus four more relations among the parity-violating coefficients. This leaves eleven transport coefficients (one thermodynamic and ten non-equilibrium) for a conducting parity-invariant system in magnetic field in 3+1 dimensions. In a conformal theory, the tracelessness condition 2 will in addition impose The constitutive relations may be simplified further if we note that the shear tensor can be decomposed with respect to the magnetic field as Here σ µν whereσ µν ⊥ is transverse to both u µ and B µ , symmetric, and traceless. For completeness, let us summarize the constitutive relations for a parity-invariant theory in the thermodynamic frame. Defining M Ω ≡ M 4 , the energy-momentum tensor is given by eq. (3.1) with the following coefficients: In a conformal theory subject to external fields g µν and A µ , the trace of the energy-momentum tensor receives an anomalous contribution T µ µ = κF 2 + O(∂ 4 ), where κ is a theory-dependent constant that counts the number of charged degrees of freedom, and the terms O(∂ 4 ) are due to curvature invariants. It was shown in ref. [19] that the conformal anomaly may be captured by a certain local term in the hydrostatic generating functional, which for our purposes amounts to a term in p(T, µ, B 2 ) proportional to κ. and the current is given by eq. (3.2) with the following coefficients: The current is written in terms of the magnetic polarization vector while the electric polarization vector vanishes at leading order in a parity-invariant system. The comma subscript denotes the derivative with respect to the argument that follows. Note that we are keeping O(∂ 2 ) thermodynamic terms in the constitutive relations (coming from the variation of M 4 s (1) 4 ) that are needed to ensure that the conservation laws (2.2) are satisfied identically for time-independent background fields. In writing down the constitutive relations (3.11), (3.12), we have relabeled the non-equilibrium transport coefficients as The coefficients σ ⊥ , σ are the transverse and longitudinal conductivities, and η ⊥ , η are the transverse and longitudinal shear viscosities. The coefficients ζ 1 , ζ 2 , η 1 and η 2 may all be called "bulk viscosities", of which only three are independent due to the Onsager relation. The coefficientsη ⊥ ,η are the two Hall viscosities, andσ is the Hall conductivity. 3 When the external electromagnetic field vanishes, the system becomes isotropic, and we expect to recover the constitutive relations of the standard isotropic hydrodynamics, with shear viscosity η, bulk viscosity ζ, and electrical conductivity σ. Thus as B → 0 we expect

Eigenmodes
As a simple application of the hydrodynamic equations (2.2) together with the constitutive relations (3.11), (3.12), one can study the eigenmodes of small oscillations about the thermal equilibrium state. We set the external sources to zero, and linearize the hydrodynamic equations near the flat-space equilibrium state with constant T = T 0 , µ = µ 0 , u α = (1, 0), and B α = (0, 0, 0, B 0 ). Taking the fluctuating hydrodynamic variables proportional to exp(−iωt+ ik·x), the source-free system admits five eigenmodes, two gapped (ω(k→0) = 0), and three gapless (ω(k→0) = 0). The frequencies of the gapped eigenmodes are where w 0 ≡ ǫ 0 + p 0 is the equilibrium enthalpy density, and we have taken As the imaginary part of the eigenfrequency must be negative for stability, this implies σ ⊥ > 0. The mode has a circular polarization (at k = 0), with δu x and δu y oscillating with a π/2 phase difference. The analogous mode in 2+1 dimensional hydrodynamics was christened the hydrodynamic cyclotron mode in ref. [12], which also explored its implications for transport near two-dimensional quantum critical points.
For momenta k B 0 , the three gapless eigenmodes are the two sound waves, and one diffusive mode. The eigenfrequencies in the small momentum limit are where v s is the speed of sound. As in ref. [18], we can write the coefficients in terms of the elements of the susceptibility matrix in the grand canonical ensemble. The non-zero elements of the 3 × 3 susceptibility matrix are χ 11 = T (∂ǫ/∂T ) µ/T , χ 13 = χ 31 = (∂ǫ/∂µ) T , The positivity of the diffusion constant implies σ > 0. The speed of sound squared expressed in terms of the elements of the susceptibility matrix is given by and the damping coefficient is The expression for v s and D in terms of the thermodynamic functions formally look the same as in hydrodynamics without external O(1) magnetic fields [18]. All of v s , Γ s, , and D depend on B 0 through p = p(T, µ, B 2 ) and the transport coefficients. For momenta k ⊥ B 0 , the three gapless eigenmodes include two diffusive modes, and one "subdiffusive" mode with a quartic dispersion relation, The transverse diffusion constant is determined by the transverse resistivity. We define the 2 × 2 conductivity matrix in the plane transverse to B 0 as σ ab ≡ σ ⊥ δ ab + n 0 |B 0 | +σ ǫ ab , and the corresponding resistivity matrix as ρ ab ≡ (σ −1 ) ab = ρ ⊥ δ ab +ρ ⊥ ǫ ab , which defines ρ ⊥ andρ ⊥ . The transverse diffusion constant is then Stability of the equilibrium state now implies η ⊥ > 0, η > 0. For modes propagating at an angle θ with respect to B 0 , the gapless modes include sound waves (unless θ = π/2), and a diffusive mode. For a fixed value of θ, the small-momentum The coefficient D c in the cyclotron mode eigenfrequency (3.14) at small B 0 is Note that the limits θ → π/2 and k → 0 in the eigenfrequencies do not commute.

Entropy production
The simple flat-space eigenfrequency analysis in the previous subsection imposes certain constraints on non-equilibrium transport coefficients. In order to find more general constraints, one method is to impose a local version of the second law of thermodynamics: the existence of a local entropy current with positive semi-definite divergence for every non-equilibrium configuration consistent with the hydrodynamic equations. We will not attempt to construct the most general entropy current from scratch. Rather, we will use the result of [7,8] saying that the constraints on transport coefficients derived from the entropy current are the same as those derived from the equilibrium generating functional, plus the inequality constraints on dissipative transport coefficients. We take the entropy current to be where the canonical part of the entropy current is and S µ eq. is found from the equilibrium partition function, as described in [7,8]. The constraints on transport coefficients follow by demanding ∇ µ S µ 0. Using conservation laws (2.2), the divergence of the canonical entropy current is The S µ eq. part of the entropy current is explicitly built to cancel out the part of ∇ µ S µ canon that arises from the equilibrium terms in the constitutive relations, i.e. the terms in T µν and J µ derived from the equilibrium generating functional. In fact, ref. [8] has already found S µ eq. in the case when the generating functional contains a contribution proportional to B·Ω. We thus focus on non-equilibrium terms, and write the thermodynamic frame constitutive relations (3.7) as T µν = T µν eq. + T µν non-eq. and J µ = J µ eq. + J µ non-eq. . The divergence of the entropy current is then Using the constitutive relations (3.11), (3.12), this leads to where again S 3 ≡ ∇·u and together with the condition that the quadratic form made out of S 3 , S 4 in the second line of eq. (3.18) is non-negative, which implies The coefficientsη ⊥ ,η , andσ do not contribute to entropy production, and are not constrained by the above analysis. Thus,η ⊥ ,η , andσ are non-equilibrium non-dissipative coefficients.

Kubo formulas
When the microscopic system is time-reversal invariant (i.e. the only source of time-reversal breaking is due to the external magnetic field), transport coefficients can be further constrained by the Onsager relations. The retarded two-point functions of operators O a and O b in a time-reversal invariant theory in equilibrium obey where ǫ a and ǫ b are time-reversal eigenvalues of the operators O a and O b . We take our operators to be various components of T µν and J µ , and evaluate the retarded two-point functions by varying one-point functions in the presence of the external source with respect to the source. Namely, we solve the hydrodynamic equations in the presence of fluctuating external sources δA, δg (proportional to exp(−iωt+ik·x)) to find δT [A, g], δµ[A, g], δu α [A, g], and then vary the resulting hydrodynamic expressions T µν [A, g] and J µ [A, g] with respect to g αβ , A α to find the retarded functions. Specifically, where the subscript "on-shell" signifies that the corresponding hydrodynamic T µν [A, g] and J µ [A, g] are evaluated on the solutions to (2.2), and the sources δA, δg are set to zero after the variation is taken. The expressions (3.21) are to be understood as etc. This provides a direct method to evaluate the retarded functions, and allows both to check the Onsager relations and to derive Kubo formulas for transport coefficients. 4 The constraint on transport coefficients we find by demanding that eq. (3.20) holds is 5 For the rest of the paper, we will assume that (3.22) holds, which leaves us with ten nonequilibrium transport coefficients for a parity-invariant microscopic system. Using eq. (3.22) to eliminate ζ 2 , the inequality constraint in eq. (3.19c) turns into We next list the expressions for transport coefficients in terms of retarded functions evaluated in flat-space equilibrium with external magnetic field in the z direction, as in sec. 3.5. In the limit k → 0 first, ω → 0 second we find the following Kubo formulas. The two-point function of the longitudinal current J z gives the longitudinal conductivity, while the two-point functions of the transverse currents J x , J y give the transverse resistivities, where the resistivities ρ ⊥ andρ ⊥ were defined below eq. (3.16). Alternatively, the resistivities can be found from correlation functions of momentum density, The shear viscosities are given by while the "bulk" viscosities may be expressed as . Correlation functions at nonzero momentum may be obtained in a straightforward way from the variational procedure described earlier.

Inequality constraints on transport coefficients
Finally, let us show that the inequality constraints on transport coefficients derived from demanding that the entropy production is non-negative can also be obtained from hydrodynamic correlation functions, without using the entropy current. The argument is based on the fact that the imaginary part of the retarded function G OO (ω, k) must be positive for any Hermitean operator O and ω > 0, Now consider the operator O = aO 1 + bO 2 , with real coefficients a and b, and Hermitean for ω 0. This quadratic form in a, b must be non-negative for all a, b which implies The two terms in the right-hand side of (3.28) can be related by the Onsager relation (3.20).

Hydrodynamics with dynamical electromagnetic fields 4.1 Dynamical gauge field
We now move on to systems where the gauge field A µ is dynamical rather than external, which will lead us to MHD. In external metric g, the (microscopic) generating functional is where S is the action. Let us couple the gauge field to an external conserved current J µ ext . We do this so that the new generating functional is and W ≡ −i ln Z. The new field ϕ is a Lagrange multiplier which shifts under gauge transformations and ensures that the external current is conserved. We define the energymomentum tensor and the current by the variation of the action: In what follows, we will omit the angular brackets, writing the (non)-conservation of the energy-momentum tensor simply as In the standard hydrodynamic approach, T µν and F µν will then be taken as dynamical variables in the classical hydrodynamic theory. Note that the sign in the right-hand side of eq. (4.2) is opposite compared to eq. (2.2a), owing to the fact that the current, rather than the gauge field, is now external. In order to proceed with hydrodynamics, we need to specify a) the constitutive relations for the energy-momentum tensor to be used in eq. (4.2), and b) the equations which determine the evolution of the dynamical gauge field F µν .

Maxwell's equations in matter
Classical equations specifying the dynamics of electric and magnetic fields are usually referred to as Maxwell's equations in matter. While we don't have a recipe of deriving them in a most general form in a model-independent way, a useful starting point is provided by matter in thermal equilibrium. Maxwell's equations for equilibrium matter may be then amended to include the non-equilibrium and dissipative effects, such as the electrical conductivity. To this end, as advocated in [20], we take the static generating functional W s [g, A] to be the effective action for gauge fields in equilibrium, where F is a local gauge-invariant function of the sources g µν and A µ , and we have ignored the surface terms. To leading order in the derivative expansion, F is simply the pressure. We can always write F = − 1 4 F µν F µν + F m , where the vacuum action is − 1 4 F µν F µν = 1 2 (E 2 − B 2 ), and F m is the "matter" contribution. The isolation of the vacuum term is arbitrary, but it will allow us to make contact with the textbook form of Maxwell's equations in matter. Our (equilibrium) effective theory is then given by the partition function (4.1), with S replaced by S eff , and the total action is The current derived by varying the total action with respect to A µ is J µ tot = J µ + J µ ext , or where the polarization tensor M µν m is defined by δ F d 4 x √ −g F m = 1 2 d 4 x √ −g M µν m δF µν , and the density of "free" charges is n ≡ ∂F m /∂µ. The equation of motion for the gauge field follows from δ A S tot = 0, or equivalently J µ tot = 0, and becomes where H µν ≡ F µν − M µν m . This is the desired equation that must be satisfied by electromagnetic fields in equilibrium. Following the standard hydrodynamic lore and assuming that eq. (4.4) also holds for small departures away from equilibrium, one obtains hydrodynamics of "perfect fluids", now with dynamical electric and magnetic fields. For these perfect fluids, equations (4.4) have to be solved together with the stress tensor (non)-conservation (4.2), where T µν is derived from the effective action (4.3).
In fact, eq. (4.4) is nothing but the standard Maxwell's equations in matter. The polarization tensor M µν m defines electric and magnetic polarization vectors P µ and M µ through the decomposition The antisymmetric tensor H µν can be decomposed in the same way as the field strength F µν , which defines D µ ≡ H µν u ν and H µ ≡ 1 2 ǫ µναβ u ν H αβ , so that It is then clear that eq.

Hydrodynamics
We take the MHD equations to be as follows: The last equation is the electromagnetic "Bianchi identity", expressing the fact that the electric and magnetic fields are derived from the vector potential A µ . The second equation (Maxwell's equations in matter) can be rewritten as ∇ ν (F µν −M µν m ) = J µ free +J µ ext which defines J µ free , the current of "free charges". While eqs. (4.7a) and (4.7c) are true microscopically, the Maxwell's equations in matter (4.7b) are written based on the above intuition of the equilibrium effective action. Note that ∇ µ J µ free = 0 is a consequence of (4.7b), and is not an independent equation. The hydrodynamic variables are T , u α , µ, as well as the electric and magnetic fields which satisfy u α E α = 0, u α B α = 0. Hydrodynamic equations (4.7) must be supplemented by constitutive relations, which express T µν , J µ (or J µ free and M µν m ) in terms of the hydrodynamic variables. These constitutive relations will contain equilibrium contributions coming from the equilibrium effective action (4.3). In addition, the constitutive relations will contain non-equilibrium contributions, such as the electrical conductivity and the shear viscosity.
Taking the divergence of eq. (4.7b) and using J µ ext = −J µ gives which shows that the variables T , u α , and µ satisfy exactly the same equations (2.2) as they did in the theory with a non-dynamical, external A µ . Thus in order to "solve" the MHD theory (4.7) one can i) solve the hydrodynamic equations with an external gauge field (4.7) to find . MHD correlation functions may then be obtained through variations with respect to the external sources J λ ext and g µν . An equivalent way to understand the classical effective theory (4.7) is to promote the real-time generating functional to the non-equilibrium effective action [20], i.e. to write where W r [A, g] is low-energy, real-time generating functional for retarded correlation functions in the theory with a non-dynamical A µ . The functional W r [g, A] is non-local due to the gapless low-energy degrees of freedom (sound waves etc). However, for the purposes of MHD we do not need the actual generating functional, but only the equations of motion for the effective action S tot . These equations of motion are J µ [A, g] + J µ ext = 0, where J µ [A, g] is the on-shell current in the theory with a non-dynamical A µ . One can then solve the theory as described in the previous paragraph.
We will thus adopt the simplest hydrodynamic effective theory (4.7) where the constitutive relations for T µν and J µ are the same as in the case of external non-dynamical electromagnetic fields. Under this "mean-field" assumption, transport coefficients which are naively independent would still be related by the conditions originating from the static generating functional.
Further, any solution T [A, g], u α [A, g], µ[A, g] to the MHD equations is also a solution to the hydrodynamic equations (2.2) in the theory with a non-dynamical A µ . Thus the entropy current with a non-negative divergence on the solutions to (2.2) will also have non-negative divergence when evaluated on the solutions to the MHD equations (4.7). This means that the entropy current in MHD may be taken the same as the entropy current in the theory with a non-dynamical gauge field [20], and we do not need to perform a separate entropy current analysis beyond what was already done in sec. 3.
To sum up, with the MHD scaling B ∼ O(1), E ∼ O(∂), the equilibrium effective action is given by eq. (2.7), For a parity-invariant theory, only the M 4 term in the sum contributes. The constitutive relations for the energy-momentum tensor and the current were already found in the previous section, where now we have p(T, µ, B 2 ) = − 1 2 B 2 +p m (T, µ, B 2 ). The energy-momentum tensor appearing in eq. (4.7) and the current J µ satisfying J µ + J µ ext = 0 take the form (3.1), (3.2), and the constitutive relations for a parity-invariant theory in the thermodynamic frame are given by Eqs. (3.11), (3.12).
We will find it useful to modify the above effective theory by giving dynamics to the electric field. To do so, we add an O(∂ 2 ) term 1 2 ε e E 2 to the effective action (4.8), where ε e is the electric permittivity which we take constant. This term is one of the many O(∂ 2 ) terms, and we add it as a "ultraviolet regulator" which improves the high-frequency behaviour of the theory. When studying the near-equilibrium eigenmodes of the system, this term will affect the frequency gaps, but not the leading-order dispersion relations of the gapless modes. With this new term, the following contributions have to be added to the constitutive relations (3.11), (3.12): The current J µ El. contains the kinetic term for the electric field in Maxwell's equations, as well as the "bound" current due to electric polarization.

Eigenmodes
As a simple application of the above MHD theory, one can study the eigenmodes of small oscillations about the thermal equilibrium state. As we did earlier, we set the external sources to zero, and linearize the hydrodynamic equations near the flat-space equilibrium state with constant T = T 0 , µ = µ 0 , u α = (1, 0), and B α = (0, 0, 0, B 0 ). For simplicity, we will take the magnetic permeability µ m constant, though it is straightforward to find how the eigenfrequencies below are modified for non-constant µ m = µ m (T, µ, B 2 ).

Neutral state
We begin with the neutral state at µ 0 = 0 and n 0 = 0. The system admits nine eigenmodes, three gapped, and six gapless.
Let us start with the familiar case of vanishing magnetic field in equilibrium. The system is then isotropic, with shear viscosity η, bulk viscosity ζ, and conductivity σ ≡ σ ⊥ = σ . The fluctuations of δT , δu i decouple from the fluctuations of δµ, δE i , δB i . The eigenmodes include two transverse shear modes with eigenfrequency ω = −iηk 2 /(ǫ 0 +p 0 ), and longitudinal sound waves with v 2 s = ∂p/∂ǫ and Γ s = ( 4 3 η + ζ)/(ǫ 0 + p 0 ). In addition, there is a longitudinal charge diffusion mode which becomes gapped because of non-zero electrical conductivity, Thus, charge fluctuations in a neutral conducting medium do not diffuse. Instead, what diffuses are the transverse magnetic and electric fields: there are two sets of transverse conductor modes whose eigenfrequencies are determined by Recall that ε e is the electric permittivity and µ m = 1/(1−2∂p m /∂B 2 ) is the magnetic permeability, so √ ε e µ m is the elementary index of refraction. The conductor modes have the following frequencies at small momenta: The gapless conductor mode is responsible for the skin effect in metals. We now turn on non-zero magnetic field and consider modes propagating at an angle θ with respect to B 0 . Thermal and mechanical fluctuations now no longer decouple from electromagnetic fluctuations. There is one longitudinal gapped mode, and two transverse gapped modes, In writing down the transverse eigenfrequencies, we have assumed B 2 0 ≪ ǫ 0 + p 0 . All six gapless modes have linear dispersion relation at small momenta. Two of the gapless modes are the Alfvén waves, whose speed and damping are determined by where ρ ≡ 1/σ , and ρ ⊥ was defined below eq. (3.16). In writing down the damping coefficient, we have taken B 2 0 ≪ ǫ 0 +p 0 , the corrections of order B 2 0 /(ǫ 0 +p 0 ) are straightforward to write down. The other four gapless modes are the two branches of magnetosonic waves, ω = ±v ms k − iΓ ms 2 k 2 , (4.10a) whose speed is determined by the quadratic equation where v 2 s = (s/T )/(∂s/∂T ) = ∂p/∂ǫ is the speed of sound at n 0 = 0. The two solutions of (4.10b) correspond to the sound-type (or "fast") branch, and the Alfvén-type (or "slow") branch. At θ = 0, the slow branch turns into a second set of Alfvén waves, while the fast branch becomes the sound wave. See e.g. ref. [21] for an early derivation of v A and v ms in relativistic MHD. The damping coefficients of the magnetosonic waves are straightforward to evaluate, but are quite lengthy to write down in general, and we will only present them in the limits of small B 0 and small θ. As B 0 → 0, the damping coefficients become slow: Γ ms = η ǫ 0 +p 0 + 1 σµ m , (4.10c) fast: Γ ms = 1 ǫ 0 +p 0 4 3 η + ζ . (4.10d) On the other hand, as θ → 0, the damping coefficients become fast: Γ ms = 1 ǫ 0 +p 0 10 3 η 1 + 2η 2 + ζ 1 . (4.10f) We have again taken B 2 0 ≪ ǫ 0 + p 0 , the corrections of order B 2 0 /(ǫ 0 +p 0 ) are straightforward to write down. At θ = 0, both polarizations of Alfvén waves have the same damping.

Charged state offset by background charge
We now consider a state with a non-zero value of µ 0 , which gives rise to a constant non-zero charge density n 0 . In order to ensure that the equilibrium state is stable, we will offset this equilibrium value of the dynamical charge density by a constant non-dynamical external background charge density −n 0 . This can be achieved by choosing the external current in the hydrodynamic equations (4.7) as J µ ext = (−n 0 , 0). In the particle language, this would correspond to a state where the excess of electrically charged particles over antiparticles (or vice versa) is compensated by a constant charge density of immobile background "ions". Even though the system is overall electrically neutral, its dynamics is not equivalent to that of the system with µ 0 = 0, n 0 = 0: for example, the fluctuation of the spatial electric current has a convective contribution n 0 δu i . More formally, when analyzing hydrodynamic modes, the limits n 0 → 0 and k → 0 do not commute. We now find six gapped modes and three gapless modes.
To get some intuition about the gapped modes, let us set all transport coefficients to zero, as well as set B 0 = 0. Then at small momenta there are two longitudinal gapped modes whose frequencies are determined by where Ω 2 p ≡ n 2 0 /[(ǫ 0 +p 0 )ε e ], and v s is the speed of sound that the charged fluid would have, if the electromagnetic fields were not dynamical, see Sec 3.5. These modes are the relativistic analogues of Langmuir oscillations, and Ω p is the relativistic "plasma frequency" which gaps out the sound waves. In addition, there are four transverse gapped modes whose frequencies are determined by These are electromagnetic waves in the fluid, gapped by the same plasma frequency Ω p as the sound waves. If we now turn on the transport coefficients, the gaps are determined by indicating the damping of plasma oscillations. At non-zero B 2 0 ≪ ǫ 0 +p 0 , the gaps will receive dependence on the magnetic field.
At B 0 = 0 the system is isotropic. The gapless modes (B 0 → 0 first, k → 0 second) include two transverse shear modes with quartic dispersion relation, and one longitudinal diffusive mode, where again w 0 ≡ T 0 s 0 + µ 0 n 0 , and the susceptibility matrix χ was defined below eq. (3.15). At non-zero B 0 , the three gapless modes all have quadratic dispersion relation at small momenta. There are two propagating waves with real frequencies where θ is the angle between k and B 0 , and one diffusive mode. For B 2 0 M Ω,µ ≪ ǫ 0 + p 0 , the diffusive frequency is (4.14) For gapless modes propagating at θ = π/2 at small momenta (θ → π/2 first, k → 0 second), we again find the diffusive mode ω = −iD ⊥ k 2 , with the same coefficient D ⊥ as in sec. 3.5.
In addition, at θ = π/2 there are two "subdiffusive" modes with quartic dispersion relation, The eigenfrequencies are noticeably different from the ones in a theory with fixed, nondynamical electromagnetic field discussed in sec. 3.5. Compared to the case of n 0 = 0 earlier in this section, one can say that non-vanishing dynamical charge density gaps out the magnetosonic waves, and turns Alfvén waves into waves whose frequency is quadratic in momentum.

Kubo formulas
We can find MHD correlation functions following the same variational procedure outlined in sec. 3.7. As the total current vanishes by the equations of motion, the objects whose correlation functions it makes sense to evaluate in MHD are the energy-momentum tensor T µν and the electromagnetic field strength tensor F µν . It is straightforward to evaluate retarded functions in flat space, in an equilibrium state with constant T = T 0 , µ = µ 0 , u α = (1, 0), and constant magnetic field. We solve the hydrodynamic equations in the presence of fluctuating external sources δJ ext , δg (proportional to exp(−iωt + ik·x)) to find δT [J ext , g], δµ[J ext , g], δu α [J ext , g], δF µν [J ext , g] and then vary the resulting hydrodynamic expressions T µν [J ext , g] and F µν [J ext , g] with respect to g αβ , J α ext to find the retarded functions. The metric variations are performed as usual, The subscript "on-shell" signifies that T µν and F µν are evaluated on the solutions to (4.7) with the constitutive relations (3.11), (3.12). Further, recall that the external current must be conserved, which can be implemented by choosing δJ 0 ext = k i δJ i ext /ω + 1 2 n 0 δg µ µ . The coupling A µ J µ ext then implies that iω δ/δJ l ext (k) produces an insertion of F 0l (−k), while ik m ǫ nml δ/δJ l ext (k) produces an insertion of 1 2 ǫ nml F lm (−k). For example, for electric field correlation functions we have and similarly for the magnetic field. 6 6 Alternatively, one can introduce an antisymmetric "polarization source" M µν ext , by taking the conserved current as J µ ext = ∇ ν M µν ext . The coupling A µ J µ ext then becomes 1 2 M µν ext F µν upon integration by parts, and correlation functions of F µν may be obtained as variations with respect to M µν ext .
Choosing the external magnetic field in the z-direction, we find the same Kubo formulas (3.25) and (3.26). The electrical resistivities may also be expressed in terms of correlation functions of the electric field. In the zero-density state with µ 0 = 0, n 0 = 0 we find 1 ω Im G F z0 F z0 (ω, k=0) = ρ , (4.15a) at small frequency, where ρ ≡ 1/σ . Similarly, for the transverse resistivities we find where again w 0 ≡ ǫ 0 +p 0 , and ρ ⊥ ,ρ ⊥ were defined below eq. (3.16). We have taken B 2 0 ≪ w 0 , otherwise there is a multiplicative factor of w 0 (w 0 −B 2 0 M Ω,µ )µ 2 m /(w 0 µ m +B 2 0 ) 2 in the righthand side of (4.15b), (4.15c). In a charged state (offset by non-dynamical −n 0 ), the correlation functions change, for example G F x0 F y0 (ω, k=0) = iω B 0 n 0 , while σ can be found from Retarded functions at non-zero momentum may be found from the above variational procedure. For example, the function G F x0 F x0 (ω, k) in a state with n 0 = 0 and with k B 0 has singularities at the eigenfrequencies of Alfvén waves for small momenta.

A dual formulation
As this paper was being completed, an interesting article [22] (abbreviated below as GHI) came out which approached magnetohydrodynamics from a different perspective. The dual electromagnetic field strength tensor J µν ≡ 1 2 ǫ µναβ F αβ was taken as a conserved current, and the constitutive relations were written down for J µν , rather than for the electric current J µ as was done in MHD historically. This "dual" construction follows the earlier work of ref. [23] which studied a similar MHD-like setup for "string fluids". The paper [22] identifies six transport coefficients in MHD, compared to eleven transport coefficients (in a paritypreserving system) found here. In this section we revisit the analysis of GHI, and show that the dual formulation allows for the same eleven transport coefficients we described earlier in Sections 3 and 4.

Constitutive relations
The conservation laws are taken as follows: These are the same equations (4.7a), (4.7c) we had earlier. The conserved external current is taken as J µ ext = 1 2 ǫ µνρσ ∂ ν Π ext ρσ , where Π ext µν may be viewed as the dual of the external polarization tensor M µν ext . The coupling A µ J µ ext then becomes 1 2 Π ext µν J µν upon integration by parts, and correlation functions of J µν may be obtained as variations with respect to Π ext µν . The tensor H in (5.1) is H = 1 2 dΠ ext , or in components H αβγ = 1 4 ∂ α Π ext βγ + (signed permutations).
In order to relate the GHI thermodynamic parameters to ours, we can compare equilibrium currents. The currents at zeroth order in derivatives are given by The subscript "d" for "dual" is used to differentiate the parameters from those used earlier in the paper. The currents can be compared with our eq. (3.11) and the dual of eq. (2.4) at zeroth order: where w m ≡ T p m,T + µp m,µ = T s + µn is the enthalpy density, and µ m = 1/(1 − 2∂p m /∂B 2 ) is the magnetic permeability. Using h 2 = 1, we can identify terms. Out of equilibrium, h µ and µ d are auxiliary dynamical variables (without a unique microscopic definition) designed to capture the dynamics of the magnetic field. The entropy density is s d = p m,T + µ T p m,µ , as follows from ε d + p d = T s d + µ d ρ d . The energy densities coincide, ε d = −p + T s + µn = ǫ, again with p = − 1 2 B 2 + p m (T, µ, B 2 ). At order O(∂), our constitutive relations can not be directly compared to those of GHI because of different hydrodynamic variables. However, we can compare the number of transport coefficients. The comparison may be done based on the entropy current argument which we review below.
In a particular hydrodynamic "frame", the one-derivative contributions to the GHI constitutive relations are given in eq. (3.4), (3.5) of ref. [22], are all transverse to both u µ and h µ , the tensor t µν d is symmetric and traceless, and the tensor s µν d is anti-symmetric. We do not write the subscript on the temperature and fluid velocity, even though the GHI's T and u µ differ from ours at O(∂). Further, GHI impose charge conjugation as a constraint on the dynamics.
The tilded vectors are defined asṼ µ = ǫ µναβ u ν h α V β , and the tilded shear tensor is as in eq. (3.8). The tensor s ρσ d has only one degree of freedom, hence it contains only one transport coefficient. The divergence of the entropy current (5.6) is then The three tilded coefficients do not contribute to entropy production in eq. (5.6) due tõ V µ V µ = 0 and σ ⊥µνσ µν ⊥ = 0, and can take any real values, Demanding that ∇ µ S µ d in eq. (5.9) is non-negative now implies together with the condition that the quadratic form in the first line of eq. (5.9) is positive semi-definite. The latter gives Thus there are eleven apriori independent non-equilibrium transport coefficients listed in Eqs. (5.8) that are consistent with non-negative entropy production, provided the constraints (5.11) are satisfied. The coefficientsr ⊥ ,η ⊥ ,η are odd under charge conjugation C, and can be eliminated if one demands C-invariance of hydrodynamics. An implicit assumption of ref. [22] amounts to choosing f 1 = −f 2 = ζ ⊥ , τ 1 = 0, τ 2 = 2ζ .

Kubo formulas
Assuming time-reversal covariance, the above transport coefficients can be further constrained by the Onsager relation (3.20). In order to find the retarded functions, we can use exactly the same variational procedure as in sec. 4.5: as well as Again, the subscript "on-shell" signifies that T µν and J µν are evaluated on the solutions to the conservation equations (5.1) with the constitutive relations (5.8). We use the above prescription to evaluate correlation functions at zero spatial momentum, which gives rise to Kubo formulas. Demanding that the correlation functions satisfy (3.20) now gives the Onsager relation We further find the following Kubo formulas for transport coefficients in the constitutive relations (5.8). The resistivities are given by the "shear viscosities" are given by , (5.14e) and the "bulk viscosities" are given by The Onsager relation (3.22) maps to the Onsager relation (5.13), as expected. The entropy current constraints (3.19) map to the entropy current constraints (5.11), as expected. Finally, the mapping of transport coefficients (5.15) can be used to compare the eigenfrequencies of small oscillations of the (dynamically) neutral state found in eq. (4.9), (4.10) to those found in ref. [22]. Using the map of thermodynamic parameters spelled out below eq. (5.3), the speed of Alfvén waves agrees with ref. [22]. The damping coefficient of Alfvén waves in eq. (4.9) agrees with ref. [22] when B 2 /µ m ≪ ǫ + p. The speed of magnetosonic waves in eq. (4.10b) agrees with ref. [22]: in order to see this, note that the assumption of constant magnetic permeability amounts to assuming that the equation of state takes the form p d = 1 2 µ m µ 2 d + F (T ), or p = − 1 2µm B 2 + F (T ), with some F (T ). In general, the speed of magnetosonic waves derived from the formalisms of sec. 4 and sec. 5 will not agree, except when B 2 /µ m ≪ (ǫ + p). One reason is that the chemical potential for the electric charge is treated as a thermodynamic variable in sec. 4, hence the magnetosonic wave speed will in general depend on the charge susceptibility (∂n/∂µ) µ=0 . This thermodynamic derivative is not present in the formalism of sec. 5. Finally, note that the transport coefficient τ 1 contributes to damping of fast magnetosonic waves, for example at θ = 0 we have Γ ms = (τ 1 + τ 2 )/(T s d ), in agreement with eq. (4.10f).

Discussion
In this paper we have presented the equations of relativistic magnetohydrodynamics, by which we mean the hydrodynamics of a conducting fluid in local thermal equilibrium, with dynamical electromagnetic fields. MHD is naturally formulated in a derivative expansion with magnetic field B ∼ O(1). Electric screening does not imply that the electric field vanishes: rather, it implies E ∼ O(∂) is subleading in the derivative expansion. We have adopted the simplest "mean-field" formulation in which the constitutive relations in the theory with dynamical electromagnetic fields are inherited from the theory with external electromagnetic fields. Our main focus was on transport coefficients. For a parity-symmetric microscopic system, we find eleven transport coefficients at one-derivative order. One transport coefficient is thermodynamic: it is a part of the equation of state in curved space, and contributes to flat-space correlations. Transport coefficients of this type in relativistic hydrodynamics were first identified in [2] where they appeared at second order in derivatives. In 2+1 dimensional hydrodynamics, thermodynamic transport coefficients can already appear at first order in derivatives [24]. Of the remaining ten transport coefficients, three are non-equilibrium and non-dissipative, and seven are non-equilibrium and dissipative. There are more transport coefficients for parity-violating fluids, as listed in sec. 3. We now comment on questions not discussed in detail in the main body of the paper. Angular momentum generated by the magnetic field.-The thermodynamic transport coefficient M Ω determines the response of equilibrium magnetic polarization to vorticity, as can be seen from eq. (3.13). One way to view M Ω is to note that a system of charged particles in external magnetic field will develop angular momentum. One can see this in the thermodynamic framework of sec. 2. For a bounded system, the equilibrium energymomentum tensor obtained by varying the equilibrium free energy (2.1), (2.7) with respect to the metric will have a boundary contribution after the variation M Ω B·δ g Ω is integrated by parts [17]. The surface momentum density Q α s = M Ω ǫ αµνρ u µ B ν n ρ (where n µ is the unit spacelike normal vector to the boundary) will give rise to angular momentum induced by the magnetic field. Consider a system at rest in flat space at constant temperature, charge density, and constant magnetic field B. The angular momentum L derived from the energy-momentum tensor only receives a boundary contribution, and one finds where V is the spatial volume. In this sense M Ω determines "angular momentum density". As the coefficient M Ω is odd under charge conjugation C, this generation of angular momentum only happens in a C-invariant theory if the equilibrium state has non-zero charge density. Similarly, for a system not subject to the magnetic field, in flat space, which rotates uniformly with small (namely |ω|R ≪ 1 where R is the size of the system) angular velocity ω, the magnetization density is m = 2M Ω ω. More generally, the susceptibility M Ω provides a macroscopic parametrization of gyromagnetic phenomena such as the Barnett and Einstein-de Haas effects.
Previous work on transport coefficients.-Papers [25,26] studied transport coefficients for relativistic fluids subject to an external magnetic field. While this does not correspond to MHD in the sense described in this paper (we define MHD as a theory in which magnetic field or its auxiliary is a dynamical degree of freedom), a fluid in external field is a fundamental building block for MHD. Parts of Refs. [25,26] overlap with our Section 3. Some of our results differ from those in Refs. [25,26]: the analysis of thermodynamics, the number of transport coefficients, constraints on transport coefficients imposed by the positivity of entropy production, and some of the Kubo formulas. The details are given in Appendix B.
Dual formulation of magneto-hydrodynamics.-In sec. 5 we compared our results with the recent "dual" formulation of MHD in ref. [22]. We found the same number of transport coefficients in the two approaches, provided the bulk viscosity missed in ref. [22] is restored, and the constraint of C-invariance imposed in ref. [22] is lifted. It would be interesting to investigate the relation between the "dual" and "conventional" formulations of MHD further, in particular with regard to the description of electric charge fluctuations.
Applicability regime.-The MHD described in this paper treats electromagnetic fields classically. This means that the electromagnetic coupling constant must be small so that quantum fluctuations of the electromagnetic field can be ignored. The applicability regime of MHD also includes B ≪ T 2 (or restoring the fundamental constants ceB ≪ (k B T ) 2 ), as is necessary to restrict the hydrodynamic degrees of freedom to those inherited from thermodynamics. We do not have a method to systematically incorporate the effects of larger magnetic fields within the MHD description of sec. 4. The classical hydrodynamic theory also ignores statistical fluctuations, which are known to invalidate classical second-order hydrodynamics in 3+1 dimensions (and classical first-order hydrodynamics in 2+1 dimensions). Understanding the effects of statistical fluctuations in magnetic field requires further work.
Transport coefficients at strong coupling.-While the small electromagnetic coupling allows one to treat magnetic fields classically, other interactions in the theory do not have to be small. For strongly interacting non-abelian gauge theories in external U(1) magnetic field, methods of gauge-gravity duality provide a window into non-equilibrium physics, both within and outside the hydrodynamic regime. Some of the hydrodynamic transport coefficients discussed in this paper were evaluated in holographic models in refs. [26,27]. The full set of transport coefficients for fluids in external magnetic field has not yet been explored holographically.
connections to our work. We thank Shira Chapman, Kristan Jensen, and Adam Ritz for helpful conversations. This work was supported in part by NSERC of Canada.

A Equilibrium T µν and J µ
The coefficients ǫ n , π n , φ n , γ n , δ n , θ n in the equilibrium energy-momentum tensor and the current (3.5) have the following expressions in terms of the five parameters M n (T, µ, B 2 ) of the generating functional (2.7). The O(∂) correction to the energy density is determined by The O(∂) correction to the spatial current is determined by the magnetic susceptibility, In this appendix we will comment on how our work relates to some earlier studies of transport coefficients, for the benefit of the reader who might want to compare different approaches. Ref. [25], abbreviated below as HSR, studied relativistic hydrodynamics of parity-invariant fluids in external non-dynamical magnetic field. HSR enumerated the transport coefficients, giving a relativistic version of the classification in the book [30], §13, and derived the Kubo formulas for transport coefficients in an operator formalism. Parts of the HSR paper overlap with our Section 3.
There are also some differences between our Section 3 and HSR. In terms of the setup, the HSR treatment neglects electric fields, while we include them and explain how to do so systematically. Related to that, the treatment of polarization effects in HSR was incomplete. A direct way to obtain the equilibrium energy-momentum tensor and the current in the presence of external fields is by varying the corresponding generating functional with respect to the metric and the gauge field, as was done for example in ref. [17]. As a result, HSR did not include the thermodynamic transport coefficient, denoted in Section 3 as M Ω , and did not distinguish between the Landau-Lifshitz and thermodynamic frames. In the Landau-Lifshitz frame, M Ω would contribute to all frame invariants in eq. (3.6) inducing O(∂) contributions to pressure, electric current, and spatial stress.

(B.2)
On the other hand, the constraints coming from the second law in ref. [25] state that all the dissipative HSR transport coefficients must be positive. We find that the constraints on dissipative transport coefficients (B.2) are in fact weaker. In other words, the constraints of ref. [25] are too restrictive: some of the dissipative transport coefficients in the HSR notation can be negative, while still satisfying (B.2), and therefore still leading to positive entropy production. Finally, there are differences between our Kubo formulas and those of HSR. In particular our Kubo formulas for conductivities transverse to the external magnetic field are markedly