Thermodynamics of polarized relativistic matter

We give the free energy of equilibrium relativistic matter subject to external gravitational and electromagnetic fields, to one-derivative order in the gradients of the external fields. The free energy allows for a straightforward derivation of bound currents and bound momenta in equilibrium. At leading order, the energy-momentum tensor admits a simple expression in terms of the polarization tensor. Beyond the leading order, electric and magnetic polarization vectors are intrinsically ambiguous. The physical effects of polarization, such as the correlation between the magneto-vortically induced surface charge and the electro-vortically induced surface current, are not ambiguous.


Introduction
We would like to understand collective macroscopic behaviour of matter subject to external fields. In the high-temperature limit this is often captured by classical hydrodynamics. The ingredients that go into writing down the hydrodynamic equations are: the identification of relevant variables (conserved densities, order parameters), the derivative expansion (small gradients near equilibrium), and symmetry constraints. The hydrodynamic equations are modified when the system is subject to external electric and magnetic fields. The latter will induce polarization (electric, magnetic, or both) in a fluid, and as a result the transport properties of the fluid will change. Our focus here will be on isotropic relativistic matter because a) electromagnetic fields are intrinsically relativistic, b) relativistic fluids have more symmetry than non-relativistic fluids, and c) relativistic fluids have been a subject of much recent attention in the literature due to their appearance in heavy-ion physics [1,2], in gravitational physics, through the holographic duality [3], and even in condensed matter physics [4,5]. The systematic description (including the derivative expansion) of polarized relativistic fluids is largely missing, and the present paper is a step in filling that gap.

JHEP07(2016)028
In order to understand the hydrodynamics of matter subject to external fields, one needs to understand its thermodynamics first. In what follows we will describe the procedure for obtaining the energy-momentum tensor and the current density for stationary equilibrium polarized matter subject to external gravitational and electromagnetic fields. We will find simple expressions for "bound" currents, including equilibrium surface currents and surface momenta.
Let us start with the standard description of equilibrium thermodynamics without external fields. In the grand canonical ensemble at temperature T 0 = 1/β 0 and chemical potential µ 0 , extensivity in the large-volume limit dictates that the logarithm of the grandcanonical partition function Z[T 0 , µ 0 ] is proportional to the d-dimensional spatial volume, where the pressure P (T 0 , µ 0 ) is constant in equilibrium [6]. The partition function Z may be computed from a Euclidean path integral with a Euclidean (imaginary) time compactified with period β 0 , see e.g. [7]. In the path integral action, the fundamental fields of the microscopic theory can then be coupled to time-independent external sources: the (Euclidean) metric g E µν and the (Euclidean) gauge field A E µ . The gauge field couples to the conserved current, whose time component is the charge density corresponding to the chemical potential. See ref. [8] for a convenient parametrization of the Euclidean sources g E and A E . The Euclidean path integral gives rise to the partition function Z = Z[T 0 , µ 0 , g E , A E ], where T 0 = 1/β 0 is the coordinate periodicity of the Euclidean time. We assume that the coupling to time-independent external sources leaves the system in equilibrium, so that no entropy is produced. The temperature and the chemical potential will be altered by the external sources and are not uniform any more. For example, the equilibrium temperature becomes T (x) = T 0 / g E 00 (x) [6]. Similarly, the chemical potential will be shifted by the time component of the external gauge field. We can write W = −i ln Z as where √ g E is the square root of the determinant of g E µν , and F is the negative of the grand canonical free energy density. In flat space and without external gauge fields, F reduces to the pressure P , and in general F is a complicated function of the spatially varying external sources. In a slight abuse of terminology, we will refer to F as the free energy density, and to W as the free energy. Varying W with respect to a time-independent source gives rise to a zero-frequency insertion in the Euclidean path integral of the operator coupled to the source. The relevant operators are the energy-momentum tensor (coupled to the metric), and the conserved current (coupled to the gauge field). Thus W is the generating functional for zero-frequency correlation functions of the energy-momentum tensor and the current in equilibrium.
The Euclidean external sources g E µν and A E µ may be "un-Wick-rotated" to Minkowski time to obtain the physical real-time external sources g µν and A µ , for example g E 00 = −g 00 , g E 0k = −ig 0k , A E 0 = −iA 0 etc. In what follows we will omit the dependence on T 0 and µ 0 ,

JHEP07(2016)028
and will denote the Euclidean generating functional with arguments continued to physical time as W [g, A], so that We may as well view d d+1 x as containing an integral over the physical time, as the argument of the integral does not depend on time anyway. For a relativistic microscopic theory without gauge and gravitational anomalies, the generating functional is both gauge-and diffeomorphism-invariant. Let us further assume that all long-range interactions are screened due to a non-zero temperature T , so that the spatial correlations are local on scales longer than the screening length. The effective description of static correlations on such long scales will then be given by W [g, A], where the density F is a local function of the external sources. For external fields that vary slowly in space, the above locality implies that F may be written as a derivative expansion in the gradients of the external fields. See ref. [8] for a study of the local generating functional in the Euclidean form, and ref. [9] for the Minkowski form. Here we will use the Minkowski form, in which the underlying gauge and diffeomorphism invariance is manifested in a more straightforward way.
In order to implement the derivative expansion in practice, one needs to postulate the derivative scaling of the external sources g µν and A µ . Physically, this amounts to deciding whether the external sources are taken as "strong" or "weak" on the scale of the spatial inhomogeneity in equilibrium. In refs. [8,9], the external sources were taken as "weak" in the sense that both g µν and A µ were assumed to be O(1) in the derivative expansion, so that both electric and magnetic fields appear at order O(∂) in the expansion. This choice of scaling makes the description of equilibrium polarization rather awkward: for example, in 3+1 dimensions, the thermodynamic response to constant homogeneous magnetic field B appears at the same order as the response to two derivatives of temperature, B 2 ∼ (∂T ) 2 .
In order to describe polarized matter in constant (or slowly varying in space) electric and magnetic fields, a different derivative counting scheme is more natural, one in which constant homogeneous electric and magnetic fields are taken to be O(1) in the derivative expansion, rather than order O(∂). This will be our goal here: to implement the derivative expansion of the free energy in the regime when the external gravitational field is still "weak" so that g µν is O(1), while the external electromagnetic field is "strong" so that the Following the general approach of ref. [9], we will obtain simple expressions for the energy-momentum tensor and the conserved current in relativistic polarized matter subject to external fields.

Thermodynamic parameters
Let us first outline the starting point. The free energy W [g, A] is a gauge-and diffeomorphism-invariant functional of A µ , g µν , and their derivatives. Being in equilibrium means that that there is a Killing vector V , such that the Lie derivative with respect JHEP07(2016)028 to V vanishes on all observables, £ V (. . . ) = 0. In suitable coordinates, V µ = (1, 0). In the grand canonical ensemble, the equilibrium state is parametrized by the temperature, velocity, and the chemical potential. Their relation to the external sources is where β 0 is a constant setting the normalization of temperature, and Λ V is a gauge parameter which ensures that µ is gauge-invariant. The constant µ 0 is absorbed into Λ V . Without external gauge fields, relations (2.1) are the covariant versions of the statement that T √ −g 00 and µ √ −g 00 are constant in equilibrium [6]. The vector u µ is the normalized (u 2 = −1) velocity of matter, and the coordinates in which V µ = (1, 0) correspond to the matter "at rest". Both T and µ are gauge invariant and transform as scalars under diffeomorphisms. For a discussion of gauge and diffeomorphism covariance of the equilibrium parameters see section 5 of ref. [10]. For a system occupying a spacetime region M with a boundary ∂M, we assume that the generating functional can be separated into bulk and boundary contributions, and we take Here the first term describes the bulk contribution, and the second term the boundary contribution. To leading order in the derivative expansion, F is the pressure, and L is the surface tension. In the bulk term, g is the determinant of g µν , and F is a function of T , u µ , µ, as well as of the sources A µ , g µν , and their derivatives. For the boundary with coordinates y a whose shape is specified by x µ (y a ), the tangent vectors are e µ a = ∂x µ /∂y a , and the projector onto the boundary is P µν = g µν − n µ n ν . The induced metric on the boundary is γ ab = e µ a e ν b g µν . In the boundary term, γ is the determinant of the induced metric, while L in addition may depend on n µ , the spacelike unit normal vector to the boundary.

Response to external sources
The energy-momentum tensor and the current are defined as where £ n is the Lie derivative along the normal, the dots denote boundary terms with higher normal derivatives of the sources. The variations are performed at fixed V µ and Λ V . Here T µν , J µ are the bulk energy-momentum tensor and the current, and T µν s , J µ The derivative expansion for the free energy density F in the generating functional (2.2) gives rise to the derivative expansion for the equilibrium T µν and J µ , as described in [9]. The boundary energy-momentum tensor and the current in (2.3) may be decomposed into the contributions tangential to and normal to the boundary, Here δg a is the pullback of P λ µ δg λν n ν to the boundary, δg n = n µ n ν δg µν , δA a is the pullback of δA µ to the boundary, and δA n = n µ δA µ . Similarly, one can vary the generating functional with respect to the field strength F µν , where again δF ab is the pullback of δF µν to the boundary, δF a is the pullback of P λ µ δF λν n ν to the boundary, and the dots denote boundary terms with higher normal derivatives of δF µν . This defines the bulk polarization tensor M µν , and the boundary polarization tensor M ab s . 1 The surface terms J s , Π a s , Π s , M a s depend on how the equilibrium is set up, and what the boundary conditions on ∂M are, as determined by the nature of the phase separation at ∂M.
In all the above variations, we assume that the region M occupied by matter is unchanged. One could also consider the response of the generating functional to changing the shape of ∂M, however this will not be needed for our purposes. See ref. [11] for a recent discussion of surface terms in the Euclidean generating functional.
The polarization tensor contains both electric and magnetic components. We define the electric field as E µ ≡ F µν u ν , the magnetic field as B ≡ − 1 2 µαβ u µ F αβ for d = 2, and B µ ≡ 1 2 µναβ u ν F αβ for d = 3. In 1+1 dimensions, we define the "magnetic field" as B ≡ 1 2 µν F µν , so that F µν = −B µν . The Levi-Civita tensor is µναβ = ε µναβ / √ −g, with ε 0123 = 1, and similarly in other dimensions. Both E µ and B µ are spacelike and orthogonal to u µ . We have the following decomposition of the field strength: 1 There is a gravitational analogue of the polarization tensor which involves varying the generating functional with respect to the connection coefficients. The energy-momentum tensor then takes the form analogous to eq. (2.14) below. See section 5 of ref. [10].

JHEP07(2016)028
The electric polarization vector p α and the magnetization vector m α (for d=3) are defined by rewriting the integrand in (2.5) as 1 2 M µν δF µν = p α δE α + m α δB α . For d=2, the variation is 1 2 M µν δF µν = p α δE α + mδB, which defines the magnetization m. The decomposition of the polarization tensor into the electric and magnetic parts is then Both p α and m α are transverse to u α .

Equilibrium relations
The equilibrium definitions (2.1) together with £ V (. . . ) = 0 give where a µ ≡ u λ ∇ λ u µ is the acceleration vector, u µ a µ = 0. These relations imply that T ∂ λ (µ/T ) − E λ vanishes in equilibrium. In other words, a system subject to an external electric field will develop a gradient of µ/T in order to compensate the applied field and ensure that the equilibrium is maintained. This has implication for derivative counting. For "weak" electric fields E λ ∼ O(∂), the gradients of T and µ are O(∂) as well. For "strong" electric fields E λ ∼ O(1), there will be an O(1) gradient of µ/T . How exactly this gradient is achieved depends on the nature of the microscopic degrees of freedom. Given that the chemical potential determines the number of charge carriers, we take "strong" electric fields to mean that both E and ∂µ are O(1), while ∂T is still O(∂), so that ∂µ µ ∂T T . In the generating functional, the derivatives of the chemical potential may then be traded for the electric field.
Similarly, the derivative of the velocity can be decomposed in equilibrium as The vorticity is Ω ≡ − µνλ u µ ∇ ν u λ for d = 2, and Ω µ ≡ µναβ u ν ∇ α u β for d = 3. This velocity decomposition implies that both the expansion ∇ µ u µ and the shear tensor . This is as it should be: out of equilibrium, the expansion would contribute to dissipation through bulk viscosity, and the shear tensor would contribute to dissipation through shear viscosity.

JHEP07(2016)028
Combined with the electromagnetic "Bianchi identity" µναβ ∇ ν F αβ = 0 in 3+1 dimensions, the velocity decomposition (2.8) implies These are the covariant versions of the familiar flat-space equilibrium relations ∇·B = 0 and ∇×E = 0. More generally, for the electric field in equilibrium we have as a consequence of £ V E α = 0 and E α u α = 0.

Polarization ambiguities
The electromagnetic Bianchi identity also implies that there is an ambiguity in the definition of the polarization tensor: in 3+1 dimensions, one can always add to the generating functional an identically vanishing term a function of the field strength and its derivatives. Such a term shifts the polarization tensor by The polarization vectors correspondingly shift as In 2+1 dimensions, we can add an identically vanishing term W ∅ = 1 2 √ −g C µαβ ∇ µ F αβ , where again C can be a function of the field strength and its derivatives. The polarization tensor then shifts by The electric polarization vector correspondingly changes as while the magnetic polarization m remains unchanged. The variational derivatives of W ∅ with respect to both g µν and A µ vanish. As a result, the energy-momentum tensor and the current (both bulk and boundary) are not affected by such unphysical shifts.

Bound charges and bound currents
For matter whose degrees of freedom carry gauge charges, it is conventional to separate the charge into the "free charge" and "bound charge" components. In the grand canonical ensemble, the chemical potential µ describes the coupling of the system to a reservoir JHEP07(2016)028 of "free charges". Demanding local charge neutrality for free charges in the bulk would amount to demanding ∂F /∂µ = 0. Doing so would eliminate the contribution of free charges to polarization. One may refer to µ-dependent contributions as coming from "free charges", and µ-independent contributions as coming from "bound charges", though such a separation is somewhat artificial. We will not impose ∂F /∂µ = 0, and will keep the contribution to polarization from both free charges and bound charges.
The current J µ admits a simple expression in terms of the polarization tensor to any order in the derivative expansion. Indeed, the free energy density F can be written as where the coefficients S (n) do not contain derivatives of the electromagnetic field strength.
The derivative of the chemical potential can be traded for the electric field according to (2.7), hence we can take S (n) = S (n) (T, µ, F αβ , . . . ) where dots denote the arguments which do not depend on the gauge field. The polarization tensor can be easily found in terms of S (n) through integration by parts. It is then clear that the current extracted from the generating functional according to (2.3) is 14) to any order in the derivative expansion, where ρ ≡ ∂F /∂µ. The first term in the righthand side is the standard equilibrium current in the absence of polarization: to leading order in the derivative expansion the free energy density F is just the pressure P , and ρ = ∂P/∂µ is the density of "free charges". The second term in the right-hand side is a total derivative of an anti-symmetric tensor. It therefore does not contribute to the conservation equation ∇ µ J µ = 0, and can be interpreted in terms of "bound" charges and "bound" currents. It is clear from the expression (2.14) that the unphysical polarization shifts (2.10) and (2.12) do not affect the current. The current can be decomposed with respect to the velocity u µ as 15) where N ≡ −u µ J µ is the charge density, and the spatial current J µ ≡ ∆ µλ J λ is transverse to u µ . For the polarization tensor of the form (2.6), the definitions (2.1) together with £ V (. . . ) = 0 lead to the following equilibrium expressions for the charge density: Consider the charge density in d = 3 spatial dimensions. The second term in the right-hand side is the familiar electrostatic bound charge density, which in flat space reduces to −∇·p. The third term is the bound charge density induced by gravity: in the static Newtonian gravitational field it becomes p·∇ϕ, where ϕ is the gravitational potential. The last term is the bound charge density induced in magnetized matter which is rotating. For a system JHEP07(2016)028 undergoing rotation with small (meaning |ω|R 1, where R is the size of the system) angular velocity ω, the last term in the right-hand side becomes −2m·ω.
Similarly, the definitions (2.1) together with £ V = 0 lead to the following equilibrium expressions for the spatial current: Consider the current density in d = 3 spatial dimensions. The first term in the right-hand side is the familiar bound current, which in flat space reduces to ∇×m. The second term is the bound current induced by the gravitational field: in the static Newtonian gravitational field it reduces to (∇ϕ)×m, where ϕ is the gravitational potential. We emphasize that the above expressions for bound charges and bound currents are simply a consequence of thermal equilibrium. Equations (2.16), (2.17) do not assume any particular microscopic model of matter, and moreover they hold to any order in the derivative expansion.

Derivative expansion
We close this section with a comment on the derivative expansion of the free energy. As a schematic example, consider the functional W [a, g] which depends on two sources a(x) and g(x) which both vary slowly in space. Assuming locality, the derivative expansion is and the boundary terms are implied.
Suppose now that a changes much faster than g, such that |a (x)/a(x)| |g (x)/g(x)|. Naively, one may think that the terms containing the derivatives of a(x) are more important than those with derivatives of g(x), and there is a separate derivative counting associated with a(x) and g(x). This is not in general so: for example, integrating the P 2 term by parts gives rise to (∂P 2 /∂a) g(x)a (x), which may be of the same order as the P 1 term. It is possible to count the derivatives of g differently from the derivatives of a if the "cross" terms (∂P 2 /∂a) are in some sense small. For example, we could introduce two counting parameters ε and γ ε and count the derivatives as a ∼ ε, g ∼ γ, while ∂P 2 /∂a, ∂P 4 /∂a, ∂P 5 /∂a, ∂P 7 /∂a are of order γ/ε.
Physically, a will be the external gauge potential A µ , and g the external metric g µν . By "strong" electromagnetic fields we will mean the fields such that this derivative counting is valid, i.e. electromagnetism is more important than gravity. For such "strong" fields, by the leading order in the derivative expansion we will mean: i) setting γ to zero, ii) isolating terms polynomial in a (x), and iii) summing those terms into a single function P (a, g, a ).

Weak electromagnetic fields
Let us start with "weak" electromagnetic fields. For the sources with A µ ∼ O(1), g µν ∼ O(1), there are only two gauge and diffeomorphism invariants at leading order in the derivative expansion in the bulk: T and µ. On the boundary, there is an extra invariant u n ≡ u µ n µ . The static generating functional to leading order in the derivative expansion is then The definitions (2.3) give where s = ∂P/∂T is the bulk entropy density, ρ = ∂P/∂µ is the bulk charge density. These are the standard expressions for the energy-momentum tensor and the current in a relativistic perfect fluid. The boundary energy-momentum tensor and the current are and the other boundary terms are Π a s = (T s s + µρ s ) u n u a , Π s = 1 2 (T s s + µρ s ) u 2 n + ε s u n , and J s = ρ s u n , where we have defined s s ≡ ∂L/∂T , ρ s ≡ ∂L/∂µ, ε s ≡ ∂L/∂u n . Again, these describe a perfect fluid on the boundary with pressure L. At leading order in the derivative expansion, both K µ s and K µν s vanish.

Strong electromagnetic fields
Let us now consider "strong" electromagnetic fields, such that F µν ∼ O(1) and g µν ∼ O(1) in the derivative expansion. To leading order, the static generating functional is The dependence on F αβ includes the dependence on electric and magnetic fields, and for the boundary part also on their normal components. 2 2 When d+1 is odd, there may be Chern-Simons terms in P . The Chern-Simons term is not gauge invariant on the boundary, so in this case L must contain an anomalous piece, whose gauge variation exactly cancels the gauge variation of the Chern-Simons term. In the application to the quantum Hall effect, the anomalous boundary piece comes from the massless 1+1 dimensional chiral modes on the boundary [12]. Upon integrating out the massless boundary modes, L will in general become a non-local function of the electric and magnetic fields. If the dynamics of the boundary modes can be described classically, they may be treated directly within the generating functional, similar to what is done in ref. [13] for superfluids.
In what follows, we will ignore the massless boundary modes, and will only explore the consequences of short-distance correlations on the boundary.

JHEP07(2016)028
The bulk current is given by (2.14), with the polarization tensor M µν = 2∂P/∂F µν . In what follows we will express M µν in terms of electric and magnetic susceptibilities.
In order to find the energy-momentum tensor, we need to be more specific about the dependence of P and L on the metric. Gauge and diffeomorphism invariance requires that P = P (s (0) ) is a function of scalars s (0) , which are made out of the electromagnetic field strength (we will use the term "scalar" for both scalars and pseudo-scalars). The superscript signifies that we are working to leading order in the derivative expansion. The number of scalars s (0) depends on the dimension. To leading order in the derivative expansion, we choose to work with the following independent scalars: Let us express the bulk energy-momentum tensor using the decomposition with respect to the velocity u µ , as is often done in relativistic fluid dynamics, Here E ≡ u µ u ν T µν is the energy density, P ≡ 1 d ∆ µν T µν is the pressure, the momen- T αβ is transverse to u µ , symmetric, and traceless. Given P as a function of the above scalars, the energy-momentum tensor can be read off from the definition (2.3). For d = 1, we have P = P (T, µ, B). This gives the following energy-momentum tensor: Here the "magnetization" density m ≡ ∂P/∂B determines the polarization tensor as For d = 2, P = P (T, µ, B, E 2 ). This gives the following energy-momentum tensor:

JHEP07(2016)028
Here again s = ∂P/∂T is the entropy density, ρ = ∂P/∂µ is the charge density, m = ∂P/∂B is the magnetization density, and χ E ≡ 2∂P/∂E 2 is the electric susceptibility. They determine the polarization tensor as The dependence P = P (T, µ, B, E 2 ) implies that the electric polarization vector is p µ = χ E E µ , and the polarization tensor (3.7) coincides with the general expression (2.6), as it should. For d = 3, P = P (T, µ, B 2 , E·B, E 2 ). This gives the following energy-momentum tensor: Here again s = ∂P/∂T , ρ = ∂P/∂µ are the entropy and charge densities, S µ = µρσλ u ρ E σ B λ is the Poynting vector, χ EE ≡ 2∂P/∂E 2 is the electric susceptibility, is the electro-magnetic susceptibility, χ BB ≡ 2∂P/∂B 2 is the magnetic susceptibility. They determine the polarization tensor as where G µν = 1 2 µναβ F αβ is the dual field strength. The dependence P = P (T, µ, B 2 , E·B, E 2 ), implies that the polarization vectors are The magneto-electric susceptibility χ BE is equal to the electro-magnetic susceptibility χ EB , and the polarization tensor (3.9) coincides with the general expression (2.6), as it should.
So far we have presented T µν in terms of the decomposition (3.3) with respect to the velocity u µ , whose coefficients E, P, Q µ , and T µν are expressed in terms of the electric and magnetic fields, and the susceptibilities. The same energy-momentum tensors (3.4), (3.6), (3.8) can be equivalently expressed in terms of the polarization tensor M µν = 2∂P/∂F µν as is the "electromagnetic correction" to the perfect fluid form. Note that P , s, and ρ in (3.10a) are functions of the electric and magnetic fields. The tensor (3.10b) is symmetric; if we set the external electric field to zero (in two or three spatial dimensions),

JHEP07(2016)028
then T µν EM reduces to its first term, and is still symmetric. Note that ∇ µ T µν EM does not equal F ν λ (−∇ ν M νλ ). The above expression for T µν EM was first derived by W. Israel [14], for a free gas of polarized relativistic particles. We emphasize that one does not need to assume any particular microscopic model of matter in order to arrive at the above energy-momentum tensor: expression (3.10) is a direct consequence of gauge and diffeomorphism invariance of the theory, to leading order in the derivative expansion.
We now turn to the boundary energy-momentum tensor and the current which follow from the generating functional (3.1). The boundary current may be expressed in terms of the boundary polarization tensor m µν ≡ 2∂L/∂F µν (keeping T , µ, and n µ fixed). Upon integrating by parts on the boundary, the definition (2.3) gives the following boundary currents:

11)
Here M a is the boundary current arising from integrating the variation of P (T, µ, F αβ ) by parts, e µ a M a = n λ M λµ . As one can see from the polarization tensor (2.6), in flat space in 3+1 dimensions the boundary current n λ M λµ reduces to a vector whose time component (surface charge density) is p·n, while the spatial part (surface bound current) is m×n. These are the familiar expressions from electro-and magneto-statics. The other term in the boundary current, (ρ s u a − ∇ b m ba ), arises due to the presence of charged degrees of freedom on the boundary described by L, and mimics the bulk current (2.14), with ρ s ≡ ∂L/∂µ. The other boundary currents, J s = ρ s u n , and K µ s = n λ m λµ , emerge from L as well.
At leading order in the derivative expansion, the only contribution to the boundary energy-momentum tensor arises from the surface tension term in (3.1). In 2+1 dimensions, L = L(T, µ, u n , B, E 2 , E n ), where u n ≡ u µ n µ , E n ≡ E µ n µ . At leading order K µν s vanishes, and the definition (2.3) gives Here S a = e α a P αµ µρσ E ρ u σ , and the coefficients are s s ≡ ∂L/∂T , ρ s ≡ ∂L/∂µ, m s ≡ ∂L/∂B, ε s ≡ ∂L/∂u n , α s,E ≡ 2∂L/∂E 2 , χ n,E ≡ ∂L/∂E n . The boundary energy-momentum tensor in 3+1 dimensions looks similar, and we won't write it down explicitly.

Next order in the derivative expansion
We now proceed to the next (first) order in the derivative expansion, taking into account O(∂) terms in the generating functional (2.2). We will take the surface tension L to be constant for simplicity, and will focus on the bulk contributions to thermodynamics. The free energy density at first order in the derivative expansion is F = P (s (0) ) + n M n (s (0) ) s (1) n , (4.1)

JHEP07(2016)028
where P is the leading-order pressure. For weak electromagnetic fields, the leading order scalars are s (0) = {T, µ}, while for strong electromagnetic fields s (0) are given by eq. (3.2). The functions M n (s (0) ) parametrize the thermodynamic response at first order, and are determined by the microscopic theory. The gauge-and diffeomorphism-invariant scalars s (1) n depend on T , u µ , µ, and the sources A µ , g µν . The number of such first-order scalars depends on the dimension, and on whether the external electromagnetic fields are weak or strong. We will enumerate the scalars s (1) n in what follows. The bulk current is still given by the general expression (2.14). At first order, the free energy density (4.1) may be equivalently rewritten as where P is O(∂), but contains no derivatives of F αβ , while the last term parametrizes the static response to inhomogeneous electromagnetic fields. The polarization tensor is then The equilibrium relation ∂ λ µ = E λ + O(∂) now implies that for strong electric fields the leading-order polarization tensor may receive contributions from subleading terms in the generating functional The second term in the right-hand side describes a contribution of free charges to polarization. We will assume for simplicity that the effects of free charges are less important than those of bound charges, in the sense that ∂M n /∂µ ∼ O(∂). (Alternatively, the effects of the free charges may be lumped into the leading-order free energy, but in this case isolating their contribution becomes less straightforward.) The boundary currents are where M a e µ a = n λ M λµ as before, and e µ a e ν b S ab = n λ S λαβ P µ α P ν β . Equation (4.4a) shows that for strong electromagnetic fields beyond the leading order in the derivative expansion, the surface current is not determined solely by the bulk bound current any more, even in the absence of charged degrees of freedom on the boundary. Similarly, the energy-momentum tensor will differ from the simple form (3.10) beyond leading order.

JHEP07(2016)028
In d = 1, we define the vorticity as Ω ≡ µν ∇ µ u ν = µν a µ u ν . In d = 3, both the magnetic field and the vorticity are vectors, and there are no scalars at order O(∂).
Focussing on d = 2, the equilibrium generating functional is given by eq. (2.2), with the free energy density F (T, µ, B, M Ω are functions of T and µ. Note that both B and Ω are pseudoscalars, hence we are describing thermodynamics of a microscopic system which intrinsically violates parity. The bulk current is given by (2.14), with the polarization tensor The bulk energy-momentum tensor can be expressed as a general decomposition (3.3), whose coefficients are P = P , , see ref. [9,15]. The boundary energy-momentum tensor can be expressed in terms of the vector µ ≡ M Ω µαβ n α u β which is tangent to the boundary, where again a e µ a = M Ω µαβ n α u β , and we have assumed that the surface tension L is constant. The other boundary momentum currents are Π a s = u n a , and Π s = 0. The vector a (the energy-momentum analogue of the boundary magnetization current) is the density of momentum flowing along the boundary in equilibrium, as is generically expected to happen in a parity-violating system. To sum up, the boundary current is determined by the magnetization M B , while the boundary momentum is determined by M Ω .

Strong electromagnetic fields: 1+1 dimensions
Now let us turn to strong electromagnetic fields, with F µν ∼ O(1), g µν ∼ O(1). In 1+1 dimensions, there are only two independent scalars in equilibrium at O(∂) in the derivative expansion, which may be taken to be The equilibrium generating functional is given by eq. (2.2), with the free energy density  3  3  3  3  3 n/a n/a 1 Table 1. Independent O(∂) invariants in 2+1 dimensions. The first row in the table is the number of the invariant, and the second row says what the invariant is. The rows labeled C, P, T indicate the eigenvalue of the invariant under charge conjugation, parity, and time reversal, respectively. Parity in 2+1 dimensions is defined as a reflection of one of the spatial coordinates. The row labeled W shows the weight w of the invariant under a local rescaling of the metric; the invariants which do not transform homogeneously are marked as "n/a".

Strong electromagnetic fields: 2+1 dimensions
In two spatial dimensions, there is a large number of O(∂) scalars. However, equilibrium relations such as (2.7), (2.8) reduce the number of independent non-zero invariants to just eight. One choice of the independent invariants is listed in The table indicates how the invariants transform under charge conjugation, parity, and time reversal. The table also indicates the weight of the invariants under a Weyl rescaling of the metric, g µν →g µν = e −2ϕ g µν , where ϕ satisfies V µ ∂ µ ϕ = 0. A quantity Φ transforms homogeneously with weight w under the Weyl rescaling if Φ →Φ = e wϕ Φ. For a review of Weyl rescaling in relativistic hydrodynamics, see ref. [3]. Temperature T , chemical potential µ, velocity u µ , and the electric field E µ all have w = 1. The factors of T 2 and T 4 in the first four invariants in table 1 are inserted in order to ensure that the invariant has a well-defined weight. For the scalars which transform homogeneously, their weight w coincides with their mass dimension. The invariants s (1) 6 and s (1) 7 do not transform homogeneously and can not appear in a conformally invariant generating functional.
The first five invariants are in general already non-zero in flat space. For the static Newtonian gravitational field with potential ϕ, we have s (1) 6 ∼ E i ∂ i ϕ, s (1) 7 ∼ ij E i ∂ j ϕ. The last invariant is the vorticity, s (1) 8 = Ω, which is non-zero if the system is rotating. The equilibrium generating functional is then given by eq. (2.2), with the free energy density (1) n .
There are eight scalar functions M n , in addition to pressure, which specify the thermodynamic response at first order. For a system whose microscopic dynamics is PT-invariant, the coefficients M 3 , M 4 , and M 7 must vanish, in order for the generating functional to be PT-invariant (none of the leading-order invariants are PT-odd). For a system whose microscopic dynamics is conformally invariant, the generating functional must be conformally invariant as well, hence the coefficients M 6 and M 7 must vanish.

JHEP07(2016)028
While eight might seem like a large number, if one were to naively write down the constitutive relations directly for T µν and M µν in terms of all available O(∂) scalars, vectors, and tensors, doing so would involve introducing many more than eight unknown O(1) scalar functions, even in equilibrium. The generating functional, on the other hand, allows one to obtain the simplest expressions for the equilibrium quantities without overcounting the parameters.
At leading order in the derivative expansion, the electric polarization vector p λ was simply proportional to the external electric field. At first order, electric polarization can also be induced by the gradients of T , B, and E 2 . One finds The susceptibility here is a function of the parameters M n of the generating functional, and the other coefficients are as follows: The magnetization is The above m and p µ give the polarization tensor according to eq. (2.6b), and thus determine the O(∂ 2 ) contributions to equilibrium bound charges and bound currents, following (2.14). Finally, we note that the polarization ambiguities of section 2.4 allow one to simplify the polarization vector p λ : adding to the free energy the W ∅ term with ∂C/∂T = −γ 4 , ∂C/∂B = −γ 5 , ∂C/∂E 2 = −γ 6 eliminates the γ 4 , γ 5 , γ 6 contributions in (4.7), and adds the term ∂C/∂µ λαβ u α E β .
The energy-momentum tensor can be read off from the definition (2.3), however the general expressions are rather cumbersome, involving thermodynamic derivatives of all eight M n 's. It is easy to derive the energy-momentum tensor when the external electric field vanishes (in a certain set of coordinates), in which case T µν is only determined by P and tensor are where we have defined M Ω ≡ M 8 and g 1 ≡ (2M Ω − T ∂M Ω ∂T − µ ∂M Ω ∂µ ), to mimic the notation in section 4.1, and F = P (T, µ, B) + M Ω (T, µ, B)Ω after we have set the electric field to zero. Even in flat space and without external electric fields, there is an equilibrium energy flux, caused by the inhomogeneous magnetic field. The magnetization m, which determines the spatial bound current according to eq. (2.17), simplifies to m = ∂P/∂B + Ω ∂M Ω /∂B. There is a surface momentum a flowing along the boundary in equilibrium, completely analogous to the expression in section 4.1, where a e µ a = M Ω (T, µ, B) µαβ n α u β may now depend on the external magnetic field. The other boundary momentum currents are Π a s = u n a , and Π s = 0, as before. The energymomentum tensors (4.8), (4.9) will receive extra contributions proportional to the external electric field when the latter is non-zero.

Strong electromagnetic fields: 3+1 dimensions
In 3+1 dimensions, there is again a large number of O(∂) scalars, but many are not independent due to equilibrium constraints such as (2.7), (2.8). I counted twenty-one independent non-zero invariants. One choice is listed in table 2, where S µ = µρσλ u ρ E σ B λ is the Poynting vector, a µ = −∂ µ T /T is the acceleration, and Ω µ = µναβ u ν ∇ α u β is the vorticity. The notation in the table is the same as in the 2+1 dimensional case. The linear combinations in s (1) 8 , s (1) 9 , and s (1) 10 are taken so that the invariant has a well-defined weight under Weyl rescaling.
The first fifteen invariants are in general non-zero already in flat space. The equilibrium generating functional is given by eq. (2.2), with the free energy density (1) n . (4.10) There are twenty-one scalar functions M n , in addition to pressure, which specify the thermodynamic response at first order. For a system whose microscopic dynamics is PTinvariant, the coefficients M 9 , . . . , M 15 , M 18 , and M 21 must vanish, in order for the generating functional to be PT-invariant (none of the leading-order invariants are PT-odd).
For a system whose microscopic dynamics is conformally invariant, the generating functional must be conformally invariant as well, hence the coefficients M 16 , M 17 , and M 18 must vanish.

JHEP07(2016)028
Analogously to what happens in 2+1 dimensions, polarization may be induced by the gradients of the applied fields. The electric polarization vector which follows from the free energy is where X αβ ≡ αβρσ u ρ B σ . The susceptibility coefficients χ EE etc and γ k are determined by thermodynamic derivatives of the coefficients M n , and can be easily read off from the free energy density (4.10). However, as explained in section 2.4, polarization vectors only make sense up to certain redefinitions. For example, by choosing the arbitrary vector C µ in (2.11) appropriately, one can eliminate χ EΩ , and trade χ ES , γ 5 , γ 6 , γ 7 , γ 8 in favor of a single contribution proportional to λνρσ u ν ∇ ρ B σ . The coefficients χ EE and χ EB suffer from similar ambiguities. The magnetic polarization vector which follows from the free energy is where Y αβ ≡ αβρσ u ρ E σ . The susceptibility coefficients χ BB etc and δ k can be easily read off from the free energy density (4.10). The ambiguities (2.11) also affect the magnetic polarization: adding to the free energy the W ∅ term with C µ = Cu µ shifts χ BE → χ BE + ∂C/∂µ (in addition to shifting δ 1 , δ 2 , δ 3 , δ 4 ). While the polarization vectors are ambiguous, the energy-momentum tensor and the current are not. As an example, consider the M 21 term in the free energy. It gives rise to polarization vectors p µ = M 21 µνρσ u ν B ρ Ω σ and m µ = −M 21 µνρσ u ν E ρ Ω σ which do not suffer from polarization ambiguities. Such contributions to p µ and m µ only come from M 21 , and therefore the magneto-vortical response of the surface charge density is correlated with the electro-vortical response of the surface current. The corresponding boundary current is where Ω n ≡ Ω·n, E n ≡ E·n, and we have assumed u n = 0. As another example, let us set the electric field to zero (in a certain set of coordinates), while keeping the magnetic field non-zero. The equilibrium bulk energy-momentum tensor is then determined by only four functions M 4 , M 15 , M 17 , and M 20 , in addition to JHEP07(2016)028 the leading-order pressure P = P (T, µ, B 2 ). The correction to the leading-order energymomentum tensor (3.8) is straightforward to derive, and we will not write it down explicitly. Both E, P, and T µν will receive derivative corrections, proportional to ∂ µ T , ∇ µ B ν , and Ω µ . In addition, the magneto-vortical term M 20 will give rise to equilibrium energy currents Q µ proportional to µνρσ u ν B ρ a σ , and µνρσ u ν B ρ ∂ σ B 2 . There is also a nonvanishing boundary energy-momentum tensor T αβ s , defined by (2.3). In the decomposition T αβ Here n µ is the unit normal vector the boundary as before, and we have omitted the surface tension L. One can see that in addition to the standard surface tension, even a uniform magnetic field generates energy density, pressure, energy current, and spatial stress on the boundary.

Summary
Let us summarize. We have presented the equilibrium free energy of isotropic relativistic matter, in the regime when external electromagnetic fields are more important than external gravitational fields. From a technical point of view, this amounts to generalizing the analysis of ref. [9] by i) performing a partial summation of electromagnetic contributions, and ii) by taking into account surface terms in the generating functional. From a physical point of view, this amounts to describing the effects of polarization. The equilibrium electric current can be expressed in terms of the polarization tensor M µν to all orders, J α = ρu α − ∇ λ M λα . The charge density and the spatial current are given by eqs. (2.16) and (2.17). In 3+1 dimensions in flat space 3 they reduce to where p is the electric polarization vector, m is the magnetic polarization vector, and ω is the angular velocity. These expressions generalize the familiar n = ρ − ∇·p and j = ∇×m JHEP07(2016)028 in electro-and magneto-statics. At leading order in the derivative expansion, the surface current is J µ s = n λ M λµ , which says that the surface charge density is p·n, and the surface spatial current is m×n.
The notion of polarization is ambiguous when the external fields vary in space. This is because polarization is defined as a response to electric and magnetic fields, which are not fundamental quantities, but are rather derived from the vector potential A µ . In particular, the electro-vortical susceptibility χ EΩ is unphysical, as well as the magneto-electric susceptibility χ BE in the presence of free charges. Nevertheless, most O(∂) contributions to polarization are not affected by this ambiguity and may be derived from the equilibrium free energy, as described in section 4. For example, in a parity-violating system, there is a contribution to the electric polarization vector p ∝ B × ω, and the contribution to the magnetic polarization vector m ∝ E × ω which do not suffer from this ambiguity. The corresponding surface charge density σ s = c n·(B × ω) and the surface current j s = c (E(ω·n) − ω(E·n)) are determined by the same coefficient c = 2M 21 .
When the external fields are non-uniform, the boundary charge and spatial current are no longer determined by polarization. This is not surprising: while the polarization vectors are ambiguous, the charge and the current are not. At first order in the derivative expansion the boundary current is where P λ ρ = δ λ ρ − n λ n ρ , and S αρσ is defined by (4.2). For the generating functional (4.10) in 3+1 dimensions, there are 15 contributions to S αρσ . As an example, consider the effect of the M 2 term for non-rotating matter at constant temperature in flat space. For the boundary with vanishing extrinsic curvature (∇ µ n ν = 0), the surface charge density is where E n is the normal component of the electric field. This describes the response of the boundary charge density to the changes of the external electric field along the boundary. The equilibrium energy-momentum tensor to leading order takes a simple form (3.10) which we repeat here: T µν = P g µν + (T s + µρ)u µ u ν + T µν EM , This expression is model-independent, and is a leading-order consequence of gauge invariance, diffeomorphism invariance, and locality (on scales longer than the screening length). Beyond the leading order in derivatives, the form of the equilibrium T µν is more complicated. Equilibrium µ-independent contributions to the current J α are usually referred to as "bound charges" and "bound currents". There exist analogous contributions to the equilibrium T αβ , which one may similarly christen "bound energy", "bound pressure", JHEP07(2016)028 "bound momentum", and "bound stress". Just like bound charges and bound currents, these live both in the bulk and on the surface. For matter subject to external magnetic field (and no electric field), there will be bulk energy currents Q ∝ B×∇T , Q ∝ B×∇B 2 . The boundary energy current Q s = χ BΩ B × n is determined by the same susceptibility χ BΩ which fixes the response of magnetization to rotation, m = χ BB B + 2χ BΩ ω + O(∂T, ∂B).
Finally, our discussion so far was restricted to the state of global equilibrium, i.e. to thermodynamics. It is straightforward to extend it to hydrodynamics of polarized relativistic matter, if one assumes that the external electromagnetic and gravitational fields are not dynamical. In order to do so, one promotes u µ and T to dynamical variables, and postulates the hydrodynamic equations in the form ∇ µ T µν = F νλ J λ , ∇ µ J µ = 0, with the leadingorder energy-momentum tensor given by eq. (5.4). Beyond the leading order, the energymomentum tensor becomes much more involved as discussed in section 4, plus the extra transport coefficients such as viscosity make their way into the hydrodynamic equations.
If the electromagnetic fields are dynamical, the conservation equations for T µν and J µ need to be supplemented by the evolution equations for the electromagnetic fields. These are usually taken to be Maxwell's equations, ∇ ν F µν = J µ . Substituting the equilibrium current (2.14) gives ∇ ν (F µν − M µν ) = ρu µ , (5.5) which is the standard covariant form of Maxwell's equations in matter, see e.g. [16]. In the right-hand side of (5.5), ρ is the density of free charges, while the effects of polarization are in the left-hand side. In the framework of (5.5), the derivative expansion in hydrodynamics can be implemented through the derivative expansion for M µν , however eq. (5.5) itself will receive corrections, e.g. due to the electrical conductivity. We plan to return to the study of hydrodynamics of polarized relativistic matter in the future.