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.

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] 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

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}}$.

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$: where our convention for transforming between the two is

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

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}}$: such that the Einstein convention applies to repeated matter space indices, and

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 as well as the symmetric matter space “metric” $g_{\mathrm{x}}^{AB}$ to the spacetime metric $g_{ab}$, via

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]

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

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

It is also straightforward to confirm that

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 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

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,

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

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.

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 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

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}$.