A Variational Approach to Resistive General Relativistic Two-Temperature Plasmas

1. Introduction

Two-temperature plasmas have been studied in astrophysical systems for nearly fifty years. Early work considered the formation of light nuclei in two temperature plasmas (the ion temperature being greater than the electrons) that could exist near relativistic astrophysical objects. Colgate [1,2] and, independently Hoyle and Fowler [3], looked at the synthesis of deuterium in a plasma (with ion temperature $T_i\sim {10}^{11}$ K) generated in shock waves produced by supernovae. Shapiro et al. [4] applied a two-temperature accretion disk model for Cygnus X-1 in order to produce the observed thermal emission temperatures of ${10}^9$ K and the observed X-ray spectrum above $8$ keV. More recently, Zhdankin et al. [5] looked at the role of extreme two-temperature plasmas in radiative relativistic turbulence, while Ohmura et al. [6] used simulations of two temperature magnetohydrodynamics to describe the propagation of semi-relativistic jets. Ryan et al. [7] have provided axisymmetric two-temperature general relativistic radiation magnetohydrodynamic simulations of the inner region of the accretion flow onto the supermassive black hole M87 while Meringolo et al. [8] have looked at two temperature plasmas in the context of special relativistic turbulence.

The literature on electron and ion plasmas shows there are many different scenarios under which two temperatures result, although whether or not the electrons are hotter than the ions is very much dependent on the particular scenario. In his classic text on plasmas and fusion reactions, Chen [9] writes that the positively charged ions can have a temperature which is different from that of its electrons even though they both have Maxwellian distributions. This is because the collision rate of the ions with themselves and the collision rate of electrons with themselves are much higher than that of electrons with ions. Kawazura et al. [10] argue that in a collisionless plasma heated through Alvenic turbulence electrons will be preferentially heated when magnetic energy density is greater than the thermal energy density, whereas it is the ions which are hotter when the energy densities are the other way around.

The problem with developing models of complex plasmas in dynamical spacetimes, particularly for numerical simulations, is the consistency of the approximations used. It is standard to develop the approximations by dropping terms based on scaling arguments. Any “inconsistencies” introduced in the process typically lead to some (often small) loss of total energy or generation of spurious heat. However, as discussed in detail below, in a relativistic context, heat will produce an effective mass which contributes to the dynamics of a given system and (at least in principle) the generation of gravitational waves. Therefore, even small inconsistencies in the model development will lead to systematic errors in the generated (potentially observable) signals.

Our purpose here is to use well-established action-based techniques [11] to construct the full suite of field equations for a consistent, resistive, two-fluid, five-constituent, two-temperature general relativistic plasma. The model involves a positively charged species flux comoving with a charge-neutral species and a separate negatively charged species flux. The positively and neutrally charged species are assumed to have the same temperature and there is a single entropy comoving with them. Because the negatively charged species is at a different temperature, it will have its own (comoving) entropy.

To see how this comes about, consider the simple case of ionized hydrogen, for which collective behavior of the electrons means they can be described as a fluid. They have well-defined fluid elements with their own four-velocities, and within these elements there will be a thermodynamic description based on, say, temperature and particle density. Clearly, this assumes that the electrons are thermalized, i.e. from a kinetic theory point-of-view their state can be described by an equilibrium distribution function (say, Maxwell). From that same kinetic theory point-of-view, we know that entropy is calculable from the distribution. All of this is also true for the protons, except that the difference in temperature would necessarily lead to a different (maybe not in form, but certainly in specific values) distribution and hence different values for the entropy. Since the electrons are at equilibrium among themselves, and likewise for the protons, the electron entropy flows along with the electrons and the proton entropy flows along with the protons; therefore, because the electrons flow relatively to the protons, there are two entropy fluxes. It is conceptually straightforward to allow for ionization/recombination, by adding an additional flux of “neutral” particles. This leads to particle flux creation rates for both of the charged particle fluxes as well as the neutral particle flux. Conservation of baryon number will of course link these two rates.

Given that the physical system considered is broad, and readers may have different backgrounds—plasma physics, astrophysics, numerical relativity, and so on—we have tried to make this presentation as self-contained as possible. For example, there is an extended discussion of the so-called $3+1$ approach to General Relativity. We have attempted to make this a basic exercise in projecting tensors into spacelike hypersurfaces, or onto the normals to these hypersurfaces. Moreover, in order to set-up the taking of the Newtonian limit (in Sec. Section 6), it is advantageous to keep G, c, the magnetic permeability ${\mu }_{\mathrm{o}}$, and $k_B$ in the equations. Of course, this involves introducing a set of conventions, which are initially somewhat arbitrary, but eventually self-consistent. The complexity of our total system, with its mixing of dynamical, electrodynamical, and thermodynamical energies, fluxes, and momenta, requires a careful, yet admittedly tedious, dimensional analysis of the field variables. The relevant dimensions of field variables will be discussed as the variables are introduced. This is also required for taking the Newtonian limit, where we need to have an internal calibration of what “small” is when we expand the field equations.

The plan of this effort is as follows: In Sec. Section 2 the field variables are introduced, as well as some of their kinematical features. In Sec. Section 3 the “matter space” [11,12] is introduced as it provides the arena in which fluid displacements are performed in the action principle. In Sec. Section 4 we give the independent pieces of the action principle and derive the field equations. In Sec. Section 5 we give an overview of the $3+1$ formalism, focusing on the geometric arguments, and then apply it to the coupled system of general relativistic plasmas and electromagnetism. The overview is for the reader who is knowledgeable about plasma physics but not particularly familiar with numerical relativity, and/or with how to take a generally covariant theory and introduce a global separation of space from time. We follow this up in Sec. Section 5.4 with a review of the arguments given in [13] for building simple models of resistivity, for both the charged and neutral current and entropy flows. This is used in Sec. Section 6 where we take the Newtonian limit. In Sec. Section 7 we offer some concluding remarks. Adding further details, in Appendix A we review total charge conservation, in Appendix B we derive the “3 + 1” form of the Einstein equations, and in Appendix C we adapt the “3 + 1" formalism to a preferred coordinate system. The conventions of Misner, Thorne, and Wheeler [14] are used throughout (although we use $a,b,c,...$ rather than Greek letters to represent spacetime indices). We assume that the metric $g_{ab}$ is dimensionless, the coordinates carry the unit of length l, and the time unit is given by $l/c$; e.g. the time-coordinate $x^0=ct$. As one might expect, the notation will quickly become a nightmare, and so notational conventions will be explained as the story develops.

2. The Plasma State and the Field Variables

The first step towards modelling a plasma system involves understanding the scales involved and the relevant variables. Perhaps the most important scale is the Debye length ${\lambda }_D$, which is given by [15]

$$\frac{1}{{\lambda }_D^2}=\sum _i\frac{n_i{q}_i^2}{{ϵ}_{\mathrm{o}}{T}_i},$$

(1)

where $n_i$ is the number density of the $i^{\mathrm{th}}$–species, $q_i$ its charge, and $T_i$ its temperature. The Debye length is the effective distance at which the influence of a single charge is no longer felt by other particles; that is, for a length-scale l, somewhere between the inter-particle separation $1/n_i^{1/3}$ and ${\lambda }_D$, polarization (or collective) effects will occur so that charges outside of the Debye sphere (area $\propto {\lambda }_D^2$) are shielded from the single charge. For scales L much bigger than ${\lambda }_D$, the system will exhibit fluid-like features, such as wave propagation.

This helps establish criteria through which we can define the plasma state: 1) the typical length-scale L for the system must be much larger than the Debye length—$L\gg {\lambda }_D$—and such that quasi-neutrality holds (${\sum }_i{q}_i{n}_i{L}^3\approx 0$);1 2) there must be a large enough number of particles in the Debye sphere that collective effects occur so that the shielding takes hold ($n_i{\lambda }_D^3\gg 1$); and 3) letting $\tau$ represent the mean collision time for the neutral particles and $1/\omega$ a time-scale for collective plasma phenomena, we have that the last criterion is $\omega \tau \gg 1$.

In a system like an accretion disc around a black hole there can be several length scales—the horizontal reach L of the disc, the size $2G{M}_{BH}/c^2$ of the black hole with total mass $M_{BH}$, and so on. A satisfactory fluid model of the matter and heat in the disc exists when the system can be broken up into a continuum of “boxes” of volume $l^3$, each of which is small enough that it can be considered as being microscopic with respect to the system as a whole ($l/L\ll 1$), and yet large enough that it contains enough particles N for which the Laws of Thermodynamics hold. In this case, intensive quantities such as chemical potential, pressure, and temperature will be well defined [16].

In the limit where l becomes infinitesimal, these conceptual boxes become the fluid elements of fluid models. As the fluid evolves, the fluid elements will trace out a continuum of worldlines in spacetime; i.e. smooth curves whose spacetime points are identified by a set of coordinates $x^a\left(\tau \right)$, with $\tau$ being the proper time along the curves. Because the fluid elements contain particles, then these curves form the basis for tracking particle flux. It is important to note that since a fluid element is infinitesimal with respect to the system as a whole, then changes in the gravitational field across it are negligible. The equivalence principle also implies that the local geometry can be treated as flat spacetime.

Particle flux is defined in the standard way as being a number of particles N passing through an area $l^2$ per some time $l/c$; i.e., particle flux magnitude is $\left(N/l^3\right)c$. We do the same for entropy flux, except to note that the entropy unit is $k_B$, which is energy e per temperature T. Assuming that we can count the amount of entropy as some number $N_s$ times $k_B$, then the entropy flux will be $N_s$ units of entropy passing through area $l^2$ per time $l/c$; i.e., entropy flux magnitude is $\left(N_s/l^3\right)c$.2

Our system consists of a neutrally charged species ($q_{\eta }=0$) with particle flux $n_{\eta }^a$ and a comoving entropy flux $s_{\overline{\eta }}^a/k_B$; a positively charged species ($q_{\mathcal{P}}>0$) with particle flux $n_{\mathcal{P}}^a$ and a comoving entropy flux $s_{\overline{\mathcal{P}}}^a/k_B$; and a negatively charged species ($q_{\mathcal{N}}=-q_{\mathcal{P}}$) with particle flux $n_{\mathcal{N}}^a$ and a comoving entropy flux $s_{\overline{\mathcal{N}}}^a/k_B$. As we will see later, associated with the particle fluxes $\{n_{\eta }^a,n_{\mathcal{P}}^a,n_{\mathcal{N}}^a\}$ are, respectively, canonically conjugate chemical potential covectors $\{{\mu }_a^{\eta },{\mu }_a^{\mathcal{P}},{\mu }_a^{\mathcal{N}}\}$ [cf.Eq. (24)] and for the entropy fluxes $\{s_{\overline{\eta }}^a/k_B,s_{\overline{\mathcal{P}}}^a/k_B,s_{\overline{\mathcal{N}}}^a/k_B\}$ there are respective canonically conjugate “temperature” covectors $\{k_B{\Theta }_a^{\overline{\eta }},k_B{\Theta }_a^{\overline{\mathcal{P}}},k_B{\Theta }_a^{\overline{\mathcal{N}}}\}$.

At this point, it is convenient to simplify the notation, by introducing constituent indices $\{\mathrm{x},\mathrm{y},\dots \}$ which will take the values $\mathrm{x}=1,2,\dots ,6$. With these, we will write generic particles fluxes $n_{\mathrm{x}}^a$ such that the first three are $\{n_1^a=n_{\eta }^a,n_2^a=n_{\mathcal{P}}^a,n_3^a=n_{\mathcal{N}}^a\}$, and the next three are $\{n_4^a=s_{\overline{\eta }}^a/k_B,n_5^a=s_{\overline{\mathcal{P}}}^a/k_B,n_6^a=s_{\overline{\mathcal{N}}}^a/k_B\}$. For the canonically conjugate covectors we will identify $\{{\mu }_a^1={\mu }_a^{\eta },{\mu }_a^2={\mu }_a^{\mathcal{P}},{\mu }_a^3={\mu }_a^{\mathcal{N}}\}$ and $\{{\mu }_a^4=k_B{\Theta }_a^{\overline{\eta }},{\mu }_a^5=k_B{\Theta }_a^{\overline{\mathcal{P}}},{\mu }_a^6=k_B{\Theta }_a^{\overline{\mathcal{N}}}\}$. In order to make direct contact with the First and Second Laws of Thermodynamics we use an energy e to assign to the combination ${\mu }_a^{\mathrm{x}}{n}_{\mathrm{x}}^a$ energy density units $e/l^3$. This implies that the ${\mu }_a^{\mathrm{x}}$ must have momentum units $e/c$. The energy e can take two distinct forms: a particle energy based on mass-energy, $e_m=m{c}^2$, for the set $\{{\mu }_a^1,{\mu }_a^2,{\mu }_a^3\}$, and a thermal energy $e_T=k_{B}T$ for the set $\{{\mu }_a^4,{\mu }_a^5,{\mu }_a^6\}$.

The density $n_{\mathrm{x}}$, with units $N/l^3$, associated with the flux $n_{\mathrm{x}}^a$ allows us to define a four-velocity field $u_{\mathrm{x}}^a=n_{\mathrm{x}}^a/n_{\mathrm{x}}$, which is normalized such that $g_{ab}{u}_{\mathrm{x}}^a{u}_{\mathrm{x}}^b=-c^2$. These flux worldlines are tied to those of the fluid elements by setting $u_{\mathrm{x}}^a=\mathrm{d}{x}_{\mathrm{x}}^a/\mathrm{d}{\tau }_{\mathrm{x}}$, where ${\tau }_{\mathrm{x}}$ is the proper time along the worldline traced out by $u_{\mathrm{x}}^a$. We see that $n_{\mathrm{x}}=-u_a^{\mathrm{x}}{n}_{\mathrm{x}}^a/c^2$ or $n_{\mathrm{x}}^2=-g_{ab}{n}_{\mathrm{x}}^a{n}_{\mathrm{x}}^b/c^2$. Note that in addition to the $n_{\mathrm{x}}^2$ we can have the mixed terms $n_{\mathrm{x}\mathrm{y}}^2=-g_{ab}{n}_{\mathrm{x}}^a{n}_{\mathrm{y}}^b/c^2=n_{\mathrm{y}\mathrm{x}}^2$, where it is to be understood that $\mathrm{x}\ne \mathrm{y}$.3 With respect to a flux’s rest-frame, i.e. the local frame which follows the worldline given by $u_{\mathrm{x}}^a$, we can define the fluid potentials ${\mu }_{\mathrm{x}}=-u_{\mathrm{x}}^a{\mu }_a^{\mathrm{x}}$. For $\mathrm{x}=1,2,3$, the ${\mu }_{\mathrm{x}}$ are chemical potentials, and for $\mathrm{x}=4,5,6$ the ${\mu }_{\mathrm{x}}$ are temperatures ${\mu }_4=T_{\overline{\eta }}$, ${\mu }_5=T_{\overline{\mathcal{P}}}=T_{\overline{\eta }}$, and ${\mu }_6=T_{\overline{\mathcal{N}}}\ne T_{\overline{\mathcal{P}}}$.

The remaining field variables are the four-vector potential $A_a$ and the spacetime metric $g_{ab}$. The metric couples all fields to the spacetime curvature (and vice versa). With $A_a$ and the charge density flux $j_{\mathrm{x}}^a=q_{\mathrm{x}}{n}_{\mathrm{x}}^a$ we can couple the charged fluids to the electromagnetic field (and vice versa). The total charge density flux is

$$j^a=\sum _{\mathrm{x}}{q}_{\mathrm{x}}{n}_{\mathrm{x}}^a=j_{\mathcal{P}}^a+j_{\mathcal{N}}^a.$$

(2)

The units of the charged current flux $j_{\mathrm{x}}^a$ are $\left(qN/l^3\right)c$. We note that MKS units are being used so that the electromagnetic coupling ${\mu }_{\mathrm{o}}$ combines with ${ϵ}_{\mathrm{o}}$ to give ${ϵ}_{\mathrm{o}}{\mu }_{\mathrm{o}}=1/c^2$. The four-potential $A_a$ has units of momentum per charge, or $e_{EM}/\left(qc\right)$, where $e_{EM}$ is a characteristic electromagnetic energy; for example, in the Debye limit case we would use $e_{EM}\sim q{\tilde{V}}_{\mathrm{eff}}$.

3. The Matter Space Approach to Dissipation

Our analysis builds on a well-established variational approach to relativistic multi-fluid dynamics [11], including dissipative aspects. The main fluid fields in the model are the fluxes $n_{\mathrm{x}}^a$. At the heart of the fluxes are the four-velocities $u_{\mathrm{x}}^a=\mathrm{d}{x}_{\mathrm{x}}^a/\mathrm{d}{\tau }_{\mathrm{x}}$. In general, the $u_{\mathrm{x}}^a$ are not surface forming, but they do form a fibration of spacetime. If the $u_{\mathrm{x}}^a$ are given, then $\mathrm{d}{x}_{\mathrm{x}}^a/\mathrm{d}{\tau }_{\mathrm{x}}=u_{\mathrm{x}}^a$ can be integrated so as to construct the $x_{\mathrm{x}}^a\left({\tau }_{\mathrm{x}}\right)$. Since $u_{\mathrm{x}}^a{u}_a^{\mathrm{x}}=-c^2$, then knowing, say, the three spatial pieces $\mathrm{d}{x}_{\mathrm{x}}^i/\mathrm{d}{\tau }_{\mathrm{x}}$, automatically determines the time piece $\mathrm{d}{x}_{\mathrm{x}}^0/\mathrm{d}{\tau }_{\mathrm{x}}$. For some given spacelike hypersurface, no two worldlines of, say, the $x^{\mathrm{th}}$-fluid, will intersect that hypersurface at the same point.

If we think of this surface in the context of an initial-value problem, then each worldline will be uniquely determined by the three spatial coordinates they have on that initial hypersurface. It is through this that the so-called “matter space”, or pull-back, approach enters the fluid dynamics. We replace the initial spacelike hypersurface, with an abstract, three-dimensional space endowed with coordinates $X_{\mathrm{x}}^A$ (having dimensions l and $A=1,2,3$). Instead of each worldline being identified with a point on the initial spacelike hypersurface, each point $x_{\mathrm{x}}^a\left(\tau \right)$ on the worldline gets mapped to the same point $X_{\mathrm{x}}^A$ in the matter space. Our goal here is to provide a sketch on how to reformulate our fluid model so that the $X_{\mathrm{x}}^A$ are the fundamental fields (see, e.g., Andersson and Comer [11] for complete details).

The first step in this reformulation is to introduce the three-form $n_{abc}^{\mathrm{x}}$, which is dual to $n_{\mathrm{x}}^a$:

$$n_{abc}^{\mathrm{x}}={ϵ}_{dabc}{n}_{\mathrm{x}}^d,n_{\mathrm{x}}^a=\frac{1}{3!}{ϵ}^{bcda}{n}_{bcd}^{\mathrm{x}},$$

(3)

where our convention for transforming between the two is

$${ϵ}^{bcda}{ϵ}_{ebcd}=3!{\delta }_e^a.$$

(4)

Likewise, we introduce the three-form ${\mu }_{\mathrm{x}}^{abc}$ which is dual to ${\mu }_a^{\mathrm{x}}$:

$${\mu }_{\mathrm{x}}^{abc}={ϵ}^{dabc}{\mu }_d^{\mathrm{x}},{\mu }_a^{\mathrm{x}}=\frac{1}{3!}{ϵ}_{bcda}{\mu }_{\mathrm{x}}^{bcd}.$$

(5)

Because the metric is dimensionless, we see that the three-forms carry the same units as their dual vectors.

We use the map associated with the coordinates $X_{\mathrm{x}}^A$ of the ${\mathrm{x}}^{\mathrm{th}}$-fluid’s matter space to “pullback” $n_{abc}^{\mathrm{x}}$ into the matter space where it is identified with the totally antisymmetric tensor $n_{ABC}^{\mathrm{x}}$:

$$n_{abc}^{\mathrm{x}}={}^{\mathrm{x}}{\mathcal{J}}_{abc}^{ABC}{n}_{ABC}^{\mathrm{x}},$$

(6)

such that the Einstein convention applies to repeated matter space indices, and

$${}^{\mathrm{x}}{\mathcal{J}}_{abc}^{ABC}=\frac{\partial X_{\mathrm{x}}^{[A}}{\partial x^a}\frac{\partial X_{\mathrm{x}}^B}{\partial x^b}\frac{\partial X_{\mathrm{x}}^{C]}}{\partial x^c}.$$

(7)

We also use the map associated with $X_{\mathrm{x}}^A$ to “push-forward” the fully antisymmetric matter space quantity ${\mu }_{\mathrm{x}}^{ABC}$ to the spacetime three-form ${\mu }_{\mathrm{x}}^{abc}$, via

$${\mu }_{\mathrm{x}}^{ABC}={}^{\mathrm{x}}{\mathcal{J}}_{abc}^{ABC}{\mu }_{\mathrm{x}}^{abc},$$

(8)

as well as the symmetric matter space “metric” $g_{\mathrm{x}}^{AB}$ to the spacetime metric $g_{ab}$, via

$$g_{\mathrm{x}}^{AB}=\frac{\partial X_{\mathrm{x}}^A}{\partial x^a}\frac{\partial X_{\mathrm{x}}^B}{\partial x^b}{g}^{ab}.$$

(9)

Because of the antisymmetry in the indices of $n_{ABC}^{\mathrm{x}}$ and ${\mu }_{\mathrm{x}}^{ABC}$ there are natural definitions for the volume-form ${ϵ}_{ABC}^{\mathrm{x}}$ and its inverse ${ϵ}_{\mathrm{x}}^{ABC}$ on the $\mathrm{x}$-matter space. These satisfy [13,16]

$${ϵ}_{DEF}^{\mathrm{x}}{ϵ}_{\mathrm{x}}^{ABC}=3!{\delta }_D^{[A}{\delta }_E^B{\delta }_F^{C]}⟹{ϵ}_{ABC}^{\mathrm{x}}{ϵ}_{\mathrm{x}}^{ABC}=3!.$$

(10)

We can normalize ${ϵ}_{ABC}^{\mathrm{x}}$ and ${ϵ}_{\mathrm{x}}^{ABC}$ using the determinant of $g_{\mathrm{x}}^{AB}$; i.e.

$${ϵ}_{123}^{\mathrm{x}}=\frac{1}{{ϵ}_{\mathrm{x}}^{123}}=\frac{1}{\sqrt{{\Delta }_{\mathrm{x}}}},$$

(11)

where

$${\Delta }_{\mathrm{x}}=\frac{1}{3!{\left({ϵ}_{123}^{\mathrm{x}}\right)}^2}{ϵ}_{ABC}^{\mathrm{x}}{ϵ}_{DEF}^{\mathrm{x}}{g}_{\mathrm{x}}^{AD}{g}_{\mathrm{x}}^{BE}{g}_{\mathrm{x}}^{CF}.$$

(12)

Now we can write

$$n_{ABC}^{\mathrm{x}}={\mathcal{N}}^{\mathrm{x}}{ϵ}_{ABC}^{\mathrm{x}},{\mathcal{N}}^{\mathrm{x}}=\frac{1}{3!}{ϵ}_{\mathrm{x}}^{ABC}{n}_{ABC}^{\mathrm{x}},$$

(13)

where it can be shown that ${\mathcal{N}}^{\mathrm{x}}=n_{\mathrm{x}}$ [16]. Similarly, we find

$${\mu }_{\mathrm{x}}^{ABC}={\mathcal{M}}_{\mathrm{x}}{ϵ}_{\mathrm{x}}^{ABC},{\mathcal{M}}_{\mathrm{x}}=\frac{1}{3!}{ϵ}_{ABC}^{\mathrm{x}}{\mu }_{\mathrm{x}}^{ABC},$$

(14)

where it can be shown that ${\mathcal{M}}_{\mathrm{x}}={\mu }_{\mathrm{x}}$.

It is also straightforward to confirm that

$$u_{\mathrm{x}}^a=\frac{1}{3!}{ϵ}^{bcda}{}^{\mathrm{x}}{\mathcal{J}}_{bcd}^{ABC}{ϵ}_{ABC}^{\mathrm{x}}.$$

(15)

From this we can verify that the $X_{\mathrm{x}}^A$ are conserved along their own worldlines (i.e. they are Lie-dragged by their $u_{\mathrm{x}}^a$); that is, using Eq. (15), we see

$$\frac{\mathrm{d}{X}_{\mathrm{x}}^A}{\mathrm{d}{\tau }_{\mathrm{x}}}=u_{\mathrm{x}}^a{\nabla }_a{X}_{\mathrm{x}}^A=\frac{1}{n_{\mathrm{x}}}\left(-\frac{1}{3!}{ϵ}^{abcd}\frac{\partial X_{\mathrm{x}}^A}{\partial x^{[a}}\frac{\partial X_{\mathrm{x}}^B}{\partial x^b}\frac{\partial X_{\mathrm{x}}^C}{\partial x^c}\frac{\partial X_{\mathrm{x}}^D}{\partial x^{d]}}\right)n_{BCD}^{\mathrm{x}}\equiv 0,$$

(16)

since the term in parentheses vanishes identically. The quantity ${\nabla }_a$ is the covariant derivative, with the dimension of inverse length $1/l$.

In general, dissipation is directly connected with the (matter and/or entropy) particle flux creation rate ${\Gamma }_{\mathrm{x}}$, which is given by

$${\Gamma }_{\mathrm{x}}={\nabla }_a{n}_{\mathrm{x}}^a.$$

(17)

When ${\Gamma }_{\mathrm{x}}=0$ there is no flux change and no dissipation. It is easy to see that there is a one-to-one, local identification of the divergence of a vector field with the exterior derivative of its associated three-form, i.e.${\nabla }_{[a}{n}_{bcd]}^{\mathrm{x}}$; namely,

$${\nabla }_a{n}_{\mathrm{x}}^a=\frac{1}{3!}{ϵ}^{abcd}{\nabla }_{[a}{n}_{bcd]}^{\mathrm{x}}.$$

(18)

Simply put, if the three-form is closed (e.g.${\nabla }_{[a}{n}_{bcd]}^{\mathrm{x}}=0$), then ${\nabla }_a{n}_{\mathrm{x}}^a=0$ and there is no dissipation; if the three-form is not closed (e.g.${\nabla }_{[a}{n}_{bcd]}^{\mathrm{x}}\ne 0$), then the divergence is not zero and dissipation can occur.

This is the lynchpin of the formalism for dissipative multi-fluid systems developed by Andersson and Comer [17], and another reason for invoking the matter space. In fact, it was shown by Celora et al. [16] that

$${\mu }_{\mathrm{x}}{\Gamma }_{\mathrm{x}}=\frac{1}{3!}{\mu }_{\mathrm{x}}^{ABC}\frac{\mathrm{d}}{\mathrm{d}{\tau }_{\mathrm{x}}}{n}_{ABC}^{\mathrm{x}}.$$

(19)

We see immediately that if $n_{ABC}^{\mathrm{x}}$ is a function of only the $X_{\mathrm{x}}^A$, then ${\Gamma }_{\mathrm{x}}=0$ because of Eq. (16). This is ideal when fluids are non-dissipative, because then their respective creation rates must vanish. However, if we allow $n_{ABC}^{\mathrm{x}}$ to also depend on $X_{\mathrm{y}}^A$ (for $\mathrm{y}\ne \mathrm{x}$), then the flux three-form is no longer closed and a system of fluid equations with resistive forms of dissipation4 result [13]. This will be shown later in Sec. Section 4.2.

4. The Action Principle and Field Equations

We now set up the action principle used to derive the resistive fluid/plasma, Maxwell, and Einstein set of field equations.5 The pull-back formalism will be used to build unconstrained variations of the fluid fluxes $\delta n_{\mathrm{x}}^a$ so that the fluid equations can be obtained. The Maxwell equations follow from variations of $A_a$, which appears in two pieces of the total action: one built from the antisymmetric Faraday tensor $F_{ab}$ defined as

$$F_{ab}={\nabla }_a{A}_b-{\nabla }_b{A}_a,$$

(20)

and the other constructed from a coupling term based on the scalar $j_{\mathrm{x}}^a{A}_a$. It is important to note that $F_{ab}$ satisfies a “Bianchi” identity

$${\nabla }_a{F}_{bc}+{\nabla }_c{F}_{ab}+{\nabla }_b{F}_{ca}=0⟹\frac{1}{2}{ϵ}^{abcd}{\nabla }_{[b}{F}_{cd]}=0,$$

(21)

The Faraday tensor has dimensions $e_{EM}/\left(qcl\right)$.

Gravity is incorporated (in the standard way) by using the Einstein-Hilbert action for the Einstein Equation and by the minimal coupling of the metric $g_{ab}$ to the fluid and electromagnetic fields. The minimal coupling arises from the $\sqrt{-g}$ term in spacetime volume integrals, where g is the determinant of the metric; the use of $g_{ab}$ in the inner product of vectors; and replacing partial derivatives with covariant derivatives. The energy-momentum-stress tensor $T_{ab}$, with energy density units $e/l^3$, is obtained in the usual way by varying the total action with respect to $g_{ab}$.

4.1. The Matter, Electromagnetic, Coupling, and Gravity Actions

The fluid action $S_M$ uses for its Lagrangian the so-called Master function $\Lambda$ [11], an energy density, which is a functional of all the $n_{\mathrm{x}}^2$ and $n_{\mathrm{x}\mathrm{y}}^2$. An arbitrary variation of $S_M$ with respect to the flux $n_{\mathrm{x}}^a$ and the metric results in

$$\begin{gathered} \begin{array}{ccc}\hfill \delta S_M& =& \delta \left({\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}\Lambda \right)\hfill \end{array} \\ \begin{array}{ccc}& =& {\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}\left[\sum _{\mathrm{x}}{\mu }_a^{\mathrm{x}}\delta n_{\mathrm{x}}^a+\frac{1}{2}\left(\Lambda g^{ab}+\sum _{\mathrm{x}}{n}_{\mathrm{x}}^a{\mu }_{\mathrm{x}}^b\right)\delta g_{ab}\right],\hfill \end{array} \end{gathered}$$

(22)

where

$${\mu }_a^{\mathrm{x}}={\mathcal{B}}^{\mathrm{x}}{n}_a^{\mathrm{x}}+\sum _{\mathrm{y}\ne \mathrm{x}}{\mathcal{A}}^{\mathrm{x}\mathrm{y}}{n}_a^{\mathrm{y}},$$

(23)

and

$$\begin{array}{cc}\hfill {\mathcal{B}}^{\mathrm{x}}& =-\frac{2}{c^2}\frac{\partial \Lambda }{\partial n_{\mathrm{x}}^2},\hfill \end{array}$$

(24a)

$$\begin{array}{cc}\hfill {\mathcal{A}}^{\mathrm{x}\mathrm{y}}& =-\frac{1}{c^2}\frac{\partial \Lambda }{\partial n_{\mathrm{x}\mathrm{y}}^2}.\hfill \end{array}$$

(24b)

The ${\mathcal{A}}^{\mathrm{x}\mathrm{y}}$, with units $\left(l^3/N\right)e/c^2$, provide the “entrainment” effect, which causes the fluid momenta to be “tilted” in the sense that ${\mu }_a^{\mathrm{x}}$ is not proportional to its corresponding flux $n_{\mathrm{x}}^a$. The implication is that one flux, say $n_{\mathcal{P}}^a$, carries along with it a fraction of the components of a different flux, say $n_{\mathcal{N}}^a$. This leads also to effective “mass” effects due to entrainment between any two particle fluxes, a particle flux and an entropy flux, or two entropy fluxes. Entropy flux acquires an effective mass6 (a carrier of inertia which scales like $k_{B}T/c^2$) through its (non-dissipative) energy/heat exchange within the system, which does work and can change the conjugate momenta of other fluxes [19]. Shatashvili et al. [20] have included electron effective masses in their two temperature plasma equations. It has been noted by Kotorashvili et al. [21] that the effective mass for a degenerate electron plasma arises from the degeneracy instead of kinematics and is fully determined by the plasma rest frame density (see [22] and references therein), whereas in a hot relativistic electron plasma the effective mass [23] is determined by the relativistic electron temperature.

Entrainment between neutrons and protons is known to be important in superfluid neutron star dynamics [24,25,26,27]. Entrainment between matter and entropy can be shown (see, for example, [19]) to lead to the Cattaneo equation [28], which is an important component of causal heat conductivity. This particle and entropy flux model can also be used to describe superfluid systems such as He${}^4$. In the Landau model of superfluidity [29], there is an ad hoc separation of the He${}^4$ atoms into a superfluid particle flux and a normal fluid particle flux, which are entrained with each other. In the entropy and particle flux approach, all of the He${}^4$ atoms are described with one particle flux, and the “normal fluid” flux is replaced with an entropy flux. A one-to-one mapping between the two models exists (see, for example, Andersson and Comer [30], and references therein), primarily because in the Landau model the normal fluid represents the excitations of atoms out of the ground state and are responsible for carrying the heat. This is important because it shows that the entrainment between the entropy and particle fluxes has physical impact, whether it is describing superfluid He${}^4$ or more general fluids with an independent heat flow. It is less clear whether entrainment between two entropies is important physically, or just a formally consistent piece of the overall mathematical construct.

4.1.1. The Electromagnetic and Coupling Actions

The Maxwell Action is

$$S_{Max}=\frac{1}{4{\mu }_{\mathrm{o}}}{\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}{F}_{ab}{F}^{ab},$$

(25)

and its variation with respect to $A_a$ and the metric $g_{ab}$ leads to

$$\delta S_{Max}=\frac{1}{{\mu }_{\mathrm{o}}}{\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}\left({\nabla }_a{F}^{ab}\right)\delta A_b-\frac{1}{8{\mu }_{\mathrm{o}}}{\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}\left(F_{cd}{F}^{cd}{g}^{ab}-4F^{ac}{F}^b{}_c\right)\delta g_{ab}.$$

(26)

The minimal coupling of the Maxwell field to the charge current densities $j_{\mathrm{x}}^a$ is obtained from the Coulomb action

$$S_C={\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}\left(\sum _{\mathrm{x}}{j}_{\mathrm{x}}^a\right)A_a,$$

(27)

whose variation with respect to $n_{\mathrm{x}}^a$, $A_a$, and $g_{ab}$ is

$$\delta S_C={\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}\sum _{\mathrm{x}}\left(j_{\mathrm{x}}^a\delta A_a+q_{\mathrm{x}}{A}_a\delta n_{\mathrm{x}}^a+\frac{1}{2}{j}_{\mathrm{x}}^a{A}_a{g}^{bc}\delta g_{bc}\right).$$

(28)

4.1.2. The Gravitational Einstein-Hilbert Action

At the heart of General Relativity is the Riemann tensor $R^c{}_{dab}$, with units of $1/l^2$. It can be inferred from the antisymmetric operation of two covariant derivatives on an arbitrary vector $v^a$; namely,

$${\nabla }_a{\nabla }_b{v}^c-{\nabla }_b{\nabla }_a{v}^c=R^c{}_{dab}{v}^d.$$

(29)

From the Riemann tensor we can obtain the Ricci tensor $R^{ab}=R^c{}_{acb}$ and, subsequently, the Ricci scalar $R=g_{ab}{R}^{ab}$.

The Einstein-Hilbert action is

$$S_{EH}=\frac{c^4}{16\pi G}{\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}R.$$

(30)

Varying it and the other bits of the total action written above with respect to the metric gives the Einstein equation; in particular, the left-hand-side of the Einstein equations comes from the variation of $S_{EH}$ with respect to $g_{ab}$, i.e.

$$\delta S_{EH}=-\frac{c^4}{16\pi G}{\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}{G}^{ab}\delta g_{ab},$$

(31)

where the Einstein tensor $G^{ab}$ is

$$G^{ab}=R^{ab}-\frac{1}{2}R{g}^{ab}.$$

(32)

4.1.3. The Total Action Variation

The variation of the total action S for the system is thus

$$\begin{gathered} \begin{array}{ccc}\hfill \delta S& =& \delta S_{EH}+\delta S_M+\delta S_{Max}+\delta S_C\hfill \end{array} \\ \begin{array}{ccc}& =& {\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}\left\{-\frac{c^4}{16\pi G}{G}^{ab}\delta g_{ab}+\sum _{\mathrm{x}}{\overline{\mu }}_a^{\mathrm{x}}\delta n_{\mathrm{x}}^a+\frac{1}{{\mu }_{\mathrm{o}}}\left({\nabla }_b{F}^{ba}+{\mu }_{\mathrm{o}}\sum _{\mathrm{x}}{j}_{\mathrm{x}}^a\right)\delta A_a\right.\hfill \end{array} \\ \begin{array}{ccc}& & \left.+\frac{1}{2}\left[\Lambda g^{ab}+\sum _{\mathrm{x}}\left(n_{\mathrm{x}}^a{\mu }_{\mathrm{x}}^b+j_{\mathrm{x}}^c{A}_c{g}^{ab}\right)-\frac{1}{4{\mu }_{\mathrm{o}}}\left(F_{cd}{F}^{cd}{g}^{ab}-4F^{ac}{F}^b{}_c\right)\right]\delta g_{ab}\right\},\hfill \end{array} \end{gathered}$$

(33)

where the electromagnetic minimal coupling has caused the fluid conjugate momentum to become

$${\overline{\mu }}_a^{\mathrm{x}}={\mu }_a^{\mathrm{x}}+q_{\mathrm{x}}{A}_a.$$

(34)

Imposing gauge invariance on the total action S (cf.Appendix A) leads to charge conservation in the form [cf.Eq. (A.6)]

$$q_{\mathcal{P}}{\Gamma }_{\mathcal{P}}+q_{\mathcal{N}}{\Gamma }_{\mathcal{N}}=0⟹{\Gamma }_{\mathcal{N}}={\Gamma }_{\mathcal{P}},$$

(35)

where ${\Gamma }_{\mathcal{P}}={\nabla }_a{n}_{\mathcal{P}}^a$ and ${\Gamma }_{\mathcal{N}}={\nabla }_a{n}_{\mathcal{N}}^a$. Of course, there is also baryon number conservation. The total baryon number flux is $n_B^a=n_{\eta }^a+n_{\mathcal{P}}^a$, and it is conserved if ${\Gamma }_B={\nabla }_a{n}_B^a=0$; therefore,

$$0={\nabla }_a{n}_B^a={\nabla }_a{n}_{\eta }^a+{\nabla }_a{n}_{\mathcal{P}}^a\equiv {\Gamma }_{\eta }+{\Gamma }_{\mathcal{P}}⟹{\Gamma }_{\eta }=-{\Gamma }_{\mathcal{P}}.$$

(36)

The field equations obtained from the full action variation above cannot be the final form, since the term proportional to $\delta n_{\mathrm{x}}^a$ implies that the momenta ${\overline{\mu }}_a^{\mathrm{x}}$ must vanish. This happens because the components of $\delta n_{\mathrm{x}}^a$ cannot all be varied independently; this is the main reason for using the pull-back formalism because it provides a set of variables, the $X_{\mathrm{x}}^A$, which can be varied independently.

4.2. From Matter Space to Spacetime Displacements and Resistivity

Even though we have as our unconstrained dynamical variables the scalars $X_{\mathrm{x}}^A$, ultimately we want the action principle to produce field equations for the fluxes $n_{\mathrm{x}}^a$. Fortunately, we can use the $X_{\mathrm{x}}^A$ this time to push-forward variations $\delta X_{\mathrm{x}}^A$ in matter space to Lagrangian displacements ${\xi }_{\mathrm{x}}^a$ of fluid element worldlines on spacetime; namely,

$$\delta X_{\mathrm{x}}^A=-\frac{\partial X_{\mathrm{x}}^A}{\partial x^a}{\xi }_{\mathrm{x}}^a,$$

(37)

where $\delta X_{\mathrm{x}}^A$ is an Eulerian variation (when the $X_{\mathrm{x}}^A$ are taken as scalars on spacetime). The minus sign comes in because we know that the $X_{\mathrm{x}}^A$ do not change along the fluid worldlines, meaning that their Lagrangian variation ${\Delta }_{\mathrm{x}}{X}_{\mathrm{x}}^A$ [11] has to vanish:

$${\Delta }_{\mathrm{x}}{X}_{\mathrm{x}}^A\equiv \delta X_{\mathrm{x}}^A+{\mathcal{L}}_{{\xi }_{\mathrm{x}}}{X}_{\mathrm{x}}^A=0,$$

(38)

where ${\mathcal{L}}_{{\xi }_{\mathrm{x}}}$ is the Lie derivative with respect to ${\xi }_{\mathrm{x}}^a$. Since ${\Delta }_{\mathrm{x}}{X}_{\mathrm{x}}^A=0$ we arrive at Eq. (40). Note that, because we have several fluxes, we will need also the mixed Lagrangian variation ${\Delta }_{\mathrm{x}}{X}_{\mathrm{y}}^A$ of the $X_{\mathrm{y}}^A$ with respect to the $\mathrm{x}$-fluid (and vice versa):

$${\Delta }_{\mathrm{x}}{X}_{\mathrm{y}}^A=\delta X_{\mathrm{y}}^A+{\mathcal{L}}_{{\xi }_{\mathrm{x}}}{X}_{\mathrm{y}}^A={\mathcal{L}}_{{\xi }_{\mathrm{x}}}{X}_{\mathrm{y}}^A-{\mathcal{L}}_{{\xi }_{\mathrm{y}}}{X}_{\mathrm{y}}^A=\left({\xi }_{\mathrm{x}}^a-{\xi }_{\mathrm{y}}^a\right)\frac{\partial X_{\mathrm{y}}^A}{\partial x^a}.$$

(39)

The displacements of the matter space fluid elements will lead to the variation $\delta n_{ABC}^{\mathrm{x}}$, which, in turn, will induce the variation of $n_{abc}^{\mathrm{x}}$. The Lagrangian variation of $n_{abc}^{\mathrm{x}}$, in general, is

$${\Delta }_{\mathrm{x}}{n}_{abc}^{\mathrm{x}}={}^{\mathrm{x}}{\mathcal{J}}_{abc}^{ABC}{\Delta }_{\mathrm{x}}{n}_{ABC}^{\mathrm{x}},$$

(40)

and thus

$$\delta n_{abc}^{\mathrm{x}}=-{\mathcal{L}}_{{\xi }_{\mathrm{x}}}{n}_{abc}^{\mathrm{x}}+{}^{\mathrm{x}}{\mathcal{J}}_{abc}^{ABC}{\Delta }_{\mathrm{x}}{n}_{ABC}^{\mathrm{x}},$$

(41)

where the Lie derivative of $n_{abc}^{\mathrm{x}}$ along the ${\xi }_{\mathrm{x}}^a$ is

$${\mathcal{L}}_{{\xi }_{\mathrm{x}}}{n}_{abc}^{\mathrm{x}}={\xi }_{\mathrm{x}}^d\frac{\partial n_{abc}^{\mathrm{x}}}{\partial x^d}+n_{dbc}^{\mathrm{x}}\frac{\partial {\xi }_{\mathrm{x}}^d}{\partial x^a}+n_{adc}^{\mathrm{x}}\frac{\partial {\xi }_{\mathrm{x}}^d}{\partial x^b}+n_{abd}^{\mathrm{x}}\frac{\partial {\xi }_{\mathrm{x}}^d}{\partial x^c}.$$

(42)

The resistive form of dissipation is due to the presence of $X_{\mathrm{y}}^A$ in $n_{ABC}^{\mathrm{x}}$. Applying the definitions above, we see

$${\Delta }_{\mathrm{x}}{n}_{ABC}^{\mathrm{x}}=\sum _{\mathrm{y}\ne \mathrm{x}}\frac{\partial n_{ABC}^{\mathrm{x}}}{\partial X_{\mathrm{y}}^D}{\Delta }_{\mathrm{x}}{X}_{\mathrm{y}}^D=\sum _{\mathrm{y}\ne \mathrm{x}}\frac{\partial n_{ABC}^{\mathrm{x}}}{\partial X_{\mathrm{y}}^D}\frac{\partial X_{\mathrm{y}}^D}{\partial x^a}\left({\xi }_{\mathrm{x}}^a-{\xi }_{\mathrm{y}}^a\right).$$

(43)

The sum is over $\mathrm{y}\ne \mathrm{x}$ because ${\Delta }_{\mathrm{x}}{X}_{\mathrm{x}}^A\equiv 0$.

Using the facts that

$${\Delta }_{\mathrm{x}}{g}^{ab}=\delta g^{ab}-2{\nabla }^{(a}{\xi }_{\mathrm{x}}^{b)},$$

(44)

$$\delta {ϵ}^{abcd}=-\frac{1}{2}{ϵ}^{abcd}{g}^{ef}\delta g_{ef},$$

(45)

and

$${ϵ}^{bcda}{\mathcal{L}}_{{\xi }_{\mathrm{x}}}{n}_{bcd}^{\mathrm{x}}=3!\left({\xi }_{\mathrm{x}}^b{\nabla }_b{n}_{\mathrm{x}}^a-n_{\mathrm{x}}^b{\nabla }_b{\xi }_{\mathrm{x}}^a+n_{\mathrm{x}}^a{\nabla }_b{\xi }_{\mathrm{x}}^b\right),$$

(46)

we find

$$\begin{gathered} \begin{array}{ccc}\hfill \delta n_{\mathrm{x}}^a& =& \delta \left(\frac{1}{3!}{ϵ}^{bcda}{n}_{bcd}^{\mathrm{x}}\right)\hfill \end{array} \\ \begin{array}{ccc}& =& n_{\mathrm{x}}^b{\nabla }_b{\xi }_{\mathrm{x}}^a-{\xi }_{\mathrm{x}}^b{\nabla }_b{n}_{\mathrm{x}}^a-n_{\mathrm{x}}^a\left({\nabla }_b{\xi }_{\mathrm{x}}^b+\frac{1}{2}{g}^{bc}\delta g_{bc}\right)\hfill \end{array} \\ \begin{array}{ccc}& & +\frac{1}{n_{\mathrm{x}}}{n}_{\mathrm{x}}^a\sum _{\mathrm{y}\ne \mathrm{x}}\left(\frac{1}{{\overline{\mu }}_{\mathrm{x}}}{\mathcal{R}}_b^{\mathrm{x}\mathrm{y}}\right)\left({\xi }_{\mathrm{x}}^b-{\xi }_{\mathrm{y}}^b\right),\hfill \end{array} \end{gathered}$$

(47)

where

$$\frac{1}{{\overline{\mu }}_{\mathrm{x}}}{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}\equiv \frac{1}{3!}{ϵ}_{\mathrm{x}}^{ABC}\frac{\partial n_{ABC}^{\mathrm{x}}}{\partial X_{\mathrm{y}}^D}\frac{\partial X_{\mathrm{y}}^D}{\partial x^a}.$$

(48)

The coefficient ${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}$ satisfies the identity

$$u_{\mathrm{y}}^a{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}\equiv 0⟹{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}=\left({\delta }_a^b+u_{\mathrm{y}}^b{u}_a^{\mathrm{y}}/c^2\right){\mathcal{R}}_b^{\mathrm{x}\mathrm{y}}.$$

(49)

This says that ${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}$ has only three degrees of freedom; i.e., $u_{\mathrm{y}}^a$ is timelike and therefore ${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}$ has only the spacelike components with respect to the $u_{\mathrm{y}}^a$.

We will see in the next subSection 4.3, where the equations of motion are derived, that there is a total “resistivity” current $R_a^{\mathrm{x}}$ which is given by

$$R_a^{\mathrm{x}}=\sum _{\mathrm{y}\ne \mathrm{x}}\left({\mathcal{R}}_a^{\mathrm{y}\mathrm{x}}-{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}\right),$$

(50)

and satisfies the identity

$$\sum _{\mathrm{x}}{R}_a^{\mathrm{x}}\equiv 0.$$

(51)

This identity is important because it guarantees that the energy-momentum-stress tensor $T_{ab}$ is divergenceless, i.e. ${\nabla }_b{T}^{ba}={\nabla }_b{T}^{ab}=0$ (a consequence of diffeomorphism invariance [14]).

4.3. The Field Equations

We now have everything we need to derive the full suite of field equations. Let us begin by returning to the flux variations of the total action given in Eq. (33). The fact that we are summing over all constituents leads to

$$\sum _{\mathrm{x}}\sum _{\mathrm{y}\ne \mathrm{x}}{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}\left({\xi }_{\mathrm{x}}^a-{\xi }_{\mathrm{y}}^a\right)=-\sum _{\mathrm{x}}{R}_a^{\mathrm{x}}{\xi }_{\mathrm{x}}^a,$$

(52)

so that the variation of the total action for the system is

$$\begin{gathered} \begin{array}{ccc}\hfill \delta S& =& {\int }_{\mathcal{M}}{\mathrm{d}}^{4}x\sqrt{-g}\left[-\sum _{\mathrm{x}}\left(f_a^{\mathrm{x}}+{\Gamma }_{\mathrm{x}}{\overline{\mu }}_a^{\mathrm{x}}-R_a^{\mathrm{x}}\right){\xi }_{\mathrm{x}}^a-\frac{1}{{\mu }_{\mathrm{o}}}\left({\nabla }_b{F}^{ab}-{\mu }_{\mathrm{o}}\sum _{\mathrm{x}}{j}_{\mathrm{x}}^a\right)\delta A_a\right.\hfill \end{array} \\ \begin{array}{ccc}& & \left.-\frac{1}{2}\left(\frac{c^4}{8\pi G}{G}^{ab}-T^{ab}\right)\delta g_{ab}\right].\hfill \end{array} \end{gathered}$$

(53)

where

$$f_a^{\mathrm{x}}=2n_{\mathrm{x}}^b{\nabla }_{[b}{\overline{\mu }}_{a]}^{\mathrm{x}}=2n_{\mathrm{x}}^b{\nabla }_{[b}{\mu }_{a]}^{\mathrm{x}}+q_{\mathrm{x}}{n}_{\mathrm{x}}^b{F}_{ba},$$

(54)

$$\Psi =\Lambda -\sum _{\mathrm{x}}{\mu }_c^{\mathrm{x}}{n}_{\mathrm{x}}^c,$$

(55)

and

$$T^{ab}=\Psi g^{ab}+\sum _{\mathrm{x}}{n}_{\mathrm{x}}^a{\mu }_{\mathrm{x}}^b-\frac{1}{4{\mu }_{\mathrm{o}}}\left(F_{cd}{F}^{cd}{g}^{ab}-4F^{ac}{F}^b{}_c\right).$$

(56)

It is worth noting here that the generalized pressure $\Psi$ takes the form of a Legendre transformation of $\Lambda$, which switches the roles of $n_{\mathrm{x}}^a$ and ${\mu }_a^{\mathrm{x}}$, making the latter the independent degree of freedom; i.e.

$$\delta \Psi =-\sum _{\mathrm{x}}{n}_{\mathrm{x}}^c\delta {\mu }_c^{\mathrm{x}}.$$

(57)

This will be especially useful later when we write down the Newtonian fluid/plasma field equations.

Now that the action variation is in place, we can invoke our chosen constraint that a given particle flux and its corresponding entropy flux flow together. We also restrict (by assumption!) the neutral and positively charged species to flow together. The net result is that there are only two matter spaces where $X_1^A=X_2^A=X_4^A=X_5^A\equiv X_{\mathcal{P}}^A$ and $X_3^A=X_6^A\equiv X_{\mathcal{N}}^A$. This also implies there are only two independent Lagrangian displacements: ${\xi }_1^a={\xi }_2^a={\xi }_4^a={\xi }_5^a\equiv {\xi }_{\mathcal{P}}^a$ and ${\xi }_3^a={\xi }_6^a\equiv {\xi }_{\mathcal{N}}^a$. Likewise, there are only two independent four-velocities: $u_1^a=u_2^a=u_4^a=u_5^a\equiv u_{\mathcal{P}}^a$ and $u_3^a=u_6^a\equiv u_{\mathcal{N}}^a$. We also note that $q_1=q_4=q_5=q_6=0$ and $q_2=-q_3=-q_{\mathcal{N}}$.

In order to get the field equations we employ the action principle, which states that when $\delta S=0$ for arbitrary values for the variations ${\xi }_{\mathrm{x}}^a$, $\delta A_a$, and $\delta g_{ab}$, then the coefficients multiplying them in $\delta S$ must vanish. From the coefficient of ${\xi }_{\mathcal{P}}^a$, we get a single Euler equation for the neutrally and positively charged species, which is

$$\sum _{\mathrm{x}=\{1,2,4,5\}}\left[f_a^{\mathrm{x}}+{\Gamma }_{\mathrm{x}}{\mu }_a^{\mathrm{x}}-\left(R_a^{\mathrm{x}}-q_{\mathrm{x}}{\Gamma }_{\mathrm{x}}{A}_a\right)\right]=0,$$

(58)

and from ${\xi }_{\mathcal{N}}^a$ a single Euler equation for the negative species, which is

$$\sum _{\mathrm{x}=\{3,6\}}\left[f_a^{\mathrm{x}}+{\Gamma }_{\mathrm{x}}{\mu }_a^{\mathrm{x}}-\left(R_a^{\mathrm{x}}-q_{\mathrm{x}}{\Gamma }_{\mathrm{x}}{A}_a\right)\right]=0.$$

(59)

Coming from the coefficient of $\delta A_a$ are the Maxwell equations [which must also include Eq. (21)],

$${\nabla }_b{F}^{ab}={\nabla }_b\left({\nabla }^a{A}^b-{\nabla }^b{A}^a\right)={\mu }_{\mathrm{o}}\sum _{\mathrm{x}=\{2,3\}}{j}_{\mathrm{x}}^a,$$

(60)

and from $\delta g_{ab}$ we get the Einstein equation; i.e.

$$G^{ab}=\frac{8\pi G}{c^4}{T}^{ab}.$$

(61)

An equivalent form of the Einstein equation, which will be used in Sec. Section 5, is

$$R^{ab}=\frac{8\pi G}{c^4}\left(T^{ab}-\frac{1}{2}T{g}^{ab}\right),$$

(62)

where $T=g^{ab}{T}_{ab}$.

From the process of creating the two Euler equations (58) and (59), we find that the set of resistivity vectors $R_a^{\mathrm{x}}$ is reduced from six members down to two, which we denote by $R_a^{\mathcal{P}}$ and $R_a^{\mathcal{N}}$. If we take into account that $X_1^A=X_2^A=X_4^A=X_5^A$ and $X_3^A=X_6^A\equiv X_{\mathcal{N}}^A$, then we see that Eq. (48) implies

$$\begin{array}{cc}\hfill {\mathcal{R}}_a^{12}& ={\mathcal{R}}_a^{14}={\mathcal{R}}_a^{15},{\mathcal{R}}_a^{13}={\mathcal{R}}_a^{16},\hfill \end{array}$$

(63a)

$$\begin{array}{cc}\hfill {\mathcal{R}}_a^{21}& ={\mathcal{R}}_a^{24}={\mathcal{R}}_a^{25},{\mathcal{R}}_a^{23}={\mathcal{R}}_a^{26},\hfill \end{array}$$

(63b)

$$\begin{array}{cc}\hfill {\mathcal{R}}_a^{31}& ={\mathcal{R}}_a^{32}={\mathcal{R}}_a^{34}={\mathcal{R}}_a^{35},\hfill \end{array}$$

(63c)

$$\begin{array}{cc}\hfill {\mathcal{R}}_a^{41}& ={\mathcal{R}}_a^{42}={\mathcal{R}}_a^{45},{\mathcal{R}}_a^{43}={\mathcal{R}}_a^{46},\hfill \end{array}$$

(63d)

$$\begin{array}{cc}\hfill {\mathcal{R}}_a^{51}& ={\mathcal{R}}_a^{52}={\mathcal{R}}_a^{54},{\mathcal{R}}_a^{53}={\mathcal{R}}_a^{56},\hfill \end{array}$$

(63e)

$$\begin{array}{cc}\hfill {\mathcal{R}}_a^{61}& ={\mathcal{R}}_a^{62}={\mathcal{R}}_a^{64}={\mathcal{R}}_a^{65}.\hfill \end{array}$$

(63f)

Inserting these into the definition of $R_a^{\mathrm{x}}$ in Eq. (50) leads to

$$R_a^{\mathcal{P}}=R_a^1+R_a^2+R_a^4+R_a^5=4\left({\mathcal{R}}_a^{31}+{\mathcal{R}}_a^{61}\right)-2\left({\mathcal{R}}_a^{13}+{\mathcal{R}}_a^{23}+{\mathcal{R}}_a^{43}+{\mathcal{R}}_a^{53}\right).$$

(64)

In a similar manner, we obtain

$$R_a^{\mathcal{N}}=R_a^3+R_a^6=-4\left({\mathcal{R}}_a^{31}+{\mathcal{R}}_a^{61}\right)+2\left({\mathcal{R}}_a^{13}+{\mathcal{R}}_a^{23}+{\mathcal{R}}_a^{43}+{\mathcal{R}}_a^{53}\right)=-R_a^{\mathcal{P}},$$

(65)

so the identity in Eq. (51) becomes

$$R_a^{\mathcal{P}}+R_a^{\mathcal{N}}=0.$$

(66)

Ultimately, microphysical calculations will be required to precisely specify $R_a^{\mathcal{N}}$ (e.g. as indicated by Braginskii [31]). However, the formalism itself has already provided some structure for the resistivities $R_a^{\mathrm{x}}$, as evidenced by Eqs. (19), (35), (36), (48), (49), and (66). Recall that the main assumption is that $n_{ABC}^{\mathrm{x}}$ depends on, in principle, all of the $X_{\mathrm{x}}^A$. Because of Eq. (16), then the chain-rule implies

$$\frac{\mathrm{d}}{\mathrm{d}{\tau }_{\mathrm{x}}}{n}_{ABC}^{\mathrm{x}}=u_{\mathrm{x}}^a\sum _{\mathrm{y}\ne \mathrm{x}}\frac{\partial X_{\mathrm{y}}^D}{\partial x^a}\frac{\partial n_{ABC}^{\mathrm{x}}}{\partial X_{\mathrm{y}}^D}.$$

(67)

When we substitute this into Eq. (19), and use Eq. (49), we obtain

$${\overline{\mu }}_{\mathrm{x}}{\Gamma }_{\mathrm{x}}=-u_{\mathrm{x}}^a\sum _{\mathrm{y}\ne \mathrm{x}}\left({\mathcal{R}}_a^{\mathrm{y}\mathrm{x}}-{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}\right)=-u_{\mathrm{x}}^a{R}_a^{\mathrm{x}}.$$

(68)

4.4. Impact of Change of Gauge for $A_a$

A gauge transformation will impact the fluid equations of motion because of the change to the momentum; i.e. letting ${\overline{A}}_a=A_a+{\nabla }_a\varphi$ we find

$${\overline{\mu }}_a^{\mathrm{x}}={\mu }_a^{\mathrm{x}}+q_{\mathrm{x}}{A}_a⟶{\widehat{\mu }}_a^{\mathrm{x}}={\mu }_a^{\mathrm{x}}+q_{\mathrm{x}}{\overline{A}}_a={\overline{\mu }}_a^{\mathrm{x}}+q_{\mathrm{x}}{\nabla }_a\varphi .$$

(69)

It is important here to consider in more detail the ramifications of a change of gauge, since a natural application of the present work would be to numerical evolutions [32]. In the numerical setting, we expect to be solving for the vector potential $A_a$ as we evolve the system. This will require a choice of gauge for the vector potential, which will affect the explicit values of terms (such as the resistivity) in the equations of motion.

Clearly, $R_a^{\mathrm{x}}$ is gauge-dependent, since the quantity ${\overline{\mu }}_{\mathrm{x}}^{ABC}$ in ${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}$ [cf.Eq. (48)] depends on $A_a$. Letting ${\overline{R}}_a^{\mathrm{x}}$ denote the particle resistivity in the new gauge, we find

$$\begin{gathered} \begin{array}{ccc}\hfill {\overline{R}}_a^{\mathrm{x}}& =& \sum _{\mathrm{y}\ne \mathrm{x}}\left({\overline{\mathcal{R}}}_a^{\mathrm{y}\mathrm{x}}-{\overline{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}\right)\hfill \end{array} \\ \begin{array}{ccc}& =& \sum _{\mathrm{y}\ne \mathrm{x}}\frac{1}{3!}{ϵ}^{ebcd}\left[\left({\overline{\mu }}_e^{\mathrm{y}}+q_{\mathrm{y}}{\nabla }_e\varphi \right){}^{\mathrm{y}}{\mathcal{J}}_{bcd}^{ABC}\frac{\partial n_{ABC}^{\mathrm{y}}}{\partial X_{\mathrm{x}}^D}\frac{\partial X_{\mathrm{x}}^D}{\partial x^a}\right.\hfill \end{array} \\ \begin{array}{ccc}& & \left.-\left({\overline{\mu }}_e^{\mathrm{x}}+q_{\mathrm{x}}{\nabla }_e\varphi \right){}^{\mathrm{x}}{\mathcal{J}}_{bcd}^{ABC}\frac{\partial n_{ABC}^{\mathrm{x}}}{\partial X_{\mathrm{y}}^D}\frac{\partial X_{\mathrm{y}}^D}{\partial x^a}\right]\hfill \end{array} \\ \begin{array}{ccc}& =& R_a^{\mathrm{x}}+G_a^{\mathrm{x}},\hfill \end{array} \end{gathered}$$

(70)

where

$$G_a^{\mathrm{x}}=\sum _{\mathrm{y}\ne \mathrm{x}}\left({\mathcal{G}}_a^{\mathrm{y}\mathrm{x}}-{\mathcal{G}}_a^{\mathrm{x}\mathrm{y}}\right),{\mathcal{G}}_a^{\mathrm{x}\mathrm{y}}=\frac{1}{3!}{ϵ}^{ebcd}{q}_{\mathrm{x}}\left({}^{\mathrm{x}}{\mathcal{J}}_{bcd}^{ABC}\frac{\partial n_{ABC}^{\mathrm{x}}}{\partial X_{\mathrm{y}}^D}\frac{\partial X_{\mathrm{y}}^D}{\partial x^a}\right){\nabla }_e\varphi .$$

(71)

Note that

$$\sum _{\mathrm{x}}{R}_a^{\mathrm{x}}=\sum _{\mathrm{x}}{G}_a^{\mathrm{x}}=0⟹\sum _{\mathrm{x}}{\overline{R}}_a^{\mathrm{x}}=\sum _{\mathrm{x}}{R}_a^{\mathrm{x}}+\sum _{\mathrm{x}}{G}_a^{\mathrm{x}}=0.$$

(72)

Using Eqs. (10), (15), and (48), we can re-write ${\mathcal{G}}_a^{\mathrm{x}\mathrm{y}}$ as

$$\begin{gathered} \begin{array}{ccc}\hfill {\mathcal{G}}_a^{\mathrm{x}\mathrm{y}}& =& \frac{1}{3!}{q}_{\mathrm{x}}\left({ϵ}^{ebcd}{}^{\mathrm{x}}{\mathcal{J}}_{bcd}^{ABC}{\delta }_A^{[E}{\delta }_B^F{\delta }_C^{G]}\frac{\partial n_{EFG}^{\mathrm{x}}}{\partial X_{\mathrm{y}}^D}\frac{\partial X_{\mathrm{y}}^D}{\partial x^a}\right){\nabla }_e\varphi \hfill \end{array} \\ \begin{array}{ccc}& =& -\frac{q_{\mathrm{x}}}{{\overline{\mu }}_{\mathrm{x}}}\left(\frac{1}{3!}{\overline{\mu }}_{\mathrm{x}}{ϵ}_{\mathrm{x}}^{EFG}\frac{\partial n_{EFG}^{\mathrm{x}}}{\partial X_{\mathrm{y}}^D}\frac{\partial X_{\mathrm{y}}^D}{\partial x^a}\right)\left(\frac{1}{3!}{ϵ}^{bcde}{}^{\mathrm{x}}{\mathcal{J}}_{bcd}^{ABC}{ϵ}_{ABC}^{\mathrm{x}}\right){\nabla }_e\varphi \hfill \end{array} \\ \begin{array}{ccc}& =& -q_{\mathrm{x}}\left(u_{\mathrm{x}}^b{\nabla }_b\varphi \right)\left(\frac{1}{{\overline{\mu }}_{\mathrm{x}}}{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}\right),\hfill \end{array} \end{gathered}$$

(73)

which implies

$$G_a^{\mathrm{x}}=-\sum _{\mathrm{y}\ne \mathrm{x}}\left[q_{\mathrm{y}}\left(u_{\mathrm{y}}^b{\nabla }_b\varphi \right)\left(\frac{1}{{\overline{\mu }}_{\mathrm{y}}}{\mathcal{R}}_a^{\mathrm{y}\mathrm{x}}\right)-q_{\mathrm{x}}\left(u_{\mathrm{x}}^b{\nabla }_b\varphi \right)\left(\frac{1}{{\overline{\mu }}_{\mathrm{x}}}{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}\right)\right].$$

(74)

When the sums in Eqs. (58) and (59) are performed, we see that the gauge-dependent part of each of the fluid equations of motion is

$$\begin{gathered} \begin{array}{ccc}\hfill {\overline{R}}_a^{\mathcal{N}}-q_{\mathcal{N}}{\Gamma }_3{\overline{A}}_a& =& R_a^{\mathcal{N}}-q_{\mathcal{N}}{\Gamma }_3{A}_a-4q_{\mathcal{N}}\frac{1}{{\overline{\mu }}_3}\left(u_{\mathcal{N}}^b{\mathcal{R}}_b^{31}\right){\nabla }_a\varphi \hfill \end{array} \\ \begin{array}{ccc}& & +2q_{\mathcal{N}}\left[\frac{1}{{\overline{\mu }}_2}{u}_{\mathcal{P}}^b{\mathcal{R}}_a^{23}-\frac{1}{{\overline{\mu }}_3}{u}_{\mathcal{N}}^b\left(2{\mathcal{R}}_a^{31}+{\mathcal{R}}_a^{36}\right)\right]{\nabla }_b\varphi .\hfill \end{array} \end{gathered}$$

(75)

Clearly, Eqs. (58) and (59) are modified under a gauge transformation. This was expected. The point is that we have shown how the transformation enters the field equations and therefore we can still evolve the system regardless of the choice of gauge.

It is a different story if we look at the projection of Eq. (58) along $u_{\mathcal{P}}^a$ and Eq. (59) along $u_{\mathcal{N}}^a$. Clearly, $u_{\mathcal{P}}^a{f}_a^{\mathrm{x}}=0$ for Eq. (58) and $u_{\mathcal{N}}^a{f}_a^{\mathrm{x}}=0$ for Eq. (59), leaving two equations having linear combinations of creation rates ${\Gamma }_{\mathrm{x}}$, combined with the resistivity and the gauge-dependent terms. The creation rates must be gauge invariant. Fortunately, if we use Eq. (49), and project Eq. (75) along $u_{\mathcal{P}}^a$ and then along $u_{\mathcal{N}}^a$, we get

$$\begin{array}{cc}\hfill u_{\mathcal{P}}^a\left({\overline{R}}_a^{\mathcal{N}}-q_{\mathcal{N}}{\Gamma }_3{\overline{A}}_a\right)& =u_{\mathcal{P}}^a\left(R_a^{\mathcal{N}}-q_{\mathcal{N}}{\Gamma }_3{A}_a\right),\hfill \end{array}$$

(76a)

$$\begin{array}{cc}\hfill u_{\mathcal{N}}^a\left({\overline{R}}_a^{\mathcal{N}}-q_{\mathcal{N}}{\Gamma }_3{\overline{A}}_a\right)& =u_{\mathcal{N}}^a\left(R_a^{\mathcal{N}}-q_{\mathcal{N}}{\Gamma }_3{A}_a\right),\hfill \end{array}$$

(76b)

thus verifying that the ${\Gamma }_{\mathrm{x}}$ are gauge invariant. This was also noted in [13] and is a result of starting with an action with well-defined couplings. The formalism itself takes care of gauge issues through internal consistency.

5. $\mathbf{3}+\mathbf{1}$ Formulation

Having derived the equations of motion for the plasma system, we want to make contact with applications and known results in the non-relativistic limit. In order to do this, we work out the $3+1$ form of the field equations, keeping the speed of light explicit. This makes taking the Newtonian limit a simple power counting exercise and also sets the problem up for foliation-based numerical simulations. Our approach to the $3+1$ problem follows the set of notes by Gourgoulhon [33].

5.1. The $3+1$ Setup

We begin by restricting our formalism to a special class of manifolds—globally hyperbolic. These manifolds contain a family of causal curves, which are such that every vector tangent to them is timelike or null. They also contain a Cauchy surface, which is a spacelike hypersurface that is intersected exactly once by every inextendible causal curve in the manifold. It can be shown that, on these manifolds with coordinates ${\overline{x}}^a$, a scalar “time” function $t\left({\overline{x}}^a\right)$ exists such that its level (“constant time”) hypersurfaces can be smoothly stacked on top of each other to form a foliation of the spacetime.

A normal at a point on a constant-time hypersurface is obtained in the standard way by taking the gradient of the time function, i.e. ${\nabla }_{a}t$, and then evaluating the gradient at the point under consideration. A unit normal $u^a$ ($u^a{u}_a=-c^2$) at each point is created by introducing the so-called lapse function N, which is a speed, as a normalization factor for ${\nabla }_{a}t$; that is,

$$u^a=-cN{\nabla }^{a}t.$$

(77)

If we build an initial slice of the foliation by solving $t\left({\overline{x}}^a\right)=t_o=$ constant, the next one, say for $t=t_o+\delta t$, will consist of the set of points obtained by moving the same, “small” proper distance in the $u^a$ direction. The $u^a$ will merge together from slice-to-slice to become tangents to worldlines. The acceleration $a_a$ of an observer following one of these worldlines is

$$a_a=u^b{\nabla }_b{u}_a\equiv \frac{\mathrm{D}{u}_a}{\mathrm{d}t},$$

(78)

which introduces our notion of time-derivative.

So far, we have a mechanism for stacking the spacelike hypersurfaces, but nothing for how they “slip” past each other. To take care of that we introduce a “flow-of-time” vector $t^a$ (with the units of speed) which joins spatial points ${\left.{\overline{x}}^i\right|}_{t_o}$ on the hypersurface $t=t_o$ to spatial points ${\left.{\overline{x}}^i\right|}_{t_o+\delta t}$ on the next hypersurface $t=t_o+\delta t$ such that ${\left.{\overline{x}}^i\right|}_{t_o}={\left.{\overline{x}}^i\right|}_{t_o+\delta t}$; in words, it is the observers following $t^a$ and not $u^a$ who are “at rest” with respect to the foliation slices. We normalize $t^a$ by setting

$$t^a{\nabla }_{a}t=1.$$

(79)

We can use $u^a/c$ in two ways to decompose $t^a$ into pieces perpendicular and parallel to the foliation slices; namely,

$$t^a=N\left(u^a/c\right)+N^a,N^a={\perp }_b^a{t}^b,{\perp }_b^a={\delta }_b^a+u^a{u}_b/c^2,$$

(80)

where $N^a$ is the so-called shift vector (with speed units). The tensor ${\perp }_a^b$ is the (idempotent) operator that provides the parallel (spacelike) projection and $u^a/c$ provides the perpendicular (timelike) projection. Since ${\perp }_b^a{u}^b=0$ the shift vector satisfies $\left(u_a/c\right)N^a=0$ and therefore has no perpendicular component.

Each slice of the foliation is, in principle, a curved space. The curvature information is contained in an induced three-metric $h_{ab}$ given by

$$h_{ab}={\perp }_a^c{\perp }_b^d{g}_{cd}=g_{ab}+u_a{u}_b/c^2.$$

(81)

Our notion of spatial covariant derivative $D_a$ is generated by the action of ${\perp }_b^a$ on the covariant derivative of an arbitrary vector ${\tilde{v}}^a={\perp }_b^a{v}^b$; namely,

$$D_a{\tilde{v}}^b={\perp }_a^c{\perp }_d^b{\nabla }_c{\tilde{v}}^d.$$

(82)

The three-metric $h_{ab}$ is compatible with $D_a$; i.e. $D_a{h}_{bc}=0$. The intrinsic curvature of slices of the foliation, ${}^{\left(3\right)}{R}^c{}_{dab}$, can be inferred from

$$D_a{D}_b{\tilde{v}}^c-D_b{D}_a{\tilde{v}}^c={}^{\left(3\right)}{R}^c{}_{dab}{\tilde{v}}^d.$$

(83)

The acceleration can be shown [by inserting Eq. (77) into (78)] to have the alternative form

$$a_a=c^2{D}_{a}ln\left(N/c\right).$$

(84)

Because the three-dimensional slices of the foliation are embedded in four-dimensional spacetime, they have an extrinsic curvature $K_{ab}$ (with inverse time dimensions) given by

$$K_{ab}=-\frac{1}{2}{\mathcal{L}}_u{h}_{ab}=-\frac{1}{2}\left({\perp }_b^c{\nabla }_c{u}_a+{\perp }_a^c{\nabla }_c{u}_b\right).$$

(85)

It is easy to show that the trace of the extrinsic curvature, which is $K=g^{ab}{K}_{ab}$, becomes

$$K=-{\nabla }_a{u}^a\equiv -\Theta .$$

(86)

When we develop the $3+1$ form of the field equations it will be found that the covariant derivative of $u^a$ enters repeatedly. A couple of important “tools” for dealing with this can be obtained by applying the well-known decomposition

$${\nabla }_a{u}_b={\sigma }_{ab}+\frac{1}{3}\Theta h_{ab}+{\varpi }_{ab}-a_b{u}_a/c^2=-K_{ab}+{\varpi }_{ab}-u_a{a}_b/c^2,$$

(87)

where

$$\begin{array}{cc}\hfill {\sigma }_{ab}& =\frac{1}{2}\left({\perp }_b^c{\nabla }_c{u}_a+{\perp }_a^c{\nabla }_c{u}_b\right)-\frac{1}{3}\Theta h_{ab}=-\left(K_{ab}-\frac{1}{3}K{h}_{ab}\right),\hfill \end{array}$$

(88a)

$$\begin{array}{cc}\hfill {\varpi }_{ab}& =\frac{1}{2}\left({\perp }_b^c{\nabla }_c{u}_a-{\perp }_a^c{\nabla }_c{u}_b\right).\hfill \end{array}$$

(88b)

The most useful formula is a consequence of the fact that $u^a$ is surface forming: This implies ${\varpi }_{ab}=0$, and so therefore

$${\nabla }_a{u}_b=-K_{ab}-u_a{a}_b/c^2.$$

(89)

From this we can immediately show

$${\nabla }_c{\perp }_a^b=-2g^{bd}\left[u_c{u}_{(a}{a}_{d)}/c^4+u_{(a}{K}_{d)c}/c^2\right].$$

(90)

5.2. Field Decompositions

We have just seen how the metric can be re-framed in terms of the lapse N, the shift-vector $N^a$, and the three-metric $h_{ab}$. Now we need to produce the similar re-framing for the remaining field variables $n_{\mathrm{x}}^a$ and $A_a$.

Using the projection operators $u^a/c$ and ${\perp }_b^a$, and taking into account the dimensional analysis of the flux earlier, the $3+1$ forms of the fluxes must be

$$n_{\mathrm{x}}^a={\tilde{n}}_{\mathrm{x}}{u}^a+{\tilde{n}}_{\mathrm{x}}^a,{\tilde{n}}_{\mathrm{x}}=-\left(u_a/c^2\right)n_{\mathrm{x}}^a,{\tilde{n}}_{\mathrm{x}}^a={\perp }_b^a{n}_{\mathrm{x}}^b.$$

(91)

From the definition of the four-velocity $u_{\mathrm{x}}^a=n_{\mathrm{x}}^a/n_{\mathrm{x}}$ we can infer

$$u_{\mathrm{x}}^a=\frac{{\tilde{n}}_{\mathrm{x}}}{n_{\mathrm{x}}}\left(u^a+{\tilde{u}}_{\mathrm{x}}^a\right),{\tilde{u}}_{\mathrm{x}}^a=\frac{{\tilde{n}}_{\mathrm{x}}^a}{{\tilde{n}}_{\mathrm{x}}},$$

(92)

and can therefore show

$$\frac{{\tilde{n}}_{\mathrm{x}}}{n_{\mathrm{x}}}={\tilde{\gamma }}_{\mathrm{x}},{\tilde{\gamma }}_{\mathrm{x}}=\frac{1}{\sqrt{1-{\tilde{u}}_a^{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^a/c^2}}.$$

(93)

Because $u_{\mathrm{x}}^a{u}_a=-{\tilde{\gamma }}_{\mathrm{x}}$ and $u_{\eta }^a{u}_a=u_{\mathcal{P}}^a{u}_a$ we have ${\tilde{\gamma }}_1={\tilde{\gamma }}_2={\tilde{\gamma }}_4={\tilde{\gamma }}_5\equiv {\tilde{\gamma }}_{\mathcal{P}}$ and consequently ${\tilde{u}}_1^a={\tilde{u}}_2^a={\tilde{u}}_4^a={\tilde{u}}_5^a\equiv {\tilde{u}}_{\mathcal{P}}^a$. Similarly, we have ${\tilde{\gamma }}_3={\tilde{\gamma }}_6\equiv {\tilde{\gamma }}_{\mathcal{N}}$ and ${\tilde{u}}_3^a={\tilde{u}}_6^a\equiv {\tilde{u}}_{\mathcal{N}}^a$.

For the chemical potential covector ${\mu }_a^{\mathrm{x}}$, the dimensional analysis leads to a slightly different form for the decompositions:

$${\mu }_a^{\mathrm{x}}={\tilde{\mu }}_{\mathrm{x}}{u}_a/c^2+{\tilde{\mu }}_a^{\mathrm{x}},{\tilde{\mu }}_{\mathrm{x}}=-u^a{\mu }_a^{\mathrm{x}},{\tilde{\mu }}_a^{\mathrm{x}}={\perp }_a^b{\mu }_b^{\mathrm{x}}.$$

(94)

If we substitute into the spatial part of this the initial result for ${\mu }_a^{\mathrm{x}}$, i.e. Eq. (24), we find a form more amenable for the Newtonian limit, which is

$${\tilde{\mu }}_a^{\mathrm{x}}=\frac{{\tilde{\mu }}_{\mathrm{x}}}{c^2}{\tilde{u}}_a^{\mathrm{x}}+\sum _{\mathrm{y}\ne \mathrm{x}}{\mathcal{A}}^{\mathrm{x}\mathrm{y}}{\tilde{n}}_{\mathrm{y}}{\tilde{w}}_a^{\mathrm{y}\mathrm{x}},$$

(95)

where

$${\tilde{w}}_{\mathrm{y}\mathrm{x}}^a={\tilde{u}}_{\mathrm{y}}^a-{\tilde{u}}_{\mathrm{x}}^a.$$

(96)

As an effect of the tilting of the momenta, the chemical potentials in the fluid rest frames are related with those of the foliation in more complicated ways, which are

$${\mu }_{\mathrm{x}}={\tilde{\gamma }}_{\mathrm{x}}\left({\tilde{\mu }}_{\mathrm{x}}-{\tilde{u}}_{\mathrm{x}}^a{\tilde{\mu }}_a^{\mathrm{x}}\right).$$

(97)

By direct substitution of the decompositions just above into Eq. (61), the generalized pressure $\Psi$ becomes

$$\Psi =\Lambda +\sum _{\mathrm{x}}\left({\tilde{\mu }}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}-{\tilde{\mu }}_a^{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}^a\right),$$

(98)

and the fluid/plasma part of $T^{ab}$ is

$$\begin{gathered} \begin{array}{ccc}\hfill \Psi g^{ab}+\sum _{\mathrm{x}}{n}_{\mathrm{x}}^a{\mu }_{\mathrm{x}}^b& =& \left(-\Lambda +\sum _{\mathrm{x}}{\tilde{\mu }}_c^{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}^c\right)u^a{u}^b/c^2+\sum _{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}{\tilde{\mu }}_{\mathrm{x}}^b{u}^a+\sum _{\mathrm{x}}\frac{{\tilde{\mu }}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}}{c^2}{\tilde{u}}_{\mathrm{x}}^a{u}^b\hfill \end{array} \\ \begin{array}{ccc}& & +\Psi h^{ab}+\sum _{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}^a{\tilde{\mu }}_{\mathrm{x}}^b.\hfill \end{array} \end{gathered}$$

(99)

The charge current flux $j_{\mathrm{x}}^a$ is

$$j_{\mathrm{x}}^a={\tilde{\sigma }}_{\mathrm{x}}{u}^a+{\tilde{j}}_{\mathrm{x}}^a,{\tilde{\sigma }}_{\mathrm{x}}=-\left(u_a/c^2\right)j_{\mathrm{x}}^a,{\tilde{j}}_{\mathrm{x}}^a={\perp }_b^a{j}_{\mathrm{x}}^b,$$

(100)

and the four-potential $A_a$ is

$$A_a=\tilde{V}{u}_a/c^2+{\tilde{A}}_a,\tilde{V}=-u^a{A}_a,{\tilde{A}}_a={\perp }_a^b{A}_b,$$

(101)

where we have introduced the scalar potential $\tilde{V}$ (with the standard energy per charge units) and the three-vector potential ${\tilde{A}}_a$. Inserting this into the Faraday tensor, and using Eq. (89) for the covariant derivative, we find

$$\begin{gathered} \begin{array}{ccc}\hfill F_{ab}& =& \frac{1}{c^2}\left(u_b{D}_a-u_a{D}_b\right)\tilde{V}+\frac{1}{c^2}\tilde{V}\left(a_a{u}_b-a_b{u}_a\right)+\frac{1}{c^2}{u}_b\frac{\mathrm{D}{\tilde{A}}_a}{dt}-\frac{1}{c^2}{u}_a\frac{\mathrm{D}{\tilde{A}}_b}{dt}\hfill \end{array} \\ \begin{array}{ccc}& & -\frac{1}{c^2}\left(u_b{\delta }_a^d-u_a{\delta }_b^d\right){\tilde{A}}^c{K}_{dc}+D_a{\tilde{A}}_b-D_b{\tilde{A}}_a.\hfill \end{array} \end{gathered}$$

(102)

The electric ${\tilde{E}}_a$ and magnetic ${\tilde{B}}_a$ fields are defined as

$$\begin{array}{cc}\hfill {\tilde{E}}_a& =-\frac{u^b}{c}{F}_{ba}=-D_a\tilde{V}-{\perp }_a^b\frac{\mathrm{D}{\tilde{A}}_b}{\mathrm{d}t}-\tilde{V}{a}_a/c^2+{\tilde{A}}^b{K}_{ab},\hfill \end{array}$$

(103a)

$$\begin{array}{cc}\hfill {\tilde{B}}_a& =\frac{1}{2}{\widetilde{ϵ}}_a{}^{bc}{F}_{bc}={\widetilde{ϵ}}_a{}^{bc}{D}_b{\tilde{A}}_c,{\widetilde{ϵ}}_{abc}=\frac{u^d}{c}{ϵ}_{dabc},\hfill \end{array}$$

(103b)

which implies

$$F_{ab}=\frac{2}{c^2}{u}_{[a}{\tilde{E}}_{b]}+{\widetilde{ϵ}}_{abc}{\tilde{B}}^c.$$

(104)

Finally, the electromagnetic contribution to $T^{ab}$ is

$$\begin{gathered} \begin{array}{ccc}\hfill -\frac{1}{4{\mu }_{\mathrm{o}}}\left(F_{cd}{F}^{cd}{g}^{ab}-4F^{ac}{F}^b{}_c\right)& =& \frac{1}{2c^2{\mu }_{\mathrm{o}}}\left({\tilde{E}}^2+c^2{\tilde{B}}^2\right)u^a{u}^b+\frac{1}{c^2{\mu }_{\mathrm{o}}}\left(u^a{\widetilde{ϵ}}^{bcd}+u^b{\widetilde{ϵ}}^{acd}\right){\tilde{E}}_c{\tilde{B}}_d\hfill \end{array} \\ \begin{array}{ccc}& & -\frac{1}{c^2{\mu }_{\mathrm{o}}}\left[{\tilde{E}}^a{\tilde{E}}^b+c^2{\tilde{B}}^a{\tilde{B}}^b-\frac{1}{2}\left({\tilde{E}}^2+c^2{\tilde{B}}^2\right)h^{ab}\right].\hfill \end{array} \end{gathered}$$

(105)

We end this subsection by pointing out that Eqs. (99) and (105) shows that $T_{ab}$ naturally separates into “time-time”, “time-space”, and “space-space” pieces. Respectively, these give the total mass-energy density E, the total momentum density $P_a$, and the total stress $S_{ab}$:

$$\begin{array}{cc}\hfill E& =\frac{1}{c^2}{u}^a{u}^b{T}_{ab},\hfill \end{array}$$

(106a)

$$\begin{array}{cc}\hfill P^a& =-\frac{1}{c}{u}_b{\perp }_c^a{T}^{bc}=-\frac{1}{c}{u}_b{\perp }_c^a{T}^{cb},\hfill \end{array}$$

(106b)

$$\begin{array}{cc}\hfill S^{ab}& ={\perp }_c^a{\perp }_d^b{T}^{cd},S=h_{ab}{S}^{ab}.\hfill \end{array}$$

(106c)

The terms in Eqs. (99) and (105) combine to give

$$\begin{array}{cc}\hfill E& =-\Lambda +\sum _{\mathrm{x}}{\tilde{\mu }}_a^{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}^a+\frac{1}{2c^2{\mu }_{\mathrm{o}}}\left({\tilde{E}}^2+c^2{\tilde{B}}^2\right),\hfill \end{array}$$

(107a)

$$\begin{gathered} \begin{array}{cc}\hfill P^a& =\frac{1}{c}\sum _{\mathrm{x}}{\tilde{\mu }}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^a+\frac{1}{c{\mu }_{\mathrm{o}}}{\widetilde{ϵ}}^{abc}{\tilde{E}}_b{\tilde{B}}_c\hfill \end{array} \\ \begin{array}{cc}\hfill & =c\sum _{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}{\tilde{\mu }}_{\mathrm{x}}^a+\frac{1}{c{\mu }_{\mathrm{o}}}{\widetilde{ϵ}}^{abc}{\tilde{E}}_b{\tilde{B}}_c,\hfill \end{array} \end{gathered}$$

(107b)

$$\begin{array}{cc}\hfill S^{ab}& =\Psi h^{ab}+\sum _{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}^a{\tilde{\mu }}_{\mathrm{x}}^b-\frac{1}{c^2{\mu }_{\mathrm{o}}}\left[{\tilde{E}}^a{\tilde{E}}^b+c^2{\tilde{B}}^a{\tilde{B}}^b-\frac{1}{2}\left({\tilde{E}}^2+c^2{\tilde{B}}^2\right)h^{ab}\right],\hfill \end{array}$$

(107c)

$$\begin{array}{cc}\hfill S& =3\Psi +\sum _{\mathrm{x}}{\tilde{\mu }}_a^{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}^a+\frac{1}{2c^2{\mu }_{\mathrm{o}}}\left({\tilde{E}}^2+c^2{\tilde{B}}^2\right).\hfill \end{array}$$

(107d)

5.3. The $3+1$ Field Equations

The logic of rewriting the Einstein, fluid/plasma, and electromagnetic field equations in their $3+1$ forms is the same as for the field variables — project free indices perpendicular to the foliation slices using the operator $u^a/c$ and project free indices parallel to the foliation slices using ${\perp }_b^a$, and then make substitutions of the decomposed quantities derived in the previous section. The main complication is that the field equations have derivatives, and we will need to replace everywhere covariant derivatives ${\nabla }_a$ with their $3+1$ counterparts $\mathrm{D}/\mathrm{d}t$ and $D_a$.

We will start with the Einstein equations as given in Eq. (62). The projections of the Ricci tensor $R_{ab}$ are performed in Appendix B. When these and the terms E, $P^a$, and $S^{ab}$ are substituted back into Eq. (62) we get the Hamiltonian Constraint

$${}^{\left(3\right)}R+\frac{1}{c^2}{K}^2-\frac{1}{c^2}{K}^{ab}{K}_{ab}=\frac{16\pi G}{c^4}E,$$

(108)

the Momentum Constraint

$$D_b\left(K^b{}_a-K{\delta }_a^b\right)=\frac{8\pi G}{c^3}{P}_a,$$

(109)

and finally an evolution equation

$$-\frac{1}{c^2}{\mathcal{L}}_u{K}_{ab}-\frac{1}{N}{D}_a{D}_{b}N+{}^{\left(3\right)}{R}_{ab}+\frac{1}{c^2}K{K}_{ab}-\frac{2}{c^2}{K}_{ac}{K}^c{}_b=\frac{8\pi G}{c^4}\left[S_{ab}-\frac{1}{2}\left(S-E\right)h_{ab}\right].$$

(110)

For the fluid/plasma equations, the results are long, and so it is better to break them up into individual pieces, and present them instead:

$$\begin{array}{cc}\hfill u^a{f}_a^{\mathrm{x}}& =-{\tilde{n}}_{\mathrm{x}}^a\left(D_a{\tilde{\mu }}_{\mathrm{x}}+\frac{\mathrm{D}{\tilde{\mu }}_a^{\mathrm{x}}}{\mathrm{d}t}\right)-\frac{1}{c^2}{\tilde{\mu }}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}^a{a}_a+K_{ab}{\tilde{n}}_{\mathrm{x}}^a{\tilde{\mu }}_{\mathrm{x}}^b+{\tilde{j}}_{\mathrm{x}}^a{\tilde{E}}_a,\hfill \end{array}$$

(111a)

$$\begin{gathered} \begin{array}{cc}\hfill {\perp }_a^b{f}_b^{\mathrm{x}}& ={\tilde{n}}_{\mathrm{x}}\left({\perp }_a^b\frac{\mathrm{D}{\tilde{\mu }}_b^{\mathrm{x}}}{\mathrm{d}t}+{\tilde{u}}_{\mathrm{x}}^b{D}_b{\tilde{\mu }}_a^{\mathrm{x}}\right)+{\tilde{n}}_{\mathrm{x}}{D}_a{\tilde{\mu }}_{\mathrm{x}}-{\tilde{n}}_{\mathrm{x}}^b{D}_a{\tilde{\mu }}_b^{\mathrm{x}}\hfill \end{array} \\ \begin{array}{cc}\hfill & +\frac{{\tilde{\mu }}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}}{c^2}{a}_a-{\tilde{n}}_{\mathrm{x}}{\tilde{\mu }}_{\mathrm{x}}^b{K}_{ba}-\left({\tilde{\sigma }}_{\mathrm{x}}{\tilde{E}}_a+{\widetilde{ϵ}}_{abc}{\tilde{j}}_{\mathrm{x}}^b{\tilde{B}}^c\right),\hfill \end{array} \end{gathered}$$

(111b)

$$\begin{array}{cc}\hfill \left(-u^a{\mu }_a^{\mathrm{x}}\right){\Gamma }_{\mathrm{x}}& ={\tilde{\mu }}_{\mathrm{x}}\left(D_a{\tilde{n}}_{\mathrm{x}}^a+\frac{\mathrm{D}{\tilde{n}}_{\mathrm{x}}}{\mathrm{d}t}\right)-K{\tilde{\mu }}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}+\frac{1}{c^2}{\tilde{\mu }}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}^a{a}_a,\hfill \end{array}$$

(111c)

$$\begin{array}{cc}\hfill \left({\perp }_a^b{\mu }_b^{\mathrm{x}}\right){\Gamma }_{\mathrm{x}}& =\left(D_b{\tilde{n}}_{\mathrm{x}}^b+\frac{\mathrm{D}{\tilde{n}}_{\mathrm{x}}}{\mathrm{d}t}\right){\tilde{\mu }}_a^{\mathrm{x}}-K{\tilde{n}}_{\mathrm{x}}{\tilde{\mu }}_a^{\mathrm{x}}+\frac{1}{c^2}{\tilde{n}}_{\mathrm{x}}^b{a}_b{\tilde{\mu }}_a^{\mathrm{x}},\hfill \end{array}$$

(111d)

$$\begin{array}{cc}\hfill u^a\left(R_a^{\mathrm{x}}-q_{\mathrm{x}}{\Gamma }_{\mathrm{x}}{A}_a\right)& =u^a{R}_a^{\mathrm{x}}+q_{\mathrm{x}}{\Gamma }_{\mathrm{x}}\tilde{V},\hfill \end{array}$$

(111e)

$$\begin{array}{cc}\hfill {\perp }_a^b\left(R_b^{\mathrm{x}}-q_{\mathrm{x}}{\Gamma }_{\mathrm{x}}{A}_b\right)& ={\perp }_a^b{R}_b^{\mathrm{x}}-q_{\mathrm{x}}{\Gamma }_{\mathrm{x}}{\tilde{A}}_a.\hfill \end{array}$$

(111f)

We will present a more detailed look at $u^a{R}_a^{\mathrm{x}}$ and ${\perp }_a^b{R}_b^{\mathrm{x}}$ later in Sec. Section 5.4.

Lastly, we have to evaluate the following projections of the Maxwell equations:

$$\begin{array}{cc}\hfill u_a{\nabla }_b{F}^{ab}& ={\mu }_{\mathrm{o}}{c}^2\sum _{\mathrm{x}=\{2,3\}}{u}_a{j}_{\mathrm{x}}^a,\hfill \end{array}$$

(112a)

$$\begin{array}{cc}\hfill {\perp }_c^a{\nabla }_b{F}^{cb}& ={\mu }_{\mathrm{o}}\sum _{\mathrm{x}=\{2,3\}}{\perp }_b^a{j}_{\mathrm{x}}^b,\hfill \end{array}$$

(112b)

$$\begin{array}{cc}\hfill u_a{ϵ}^{abcd}{\nabla }_{[b}{F}_{cd]}& =0,\hfill \end{array}$$

(112c)

$$\begin{array}{cc}\hfill {\perp }_e^a{ϵ}^{ebcd}{\nabla }_{[b}{F}_{cd]}& =0.\hfill \end{array}$$

(112d)

Before applying the projections, it is convenient to do a little preparatory work: take the covariant derivative of Eq. (104), and use Eq. (89) to get

$$\begin{gathered} \begin{array}{ccc}\hfill {\nabla }_a{F}_{bc}& =& \frac{1}{c^2}{u}_a\left(\frac{2}{c^2}{a}_{[c}{\tilde{E}}_{b]}+\frac{2}{c^2}{u}_{[c}{\perp }_{b]}^d\frac{\mathrm{D}{\tilde{E}}_d}{\mathrm{d}t}-\frac{1}{c}{ϵ}_{bc}{}^{de}{a}_d{\tilde{B}}_e-{\widetilde{ϵ}}_{bcd}\frac{\mathrm{D}{\tilde{B}}^d}{\mathrm{d}t}\right)\hfill \end{array} \\ \begin{array}{ccc}& & +\frac{2}{c^2}{u}_{[b}{D}_{\left|a\right|}{\tilde{E}}_{c]}+\frac{2}{c^2}{K}_{a[c}{\tilde{E}}_{b]}-\frac{1}{c}{ϵ}_{bc}{}^{de}{K}_{ad}{\tilde{B}}_e+{\widetilde{ϵ}}_{bcd}{D}_a{\tilde{B}}^d,\hfill \end{array} \end{gathered}$$

(113)

which, in turn, gives

$${\nabla }^b{F}_{ab}=\frac{1}{c^2}{u}_a{D}_b{\tilde{E}}^b-\frac{1}{c^2}{\perp }_a^b\frac{\mathrm{D}{\tilde{E}}_b}{\mathrm{d}t}+{\widetilde{ϵ}}_{abc}\left(D^b{\tilde{B}}^c+\frac{1}{c^2}{a}^b{\tilde{B}}^c\right)-\frac{1}{c^2}\left(K_{ab}-K{h}_{ab}\right){\tilde{E}}^b,$$

(114)

and

$$\begin{gathered} \begin{array}{ccc}\hfill \frac{1}{2}{ϵ}_a{}^{bcd}{\nabla }_b{F}_{cd}& =& -u_a{D}_b{\tilde{B}}^b+\frac{1}{c^2}{\perp }_a^b\frac{\mathrm{D}{\tilde{B}}_b}{\mathrm{d}t}+\frac{1}{c^2}{\widetilde{ϵ}}_{abc}\left(D^b{\tilde{E}}^c+\frac{1}{c^2}{a}^b{\tilde{E}}^c\right)\hfill \end{array} \\ \begin{array}{ccc}& & +\frac{1}{c^2}\left(K_{ab}-K{h}_{ab}\right){\tilde{B}}^b.\hfill \end{array} \end{gathered}$$

(115)

Therefore, the $u^a/c$ and ${\perp }_a^b$ projections of the Maxwell equations and the continuity equation are [34]

$$\begin{array}{cc}\hfill D_a{\tilde{E}}^a& ={\mu }_{\mathrm{o}}{c}^2\sum _{\mathrm{x}=\{2,3\}}{\tilde{\sigma }}_{\mathrm{x}},\hfill \end{array}$$

(116a)

$$\begin{array}{cc}\hfill {\widetilde{ϵ}}^a{}_{bc}\left(D^b{\tilde{B}}^c+\frac{1}{c^2}{a}^b{\tilde{B}}^c\right)& ={\mu }_{\mathrm{o}}\sum _{\mathrm{x}=\{2,3\}}{\tilde{j}}_{\mathrm{x}}^a+\frac{1}{c^2}{\perp }_b^a\frac{\mathrm{D}{\tilde{E}}^b}{\mathrm{d}t}+\frac{1}{c^2}\left(K^{ab}-K{h}^{ab}\right){\tilde{E}}_b,\hfill \end{array}$$

(116b)

$$\begin{array}{cc}\hfill D_b{\tilde{B}}^b& =0,\hfill \end{array}$$

(116c)

$$\begin{array}{cc}\hfill {\widetilde{ϵ}}_{abc}\left(D^b{\tilde{E}}^c+\frac{1}{c^2}{a}^b{\tilde{E}}^c\right)& =-{\perp }_a^b\frac{\mathrm{D}{\tilde{B}}_b}{\mathrm{d}t}-\left(K_{ab}-K{h}_{ab}\right){\tilde{B}}^b,\hfill \end{array}$$

(116d)

$$\begin{array}{cc}\hfill \sum _{\mathrm{x}=\{2,3\}}{\nabla }_a{j}_{\mathrm{x}}^a& =\sum _{\mathrm{x}=\{2,3\}}\left(D_a{\tilde{j}}_{\mathrm{x}}^a+\frac{\mathrm{D}{\tilde{\sigma }}_{\mathrm{x}}}{\mathrm{d}t}-K{\tilde{\sigma }}_{\mathrm{x}}+\frac{1}{c^2}{\tilde{j}}_{\mathrm{x}}^a{a}_a\right)=0.\hfill \end{array}$$

(116e)

5.4. Resistivities and Dissipation in the $3+1$ Formalism

We have now finished our development of the $3+1$ form of the full suite of field equations. This has been accomplished without having to make detailed statements about the specific dependence of $n_{ABC}^{\mathrm{x}}$ on $\{X_{\mathcal{P}}^A,X_{\mathcal{N}}^A\}$ nor, in turn, the specific dependence of $\Lambda$ on $n_{ABC}^{\mathrm{x}}$. In fact, we have taken the point of view that each of these are “known” a priori, meaning that once a specific application is considered the relevant forms and dependencies can, at least in principle, be constructed based on the relevant microphysics of the system. However, even without such an analysis, the action-based formalism has taken us a long way. This has been pointed out already by Andersson et al. [13]. They used this as a basic platform upon which resistivities could be built phenomenologically. Our purpose now is to give a review of the salient points, and then to apply them to the two-temperature extended system considered here.

We start by applying Eq. (49) to the $3+1$ decomposition of ${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}$, which is

$${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}={\tilde{\mathcal{R}}}^{\mathrm{x}\mathrm{y}}{u}_a/c^2+{\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}},{\tilde{\mathcal{R}}}^{\mathrm{x}\mathrm{y}}=-u^a{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}},{\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}={\perp }_a^b{\mathcal{R}}_b^{\mathrm{x}\mathrm{y}}.$$

(117)

By imposing Eq. (49) we find that ${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}$ becomes

$${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}=\left({\delta }_a^b+u_a{\tilde{u}}_{\mathrm{y}}^b/c^2\right){\tilde{\mathcal{R}}}_b^{\mathrm{x}\mathrm{y}},$$

(118)

and the resistivity $R_a^{\mathrm{x}}$ is

$$R_a^{\mathrm{x}}=\sum _{\mathrm{y}\ne \mathrm{x}}\left[\left({\tilde{u}}_{\mathrm{x}}^b{\tilde{\mathcal{R}}}_b^{\mathrm{y}\mathrm{x}}-{\tilde{u}}_{\mathrm{y}}^b{\tilde{\mathcal{R}}}_b^{\mathrm{x}\mathrm{y}}\right)u_a/c^2+{\tilde{\mathcal{R}}}_a^{\mathrm{y}\mathrm{x}}-{\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}\right].$$

(119)

Inserting this modified form for $R_a^{\mathrm{x}}$ into Eq. (68), we determine that the creation rate becomes

$${\Gamma }_{\mathrm{x}}=\frac{{\tilde{\gamma }}_{\mathrm{x}}}{{\overline{\mu }}_{\mathrm{x}}}\sum _{\mathrm{y}\ne \mathrm{x}}{\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}{\tilde{w}}_{\mathrm{x}\mathrm{y}}^a.$$

(120)

To make further progress, we impose three physical constraints — charge conservation, baryon number conservation, and the Second Law of Thermodynamics. The conservation of charge [cf.Eq. (35)] leads to

$$0=\sum _{\mathrm{x}=2,3}{e}^{\mathrm{x}}{\Gamma }_{\mathrm{x}}=\sum _{\mathrm{x}=2,3}\frac{q_{\mathrm{x}}{\tilde{\gamma }}_{\mathrm{x}}}{{\overline{\mu }}_{\mathrm{x}}}\sum _{\mathrm{y}\ne \mathrm{x}}{\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}{\tilde{w}}_{\mathrm{x}\mathrm{y}}^a,$$

(121)

while baryon number conservation [cf.Eq. (36)] says

$$0=\sum _{\mathrm{x}=1,2}{\Gamma }_{\mathrm{x}}=\sum _{\mathrm{x}=1,2}\frac{{\tilde{\gamma }}_{\mathrm{x}}}{{\overline{\mu }}_{\mathrm{x}}}\sum _{\mathrm{y}\ne \mathrm{x}}{\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}{\tilde{w}}_{\mathrm{x}\mathrm{y}}^a.$$

(122)

The Second Law of Thermodynamics gives the inequality

$$\sum _{\mathrm{x}=4,5,6}{\Gamma }_{\mathrm{x}}=\sum _{\mathrm{x}=4,5,6}\frac{{\tilde{\gamma }}_{\mathrm{x}}}{{\overline{\mu }}_{\mathrm{x}}}\sum _{\mathrm{y}\ne \mathrm{x}}{\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}{\tilde{w}}_{\mathrm{x}\mathrm{y}}^a\ge 0.$$

(123)

In order to satisfy these, we need to be more specific about the terms, meaning that we will now make an ansatz about the form of the resistivity and flux creation rates, but in a manner which is consistent with the overall formalism.

Onsager [35] (see also [36,37]) developed an approach that relies on the notions of thermodynamic fluxes and forces. In our case, the thermodynamic fluxes are the ${\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}$, and the thermodynamic forces are the ${\tilde{w}}_{\mathrm{x}\mathrm{y}}^a$. The key step is to combine the fluxes and forces in such a way that they tend to drive the system towards a dynamical equilibrium where the relative flows are zero and a thermodynamical equilibrium where ${\Gamma }_{\mathrm{x}}\to 0$ all the while maintaining the inequality of Eq. (123).

We begin with an obvious choice for the ${\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}$, which is to write

$${\tilde{\mathcal{R}}}_a^{\mathrm{x}\mathrm{y}}={\tilde{r}}_{\mathrm{x}\mathrm{y}}{\tilde{w}}_a^{\mathrm{x}\mathrm{y}}⟹{\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}={\tilde{r}}_{\mathrm{x}\mathrm{y}}\left({\delta }_a^b+u_a{\tilde{u}}_{\mathrm{y}}^b/c^2\right){\tilde{w}}_b^{\mathrm{x}\mathrm{y}}.$$

(124)

This causes the sum for the total entropy creation rate to be over the set of positive-definite terms given by $w_{\mathrm{x}\mathrm{y}}^a{w}_a^{\mathrm{x}\mathrm{y}}$. Because the relation for ${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}$ is linear in the ${\tilde{r}}_{\mathrm{x}\mathrm{y}}$, then we can reduce the number of ${\tilde{r}}_{\mathrm{x}\mathrm{y}}$ by imposing that (in their indices) they satisfy the same equalities that the ${\mathcal{R}}_a^{\mathrm{x}\mathrm{y}}$ do in Eqs. (63a)–(63f). Noting that

$$\begin{array}{cc}\hfill {\Gamma }_1& =2\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\overline{\mu }}_1}{\tilde{r}}_{13}{\tilde{w}}_a^{\mathcal{P}\mathcal{N}}{\tilde{w}}_{\mathcal{P}\mathcal{N}}^a,\hfill \end{array}$$

(125a)

$$\begin{array}{cc}\hfill {\Gamma }_2& =2\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\overline{\mu }}_2}{\tilde{r}}_{23}{\tilde{w}}_a^{\mathcal{P}\mathcal{N}}{\tilde{w}}_{\mathcal{P}\mathcal{N}}^a,\hfill \end{array}$$

(125b)

$$\begin{array}{cc}\hfill {\Gamma }_3& =4\frac{{\tilde{\gamma }}_{\mathcal{N}}}{{\overline{\mu }}_3}{\tilde{r}}_{31}{\tilde{w}}_a^{\mathcal{P}\mathcal{N}}{\tilde{w}}_{\mathcal{P}\mathcal{N}}^a,\hfill \end{array}$$

(125c)

$$\begin{array}{cc}\hfill {\Gamma }_4& =2\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\overline{\mu }}_4}{\tilde{r}}_{43}{\tilde{w}}_a^{\mathcal{P}\mathcal{N}}{\tilde{w}}_{\mathcal{P}\mathcal{N}}^a,\hfill \end{array}$$

(125d)

$$\begin{array}{cc}\hfill {\Gamma }_5& =2\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\overline{\mu }}_5}{\tilde{r}}_{53}{\tilde{w}}_a^{\mathcal{P}\mathcal{N}}{\tilde{w}}_{\mathcal{P}\mathcal{N}}^a,\hfill \end{array}$$

(125e)

$$\begin{array}{cc}\hfill {\Gamma }_6& =4\frac{{\tilde{\gamma }}_{\mathcal{N}}}{{\overline{\mu }}_6}{\tilde{r}}_{61}{\tilde{w}}_a^{\mathcal{P}\mathcal{N}}{\tilde{w}}_{\mathcal{P}\mathcal{N}}^a,\hfill \end{array}$$

(125f)

we can reduce again the number of ${\tilde{r}}_{\mathrm{x}\mathrm{y}}$, by imposing charge [cf.Eq. (35)] and baryon number [cf.Eq. (36)] conservation, since they imply

$$\begin{array}{cc}\hfill {\tilde{r}}_{23}& =-\frac{{\overline{\mu }}_2}{{\overline{\mu }}_1}{\tilde{r}}_{13},\hfill \end{array}$$

(126a)

$$\begin{array}{cc}\hfill {\tilde{r}}_{31}& =-\frac{1}{2}\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\tilde{\gamma }}_{\mathcal{N}}}\frac{{\overline{\mu }}_3}{{\overline{\mu }}_1}{\tilde{r}}_{13}.\hfill \end{array}$$

(126b)

The Second Law of Thermodynamics [cf.Eq. (123)] implies that the coefficients must satisfy

$$\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\overline{\mu }}_4}{\tilde{r}}_{43}+\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\overline{\mu }}_5}{\tilde{r}}_{53}+2\frac{{\tilde{\gamma }}_{\mathcal{N}}}{{\overline{\mu }}_6}{\tilde{r}}_{61}\ge 0.$$

(127)

The independent resistivity vector takes the final form

$$R_a^{\mathcal{N}}=2\left[\left(2{\tilde{r}}_{61}-\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\tilde{\gamma }}_{\mathcal{N}}}\frac{{\overline{\mu }}_3}{{\overline{\mu }}_1}{\tilde{r}}_{13}\right)\frac{{\tilde{w}}_{\mathcal{P}\mathcal{N}}^b{\tilde{w}}_b^{\mathcal{P}\mathcal{N}}}{c^2}+{\tilde{r}}_{\mathcal{N}}\frac{{\tilde{u}}_{\mathcal{N}}^b{\tilde{w}}_b^{\mathcal{P}\mathcal{N}}}{c^2}\right]u_a+2{\tilde{r}}_{\mathcal{N}}{\tilde{w}}_a^{\mathcal{P}\mathcal{N}},$$

(128)

where

$${\tilde{r}}_{\mathcal{N}}=2{\tilde{r}}_{61}+\left(1-\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\tilde{\gamma }}_{\mathcal{N}}}\frac{{\overline{\mu }}_3}{{\overline{\mu }}_1}-\frac{{\overline{\mu }}_2}{{\overline{\mu }}_1}\right){\tilde{r}}_{13}+{\tilde{r}}_{43}+{\tilde{r}}_{53}.$$

(129)

We see that our final model requires the four coefficients $\{{\tilde{r}}_{13},{\tilde{r}}_{43},{\tilde{r}}_{53},{\tilde{r}}_{61}\}$ to completely determine the creation rates ${\Gamma }_{\mathrm{x}}$ and the independent resistivity $R_a^{\mathcal{N}}$. Notably, as ${\tilde{w}}_{\mathcal{P}\mathcal{N}}^a\to 0$ (all the fluids are comoving) then $R_a^{\mathcal{N}}\to 0$ and ${\Gamma }_{\mathrm{x}}\to 0$. Any further development of this model would require microscopic modeling of specific systems to determine the four coefficients.

6. The “Newtonian” Limit

In order to make contact with existing work on two-temperature plasmas, which is mainly in the Newtonian setting, we will now work out the “Newtonian limit” of our equations. Poisson and Will [38] point out that when gravity is formulated as a metric theory, then the limit we are imposing is to be understood as the first-order correction to flat spacetime, which is not, a priori, the same thing as Newtonian gravity, which is based on forces and action-at-a-distance.

Our definition of the “Newtonian limit” includes the following criteria: a) The particles are moving much slower than the speed of light c; b) the gravitational field is “weak”, meaning it is a linear perturbation away from flat spacetime ($R^c{}_{dab}=0$); and, c) the gravitational field is static. The latter two criteria will be imposed by an expansion of N, $N^a$, and $h_{ab}$ away from flat spacetime. Some of this work is presented in Appendix C, where we have taken the $3+1$ formulas, and adapted them to a coordinate system such that the time coordinate ${\overline{x}}^0=ct$, where recall $t\left({\overline{x}}^a\right)$ is the scalar field from which the spacelike hypersurfaces of the foliation are constructed.

It is still an open question as to whether or not Newtonian gravity is a subset of this limit of General Relativity, or if it is all inclusive. Philosophical issues aside, we take a practical point-of-view, which is to impose the criteria a), b), and c) above on the field equations and thereby extract the terms which formally survive the limit. It then becomes a question of the particular physical scenario to which the field equations are being applied as to whether or not all of the remaining terms are required.

6.1. The Metric Expansion and Linear Corrections to Flat Spacetime

In order to take the Newtonian gravitational limit of the Einstein equations, we will need to analyze the left- and right-hand-sides separately. Here we will be setting up the left-hand-sides of the Hamiltonian and Momentum constraints — Eqs. (108) and (109), respectively — and the evolution Eq. (110). We simplify the equations by taking the ${\overline{x}}^i$ to be Cartesian coordinates.

A linear expansion of the metric away from flat spacetime takes the form

$$g_{ab}={\eta }_{ab}+\delta g_{ab},$$

(130)

where ${\eta }_{ab}=\mathrm{diag}\left[-1,1,1,1\right]$ is the Minkowski metric and the components of $\delta g_{ab}$ are taken to be small, meaning that we ignore any terms of the form $\delta g_{ab}\delta g_{cd}$, $\delta g_{ab}{\nabla }_c\delta g_{de}$, and so on. The flat-spacetime pieces of the metric are $N=c$, $N^i=0$, and $h_{ij}={\delta }_{ij}=\mathrm{diag}\left[1,1,1\right]$. The flat spacetime plus linear perturbations metric pieces are

$$\begin{array}{cc}\hfill N& =c+\delta N,\hfill \end{array}$$

(131a)

$$\begin{array}{cc}\hfill N^i& =\delta N^i,\delta N_i={\delta }_{ij}\delta N^i,\hfill \end{array}$$

(131b)

$$\begin{array}{cc}\hfill h_{ij}& ={\delta }_{ij}+\delta h_{ij}.\hfill \end{array}$$

(131c)

These expansions will be inserted into the left-hand-sides of Eqs. (108), (109), and (110), keeping only the first-order terms.

But before we take that step, it is important to note that the Einstein equations have a “gauge” symmetry that basically comes from their coordinate invariance (or, more formally, diffeomorphism invariance). We employ that here by using the harmonic gauge, which takes the form

$${\partial }_b\left(\delta g^{ba}-\frac{1}{2}{\eta }^{ba}{\eta }^{cd}\delta g_{cd}\right)=-\left({\eta }^{ac}{\eta }^{jd}+\frac{1}{2}{\eta }^{aj}{\eta }^{cd}\right){\partial }_j\delta g_{cd}=0,$$

(132)

where we have used

$$\delta g^{ab}=-{\eta }^{ac}{\eta }^{bd}\delta g_{cd}.$$

(133)

In terms of the $3+1$ decomposition, we have

$$\begin{array}{ccc}\hfill \delta g_{ab}& =& \left[\begin{array}{cc}-\frac{2}{c}\delta N& \frac{1}{c}\delta N_i\\ \frac{1}{c}\delta N_i& \delta h_{ij}\end{array}\right],\hfill \end{array}$$

(134)

and so the gauge condition leads to

$$\begin{array}{cc}\hfill 0& ={\partial }_i\delta N^i,\hfill \end{array}$$

(135a)

$$\begin{array}{cc}\hfill 0& ={\partial }^j\delta h_{ij}+{\partial }_i\left(\frac{1}{c}\delta N+\frac{1}{2}\delta h\right),\hfill \end{array}$$

(135b)

where $\delta h={\delta }^{ij}\delta h_{ij}$. The unit normal to the hypersurfaces $u^a$, the acceleration $a_a$, and non-zero components of the projection operator ${\perp }_a^b$ become, respectively,

$$\begin{array}{cc}\hfill u^a& =\left(c-\delta N,\delta N^i\right),u_a=\left(-c-\delta N,0,0,0\right),\hfill \end{array}$$

(136a)

$$\begin{array}{cc}\hfill a^a& =\left(0,c{\partial }^i\delta N\right),a_a=\left(0,c{\partial }_i\delta N\right),\hfill \end{array}$$

(136b)

$$\begin{array}{cc}\hfill {\perp }_0^i& =\frac{1}{c}\delta N^i,{\perp }_i^j={\delta }_i^j.\hfill \end{array}$$

(136c)

In order to build ${}^{\left(3\right)}{R}_{lj}$, we need to know the ${}^{\left(3\right)}{\Gamma }_{jk}^i$. Taking Eq. (C.10), and substituting in the expansions above, while keeping only the linear terms, we find

$${}^{\left(3\right)}{\Gamma }_{jk}^i=\frac{1}{2}{\delta }^{il}\left({\partial }_j\delta h_{lk}+{\partial }_k\delta h_{lj}-{\partial }_l\delta h_{jk}\right).$$

(137)

The gauge choice leads to $K=0$, but there remain linear-order $K_{ij}$ terms, which are

$$K_{ij}=\frac{1}{2}\left({\partial }_i\delta N_j+{\partial }_j\delta N_i\right).$$

(138)

We find that the linearized forms for ${}^{\left(3\right)}{R}_{lj}$ and ${}^{\left(3\right)}R$ are

$$\begin{array}{cc}\hfill {}^{\left(3\right)}{R}_{ij}& ={}^{\left(3\right)}{\Gamma }_{ij,k}^k-{}^{\left(3\right)}{\Gamma }_{ik,j}^k=-{\partial }_i{\partial }_j\left(\frac{1}{c}\delta N+\delta h\right)-\frac{1}{2}{\partial }_k{\partial }^k\delta h_{ij},\hfill \end{array}$$

(139a)

$$\begin{array}{cc}\hfill {}^{\left(3\right)}R& ={\delta }^{ij}{}^{\left(3\right)}{R}_{ij}=-{\partial }_i{\partial }^i\left(\frac{1}{c}\delta N+\frac{3}{2}\delta h\right),\hfill \end{array}$$

(139b)

The left-hand-sides of Eqs. (108), (109), and (110), respectively, now become

$$\begin{array}{cc}\hfill {}^{\left(3\right)}R+\frac{1}{c^2}{K}^2-\frac{1}{c^2}{K}^{ab}{K}_{ab}& =-{\partial }_i{\partial }^i\left(\frac{1}{c}\delta N+\frac{3}{2}\delta h\right),\hfill \end{array}$$

(140a)

$$\begin{array}{cc}\hfill D_j\left(K^j{}_i-K{\delta }_i^j\right)& =\frac{1}{2}{\partial }_j{\partial }^j\delta N_i,\hfill \end{array}$$

(140b)

$$\begin{gathered} \begin{array}{cc}\hfill -\frac{1}{c^2}{\mathcal{L}}_u{K}_{ij}-\frac{1}{N}{D}_i{D}_{j}N& +{}^{\left(3\right)}{R}_{ij}+\frac{1}{c^2}K{K}_{ij}-\frac{2}{c^2}{K}_{ik}{K}^k{}_j\hfill \end{array} \\ \begin{array}{cc}\hfill & =-{\partial }_i{\partial }_j\left(\frac{2}{c}\delta N+\delta h\right)-\frac{1}{2}{\partial }_k{\partial }^k\delta h_{ij}.\hfill \end{array} \end{gathered}$$

(140c)

6.2. Newtonian Limit of the Fluid/Plasma and $3+1$ Energy-Momentum-Stress Tensor Components

The main approximations for the flux variables are that their relative speeds ${\tilde{u}}^a$ must be much less than the speed of light—we neglect terms of order $\mathrm{O}\left({\tilde{u}}_{\mathrm{x}}^2/c^2\right)$ and higher—and energies that scale with $c^2$ (such as the rest-mass energy densities $m_{\mathrm{x}}{c}^2{n}_{\mathrm{x}}$) are much bigger than other energy densities. The typical leading-order terms in the Master function $\Lambda$ are the rest-mass energy densities, and so it is convenient to re-fashion $\Lambda$ as a sum of $m_{\mathrm{x}}{c}^2{n}_{\mathrm{x}}$ and an “internal energy” density $\mathcal{U}$ (having the same functional dependence as $\Lambda$):

$$-\Lambda =\sum _{\mathrm{x}}{m}_{\mathrm{x}}{c}^2{n}_{\mathrm{x}}+\mathcal{U}.$$

(141)

We assume that entropy has zero rest-mass, but because of entrainment, it does have an effective mass with a leading-order term proportional to $c^2$ and it enters the field equations through its inclusion in $\mathcal{U}$ [cf. Eq. ()].

We need to first consider the Newtonian limit of the momenta, as given in Eq. (94), but with the ${\mathcal{B}}^{\mathrm{x}}$ and ${\mathcal{A}}^{\mathrm{x}\mathrm{y}}$ computed using the rewritten $\Lambda$ of Eq. (141). We will also reintroduce the notation that splits the particle number fluxes into the matter $n_{\mathrm{x}}^a$ and the entropy $s_{\overline{\mathrm{x}}}^a$ pieces, and the momenta into ${\mu }_a^{\mathrm{x}}$ and ${\Theta }_a^{\overline{\mathrm{x}}}$. Here, the constituent indices for the matter are without a bar and range over $\mathrm{x},\mathrm{y},\dots =\{\eta ,\mathcal{P},\mathcal{N}\}$, whereas the indices with a bar are for the thermal pieces and range over $\overline{\mathrm{x}},\overline{\mathrm{y}},\dots =\{\overline{\eta },\overline{\mathcal{P}},\overline{\mathcal{N}}\}$. In order to generate the momentum coefficients, we have five different sets of scalars which can appear in the $\Lambda$: the first two are $n_{\mathrm{x}}^2=-g_{ab}{n}_{\mathrm{x}}^a{n}_{\mathrm{x}}^b/c^2$ and $n_{\mathrm{x}\mathrm{y}}^2=-g_{ab}{n}_{\mathrm{x}}^a{n}_{\mathrm{y}}^b/c^2=n_{\mathrm{y}\mathrm{x}}^2$, for which $\mathrm{y}\ne \mathrm{x}$; the next two are $s_{\overline{\mathrm{x}}}^2=-g_{ab}{s}_{\overline{\mathrm{x}}}^a{s}_{\overline{\mathrm{x}}}^b/c^2$, $s_{\overline{\mathrm{x}}\overline{\mathrm{y}}}^2=-g_{ab}{s}_{\overline{\mathrm{x}}}^a{s}_{\overline{\mathrm{y}}}^b/c^2=s_{\overline{\mathrm{y}}\overline{\mathrm{x}}}^2$, for which $\overline{\mathrm{y}}\ne \overline{\mathrm{x}}$; and the last is the mixed term $m_{\mathrm{x}\overline{\mathrm{y}}}^2=-g_{ab}{n}_{\mathrm{x}}^a{s}_{\overline{\mathrm{y}}}^b/c^2=m_{\overline{\mathrm{y}}\mathrm{x}}^2$.

A variation of the re-formulated $\Lambda$ yields the coefficients

$$\begin{array}{cc}\hfill {\mathcal{B}}^{\mathrm{x}}& =\frac{m_{\mathrm{x}}}{n_{\mathrm{x}}}+\frac{1}{c^2{n}_{\mathrm{x}}}\frac{\partial \mathcal{U}}{\partial n_{\mathrm{x}}},\hfill \end{array}$$

(142a)

$$\begin{array}{cc}\hfill {\mathcal{S}}^{\overline{\mathrm{x}}}& =\frac{1}{c^2{s}_{\overline{\mathrm{x}}}}\frac{\partial \mathcal{U}}{\partial s_{\overline{\mathrm{x}}}},\hfill \end{array}$$

(142b)

$$\begin{array}{cc}\hfill {\mathcal{B}}^{\mathrm{x}\mathrm{y}}& =\frac{1}{c^2}\frac{\partial \mathcal{U}}{\partial n_{\mathrm{x}\mathrm{y}}^2},\hfill \end{array}$$

(142c)

$$\begin{array}{cc}\hfill {\mathcal{S}}^{\overline{\mathrm{x}}\overline{\mathrm{y}}}& =\frac{1}{c^2}\frac{\partial \mathcal{U}}{\partial s_{\overline{\mathrm{x}}\overline{\mathrm{y}}}^2},\hfill \end{array}$$

(142d)

$$\begin{array}{cc}\hfill {\mathcal{M}}^{\mathrm{x}\overline{\mathrm{y}}}& =\frac{1}{c^2}\frac{\partial \mathcal{U}}{\partial m_{\mathrm{x}\overline{\mathrm{y}}}^2},\hfill \end{array}$$

(142e)

which combine together to give

$$\begin{array}{cc}\hfill {\mu }_a^{\mathrm{x}}& =\frac{1}{n_{\mathrm{x}}}\left(m_{\mathrm{x}}+\frac{1}{c^2}\frac{\partial \mathcal{U}}{\partial n_{\mathrm{x}}}\right)n_a^{\mathrm{x}}+\sum _{\mathrm{y}\ne \mathrm{x}}{\mathcal{B}}^{\mathrm{x}\mathrm{y}}{n}_a^{\mathrm{y}}+\sum _{\overline{\mathrm{y}}}{\mathcal{M}}^{\mathrm{x}\overline{\mathrm{y}}}{s}_a^{\overline{\mathrm{y}}},\hfill \end{array}$$

(143a)

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{x}}& =\left(m_{\mathrm{x}}+\sum _{\mathrm{y}\ne \mathrm{x}}{\mathcal{B}}^{\mathrm{x}\mathrm{y}}\frac{n_{\mathrm{x}\mathrm{y}}^2}{n_{\mathrm{x}}}+\sum _{\overline{\mathrm{y}}}{\mathcal{M}}^{\mathrm{x}\overline{\mathrm{y}}}\frac{m_{\mathrm{x}\overline{\mathrm{y}}}^2}{s_{\overline{\mathrm{y}}}}\right)c^2+\frac{\partial \mathcal{U}}{\partial n_{\mathrm{x}}},\hfill \end{array}$$

(143b)

$$\begin{array}{cc}\hfill {\Theta }_a^{\overline{\mathrm{x}}}& ={\mathcal{S}}^{\overline{\mathrm{x}}}{s}_a^{\overline{\mathrm{x}}}+\sum _{\overline{\mathrm{y}}\ne \overline{\mathrm{x}}}{\mathcal{S}}^{\overline{\mathrm{x}}\overline{\mathrm{y}}}{s}_a^{\overline{\mathrm{y}}}+\sum _{\mathrm{y}}{\mathcal{M}}^{\overline{\mathrm{x}}\mathrm{y}}{n}_a^{\mathrm{y}},\hfill \end{array}$$

(143c)

$$\begin{array}{cc}\hfill T_{\overline{\mathrm{x}}}& =\left(\sum _{\overline{\mathrm{y}}\ne \overline{\mathrm{x}}}{\mathcal{S}}^{\overline{\mathrm{x}}\overline{\mathrm{y}}}\frac{s_{\overline{\mathrm{x}}\overline{\mathrm{y}}}^2}{s_{\overline{\mathrm{x}}}}+\sum _{\mathrm{y}}{\mathcal{M}}^{\overline{\mathrm{x}}\mathrm{y}}\frac{m_{\overline{\mathrm{x}}\mathrm{y}}^2}{s_{\overline{\mathrm{x}}}}\right)c^2+\frac{\partial \mathcal{U}}{\partial s_{\overline{\mathrm{x}}}}.\hfill \end{array}$$

(143d)

In $3+1$ form we have

$$\begin{array}{cc}\hfill {\tilde{\mu }}_a^{\mathrm{x}}& =\frac{{\tilde{\mu }}_{\mathrm{x}}}{c^2}{\tilde{u}}_a^{\mathrm{x}}+\sum _{\mathrm{y}\ne \mathrm{x}}{\mathcal{B}}^{\mathrm{x}\mathrm{y}}{\tilde{n}}_{\mathrm{y}}{\tilde{w}}_a^{\mathrm{y}\mathrm{x}}+\sum _{\overline{\mathrm{y}}}{\mathcal{M}}^{\mathrm{x}\overline{\mathrm{y}}}{\tilde{s}}_{\overline{\mathrm{y}}}{\tilde{w}}_a^{\overline{\mathrm{y}}\mathrm{x}},\hfill \end{array}$$

(144a)

$$\begin{array}{cc}\hfill {\tilde{\mu }}_{\mathrm{x}}& ={\tilde{m}}_{\mathrm{x}}{c}^2+\frac{\partial \mathcal{U}}{\partial n_{\mathrm{x}}}{\tilde{\gamma }}_{\mathrm{x}},\hfill \end{array}$$

(144b)

$$\begin{array}{cc}\hfill {\tilde{\Theta }}_a^{\mathrm{x}}& =\frac{{\tilde{T}}_{\overline{\mathrm{x}}}}{c^2}{\tilde{u}}_a^{\overline{\mathrm{x}}}+\sum _{\overline{\mathrm{y}}\ne \overline{\mathrm{x}}}{\mathcal{S}}^{\overline{\mathrm{x}}\overline{\mathrm{y}}}{\tilde{s}}_{\overline{\mathrm{y}}}{\tilde{w}}_a^{\overline{\mathrm{y}}\overline{\mathrm{x}}}+\sum _{\mathrm{y}}{\mathcal{M}}^{\overline{\mathrm{x}}\mathrm{y}}{\tilde{n}}_{\mathrm{y}}{\tilde{w}}_a^{\mathrm{y}\overline{\mathrm{x}}},\hfill \end{array}$$

(144c)

$$\begin{array}{cc}\hfill {\tilde{T}}_{\overline{\mathrm{x}}}& ={\tilde{m}}_{\overline{\mathrm{x}}}{c}^2/k_B+\frac{\partial \mathcal{U}}{\partial s_{\overline{\mathrm{x}}}},\hfill \end{array}$$

(144d)

where we have defined

$$\begin{array}{cc}\hfill {\tilde{m}}_{\mathrm{x}}& =m_{\mathrm{x}}{\tilde{\gamma }}_{\mathrm{x}}+\sum _{\mathrm{y}\ne \mathrm{x}}{\mathcal{B}}^{\mathrm{x}\mathrm{y}}{\tilde{n}}_{\mathrm{y}}+\sum _{\overline{\mathrm{y}}}{\mathcal{M}}^{\mathrm{x}\overline{\mathrm{y}}}{\tilde{s}}_{\overline{\mathrm{y}}},\hfill \end{array}$$

(145a)

$$\begin{array}{cc}\hfill {\tilde{m}}_{\overline{\mathrm{x}}}& =k_B\left(\sum _{\overline{\mathrm{y}}\ne \overline{\mathrm{x}}}{\mathcal{S}}^{\overline{\mathrm{x}}\overline{\mathrm{y}}}{\tilde{s}}_{\overline{\mathrm{y}}}+\sum _{\mathrm{y}}{\mathcal{M}}^{\overline{\mathrm{x}}\mathrm{y}}{\tilde{n}}_{\mathrm{y}}\right).\hfill \end{array}$$

(145b)

We can get a handle on the lowest order impact of the condition ${\tilde{u}}_a^{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^a\ll c^2$ by expanding the parameters $n_{\mathrm{x}}^2$, $n_{\mathrm{x}\mathrm{y}}^2$, $s_{\overline{\mathrm{x}}}^2$, $s_{\overline{\mathrm{x}}\overline{\mathrm{y}}}^2$, and $m_{\mathrm{x}\overline{\mathrm{y}}}^2$:

$$\begin{array}{cc}\hfill n_{\mathrm{x}}^2& ={\tilde{n}}_{\mathrm{x}}^2{\tilde{\gamma }}_{\mathrm{x}}^{-2}\approx {\tilde{n}}_{\mathrm{x}}^2\left(1-{\tilde{u}}_i^{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^i/c^2\right),\hfill \end{array}$$

(146a)

$$\begin{array}{cc}\hfill n_{\mathrm{x}\mathrm{y}}^2& =-\frac{1}{c^2}{n}_{\mathrm{x}}{n}_{\mathrm{y}}{g}_{ab}{u}_{\mathrm{x}}^a{u}_{\mathrm{y}}^b\approx {\tilde{n}}_{\mathrm{x}}{\tilde{n}}_{\mathrm{y}}\left(1-{\tilde{u}}_k^{\mathrm{x}}{\tilde{u}}_{\mathrm{y}}^k/c^2\right),\hfill \end{array}$$

(146b)

$$\begin{array}{cc}\hfill s_{\overline{\mathrm{x}}}^2& ={\tilde{s}}_{\overline{\mathrm{x}}}^2{\tilde{\gamma }}_{\overline{\mathrm{x}}}^{-2}\approx {\tilde{s}}_{\overline{\mathrm{x}}}^2\left(1-{\tilde{u}}_i^{\overline{\mathrm{x}}}{\tilde{u}}_{\overline{\mathrm{x}}}^i/c^2\right),\hfill \end{array}$$

(146c)

$$\begin{array}{cc}\hfill s_{\overline{\mathrm{x}}\overline{\mathrm{y}}}^2& =-\frac{1}{c^2}{s}_{\overline{\mathrm{x}}}{s}_{\overline{\mathrm{y}}}{g}_{ab}{u}_{\overline{\mathrm{x}}}^a{u}_{\overline{\mathrm{y}}}^b\approx {\tilde{s}}_{\overline{\mathrm{x}}}{\tilde{s}}_{\overline{\mathrm{y}}}\left(1-{\tilde{u}}_k^{\overline{\mathrm{x}}}{\tilde{u}}_{\overline{\mathrm{y}}}^k/c^2\right),\hfill \end{array}$$

(146d)

$$\begin{array}{cc}\hfill m_{\mathrm{x}\overline{\mathrm{y}}}^2& =-\frac{1}{c^2}{n}_{\mathrm{x}}{s}_{\overline{\mathrm{y}}}{g}_{ab}{u}_{\mathrm{x}}^a{u}_{\overline{\mathrm{y}}}^b\approx {\tilde{n}}_{\mathrm{x}}{\tilde{s}}_{\overline{\mathrm{y}}}\left(1-{\tilde{u}}_k^{\mathrm{x}}{\tilde{u}}_{\overline{\mathrm{y}}}^k/c^2\right).\hfill \end{array}$$

(146e)

We see from this that the differences ${\tilde{n}}_{\mathrm{x}}^2-n_{\mathrm{x}}^2$, ${\tilde{n}}_{\mathrm{x}}{\tilde{n}}_{\mathrm{y}}-n_{\mathrm{x}\mathrm{y}}^2$, etc. are small. The expansion of $\Lambda$ gives

$$\begin{gathered} \begin{array}{ccc}\hfill -\Lambda & \approx & \sum _{\mathrm{x}}{m}_{\mathrm{x}}{c}^2{\tilde{n}}_{\mathrm{x}}+{\mathcal{U}}_{\mathrm{o}}\left({\tilde{n}}_{\mathrm{x}}^2,{\tilde{s}}_{\overline{\mathrm{x}}}^2\right)\hfill \end{array} \\ \begin{array}{ccc}& & -\frac{1}{2}\sum _{\mathrm{x}}\left[\frac{1}{{\tilde{n}}_{\mathrm{x}}}\left(m_{\mathrm{x}}+{\left.\frac{\partial \left(\mathcal{U}/c^2\right)}{\partial n_x}\right|}_{\mathrm{o}}\right){\tilde{n}}_i^{\mathrm{x}}+\sum _{\mathrm{y}\ne \mathrm{x}}{\mathcal{B}}_{\mathrm{o}}^{\mathrm{x}\mathrm{y}}{\tilde{n}}_i^{\mathrm{y}}+\sum _{\overline{\mathrm{y}}}{\mathcal{M}}_{\mathrm{o}}^{\mathrm{x}\overline{\mathrm{y}}}{\tilde{s}}_i^{\overline{\mathrm{y}}}\right]{\tilde{n}}_{\mathrm{x}}^i\hfill \end{array} \\ \begin{array}{ccc}& & -\frac{1}{2}\sum _{\overline{\mathrm{x}}}\left[\frac{1}{{\tilde{s}}_{\mathrm{x}}}{\left.\frac{\partial \left(\mathcal{U}/c^2\right)}{\partial s_{\overline{\mathrm{x}}}}\right|}_{\mathrm{o}}{\tilde{s}}_i^{\overline{\mathrm{x}}}+\sum _{\overline{\mathrm{y}}\ne \overline{\mathrm{x}}}{\mathcal{S}}_{\mathrm{o}}^{\overline{\mathrm{x}}\overline{\mathrm{y}}}{\tilde{s}}_i^{\overline{\mathrm{y}}}+\sum _{\mathrm{y}}{\mathcal{M}}_{\mathrm{o}}^{\overline{\mathrm{x}}\mathrm{y}}{\tilde{n}}_i^{\mathrm{y}}\right]{\tilde{s}}_{\overline{\mathrm{x}}}^i,\hfill \end{array} \end{gathered}$$

(147)

where the “$\mathrm{o}$” subscript means the quantity is evaluated for the ratio ${\tilde{u}}_i^{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^i/c^2\to 0$. Because of effective mass effects, the combination $\mathcal{U}/c^2$ as it appears in, say, ${\mathcal{B}}_{\mathrm{o}}^{\mathrm{x}\mathrm{y}}$ is not necessarily small.

The limiting form of the generalized pressure ${\Psi }_{\mathrm{o}}$ [cf.Eq. (61)] is

$${\Psi }_{\mathrm{o}}=-{\mathcal{U}}_{\mathrm{o}}+\sum _{\mathrm{x}}\left({\left.{\tilde{\mu }}_{\mathrm{x}}\right|}_{\mathrm{o}}-m_{\mathrm{x}}{c}^2\right){\tilde{n}}_{\mathrm{x}}+\sum _{\overline{\mathrm{x}}}{\left.{\tilde{T}}_{\overline{\mathrm{x}}}\right|}_{\mathrm{o}}{\tilde{s}}_{\overline{\mathrm{x}}},$$

(148)

and the $3+1$ total energy density E, momentum $P^a$, and stress $S_{ab}$ tend towards the values

$$\begin{array}{cc}\hfill E_{\mathrm{o}}& \to \sum _{\mathrm{x}}{m}_{\mathrm{x}}{c}^2{\tilde{n}}_{\mathrm{x}},\hfill \end{array}$$

(149a)

$$\begin{array}{cc}\hfill P_{\mathrm{o}}^i& \to \sum _{\mathrm{x}}{\tilde{\mu }}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}\left({\tilde{u}}_{\mathrm{x}}^i/c\right)+\sum _{\overline{\mathrm{x}}}{\tilde{T}}_{\overline{\mathrm{x}}}{\tilde{s}}_{\overline{\mathrm{x}}}\left({\tilde{u}}_{\overline{\mathrm{x}}}^i/c\right)+\frac{1}{{\mu }_{\mathrm{o}}}\left({\widetilde{ϵ}}^{ijk}{\tilde{E}}_j{\tilde{B}}_k/c\right)\to 0,\hfill \end{array}$$

(149b)

$$\begin{gathered} \begin{array}{cc}\hfill S_{\mathrm{o}}^{ij}& \to \sum _{\mathrm{x}}\left({\tilde{m}}_{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^j+\sum _{\mathrm{y}\ne \mathrm{x}}{\mathcal{B}}_{\mathrm{o}}^{\mathrm{x}\mathrm{y}}{\tilde{n}}_{\mathrm{y}}{\tilde{w}}_{\mathrm{y}\mathrm{x}}^j+\sum _{\overline{\mathrm{y}}}{\mathcal{M}}_{\mathrm{o}}^{\mathrm{x}\overline{\mathrm{y}}}{\tilde{s}}_{\overline{\mathrm{y}}}{\tilde{w}}_{\overline{\mathrm{y}}\mathrm{x}}^j\right){\tilde{n}}_{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^i\hfill \end{array} \\ \begin{array}{cc}\hfill & +\sum _{\overline{\mathrm{x}}}\left({\tilde{m}}_{\overline{\mathrm{x}}}{\tilde{u}}_{\overline{\mathrm{x}}}^j+\sum _{\overline{\mathrm{y}}\ne \overline{\mathrm{x}}}{\mathcal{S}}_{\mathrm{o}}^{\overline{\mathrm{x}}\overline{\mathrm{y}}}{\tilde{s}}_{\overline{\mathrm{y}}}{\tilde{w}}_{\overline{\mathrm{y}}\overline{\mathrm{x}}}^j+\sum _y{\mathcal{M}}_{\mathrm{o}}^{\overline{\mathrm{x}}\mathrm{y}}{\tilde{n}}_{\mathrm{y}}{\tilde{w}}_{\mathrm{y}\overline{\mathrm{x}}}^j\right){\tilde{s}}_{\overline{\mathrm{x}}}{\tilde{u}}_{\overline{\mathrm{x}}}^i\hfill \end{array} \\ \begin{array}{cc}\hfill & +\left({\Psi }_{\mathrm{o}}+\frac{1}{2{\mu }_{\mathrm{o}}}{\tilde{B}}^2\right)h^{ij}-\frac{1}{{\mu }_{\mathrm{o}}}{\tilde{B}}^i{\tilde{B}}^j.\hfill \end{array} \end{gathered}$$

(149c)

We have assumed that the so-called “$\overrightarrow{E}\times \overrightarrow{B}$” drift velocity for plasmas, i.e., ${\overrightarrow{v}}_{dr}=\overrightarrow{E}\times \overrightarrow{B}/{\left|\overrightarrow{B}\right|}^2$ must be small with respect to c. This leads to the constraint that $|{\overrightarrow{v}}_{dr}|\sim |\overrightarrow{E}|/|\overrightarrow{B}|\ll c$. We have assumed also that $\{{\mathcal{U}}_{\mathrm{o}},{\tilde{B}}^i{\tilde{B}}^j/{\mu }_{\mathrm{o}}\}\ll {\tilde{m}}_{\mathrm{x}}{c}^2{n}_{\mathrm{x}}$.

6.3. The Field Equations

To obtain the limiting form of the Einstein equation, we first work out the leading-order of the right-hand-sides of Eqs. (108), (109), and (110):

$$\begin{array}{cc}\hfill \frac{16\pi G}{c^4}{E}_{\mathrm{o}}& \to \frac{16\pi G}{c^2}\sum _{\mathrm{x}}{m}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}},\hfill \end{array}$$

(150a)

$$\begin{array}{cc}\hfill \frac{8\pi G}{c^3}{P}_{\mathrm{o}}^i& \to 0,\hfill \end{array}$$

(150b)

$$\begin{array}{cc}\hfill \frac{8\pi G}{c^4}{S}_{\mathrm{o}}^{ij}& \to 0,\hfill \end{array}$$

(150c)

$$\begin{array}{cc}\hfill \frac{8\pi G}{c^4}{S}_{\mathrm{o}}& \to 0.\hfill \end{array}$$

(150d)

Here, a factor of $1/c^2$ combines with the velocity terms ${\tilde{u}}_{\mathrm{x}}^i{\tilde{u}}_{\mathrm{y}}^j$ to drive to zero the stress terms $S_{\mathrm{o}}^{ij}$ and $S_{\mathrm{o}}$; the same factor drives ${\Psi }_{\mathrm{o}}/c^2\to 0$. The limiting forms of the Einstein equation components are

$$\begin{array}{cc}\hfill -{\partial }_i{\partial }^i\left(\frac{1}{c}\delta N+\frac{3}{2}\delta h\right)& \approx \frac{16\pi G}{c^2}\sum _{\mathrm{x}}{m}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}},\hfill \end{array}$$

(151a)

$$\begin{array}{cc}\hfill \frac{1}{2}{\partial }_j{\partial }^j\delta N_i& \approx 0,\hfill \end{array}$$

(151b)

$$\begin{array}{cc}\hfill -{\partial }_i{\partial }_j\left(\frac{2}{c}\delta N+\delta h\right)-\frac{1}{2}{\partial }_k{\partial }^k\delta h_{ij}& \approx \frac{4\pi G}{c^2}\left(\sum _{\mathrm{x}}{m}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}\right)h_{ij}.\hfill \end{array}$$

(151c)

If we take the trace of Eq. (151c), we can solve for ${\partial }_i{\partial }^i\delta h$. Substituting this into Eq. (151a) gives

$${\partial }_i{\partial }^i\Phi \approx 4\pi G\sum _{\mathrm{x}}{m}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}},$$

(152)

where $c\delta N\equiv \Phi$ is the standard gravitational potential. As a check of this identification we note that the geodesic equation — $u_p^b{\nabla }_b{u}_p^a=0$, where $u_p^a$ is a point particle four-velocity — gives in this limit

$$\frac{{\mathrm{d}}^2{x}^i}{\mathrm{d}{t}^2}\approx -c^2{\Gamma }_{00}^i=-c^2{\partial }^i\left(\frac{1}{c}\delta N\right)=-{\partial }^i\Phi =a^i,$$

(153)

where the last equality follows from Eq. (136b).

Using again the trace Eq. (151c), but substituting it into Eq. (152), then we find (to consistent order in c)

$${\partial }_i{\partial }^i\delta h\approx 0,$$

(154)

which then implies

$${\partial }_k{\partial }^k\delta h_{ij}\approx 0.$$

(155)

In this Newtonian context, we assume our system has compact support and is such that an asymptotically flat infinity exists for which $\delta N_i\to 0$ and $\delta h_{ij}\to 0$. Given that they both satisfy the Laplace equation it is consistent to have $\delta N_i=0$ and $\delta h_{ij}=0$ everywhere.

With this, we can implement now the limit of the fluid/plasma equations. Taking into account the fact that ${\tilde{\mu }}_{\mathrm{x}}/c^2$ and ${\tilde{T}}_{\overline{\mathrm{x}}}/c^2$ can have non-zero terms in the limit ${\tilde{u}}_{\mathrm{x}}^i/c\to 0$, then the individual pieces of the fluid/plasma equations in Eqs. (111a) — (111f) and the projections of the final form of $R_a^{\mathcal{N}}$ given in Eq. (128) become

$$\begin{array}{cc}\hfill u^a{f}_a^{\mathrm{x}}& =-{\tilde{n}}_{\mathrm{x}}^i\left(\frac{\partial {\tilde{\mu }}_i^{\mathrm{x}}}{\partial t}+{\partial }_i{\tilde{\mu }}_{\mathrm{x}}\right)-{\tilde{m}}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^i{a}_i+{\tilde{j}}_{\mathrm{x}}^i{\tilde{E}}_i,\hfill \end{array}$$

(156a)

$$\begin{array}{cc}\hfill u^a{f}_a^{\overline{\mathrm{x}}}& =-{\tilde{s}}_{\overline{\mathrm{x}}}^i\left(\frac{\partial {\tilde{\Theta }}_i^{\overline{\mathrm{x}}}}{\partial t}+{\partial }_i{\tilde{T}}_{\overline{\mathrm{x}}}\right)-{\tilde{m}}_{\overline{\mathrm{x}}}\left({\tilde{s}}_{\overline{\mathrm{x}}}/k_B\right){\tilde{u}}_{\overline{\mathrm{x}}}^i{a}_i,\hfill \end{array}$$

(156b)

$$\begin{gathered} \begin{array}{cc}\hfill {\perp }_i^j{f}_j^{\mathrm{x}}& ={\tilde{n}}_{\mathrm{x}}\left(\frac{\partial }{\partial t}+{\tilde{u}}_{\mathrm{x}}^j{\partial }_j\right){\tilde{\mu }}_i^{\mathrm{x}}-{\tilde{n}}_{\mathrm{x}}^j{\partial }_i{\tilde{\mu }}_j^{\mathrm{x}}+{\tilde{n}}_{\mathrm{x}}{\partial }_i{\tilde{\mu }}_{\mathrm{x}}+{\tilde{m}}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}{a}_i\hfill \end{array} \\ \begin{array}{cc}\hfill & -\left({\tilde{\sigma }}_{\mathrm{x}}{\tilde{E}}_i+{\widetilde{ϵ}}_{ijk}{\tilde{j}}_{\mathrm{x}}^j{\tilde{B}}^k\right),\hfill \end{array} \end{gathered}$$

(156c)

$$\begin{array}{cc}\hfill {\perp }_i^j{f}_j^{\overline{\mathrm{x}}}& ={\tilde{s}}_{\overline{\mathrm{x}}}\left(\frac{\partial }{\partial t}+{\tilde{u}}_{\overline{\mathrm{x}}}^j{\partial }_j\right){\tilde{\Theta }}_i^{\overline{\mathrm{x}}}-{\tilde{s}}_{\overline{\mathrm{x}}}^j{\partial }_i{\tilde{\Theta }}_j^{\overline{\mathrm{x}}}+{\tilde{s}}_{\overline{\mathrm{x}}}{\partial }_i{\tilde{T}}_{\overline{\mathrm{x}}}+{\tilde{m}}_{\overline{\mathrm{x}}}\left({\tilde{s}}_{\overline{\mathrm{x}}}/k_B\right)a_i,\hfill \end{array}$$

(156d)

$$\begin{array}{cc}\hfill \left(-u^a{\mu }_a^{\mathrm{x}}\right){\Gamma }_{\mathrm{x}}& ={\tilde{\mu }}_{\mathrm{x}}\left(\frac{\partial {\tilde{n}}_{\mathrm{x}}}{\partial t}+{\partial }_i{\tilde{n}}_{\mathrm{x}}^i\right)+{\tilde{m}}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^i{a}_i,\hfill \end{array}$$

(156e)

$$\begin{array}{cc}\hfill \left(-u^a{\Theta }_a^{\overline{\mathrm{x}}}\right){\Gamma }_{\overline{\mathrm{x}}}& ={\tilde{T}}_{\overline{\mathrm{x}}}\left(\frac{\partial {\tilde{s}}_{\overline{\mathrm{x}}}}{\partial t}+{\partial }_i{\tilde{s}}_{\overline{\mathrm{x}}}^i\right)+{\tilde{m}}_{\overline{\mathrm{x}}}\left({\tilde{s}}_{\overline{\mathrm{x}}}/k_B\right){\tilde{u}}_{\overline{\mathrm{x}}}^i{a}_i,\hfill \end{array}$$

(156f)

$$\begin{array}{cc}\hfill \left({\perp }_i^j{\mu }_j^{\mathrm{x}}\right){\Gamma }_{\mathrm{x}}& =\left(\frac{\partial {\tilde{n}}_{\mathrm{x}}}{\partial t}+{\partial }_j{\tilde{n}}_{\mathrm{x}}^j\right){\tilde{\mu }}_i^{\mathrm{x}}+\frac{1}{c^2}{\tilde{n}}_{\mathrm{x}}^j{a}_j{\tilde{\mu }}_i^{\mathrm{x}},\hfill \end{array}$$

(156g)

$$\begin{array}{cc}\hfill \left({\perp }_i^j{\Theta }_j^{\overline{\mathrm{x}}}\right){\Gamma }_{\overline{\mathrm{x}}}& =\left(\frac{\partial {\tilde{s}}_{\overline{\mathrm{x}}}}{\partial t}+{\partial }_j{\tilde{s}}_{\overline{\mathrm{x}}}^j\right){\tilde{\Theta }}_i^{\overline{\mathrm{x}}}+\frac{1}{c^2}{\tilde{s}}_{\overline{\mathrm{x}}}^j{a}_j{\tilde{\Theta }}_i^{\overline{\mathrm{x}}},\hfill \end{array}$$

(156h)

$$\begin{array}{cc}\hfill u^a\left(R_a^{\mathcal{N}}-q_{\mathcal{N}}{\Gamma }_{\mathcal{N}}{A}_a\right)& =-2\left[\left(2{\tilde{r}}_{61}-\frac{{\tilde{\gamma }}_{\mathcal{P}}}{{\tilde{\gamma }}_{\mathcal{N}}}\frac{{\overline{\mu }}_3}{{\overline{\mu }}_1}{\tilde{r}}_{13}\right){\tilde{w}}_{\mathcal{P}\mathcal{N}}^i+{\tilde{r}}_{\mathcal{N}}{\tilde{u}}_{\mathcal{N}}^i\right]{\tilde{w}}_i^{\mathcal{P}\mathcal{N}}+q_{\mathcal{N}}\tilde{V}{\Gamma }_{\mathcal{N}},\hfill \end{array}$$

(156i)

$$\begin{array}{cc}\hfill {\perp }_i^j\left(R_j^{\mathcal{N}}-q_{\mathcal{N}}{\Gamma }_{\mathcal{N}}{A}_j\right)& ={\tilde{r}}_{\mathcal{N}}{\tilde{w}}_i^{\mathcal{P}\mathcal{N}}-q_{\mathcal{N}}{\Gamma }_{\mathcal{N}}{\tilde{A}}_i.\hfill \end{array}$$

(156j)

The Maxwell equations and the continuity equation take the expected form of

$$\begin{array}{cc}\hfill {\partial }_i{\tilde{E}}^i& =c^2{\mu }_{\mathrm{o}}{q}_{\mathcal{P}}\left({\tilde{n}}_{\mathcal{P}}-{\tilde{n}}_{\mathcal{N}}\right),\hfill \end{array}$$

(157a)

$$\begin{array}{cc}\hfill {\widetilde{ϵ}}^i{}_{jk}{\partial }^j{\tilde{B}}^k& ={\mu }_{\mathrm{o}}{q}_{\mathcal{P}}\left({\tilde{n}}_{\mathcal{P}}^i-{\tilde{n}}_{\mathcal{N}}^i\right)+\frac{\partial {\tilde{E}}^i}{\partial t},\hfill \end{array}$$

(157b)

$$\begin{array}{cc}\hfill {\partial }_i{\tilde{B}}^i& =0,\hfill \end{array}$$

(157c)

$$\begin{array}{cc}\hfill {\widetilde{ϵ}}_{ijk}{\partial }^j{\tilde{E}}^k& =-\frac{\partial {\tilde{B}}_i}{\partial t},\hfill \end{array}$$

(157d)

$$\begin{array}{cc}\hfill \sum _{\mathrm{x}}\left({\partial }_i{\tilde{j}}_{\mathrm{x}}^i+\frac{\partial {\tilde{\sigma }}_{\mathrm{x}}}{\partial t}\right)& =0.\hfill \end{array}$$

(157e)

6.4. The Final Fluid/Plasma Newtonian Equations

Now we will write the final set of field equations so that we can point to some differences with those of the extant literature (such as [39,40]). We have clearly recovered the Newton equation for gravity and the Maxwell equations. The last thing is to collect all the fluid/plasma pieces to write the final form of their equations. To get the spirit of their role, we will assume that the gravitational and electromagnetic terms are known.

In total, we have to determine the six components ${\tilde{u}}_{\overline{\eta }}^i={\tilde{u}}_{\overline{\mathcal{P}}}^i={\tilde{u}}_{\eta }^i={\tilde{u}}_{\mathcal{P}}^i$ and ${\tilde{u}}_{\overline{\mathcal{N}}}^i={\tilde{u}}_{\mathcal{N}}^i$, as well as the six scalars $\left\{{\tilde{n}}_{\mathrm{x}},{\tilde{s}}_{\overline{\mathrm{x}}}\right\}$. Once the components $\left\{{\tilde{u}}_{\mathcal{P}}^i,{\tilde{u}}_{\mathcal{N}}^i\right\}$ are known, then we can use the divergence formulas in Eqs. (125a)–(125f), taken in combination with Eqs. (156e) and (156f), to determine the six scalars. Likewise, we can use the non-relativistic limit of the Euler Eqs. (58) and (59) to determine $\left\{{\tilde{u}}_{\mathcal{P}}^i,{\tilde{u}}_{\mathcal{N}}^i\right\}$ if the six scalars are known.

Using the sum of the non-relativistic forms of Eqs. (58) and (59) as the first Euler equation and keeping the non-relativistic form of Eq. (59) as the second, we find

$$\begin{gathered} \begin{array}{cc}\hfill 0& =\sum _{\mathrm{x}}\left[{\tilde{n}}_{\mathrm{x}}^i\left(\frac{\partial {\tilde{\mu }}_i^{\mathrm{x}}}{\partial t}+{\partial }_i{\tilde{\mu }}_{\mathrm{x}}\right)+{\Gamma }_{\mathrm{x}}{\tilde{\mu }}_{\mathrm{x}}\right]+\sum _{\overline{\mathrm{x}}}\left[{\tilde{s}}_{\overline{\mathrm{x}}}^i\left(\frac{\partial {\tilde{\Theta }}_i^{\overline{\mathrm{x}}}}{\partial t}+{\partial }_i{\tilde{T}}_{\overline{\mathrm{x}}}\right)+{\Gamma }_{\overline{\mathrm{x}}}{\tilde{T}}_{\overline{\mathrm{x}}}\right]\hfill \end{array} \\ \begin{array}{cc}\hfill & +\left[\sum _{\mathrm{x}}{\tilde{m}}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}{\tilde{u}}_{\mathrm{x}}^i+\sum _{\overline{\mathrm{x}}}{\tilde{m}}_{\overline{\mathrm{x}}}\left({\tilde{s}}_{\overline{\mathrm{x}}}/k_B\right){\tilde{u}}_{\overline{\mathrm{x}}}^i\right]a_i+\sum _{\mathrm{x}}{\tilde{j}}_{\mathrm{x}}^i{\tilde{E}}_i,\hfill \end{array} \end{gathered}$$

(158a)

$$\begin{gathered} \begin{array}{cc}\hfill 0& ={\tilde{n}}_{\mathcal{N}}^i\left(\frac{\partial {\tilde{\mu }}_i^{\mathcal{N}}}{\partial t}+{\partial }_i{\tilde{\mu }}_{\mathcal{N}}\right)+\left({\tilde{\mu }}_{\mathcal{N}}-{\overline{\mu }}_{\mathcal{N}}+q_{\mathcal{N}}\tilde{V}\right){\Gamma }_{\mathcal{N}}+{\tilde{s}}_{\overline{\mathcal{N}}}^i\left(\frac{\partial {\tilde{\Theta }}_i^{\overline{\mathcal{N}}}}{\partial t}+{\partial }_i{\tilde{T}}_{\overline{\mathcal{N}}}\right)\hfill \end{array} \\ \begin{array}{cc}\hfill & +\left({\tilde{T}}_{\overline{\mathcal{N}}}-T_{\overline{\mathcal{N}}}\right){\Gamma }_{\overline{\mathcal{N}}}+\left[{\tilde{m}}_{\mathcal{N}}{\tilde{n}}_{\mathcal{N}}+{\tilde{m}}_{\overline{\mathcal{N}}}\left({\tilde{s}}_{\overline{\mathcal{N}}}/k_B\right)\right]{\tilde{u}}_{\mathcal{N}}^i{a}_i+{\tilde{j}}_{\mathcal{N}}^i{\tilde{E}}_i-2{\tilde{r}}_{\mathcal{N}}{\tilde{u}}_{\mathcal{N}}^i{\tilde{w}}_i^{\mathcal{P}\mathcal{N}},\hfill \end{array} \end{gathered}$$

(158b)

$$\begin{gathered} \begin{array}{cc}\hfill 0& =\sum _{\mathrm{x}}\left[{\tilde{n}}_{\mathrm{x}}\left(\frac{\partial }{\partial t}+{\tilde{u}}_{\mathrm{x}}^j{\partial }_j\right){\tilde{\mu }}_i^{\mathrm{x}}+{\Gamma }_{\mathrm{x}}{\tilde{\mu }}_i^{\mathrm{x}}\right]+\sum _{\overline{\mathrm{x}}}\left[{\tilde{s}}_{\overline{\mathrm{x}}}\left(\frac{\partial }{\partial t}+{\tilde{u}}_{\overline{\mathrm{x}}}^j{\partial }_j\right){\tilde{\Theta }}_i^{\overline{\mathrm{x}}}+{\Gamma }_{\overline{\mathrm{x}}}{\tilde{\Theta }}_i^{\overline{\mathrm{x}}}\right]\hfill \end{array} \\ \begin{array}{cc}\hfill & +\left[\sum _{\mathrm{x}}{\tilde{m}}_{\mathrm{x}}{\tilde{n}}_{\mathrm{x}}+\sum _{\overline{\mathrm{x}}}{\tilde{m}}_{\overline{\mathrm{x}}}\left({\tilde{s}}_{\overline{\mathrm{x}}}/k_B\right)\right]a_i+{\partial }_i\Psi -\sum _{\mathrm{x}}\left({\tilde{\sigma }}_{\mathrm{x}}{\tilde{E}}_i+{\widetilde{ϵ}}_{ijk}{\tilde{j}}_{\mathrm{x}}^j{\tilde{B}}^k\right),\hfill \end{array} \end{gathered}$$

(158c)

$$\begin{gathered} \begin{array}{cc}\hfill 0& ={\tilde{n}}_{\mathcal{N}}\left(\frac{\partial }{\partial t}+{\tilde{u}}_{\mathcal{N}}^j{\partial }_j\right){\tilde{\mu }}_i^{\mathcal{N}}+{\Gamma }_{\mathcal{N}}{\tilde{\mu }}_i^{\mathcal{N}}+{\tilde{n}}_{\mathcal{N}}{\partial }_i{\tilde{\mu }}_{\mathcal{N}}-{\tilde{n}}_{\mathcal{N}}^j{\partial }_i{\tilde{\mu }}_j^{\mathcal{N}}\hfill \end{array} \\ \begin{array}{cc}\hfill & +{\tilde{s}}_{\overline{\mathcal{N}}}\left(\frac{\partial }{\partial t}+{\tilde{u}}_{\mathcal{N}}^j{\partial }_j\right){\tilde{\Theta }}_i^{\overline{\mathcal{N}}}+{\Gamma }_{\overline{\mathcal{N}}}{\tilde{\Theta }}_i^{\overline{\mathcal{N}}}+{\tilde{s}}_{\overline{\mathcal{N}}}{\partial }_i{\tilde{T}}_{\overline{\mathcal{N}}}-{\tilde{s}}_{\overline{\mathcal{N}}}^j{\partial }_i{\tilde{\Theta }}_j^{\overline{\mathcal{N}}}\hfill \end{array} \\ \begin{array}{cc}\hfill & +\left[{\tilde{m}}_{\mathcal{N}}{\tilde{n}}_{\mathcal{N}}+{\tilde{m}}_{\overline{\mathcal{N}}}\left({\tilde{s}}_{\overline{\mathcal{N}}}/k_B\right)\right]a_i-2{\tilde{r}}_{\mathcal{N}}{\tilde{w}}_i^{\mathcal{P}\mathcal{N}}+q_{\mathcal{N}}{\Gamma }_{\mathcal{N}}{\tilde{A}}_i-\left({\tilde{\sigma }}_{\mathcal{N}}{\tilde{E}}_i+{\widetilde{ϵ}}_{ijk}{\tilde{j}}_{\mathcal{N}}^j{\tilde{B}}^k\right),\hfill \end{array} \end{gathered}$$

(158d)

where we have used Eq. (57) to infer

$${\partial }_i\Psi =\sum _{\mathrm{x}}\left({\tilde{n}}_{\mathrm{x}}{\partial }_i{\tilde{\mu }}_{\mathrm{x}}-{\tilde{n}}_{\mathrm{x}}^j{\partial }_i{\tilde{\mu }}_j^{\mathrm{x}}\right)+\sum _{\overline{\mathrm{x}}}\left({\tilde{s}}_{\overline{\mathrm{x}}}{\partial }_i{\tilde{T}}_{\overline{\mathrm{x}}}-{\tilde{s}}_{\overline{\mathrm{x}}}^j{\partial }_i{\tilde{\Theta }}_j^{\overline{\mathrm{x}}}\right).$$

(159)

The obvious difference with the current literature is the impact of entrainment. We see that its effect of “tilting” the fluid momenta for the particles has survived the non-relativistic limit. Something else that survives is the entropy momentum. An unanticipated difference is the coupling of the particle ${\tilde{m}}_{\mathrm{x}}$ and thermal ${\tilde{m}}_{\overline{\mathrm{x}}}$ effective masses to gravity (via the acceleration $a_i$).

Tracing back, it is the presence of $n_{\mathrm{x}\mathrm{y}}^2$ in $\Lambda$ that leads to ${\tilde{m}}_{\mathrm{x}}$ and ${\tilde{m}}_{\overline{\mathrm{x}}}$ in the first place. Given the approach taken here, there is no a priori, generic principle for why the entrainment pieces in the gravitational couplings should be negligible; obviously they survive the $c\to \infty$ limit. In the absence of a generic principle for why it should be, say, $m_{\mathrm{x}}$ and not ${\tilde{m}}_{\mathrm{x}}$ that couples to gravity one must rely on the microscopic details of the particular system to be modelled. The difference between ${\tilde{m}}_{\mathrm{x}}$ and $m_{\mathrm{x}}$ can be assessed and then compared with the “smallness” of other approximations in the model.7

7. Conclusions and Follow-On Work

We have presented an action principle which yields, from start to finish, the field equations for a dissipative/resistive general relativistic two-fluid two-temperature plasma, with a neutrally charged component. The model is distinct from previous general relativistic formulations of the two-temperature plasma system (some of which are cited throughout the text), none of which rely on action principles, as far as we know.

Due to the very nature of action principles, the couplings between the fields are self-consistently incorporated into the full suite of field equations. For example, $T_{ab}$ follows automatically from the fields and couplings built into the total action, and its covariant divergence ${\nabla }_b{T}_a^b$ vanishes identically when the field equations are satisfied; i.e. ${\nabla }_b{T}_a^b=0$ is not itself an equation of motion, but rather an identity (as it should be because of diffeomorphism invariance). Along these same lines, we have shown how the formal inclusion of terms like $n_{\mathrm{x}\mathrm{y}}^2$ in the fluid action naturally leads to entrainment between different fluids and effective masses for particles and entropy. We have also seen that electromagnetic gauge issues are automatically accounted for by the internal consistency of the overall formalism.

Because of the fact that systems containing plasmas occur across many independent branches of physics, we made an effort to provide a, more or less, self-contained presentation. This is especially true for the $3+1$ decomposition discussion, which includes steps that are textbook material. However, while these steps are well-known in the general relativity community, they may well be new to other readers. Moreover, one of our main goals was to derive the Newtonian limit in a self-consistent way. This way we recovered field equations very much like those in the extant literature, but we also saw a new element emerge: the effective mass of entropy.

By developing the framework from the fibration picture into $3+1$ language, a step was taken towards a practical implementation of a two-temperature plasma within a general relativistic numerical simulation, as needed for neutron-star merger. There are, however, many further steps that are required. As noted in [32], as soon as an entrained multifluid system is constructed from this action approach, not all the equations of motion can be written in a conservation law form. Standard approaches for numerically evolving solutions with discontinuities, particularly the shocks forming during mergers, then do not apply. Instead, path-conservative methods are required (see, e.g., [41] for a brief review). However, these methods require a deeper understanding of the correct form of the dissipative terms appropriate to the model. Whilst the form of these terms can be deduced from the action framework, as detailed in [17], we have not taken those steps here. Furher work in this direction is required.

Moving forward there are several things that should be done: The first step would be to analyze local waves and modes of oscillation, to get a basic understanding of the stability/instability properties of the system. This would provide some insight on when the temperature difference is driven to zero or forced to diverge. Another step would be to allow for the additional terms in the fluid action that lead to bulk and shear viscosity, so as to tackle the numerical evolution issue raised above. Finally, a post-Newtonian expansion of the field equations will further unravel the role of (particle and entropy) effective masses and their coupling to the gravitational field. This may shed further light on the relevance of the entropy entrainment.