format: text

text:
# Conservative formulations of general relativistic kinetic theory

## I Introduction

The state of the art in core-collapse supernova simulations now includes energy- and angle-dependent neutrino transport Rampp and Janka (2000); Mezzacappa et al. (2001); Liebendörfer et al. (2001); Rampp and Janka (2002); Liebendörfer et al. (2002); Thompson et al. (2002); Janka et al. (2002a, b). However, experience shows that simultaneous conservation of both energy and lepton number is difficult numerically Mezzacappa and Bruenn (1993a); Liebendörfer et al. (2002). This challenge motivates us to develop new conservative formulations of relativistic kinetic theory that are specifically attuned to the need for accurate accounting of particle number and energy in numerical simulations of radiative transfer problems.

Before describing the conservative formulations of kinetic theory we seek, we explain why energy- and angle-dependent neutrino transport is necessary in supernova simulations and detail the magnitude of the challenge of energy conservation.

Sophisticated treatments of neutrino transport are necessary because the ultimate energy source of the supernova explosion—the gravitational potential energy of the stellar progenitor’s core—is eventually converted almost completely into neutrinos. Some of the gravitational potential energy is lost to escaping neutrinos during core collapse, but most of it is converted into a thermal bath of dense nuclear matter, photons, electron/positron pairs, and trapped neutrinos deep inside the nascent neutron star. Neutrinos, having the weakest interactions, are the most efficient means of cooling; they diffuse outward on a time scale of seconds towards a semi-transparent region near the surface of the newly forming neutron star, and eventually escape with about 99% of the released gravitational energy. In modeling the conversion of gravitational potential energy into neutrino fluxes, energy- and angle-dependent neutrino transport is necessary to accurately follow the transition from quasi-isotropic diffusion to forward-peaked free streaming. In this transition region energy is transferred from the neutrino radiation to the matter behind a stalled shock wave, and this energy transfer may be necessary to propel the shock through the outer layers in an explosion Colgate and White (1966); Bethe and Wilson (1985). But whether or not such neutrino heating is the proximate cause of explosion, the fact that neutrinos dominate the energetics implies that accurate neutrino transport is integral to any realistic and comprehensive study of the explosion mechanism.

The importance of accurate neutrino transport is a lesson learned from experience with supernova simulations. Parametrized studies Janka and Mueller (1996) highlight the sensitivity of explosions to neutrino luminosities and to conditions in the semi-transparent region near the nascent neutron star’s surface (see also Ref. Janka (2001))—precisely the region where neutrino energy and angle dependence must be tracked carefully. Moreover, there remains a nagging qualititative uncertainty in simulations with multidimensional hydrodynamics: Those with neutrino transport that depends on direction in configuration space but is averaged over energy and angles exhibit explosions Herant et al. (1994); Burrows et al. (1995); Fryer and Warren (2002), while those with neutrino transport that depends on energy but is averaged over angles in both configuration and momentum space do not show explosions Mezzacappa et al. (1998a, b). It may be that these differing outcomes are due to the different neutrino transport schemes, both of which are ultimately inadequate.

Moving beyond these general arguments for the necessity of accurate neutrino transport, quantitative consideration of the energetics shows how severe the requirements are on one aspect of accuracy: energy conservation. As mentioned above, virtually all of the gravitational potential energy ($`10^{53}\mathrm{erg}`$) released during collapse is eventually converted into intense neutrino fluxes lasting several seconds. However, supernova explosion energies (the kinetic energy of the ejecta) are observed to be only $`10^{51}\mathrm{erg}`$. Now, because it is difficult to argue with any rigor about the physical impact of any energy lost or gained due to numerical error, the total energy should be conserved to a precision corresponding to the phenomena of interest in the problem. Hence a simulation’s result for explosion energy accurate to, say, $`10\%`$ would require total energy conservation to an accuracy of about $`0.1\%`$. Allowing for systematic error accrual, the total energy would have to be conserved at a level of $`0.1\%/N`$ per time step, where $`N10^5`$ is the total number of time steps in the simulation.

Conservative formulations of kinetic theory would help to meet the numerical challenge of particle number and energy conservation in core-collapse supernova simulations. To give an idea of the kind of formulation of kinetic theory that we seek, we first review familiar descriptions of the dynamics of a fluid medium. The dynamics can be described in two different ways, which might be called elemental and conservative.

The elemental formulation expresses the evolution of the fluid in terms of equations of motion for the velocity and some independent set of quantities (e.g. temperature, densities of various species comprising the fluid, etc.) measured by an observer moving along with the fluid (“comoving observer”). For example, consider a spacetime with metric components $`\{g_{\mu \nu }\}`$ and metric determinant $`gdet(g_{\mu \nu })`$, containing a perfect fluid with 4-velocity components $`\{u^\mu \}`$ and comoving frame total energy density $`\rho `$, pressure $`p`$, and baryon density $`n`$. Iin the absence of radiative transfer and significant energy input from nuclear reactions, the perfect fluid evolves according to

$`(\rho +p)u^\mu \left({\displaystyle \frac{u^i}{x^\mu }}+\mathrm{\Gamma }_{}^{i}{}_{\rho \mu }{}^{}u^\rho \right)+(g^{i\mu }+u^iu^\mu ){\displaystyle \frac{p}{x^\mu }}`$ $`=`$ $`0,`$ (1)
$`u^\mu {\displaystyle \frac{\rho }{x^\mu }}+{\displaystyle \frac{(\rho +p)}{\sqrt{g}}}{\displaystyle \frac{}{x^\mu }}\left(\sqrt{g}u^\mu \right)`$ $`=`$ $`0,`$ (2)
$`u^\mu {\displaystyle \frac{n}{x^\mu }}+{\displaystyle \frac{n}{\sqrt{g}}}{\displaystyle \frac{}{x^\mu }}\left(\sqrt{g}u^\mu \right)`$ $`=`$ $`0.`$ (3)
(Greek and latin letters are spacetime and space indices respectively.) Supplementary relations between $`\rho `$, $`p`$, and $`n`$—referred to as the equation of state—are determined by the microphysics of the fluid. The name elemental denotes the fact that by writing down separate equations of motion for the velocity and comoving-frame quantities, the kinetic and “intrinsic” fluid energies—two “elements” of the system—are analytically separated.

The conservative approach expresses the evolution of the system in terms of the divergence of the stress-energy tensor $`T^{\mu \nu }`$. For a perfect fluid, $`T^{\mu \nu }=(\rho +p)u^\mu u^\nu +pg^{\mu \nu }`$, and Eqs. (1) and (2) are replaced by

$$\frac{1}{\sqrt{g}}\frac{}{x^\mu }\left(\sqrt{g}T^{\nu \mu }\right)=\mathrm{\Gamma }_{}^{\nu }{}_{\rho \mu }{}^{}T^{\rho \mu },$$
(4)
while Eq. (3) is replaced by

$$\frac{1}{\sqrt{g}}\frac{}{x^\mu }\left(\sqrt{g}nu^\mu \right)=0.$$
(5)
Volume integrals of the left-hand sides of Eqs. (4) and (5)—obtained by multiplying by the invariant spacetime volume element $`\sqrt{g}d^4x`$ and integrating—are related to surface integrals in an obvious manner. Physically, this relates the time rates of change of 4-momentum and baryon number in a volume to fluxes through a surface surrounding that volume; hence the labeling of Eqs. (4) and (5) as conservative.

The right-hand side of Eq. (4) deserves note. After discussing the reasons for this term’s existence, we comment on what it means for conservation issues.

There are several important cases where the connection coefficients $`\mathrm{\Gamma }_{}^{\nu }{}_{\rho \mu }{}^{}`$ on the right-hand side of Eq. (4) might not vanish. They are present in curved spacetime, where (at least in part) they embody gravitational forces. But even in flat spacetime, coordinates employed by accelerated observers give rise to connection coefficients. And even without spacetime curvature or accelerated reference frames, connection coefficients arise from the use of curvilinear coordinates.

What does the right-hand side of Eq. (4) mean for 4-momentum conservation? Only when it vanishes—that is, only for inertial observers in flat spacetime employing rectilinear coordinates—are the components of total 4-momentum constant in time (“conserved”). For only in this case do the coordinates reflect the translation invariance of flat spacetime, the physical origin of 4-momentum conservation. (Curved spacetime lacks translation invariance, so there is no 4-momentum conservation.) Because the presence, for whatever reason, of a source term like the right-hand side of Eq. (4) means that the 4-momentum components in such a basis are not conserved, it might more properly be called a “balance equation” than a “conservation law”. But because the volume integral of the left-hand side easily translates into a surface integral, we still call it a “conservative” formulation.

There are some special cases in which a conserved quantity associated with the time coordinate $`t`$ can be found, however. For unaccelerated observers in flat spacetime, $`\mathrm{\Gamma }_{}^{t}{}_{\rho \mu }{}^{}=0`$, even in curvilinear coordinates. This means that energy is conserved, though the 3-momentum components in cuvilinear coordinates are not. Another special case is the restriction to spherical symmetry in general relativity: Here certain coordinate choices allow the non-vanishing $`\mathrm{\Gamma }_{}^{t}{}_{\rho \mu }{}^{}`$ terms to be absorbed into the left-hand aside, leading to the identification of a conserved energy-like quantity (see e.g. Ref. Liebendörfer et al. (2001)).

Having used the familiar example of a perfect fluid to discuss what we mean by “elemental” and “conservative” formulations, we now consider kinetic theory in terms of these categories. The evolution of a particle type described by kinetic theory is often expressed as an equation of motion for the distribution function $`f`$, the ensemble-averaged density of a given particle type in phase space. (These concepts will be defined with greater precision in subsequent sections; for the present discussion it is sufficient to understand that phase space is the combination of configuration space and momentum space, and that multiplying $`f`$ by the volume of an infinitesimal cell in phase space gives the number of particles having positions and momenta within the ranges defined by that cell.) The distribution function evolves due to advection through phase space and particle interactions.

Advection through phase space gives rise to derivatives of $`f`$ with respect to the components of the position vector $`x`$ and the momentum vector $`p`$. For numerical evolution it is convenient to parametrize the distribution function in terms of $`\{x^\mu \}`$, the components of $`x`$ in a global “coordinate basis” 111In a “coordinate basis” the basis vectors are $`\stackrel{}{e}_\mu =\frac{}{x^\mu }`$. These basis vectors “commute” because the order of partial derivatives can be interchanged freely. The connection coefficients in a coordinate basis are given by $`\mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}=\frac{1}{2}g^{\mu \sigma }\left(\frac{g_{\sigma \nu }}{x^\rho }+\frac{g_{\sigma \rho }}{x^\nu }\frac{g_{\nu \rho }}{x^\sigma }\right).`$ In a “non-coordinate basis” obtained from the coordinate basis by a transformation $`\stackrel{}{e}_\mu ^{}=L_{}^{\mu }{}_{\mu ^{}}{}^{}\stackrel{}{e}_\mu `$, the basis vectors $`\stackrel{}{e}_\mu ^{}=L_{}^{\mu }{}_{\mu ^{}}{}^{}\frac{}{x^\mu }`$ do not commute if $`L_{}^{\mu }{}_{\mu ^{}}{}^{}`$ depends on position in spacetime. The connection coefficients in a non-coordinate basis have extra terms associated with the non-vanishing commutators of the basis vectors. In this paper, the connection coefficients in a non-coordinate basis are obtained from those of a coordinate basis via the transformation $`\mathrm{\Gamma }_{}^{\mu ^{}}{}_{\nu ^{}\rho ^{}}{}^{}=L_{}^{\mu ^{}}{}_{\mu }{}^{}L_{}^{\nu }{}_{\nu ^{}}{}^{}L_{}^{\rho }{}_{\rho ^{}}{}^{}\mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}+L_{}^{\mu ^{}}{}_{\mu }{}^{}L_{}^{\rho }{}_{\rho ^{}}{}^{}\frac{L_{}^{\mu }{}_{\nu ^{}}{}^{}}{x^\rho }`$. . (Here and throughout this paper, quantities defined with respect to the coordinate basis have indices without accents.) For momentum we make a different choice of basis, because interactions with a fluid are most easily described—and best handled numerically—in terms of momentum components measured by a comoving observer. The change from the coordinate basis to an orthonormal basis associated with the comoving observer (a “non-coordinate basis”) has two parts. First there is a transformation $`dx^{\overline{\mu }}=e_{}^{\overline{\mu }}{}_{\mu }{}^{}dx^\mu `$ to an (in general non-comoving) orthonormal basis. (Here and throughout this paper, quantities defined with respect to the non-comoving orthonormal basis have indices accented with a bar.) This is followed by a Lorentz transformation $`dx^{\widehat{\mu }}=\mathrm{\Lambda }_{}^{\widehat{\mu }}{}_{\overline{\mu }}{}^{}dx^{\overline{\mu }}`$ to a comoving orthonormal basis. (Here and throughout this paper, quantities defined with respect to the comoving orthonormal basis have indices accented with a hat.) Hence advection through phase space will involve derivatives of $`f`$ with respect to the coordinate basis position components $`\{x^\mu \}`$ and the comoving orthonormal basis momentum components $`\{p^{\widehat{\mu }}\}`$.

Turning from advection to particle interactions, we consider the case where the particle species are sufficiently dilute that the interactions can be described in terms of a collision integral $`[f]`$ depending only on the distribution functions $`f`$ of the individual particle species. (This approximation ignores correlations between particles; that is, the number of instances of finding particles at the same position is obtained from the product of their distribution functions.) In this case, the equation of motion for the distribution function $`f`$ is Lindquist (1966); Mezzacappa and Matzner (1989)

$$p^{\widehat{\mu }}\left(\mathrm{\Lambda }_{}^{\overline{\mu }}{}_{\widehat{\mu }}{}^{}e_{}^{\mu }{}_{\overline{\mu }}{}^{}\frac{f}{x^\mu }\mathrm{\Gamma }_{}^{\widehat{\nu }}{}_{\widehat{\rho }\widehat{\mu }}{}^{}p^{\widehat{\rho }}\frac{f}{p^{\widehat{\nu }}}\right)=[f].$$
(6)
This is the Boltzmann equation.

In terms of the categories described above in connection with fluid evolution, the Boltzmann equation is “elemental”: The most fundamental quantity—the distribution function—is the evolved variable, and volume integrals of the equation in both configuration space and momentum space are not obviously related to surface integrals. This “non-conservative” character is present even in flat spacetime and rectangular coordinates. The first term of Eq. (6) is non-conservative even in flat spacetime and rectilinear coordinates because the boost $`\mathrm{\Lambda }_{}^{\overline{\mu }}{}_{\widehat{\mu }}{}^{}`$—which depends on the coordinates $`\{x^\mu \}`$—sits outside the derivative $`f/x^\mu `$. The spatial dependence of the boost also gives rise to non-vanishing $`\mathrm{\Gamma }_{}^{\widehat{\nu }}{}_{\widehat{\rho }\widehat{\mu }}{}^{}`$, even in flat spacetime and rectilinear coordinates; therefore the second term of Eq. (6) does not vanish. This second term is non-conservative because the factor $`p^{\widehat{\rho }}`$ sits outside the derivative $`f/p^{\widehat{\nu }}`$.

While the Boltzmann equation is “elemental” or “non-conservative”, it is well known (e.g., see Refs. Lindquist (1966); Ehlers (1971); Israel (1972)) the first two momentum moments of the distribution function (integrated over a suitable invariant momentum space volume element $`dP`$) constitute a particle number current $`N^\mu `$ and particle stress-energy tensor $`T^{\mu \nu }`$,

$`N^\mu `$ $`=`$ $`{\displaystyle fp^\mu 𝑑P}={\displaystyle fp^{\widehat{\mu }}\mathrm{\Lambda }_{}^{\overline{\mu }}{}_{\widehat{\mu }}{}^{}e_{}^{\mu }{}_{\overline{\mu }}{}^{}𝑑P},`$ (7)
$`T^{\mu \nu }`$ $`=`$ $`{\displaystyle fp^\mu p^\nu 𝑑P}={\displaystyle fp^{\widehat{\mu }}\mathrm{\Lambda }_{}^{\overline{\mu }}{}_{\widehat{\mu }}{}^{}e_{}^{\mu }{}_{\overline{\mu }}{}^{}p^{\widehat{\nu }}\mathrm{\Lambda }_{}^{\overline{\nu }}{}_{\widehat{\nu }}{}^{}e_{}^{\nu }{}_{\overline{\nu }}{}^{}𝑑P},`$ (8)
which (for electrically neutral particles) obey the balance equations

$`{\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{}{x^\mu }}\left(\sqrt{g}N^\mu \right)`$ $`=`$ $`{\displaystyle [f]𝑑P},`$ (9)
$`{\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{}{x^\mu }}\left(\sqrt{g}T^{\nu \mu }\right)`$ $`=`$ $`\mathrm{\Gamma }_{}^{\nu }{}_{\rho \mu }{}^{}T^{\rho \mu }+{\displaystyle [f]p^\nu 𝑑P}.`$ (10)
While this result is often stated in the literature (e.g., see Refs. Lindquist (1966); Ehlers (1971); Israel (1972)), because of the non-conservative character of Eq. (6) it is not obvious how its momentum moments give rise to Eqs. (9) and (10). Equation (6) contains factors $`\mathrm{\Lambda }_{}^{\overline{\mu }}{}_{\widehat{\mu }}{}^{}e_{}^{\mu }{}_{\overline{\mu }}{}^{}`$ outside $`f/x^\mu `$; but from Eqs. (7)-(10) we see that these factors have come inside the derivative with respect to $`x^\mu `$. What happens to the spacetime derivatives of $`\mathrm{\Lambda }_{}^{\overline{\mu }}{}_{\widehat{\mu }}{}^{}e_{}^{\mu }{}_{\overline{\mu }}{}^{}`$ that are generated in taking these factors inside the derivative? Furthermore, according to Eqs. (7) and (8), the quantities $`N^\mu `$ and $`T^{\mu \nu }`$ involve no momentum derivatives of $`f`$. But the second term of Eq. (6) contains a factor $`p^{\widehat{\rho }}(f/p^{\widehat{\nu }})`$. Because of the momentum factor outside the momentum derivative it is not obvious how this term will contribute to Eqs. (9) and (10) when integrated over momentum space.

How, then, is the connection between Eq. (6) and Eqs. (9) and (10) established in detail? The reviews of relativistic kinetic theory by Lindquist Lindquist (1966) and Israel Israel (1972) do not provide detailed proofs. As discussed in Sec. III, an explicit proof by Ehlers Ehlers (1971) relies on the fact that $`N^\mu `$ and $`T^{\mu \nu }`$ are momentum-integrated quantities, and no direct insight is gained into what happens to the momentum derivatives of $`f`$ in the integration over momentum space.

For those interested in computer models of radiative transfer problems, these are not idle academic questions; they are issues that must be faced in order to build simulations capable of making meaningful scientific statements about the core-collapse supernova explosion mechanism. Experience with supernova simulations in spherical symmetry shows that understanding the detailed connection between the Boltzmann equation and the particle number and 4-momentum balance equations has important consequences for how well these quantities are conserved in a simulation Mezzacappa and Bruenn (1993b); Liebendörfer et al. (2002). While a discretization of the Boltzmann equation is a natural numerical method of evolving the neutrino species, naïve differencings of the various terms in Eq. (6) generally will not be consistent with a straightfoward differencing of Eqs. (9) and (10), leading to unacceptable numerical errors in particle number and energy conservation. In Lagrangian (or “comoving”) coordinates in spherical symmetry, Mezzacappa and Bruenn Mezzacappa and Bruenn (1993b) derive a conservative formulation of the Boltzmann equation transparently related to particle number balance as expressed in Eq. (9). They also devise methods of handling momentum derivatives that are consistent with both number and energy conservation Mezzacappa and Bruenn (1993a). Liebendörfer et al. Liebendörfer et al. (2002) went a step further in this spherically symmetric case, deducing the connnection between the Mezzacappa and Bruenn “number conservative” Boltzmann equation in comoving coordinates and energy conservation as represented by an Eulerian (or “lab frame”) version of Eq. (10). Using complicated, non-intuitive differencings of hydrodynamic and gravitational variables, they construct a numerical implementation of the Mezzacappa and Bruenn “number conservative” Boltzmann equation that is stable and faithful to its analytic connection to the lab frame version of Eq. (10) to the accuracy necessary to make solid scientific statements about the neutrino-driven explosion mechanism in spherical symmetry Mezzacappa et al. (2001); Liebendörfer et al. (2001, 2002).

Because spherically symmetric models with Boltzmann transport fail to reproduce some important observable characteristics of core-collapse supernovae (e.g. the launch of an explosion Mezzacappa et al. (2001); Liebendörfer et al. (2001, 2002); Rampp and Janka (2000, 2002); Thompson et al. (2002)), this work must be followed up in multiple spatial dimensions (see Refs. Janka et al. (2002a, b) for some early efforts). In this paper we develop—allowing for full relativity and multiple spatial dimensions—conservative formulations of kinetic theory, such that volume integrals in both configuration and momentum space are trivially related to surface integrals. These conservative expressions make transparent the connection between Eq. (6) and Eqs. (9) and (10). They can be used to deduce the term-by-term cancellations involved in this connection, thereby illuminating the complicated differencings required to achieve the cancellations numerically. Our conservative formulations also suggest new primary variables of radiation transport: particle number and energy variables describing the contribution of each comoving orthonormal basis momentum bin to the coordinate basis number and energy densities. It may be that the use of these new radiation variables could provide a simpler path to an accurate accounting of particle number and energy in simulations of radiative transfer problems.

The organization of this paper is as follows. Differential forms and exterior calculus are natural mathematical tools for handling the volume elements and integrations needed in relativistic kinetic theory. We closely follow (and slightly extend) Ehlers’ work Ehlers (1971) in reviewing their use in the description of phase space for particles of definite mass (Sec. II) and the derivation of the Boltzmann equation (Sec. III). The centerpieces of this paper—two conservative reformulations of the Boltzmann equation, which provide transparent connections to particle number and energy-momentum balance as expressed in Eqs. (9) and (10)—are presented in Sec. IV. Because differential forms and exterior calculus may be unfamiliar to those whose primary interests are radiation transport, in general, or supernova science, in particular, Secs. II-IV each will contain two subsections: one containing a general derivation, and a second that explicitly demonstrates aspects of the derivation in the familiar case of the $`O(v)`$ limit in Lagrangian coordinates in spherical symmetry. (While we review some aspects of exterior calculus in our presentation, these are mostly in endnotes, and are more in the character of reminders than a self-contained introduction. For the latter we refer the reader to Refs. Ehlers (1971); Misner et al. (1973).) A conclusion (Sec. V) summarizes our results, comments on their connection to moment formalisms, and discusses the utility of these formulations for supernova simulations. As an application of our formalism, an appendix contains $`O(v)`$ equations in flat spacetime, but in coordinates sufficiently general to represent rectangular, spherical, and cylindrical coordinate systems.

## II Phase space for particles of definite mass

In this section we consider phase space for particles of definite mass. We specify particle trajectories—the “states” whose average occupation is specified by the distribution function $`f`$—and the volume elements needed to derive macroscopic observables from $`f`$ and operators that act upon it. In discussing the general case (Subsec. II.1) we follow Ehlers Ehlers (1971) closely, and refer the reader to his work for more detailed discussions and proofs of some of the assertions made here. We extend his discussion to the use of momentum components measured in a frame comoving with the fluid that interacts with the particle species treated by kinetic theory. The specific case of Lagrangian coordinates in spherical symmetry to $`O(v)`$ is treated in Subsec. II.2.

### II.1 Phase space for particles of definite mass: the general case

A study of kinetic theory begins with consideration of the trajectories of individual particles. The worldline $`x(\lambda )`$ of a particle of mass $`m`$ with 4-momentum $`p`$ and electric charge $`e`$, moving through a spacetime with metric components $`\{g_{\mu \nu }\}`$ and electromagnetic field tensor components $`\{F_{\mu \nu }\}`$, is determined by

$`{\displaystyle \frac{dx^\mu }{d\lambda }}`$ $`=`$ $`p^\mu ,`$ (11)
$`{\displaystyle \frac{dp^\mu }{d\lambda }}`$ $`=`$ $`eF_{}^{\mu }{}_{\nu }{}^{}p^\nu \mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}p^\nu p^\rho .`$ (12)
In the coordinate basis associated with the position vector components $`\{x^\mu \}`$, the connection coefficients are

$$\mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}=\frac{1}{2}g^{\mu \sigma }\left(\frac{g_{\sigma \nu }}{x^\rho }+\frac{g_{\sigma \rho }}{x^\nu }\frac{g_{\nu \rho }}{x^\sigma }\right).$$
(13)
As mentioned in Sec. I, we seek results expressed in terms of (orthonormal) momentum components $`\{p^{\widehat{\mu }}\}`$ measured by an observer comoving with the fluid with which the particles interact. Combining the transformation $`e_{}^{\overline{\mu }}{}_{\mu }{}^{}`$ to an orthonormal tetrad and boost $`\mathrm{\Lambda }_{}^{\widehat{\mu }}{}_{\overline{\mu }}{}^{}`$ to the comoving frame into the composite transformation

$$_{}^{\widehat{\mu }}{}_{\mu }{}^{}=\mathrm{\Lambda }_{}^{\widehat{\mu }}{}_{\overline{\mu }}{}^{}e_{}^{\overline{\mu }}{}_{\mu }{}^{},$$
(14)
Eqs. (11) and (12) are replaced by

$`{\displaystyle \frac{dx^\mu }{d\lambda }}`$ $`=`$ $`_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }},`$ (15)
$`{\displaystyle \frac{dp^{\widehat{\mu }}}{d\lambda }}`$ $`=`$ $`eF_{}^{\widehat{\mu }}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{\mu }}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}.`$ (16)
The connection coefficients do not transform exactly as tensor components, but as

$$\mathrm{\Gamma }_{}^{\widehat{\mu }}{}_{\widehat{\nu }\widehat{\rho }}{}^{}=_{}^{\widehat{\mu }}{}_{\mu }{}^{}_{}^{\nu }{}_{\widehat{\nu }}{}^{}_{}^{\rho }{}_{\widehat{\rho }}{}^{}\mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}+_{}^{\widehat{\mu }}{}_{\mu }{}^{}_{}^{\rho }{}_{\widehat{\rho }}{}^{}\frac{_{}^{\mu }{}_{\widehat{\nu }}{}^{}}{x^\rho }.$$
(17)
This completes the desired specification of particle trajectories.

The set of position vectors $`x`$ and momentum vectors $`p`$ comprise an 8-dimensional manifold $`M`$, the one-particle phase space for particles of arbitrary rest masses. The positions $`x`$ are points in spacetime, and the momenta $`p`$ are points in that portion of the tangent space at $`x`$ characterized by $`p^20`$ and $`p`$ future-directed. The curves $`(x(\lambda ),p(\lambda ))`$ obtained from the particle equations of motion comprise the phase flow in $`M`$. The phase flow is generated by the Liouville vector

$$L=p^{\widehat{\mu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}\frac{}{x^\mu }+\left(eF_{}^{\widehat{\mu }}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{\mu }}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right)\frac{}{p^{\widehat{\mu }}}.$$
(18)
Here $`(\left\{\frac{}{x^\mu }\right\},\left\{\frac{}{p^{\widehat{\mu }}}\right\})`$ have been chosen as basis vectors on $`M`$.

Specification of a particle mass $`m`$ defines a hypersurface of $`M`$, a 7-dimensional manifold $`M_m`$ called the one-particle phase space for particles of definite mass $`m`$. The mass satisifies

$$m^2=g_{\mu \nu }p^\mu p^\nu =g_{\mu \nu }_{}^{\mu }{}_{\widehat{\mu }}{}^{}_{}^{\nu }{}_{\widehat{\nu }}{}^{}p^{\widehat{\mu }}p^{\widehat{\nu }}.$$
(19)
Considered as a scalar function $`m(x,p)`$ on $`M`$, the particle mass satisfies

$$L[m]=0,$$
(20)
which expresses the constancy of $`m`$ on each phase space trajectory. Hence the 7-dimensional manifold $`M_m`$ is generated by all the phase space trajectories of particles of mass $`m`$. Equation (20) shows that $`L`$ is tangent to $`M_m`$, and Eq. (19) is a constraint indicating that it is one of the momentum dimensions of $`M`$ that has been lost in going over to $`M_m`$. Hence the Liouville operator restricted to $`M_m`$ can be expressed

$$L_m=p^{\widehat{\mu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}\frac{}{x^\mu }+\left(eF_{}^{\widehat{i}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{i}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right)\frac{}{p^{\widehat{i}}},$$
(21)
where $`(\left\{\frac{}{x^\mu }\right\},\left\{\frac{}{p^{\widehat{i}}}\right\})`$ have been chosen as basis vectors in $`M_m`$. To allow for a change of coordinates in momentum space from $`\left\{p^{\widehat{i}}\right\}`$ to $`\left\{u^{\widehat{i}}\right\}`$ (to momentum space spherical coordinates, for example), we rewrite Eq. (21) as

$$L_m=p^{\widehat{\mu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}\frac{}{x^\mu }+\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right)\frac{u^{\widehat{i}}}{p^{\widehat{j}}}\frac{}{u^{\widehat{i}}},$$
(22)
where $`(\left\{\frac{}{x^\mu }\right\},\left\{\frac{}{u^{\widehat{i}}}\right\})`$ are now the chosen basis vectors on $`M_m`$.

Before discussing volume and “surface” elements on $`M_m`$, we give a brief reminder of some properties of differential forms and their exterior derivatives (e.g., see Refs. Ehlers (1971); Misner et al. (1973)). If $`\{dz^a\}`$ are basis 1-forms on some manifold, an $`m`$-form $`\psi `$ can be expanded as

$$\psi =\frac{1}{m!}\psi _{a_1\mathrm{}a_m}dz^{a_1\mathrm{}a_m},$$
(23)
where $`dz^{a_1\mathrm{}a_m}`$ denotes the wedge product $`dz^{a_1}\mathrm{}dz^{a_m}`$. The quantities $`\psi _{a_1\mathrm{}a_m}`$ are the components of $`\psi `$. The components are completely antisymmetric ($`\psi _{a_1\mathrm{}a_m}=\psi _{[a_1\mathrm{}a_m]}`$), as is the wedge product, so that $`\psi `$ can also be expressed

$$\psi =\psi _{|a_1\mathrm{}a_m|}dz^{a_1\mathrm{}a_m},$$
(24)
where the vertical bars indicate that summation is taken only over $`a_1<a_2<\mathrm{}<a_m`$. The exterior derivative takes an $`m`$-form into an $`(m+1)`$-form as follows:

$$d\psi =\frac{\psi _{|a_1\mathrm{}a_m|}}{z^{a_{m+1}}}dz^{a_{m+1}a_1\mathrm{}a_m}.$$
(25)
Two properties of the exterior derivative will be used later in this paper:

$`d(\varphi +\psi )`$ $`=`$ $`d\varphi +d\psi ,`$ (26)
$`d(\varphi \psi )`$ $`=`$ $`d\varphi \psi +(1)^p\varphi d\psi ,`$ (27)
where $`\varphi `$ is a $`p`$-form.

Now we turn to a discussion of volume and “surface” elements on $`M_m`$. The invariant volume element in spacetime is the 4-form

$$\eta =\frac{1}{4!}\sqrt{g}ϵ_{\mu \nu \rho \sigma }dx^{\mu \nu \rho \sigma },$$
(28)
where $`g`$ is the determinant of the metric, $`ϵ_{\mu \nu \rho \sigma }`$ is the completely antisymmetric symbol with $`ϵ_{0123}=+1`$, and the wedge product of basis 1-forms is abbreviated by the notation $`dx^{\mu \nu \rho \sigma }dx^\mu dx^\nu dx^\rho dx^\sigma `$. Shifting attention to the momentum space at a given spacetime point, an invariant volume element on a mass shell corresponding to mass $`m`$ in the tangent space is

$$\pi _m=\frac{1}{3!}\frac{\sqrt{g}}{|p_0|}ϵ_{0ijk}dp^{ijk}.$$
(29)
Expressed in terms of basis 1-forms conjugate to orthonormal comoving frame momentum components, it becomes

$$\pi _m=\frac{1}{E(𝐩)}dp^{\widehat{1}\widehat{2}\widehat{3}},$$
(30)
where the bold character denotes a 3-vector, and $`E(𝐩)\sqrt{|𝐩|^2+m^2}`$. After a final change of momentum coordinates from $`𝐩`$ to $`𝐮`$, the momentum space volume element can also be expressed as

$$\pi _m=\frac{1}{E(𝐩)}\left|det\left[\frac{𝐩}{𝐮}\right]\right|du^{\widehat{1}\widehat{2}\widehat{3}}.$$
(31)
A volume element on $`M_m`$ is a 7-form, which can be expressed as

$`\mathrm{\Omega }`$ $`=`$ $`\eta \pi _m`$ (32)
$`=`$ $`{\displaystyle \frac{\sqrt{g}}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|dx^{0123}du^{\widehat{1}\widehat{2}\widehat{3}}.`$ (33)
The “surface element” normal to the phase flow in $`M_m`$ is an important quantity; it is a 6-form, obtained by contracting 222The contraction of a basis vector $`\frac{}{z^a}`$ with a 1-form $`dz^b`$ is $`\frac{}{z^a}dz^b=\delta _a^b`$. In the contraction of a vector with a direct product of 1-forms, the contraction is with the first 1-form in the direct product. This means that the contraction of a vector with a wedge product of $`n`$ 1-forms gives $`n`$ terms, because the wedge product is the completely antisymmetrized sum of $`n`$ direct products. $`\mathrm{\Omega }`$ with the Liouville vector $`L_m`$:

$`\omega `$ $`=`$ $`L_m\mathrm{\Omega }`$ (34)
$`=`$ $`{\displaystyle \frac{\sqrt{g}}{3!}}p^{\widehat{\mu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}ϵ_{\mu \nu \rho \sigma }(dx^{\nu \rho \sigma }\pi _m)+`$
$`{\displaystyle \frac{1}{2!E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{i}}}{p^{\widehat{j}}}}ϵ_{\widehat{0}\widehat{i}\widehat{k}\widehat{n}}\left(du^{\widehat{k}\widehat{n}}\eta \right).`$
The first term in Eq. (34) can be written in terms of surface elements in spacetime. A surface element in spacetime with normal in the $`x^\mu `$ direction is given by the contraction of the vector $`/x^\mu `$ with the spacetime volume element $`\eta `$:

$`\sigma _\mu `$ $`=`$ $`{\displaystyle \frac{}{x^\mu }}\eta `$ (35)
$`=`$ $`{\displaystyle \frac{\sqrt{g}}{3!}}ϵ_{\mu \nu \rho \sigma }dx^{\nu \rho \sigma }.`$
In terms of these spacetime surface elements, the first term in Eq. (34) can be expressed

$$\frac{\sqrt{g}}{3!}p^{\widehat{\mu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}ϵ_{\mu \nu \rho \sigma }(dx^{\nu \rho \sigma }\pi _m)=p^{\widehat{\mu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}(\sigma _\mu \pi _m).$$
(36)
This first term of the “surface element” $`\omega `$ is all that contributes to an integration over momentum space 333The computational procedure for integrals over differential forms involves (1) forming an infinitesimal parallelpiped—a wedge product of the vectors forming the edges of the infinitesimal region—at each point in the region over which the integration is performed; (2) contracting this infinitesimal parallelpiped with the differential form at each point; (3) calculating the integral with the usual rules of elementary calculus. An infinitesimal parallelpiped in momentum space is a wedge product of displacement vectors in each of the three directions in momentum space. In the contraction of this parallelpiped with $`\omega `$, the only terms that do not vanish are those in Eq. (36); the other terms of Eq. (34) have only two momentum space 1-forms.; it corresponds to a cell in the familiar 6-dimensional Newtonian phase space. The surface element $`\omega `$ also shares an important property with a Newtonian phase space volume element: It is invariant under phase flow. This result follows from the vanishing exterior derivative of $`\omega `$,

$$d\omega =0,$$
(37)
and the use of Stokes’ theorem on a small tube of worldlines Ehlers (1971).

### II.2 Phase space for particles of definite mass: Lagrangian coordinates in spherical symmetry to $`O(v)`$

Here we specialize certain results of the previous subsection to massless, electrically neutral particles propagating through flat spacetime, employing Lagrangian coordinates in spherical symmetry and retaining only terms to $`O(v)`$.

First we introduce these coordinates. In flat spacetime and Eulerian spherical coordinates, the metric components are specified by the line element

$$ds^2=d\stackrel{~}{t}^2+dr^2+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2.$$
(38)
We seek Lagrangian coordinates $`t`$ and $`m`$ to replace the Eulerian coordinates $`\stackrel{~}{t}`$ and $`r`$. The enclosed baryon mass $`m`$ is defined by

$$m=^r4\pi r^2\rho 𝑑r,$$
(39)
where $`\rho `$ is the baryon mass density. The mass density obeys the $`O(v)`$ conservation law

$$\frac{\rho }{\stackrel{~}{t}}+\frac{1}{r^2}\frac{}{r}\left(r^2\rho v\right)=0,$$
(40)
where

$$v=\frac{dr}{d\stackrel{~}{t}}$$
(41)
is the radial velocity measured by an Eulerian observer. We choose the Lagrangian time coordinate such that the metric components are specified by a line element of the form

$$ds^2=dt^2+g_{mm}dm^2+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2.$$
(42)
The transformations between the Eulerian and Lagrangian coordinates involve eight nontrivial quantities, the elements of the Jacobian matrices

$`J_{\{\stackrel{~}{t},r\}\{t,m\}}`$ $`=`$ $`\left(\begin{array}{cc}{\displaystyle \frac{\stackrel{~}{t}}{t}}\hfill & {\displaystyle \frac{\stackrel{~}{t}}{m}}\hfill \\ {\displaystyle \frac{r}{t}}\hfill & {\displaystyle \frac{r}{m}}\hfill \end{array}\right),`$ (43)
$`J_{\{t,m\}\{\stackrel{~}{t},r\}}`$ $`=`$ $`\left(\begin{array}{cc}{\displaystyle \frac{t}{\stackrel{~}{t}}}\hfill & {\displaystyle \frac{t}{r}}\hfill \\ {\displaystyle \frac{m}{\stackrel{~}{t}}}\hfill & {\displaystyle \frac{m}{r}}\hfill \end{array}\right).`$ (44)
Two of these elements, $`m/r`$ and $`m/\stackrel{~}{t}`$, are determined by Eqs. (39) and (40). The relationship between Eqs. (38) and (42)—specifically, the requirements $`g_{tt}=1`$ and $`g_{tm}=0`$—provide two more equations. The identity $`J_{\{\stackrel{~}{t},r\}\{t,m\}}=\left(J_{\{t,m\}\{\stackrel{~}{t},r\}}\right)^1`$ provides the remaining four equations, and we find that to $`O(v)`$

$`\left(\begin{array}{cc}{\displaystyle \frac{\stackrel{~}{t}}{t}}\hfill & {\displaystyle \frac{\stackrel{~}{t}}{m}}\hfill \\ {\displaystyle \frac{r}{t}}\hfill & {\displaystyle \frac{r}{m}}\hfill \end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}1& {\displaystyle \frac{v}{4\pi r^2\rho }}\hfill \\ v& {\displaystyle \frac{1}{4\pi r^2\rho }}\hfill \end{array}\right),`$ (45)
$`\left(\begin{array}{cc}{\displaystyle \frac{t}{\stackrel{~}{t}}}\hfill & {\displaystyle \frac{t}{r}}\hfill \\ {\displaystyle \frac{m}{\stackrel{~}{t}}}\hfill & {\displaystyle \frac{m}{r}}\hfill \end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}1& v\\ 4\pi r^2\rho v& 4\pi r^2\rho \end{array}\right).`$ (46)
From these relations we see that to $`O(v)`$, the line element in our Lagrangian coordinates is

$`ds^2`$ $`=`$ $`g_{\mu \nu }dx^\mu dx^\nu `$ (47)
$`=`$ $`dt^2+\left({\displaystyle \frac{1}{4\pi r^2\rho }}\right)^2dm^2+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2,`$
whose determinant is

$$g=\left(\frac{\mathrm{sin}\theta }{4\pi \rho }\right)^2.$$
(48)
These metric components can be used to obtain results valid to $`O(v)`$ Castor (1972).

Now we consider the equations of motion that determine the particle trajectories. The rates of change, along a particle worldline, of the coordinate basis position components $`\{x^\mu \}=\{t,m,\theta ,\varphi \}`$ and comoving orthornomal basis momentum components $`\{p^{\widehat{\mu }}\}=\{p^{\widehat{0}},p^{\widehat{1}},p^{\widehat{2}},p^{\widehat{3}}\}`$ are given by Eqs. (15) and (16):

$`{\displaystyle \frac{dx^\mu }{d\lambda }}`$ $`=`$ $`_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }},`$ (49)
$`{\displaystyle \frac{dp^{\widehat{\mu }}}{d\lambda }}`$ $`=`$ $`\mathrm{\Gamma }_{}^{\widehat{\mu }}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}.`$ (50)
We now consider the quantities $`_{}^{\mu }{}_{\widehat{\mu }}{}^{}`$ and $`\mathrm{\Gamma }_{}^{\widehat{\mu }}{}_{\widehat{\nu }\widehat{\rho }}{}^{}`$ appearing in these equations. As discussed in the previous subsection, $`_{}^{\widehat{\mu }}{}_{\mu }{}^{}=\mathrm{\Lambda }_{}^{\widehat{\mu }}{}_{\overline{\mu }}{}^{}e_{}^{\overline{\mu }}{}_{\mu }{}^{}`$ is a composite transformation from the coordinate basis (denoted by unadorned indices) to an orthonormal comoving basis (denoted by hatted indices). In the present case the coordinate basis is already comoving, so that the boost $`\mathrm{\Lambda }_{}^{\widehat{\mu }}{}_{\overline{\mu }}{}^{}`$ is simply the identity transformation. Under the remaining transformation $`e_{}^{\overline{\mu }}{}_{\mu }{}^{}`$ to an orthonormal basis, the metric of Eq. (47) must become the Lorentz metric. We make an obvious choice and take the non-zero transformation elements $`_{}^{\widehat{\mu }}{}_{\mu }{}^{}`$ to be

$`_{}^{\widehat{0}}{}_{t}{}^{}`$ $`=`$ $`1,`$ (51)
$`_{}^{\widehat{1}}{}_{m}{}^{}`$ $`=`$ $`{\displaystyle \frac{1}{4\pi r^2\rho }},`$ (52)
$`_{}^{\widehat{2}}{}_{\theta }{}^{}`$ $`=`$ $`r,`$ (53)
$`_{}^{\widehat{3}}{}_{\varphi }{}^{}`$ $`=`$ $`r\mathrm{sin}\theta .`$ (54)
The elements $`_{}^{\mu }{}_{\widehat{\mu }}{}^{}`$ of the inverse transformation are obvious. Turning to the connection coefficients, we employ Eqs. (13), (47), and (45) to find the non-zero coordinate basis connection coefficients $`\mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}`$ to be

$`\mathrm{\Gamma }_{}^{t}{}_{mm}{}^{}`$ $`=`$ $`\left({\displaystyle \frac{1}{4\pi r^2\rho }}\right)^2\left({\displaystyle \frac{2v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right),`$ (55)
$`\mathrm{\Gamma }_{}^{t}{}_{\theta \theta }{}^{}`$ $`=`$ $`rv,`$ (56)
$`\mathrm{\Gamma }_{}^{t}{}_{\varphi \varphi }{}^{}`$ $`=`$ $`rv\mathrm{sin}^2\theta ,`$ (57)
$`\mathrm{\Gamma }_{}^{m}{}_{tm}{}^{}=\mathrm{\Gamma }_{}^{m}{}_{mt}{}^{}`$ $`=`$ $`\left({\displaystyle \frac{2v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right),`$ (58)
$`\mathrm{\Gamma }_{}^{m}{}_{mm}{}^{}`$ $`=`$ $`\left({\displaystyle \frac{1}{2\pi r^3\rho }}+{\displaystyle \frac{\mathrm{ln}\rho }{m}}\right),`$ (59)
$`\mathrm{\Gamma }_{}^{m}{}_{\theta \theta }{}^{}`$ $`=`$ $`4\pi r^3\rho ,`$ (60)
$`\mathrm{\Gamma }_{}^{m}{}_{\varphi \varphi }{}^{}`$ $`=`$ $`4\pi r^3\rho \mathrm{sin}^2\theta ,`$ (61)
$`\mathrm{\Gamma }_{}^{\theta }{}_{t\theta }{}^{}=\mathrm{\Gamma }_{}^{\theta }{}_{\theta t}{}^{}=\mathrm{\Gamma }_{}^{\varphi }{}_{t\varphi }{}^{}=\mathrm{\Gamma }_{}^{\varphi }{}_{\varphi t}{}^{}`$ $`=`$ $`{\displaystyle \frac{v}{r}},`$ (62)
$`\mathrm{\Gamma }_{}^{\theta }{}_{m\theta }{}^{}=\mathrm{\Gamma }_{}^{\theta }{}_{\theta m}{}^{}=\mathrm{\Gamma }_{}^{\varphi }{}_{m\varphi }{}^{}=\mathrm{\Gamma }_{}^{\varphi }{}_{\varphi m}{}^{}`$ $`=`$ $`{\displaystyle \frac{1}{4\pi r^3\rho }},`$ (63)
$`\mathrm{\Gamma }_{}^{\theta }{}_{\varphi \varphi }{}^{}`$ $`=`$ $`\mathrm{sin}\theta \mathrm{cos}\theta ,`$ (64)
$`\mathrm{\Gamma }_{}^{\varphi }{}_{\theta \varphi }{}^{}=\mathrm{\Gamma }_{}^{\varphi }{}_{\varphi \theta }{}^{}`$ $`=`$ $`\mathrm{cot}\theta .`$ (65)
Now we can obtain the connection coefficients in the orthonormal comoving frame by employing Eqs. (17), (51)-(54) and (55)-(65). We find that the non-zero $`\mathrm{\Gamma }_{}^{\widehat{\mu }}{}_{\widehat{\nu }\widehat{\rho }}{}^{}`$ are

$`\mathrm{\Gamma }_{}^{\widehat{0}}{}_{\widehat{1}\widehat{1}}{}^{}=\mathrm{\Gamma }_{}^{\widehat{1}}{}_{\widehat{1}\widehat{0}}{}^{}`$ $`=`$ $`\left({\displaystyle \frac{2v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right),`$ (66)
$`\mathrm{\Gamma }_{}^{\widehat{0}}{}_{\widehat{2}\widehat{2}}{}^{}=\mathrm{\Gamma }_{}^{\widehat{0}}{}_{\widehat{3}\widehat{3}}{}^{}=\mathrm{\Gamma }_{}^{\widehat{2}}{}_{\widehat{2}\widehat{0}}{}^{}=\mathrm{\Gamma }_{}^{\widehat{3}}{}_{\widehat{3}\widehat{0}}{}^{}`$ $`=`$ $`{\displaystyle \frac{v}{r}},`$ (67)
$`\mathrm{\Gamma }_{}^{\widehat{1}}{}_{\widehat{2}\widehat{2}}{}^{}=\mathrm{\Gamma }_{}^{\widehat{1}}{}_{\widehat{3}\widehat{3}}{}^{}=\mathrm{\Gamma }_{}^{\widehat{2}}{}_{\widehat{2}\widehat{1}}{}^{}\mathrm{\Gamma }_{}^{\widehat{3}}{}_{\widehat{3}\widehat{1}}{}^{}`$ $`=`$ $`{\displaystyle \frac{1}{r}},`$ (68)
$`\mathrm{\Gamma }_{}^{\widehat{2}}{}_{\widehat{3}\widehat{3}}{}^{}=\mathrm{\Gamma }_{}^{\widehat{3}}{}_{\widehat{3}\widehat{2}}{}^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{cot}\theta }{r}}.`$ (69)
Note that unlike the $`\mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}`$, the $`\mathrm{\Gamma }_{}^{\widehat{\mu }}{}_{\widehat{\nu }\widehat{\rho }}{}^{}`$ are not symmetric in their lower indices $`\widehat{\nu }\widehat{\rho }`$ because the comoving orthonormal basis is not a coordinate basis. With $`_{}^{\mu }{}_{\widehat{\mu }}{}^{}`$ and $`\mathrm{\Gamma }_{}^{\widehat{\mu }}{}_{\widehat{\nu }\widehat{\rho }}{}^{}`$ now explicitly specified, our consideration of the particle equations of motion is complete.

The next topic to consider from the previous subsection is the Liouville vector. This vector generates the phase flow—the set of trajectories $`(x(\lambda ),p(\lambda ))`$ obtained from the particle equations of motion, Eqs. (49) and (50). These are trajectories through the 8-dimensional manifold comprised of position vectors $`x`$ and momentum vectors $`p`$. But a 7-dimensional submanifold, the phase space for particles of definite mass, is of more immediate interest: The physical fact that the particles have definite mass—as embodied in the “mass shell” constraint of Eq. (19)—implies that there are only three independent momentum variables. We choose these to be $`\{u^{\widehat{i}}\}=\{ϵ,\mu ,\phi \}`$, related to the comoving orthonormal basis momentum components $`\{p^{\widehat{i}}\}`$ by

$`p^{\widehat{1}}`$ $`=`$ $`ϵ\mu ,`$ (70)
$`p^{\widehat{2}}`$ $`=`$ $`ϵ\sqrt{1\mu ^2}\mathrm{cos}\phi ,`$ (71)
$`p^{\widehat{3}}`$ $`=`$ $`ϵ\sqrt{1\mu ^2}\mathrm{sin}\phi .`$ (72)
In the present case of massless particles, the mass shell condition gives

$$p^{\widehat{0}}=ϵ.$$
(73)
From these expressions one can compute the Jacobian $`(p^{\widehat{i}}/u^{\widehat{j}})`$ and its inverse $`(u^{\widehat{i}}/p^{\widehat{j}})`$; this inverse is

$$\left(\begin{array}{ccc}\frac{ϵ}{p^{\widehat{1}}}\hfill & \frac{ϵ}{p^{\widehat{2}}}\hfill & \frac{ϵ}{p^{\widehat{3}}}\hfill \\ \frac{\mu }{p^{\widehat{1}}}\hfill & \frac{\mu }{p^{\widehat{2}}}\hfill & \frac{\mu }{p^{\widehat{3}}}\hfill \\ \frac{\phi }{p^{\widehat{1}}}\hfill & \frac{\phi }{p^{\widehat{2}}}\hfill & \frac{\phi }{p^{\widehat{3}}}\hfill \end{array}\right)=\left(\begin{array}{ccc}\mu & \sqrt{1\mu ^2}\mathrm{cos}\phi & \sqrt{1\mu ^2}\mathrm{sin}\phi \\ \frac{1\mu ^2}{ϵ}\hfill & \frac{\mu \sqrt{1\mu ^2}\mathrm{cos}\phi }{ϵ}\hfill & \frac{\mu \sqrt{1\mu ^2}\mathrm{sin}\phi }{ϵ}\hfill \\ 0& \frac{\mathrm{sin}\phi }{ϵ\sqrt{1\mu ^2}}\hfill & \frac{\mathrm{cos}\phi }{ϵ\sqrt{1\mu ^2}}\hfill \end{array}\right).$$
(74)
We are now ready to write down the Liouville vector on the phase space for particles of definite mass, obtained with the help of Eqs. (22), (51)-(54), (66)-(69), and (70)-(74):

$`L_m`$ $`=`$ $`ϵ{\displaystyle \frac{}{t}}+4\pi r^2\rho ϵ\mu {\displaystyle \frac{}{m}}+ϵ^2\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{}{ϵ}}+`$ (75)
$`ϵ(1\mu ^2)\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]{\displaystyle \frac{}{\mu }}+\mathrm{},`$
where “…” represents $`/\theta `$, $`/\varphi `$, and $`/\phi `$ terms whose explicit form we will not need in this spherically symmetric case.

Finally we consider volume and “surface” elements in phase space. From Eqs. (28) and (48), the spacetime volume element is

$$\eta =\frac{\mathrm{sin}\theta }{4\pi \rho }dtdmd\theta d\varphi .$$
(76)
Contraction of $`\eta `$ with the coordinate basis vectors $`/t`$, $`/m`$, $`/\theta `$, and $`/\varphi `$ produces (see Eq. (35)) the spacetime surface elements

$`\sigma _t`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}\theta }{4\pi \rho }}dmd\theta d\varphi ,`$ (77)
$`\sigma _m`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}\theta }{4\pi \rho }}dtd\theta d\varphi ,`$ (78)
$`\sigma _\theta `$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}\theta }{4\pi \rho }}dtdmd\varphi ,`$ (79)
$`\sigma _\varphi `$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}\theta }{4\pi \rho }}dtdmd\theta .`$ (80)
Turning to momentum space, from Eqs. (31) and (70)-(73) the invariant volume element in momentum space is found to be

$$\pi _m=ϵdϵd\mu d\phi .$$
(81)
We form the volume element on phase space, using Eqs. (32), (76), and (81):

$`\mathrm{\Omega }`$ $`=`$ $`\eta \pi _m`$ (82)
$`=`$ $`{\displaystyle \frac{\mathrm{sin}\theta ϵ}{4\pi \rho }}dtdmd\theta d\varphi dϵd\mu d\phi .`$ (83)
With the help of Eqs. (75) and (83), we also construct the surface element in phase space that is normal to the phase flow, the 6-form $`\omega `$ of Eq. (34):

$`\omega `$ $`=`$ $`L_m\mathrm{\Omega }`$ (86)
$`=`$ $`{\displaystyle \frac{\mathrm{sin}\theta ϵ}{4\pi \rho }}[(L_m)^tdmd\theta d\varphi dϵd\mu d\phi (L_m)^mdtd\theta d\varphi dϵd\mu d\phi +`$
$`(L_m)^ϵdtd\theta d\varphi dmd\mu d\phi (L_m)^\mu dtdmd\theta d\varphi dϵd\phi +\mathrm{}]`$
$`=`$ $`{\displaystyle \frac{\mathrm{sin}\theta ϵ^2}{4\pi \rho }}dmd\theta d\varphi dϵd\mu d\phi r^2\mathrm{sin}\theta ϵ^2\mu dtd\theta d\varphi dϵd\mu d\phi +`$
$`{\displaystyle \frac{\mathrm{sin}\theta ϵ^3}{4\pi \rho }}\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]dtd\theta d\varphi dmd\mu d\phi `$
$`{\displaystyle \frac{\mathrm{sin}\theta ϵ^2(1\mu ^2)}{4\pi \rho }}\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]dtdmd\theta d\varphi dϵd\phi +\mathrm{},`$
where “…” represents terms arising from the $`/\theta `$, $`/\varphi `$, and $`/\phi `$ terms of the Liouville vector. The validity of Eq. (37) can be verified (see Eq. (25) for the definition of the exterior derivative):

$`d\omega `$ $`=`$ $`({\displaystyle \frac{}{t}}\left({\displaystyle \frac{\mathrm{sin}\theta ϵ^2}{4\pi \rho }}\right)+{\displaystyle \frac{}{m}}\left(r^2\mathrm{sin}\theta ϵ^2\mu \right)+`$ (87)
$`{\displaystyle \frac{}{ϵ}}\left\{{\displaystyle \frac{\mathrm{sin}\theta ϵ^3}{4\pi \rho }}\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]\right\}+`$
$`{\displaystyle \frac{}{\mu }}\{{\displaystyle \frac{\mathrm{sin}\theta ϵ^2(1\mu ^2)}{4\pi \rho }}[{\displaystyle \frac{1}{r}}+\mu ({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}})]+\mathrm{}\})\times `$
$`dtdmd\theta d\varphi dϵd\mu d\phi `$
$`=`$ $`0,`$ (88)
where the expression for $`r/m`$ is taken from Eq. (45). The terms explicitly displayed in Eq. (87) do in fact sum to zero by themselves, and it can be shown that those represented by “…” do as well.

## III The distribution function, Boltzmann equation, and balance equations

In this section—again following Ehlers Ehlers (1971) for the general case—we define the distribution function $`f`$, derive the Boltzmann equation, and present the equations describing the balance of particle number and 4-momentum. In the general case we reproduce Ehlers’ derivation of the momentum-integrated number balance equation, and explain why it does not yield a conservative reformulation of the Boltzmann equation. We also present the balance equation for 4-momentum without derivation. In the special case of Lagrangian coordinates in spherical symmetry to $`O(v)`$, we verify part of the general derivation of the Boltzmann equation by explicit calculation. We also display the Boltzmann equation and the momentum-integrated number and 4-momentum balance equations pertaining to this special case, including the transformation to the lab frame needed to obtain a “conserved” energy.

### III.1 The distribution function, Boltzmann equation, and balance equations: the general case

The distribution function $`f`$ for a particle of a given type represents the density of particles in phase space. The particle type’s 7-dimensional phase space for particles of definite mass, $`M_m`$, is filled with trajectories $`(x(\lambda ),p(\lambda ))`$, or “states”. As a collection of particles evolves, the number of particles in each state changes due to collisions. If one considers a 6-dimensional hypersurface $`\mathrm{\Sigma }`$ in $`M_m`$, the ensemble-averaged number $`N[\mathrm{\Sigma }]`$ of occupied states crossing $`\mathrm{\Sigma }`$ is

$$N[\mathrm{\Sigma }]=_\mathrm{\Sigma }f\omega ,$$
(89)
where the surface element $`\omega `$ is given by Eq. (34). (Hypersurfaces like $`\mathrm{\Sigma }`$ are oriented, and the particle trajectories have a direction associated with them as well. The crossing of a hypersurface by an occupied trajectory can give a positive or negative contribution to Eq. (89), depending on their relative orientation.) With this definition, the distribution function $`f(x,p)`$ is a scalar Ehlers (1971).

An equation governing the evolution of the distribution function is obtained by considering a closed 6-dimensional hypersurface $`D`$ bounding a region $`D`$ in $`M_m`$. The net number of occupied states emerging from $`D`$ is, from Eq. (89) and the generalized Stokes’ theorem,

$$N[D]=_Df\omega =_Dd(f\omega ).$$
(90)
(We note in passing—and discuss in more detail in Sec. IV—that this expression, which relates a volume integral to a surface integral, is key to obtaining conservative formulations of kinetic theory.) Using the “product rule” of Eq. (27), the vanishing exterior derivative of $`\omega `$ (Eq. (37)), and the definition of $`\omega `$ as the contraction of the Liouville vector $`L_m`$ with the volume element $`\mathrm{\Omega }`$ (Eq. (34)), we have

$$d(f\omega )=df\omega =df(L_m\mathrm{\Omega }).$$
(91)
Any scalar function $`f`$, vector field $`L`$, and $`n`$-form field $`\mathrm{\Omega }`$ on an $`n`$-dimensional manifold obey the identity Ehlers (1971)

$$df(L\mathrm{\Omega })=L[f]\mathrm{\Omega }.$$
(92)
Hence the net number of occupied states emerging from $`D`$ is

$$N[D]=_DL_m[f]\mathrm{\Omega }.$$
(93)
(Equations (90)-(93) amount to the Liouville theorem in a relativistic context.) The domain $`D`$ in phase space consists of a region $`H`$ in spacetime, together with regions $`K_x`$ in the momentum space (the mass shell in the tangent space) at each spacetime point $`x`$. Recall also that the volume element $`\mathrm{\Omega }`$ in phase space is the product of the spacetime and momentum space volume elements, $`\mathrm{\Omega }=\eta \pi _m`$. Hence the integral over $`D`$ in Eq. (93) can be expressed as the iterated integral

$$N[D]=_H\eta \left(_{K_x}L_m[f]\pi _m\right).$$
(94)
Because $`D`$ is a closed surface, the net number of particles in trajectories emerging from $`D`$ must equal the net number of collisions in $`D`$. If correlations between particles can be neglected, the spacetime density of collisions can be expressed in terms of a collision integral $`[f]`$ that depends only on one-particle distribution functions. Therefore, if $`_{K_x}[f]\pi _m`$ denotes the spacetime density of collisions, we have

$$N[D]=_H\eta \left(_{K_x}[f]\pi _m\right).$$
(95)
Because the regions $`H`$ and $`K_x`$ are arbitrary, comparison of Eqs. (94) and (95) shows that the evolution of $`f`$ is determined by

$$L_m[f]=[f],$$
(96)
or, using Eq. (22),

$$p^{\widehat{\mu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}\frac{f}{x^\mu }+\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right)\frac{u^{\widehat{i}}}{p^{\widehat{j}}}\frac{f}{u^{\widehat{i}}}=[f].$$
(97)
This is the Boltzmann equation.

Next we consider the particle number 4-current, electromagnetic 4-current, and stress-energy tensor, and present the balance equations they satisfy. In Eq. (89), specialize the hypersurface of $`M_m`$ to be $`\mathrm{\Sigma }=G\times P_m(x)`$, where $`G`$ is an infinitesimal spacelike hypersurface in spacetime at point $`x`$, and $`P_m(x)`$ is the entire momentum space at point $`x`$ 444Other 6-dimensional hypersurfaces can be formed, which contain a 4-dimensional region in spacetime and a 2-dimensional surface in momentum space or a timelike region in spacetime together with a 3-dimensional region in momentum space; such hypersurfaces are not of use in defining the particle number 4-current and stress-energy tensor.. As explained in the endnote to the discussion following Eq. (36), on such a hypersurface the only relevant term of $`\omega `$ is its first term in Eq. (34) (given also by Eq. (36)), so that Eq. (89) becomes

$$N[G\times P_m(x)]=_G\sigma _\mu \left(_{P_m(x)}f_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}\pi _m\right).$$
(98)
This is the number of particles whose worldlines cross $`G`$ (at $`x`$); hence the 4-vector

$$N^\mu (x)=_{P_m(x)}f_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}\pi _m=_{P_m(x)}fp^\mu \pi _m$$
(99)
is the particle 4-current density. The electromagnetic 4-current is obtained by multiplying by the electric charge: $`J^\mu =eN^\mu `$. The particle 4-current satisfies the balance equation

$$\frac{1}{\sqrt{g}}\frac{}{x^\mu }\left(\sqrt{g}N^\mu \right)=[f]\pi _m,$$
(100)
and similarly for the electromagnetic 4-current. (Of course the collisions will be such that the sum of the divergences of the electromagnetic 4-currents of all particle species will vanish.) The stress-energy tensor is

$$T^{\mu \nu }=_{P_m(x)}f_{}^{\mu }{}_{\widehat{\mu }}{}^{}_{}^{\nu }{}_{\widehat{\nu }}{}^{}p^{\widehat{\mu }}p^{\widehat{\nu }}\pi _m=_{P_m(x)}fp^\mu p^\nu \pi _m,$$
(101)
which obeys the balance equation

$$\frac{1}{\sqrt{g}}\frac{}{x^\mu }\left(\sqrt{g}T^{\nu \mu }\right)=F_{}^{\nu }{}_{\mu }{}^{}J^\mu \mathrm{\Gamma }_{}^{\nu }{}_{\rho \mu }{}^{}T^{\rho \mu }+[f]_{}^{\nu }{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\pi _m.$$
(102)
To conclude this subsection, we reproduce Ehlers’ derivation of Eq. (100) for number balance, and explain why this proof does not yield a conservative reformulation of the Boltzmann equation. The proof involves forming an integral over a suitable hypersurface in phase space, evaluating that integral in two different ways, and comparing the results.

First we specialize the integral in Eq. (89) to a specific hypersurface of integration. Consider an arbitrary region in spacetime, $`\widehat{D}`$, whose boundary is $`\widehat{D}`$. Form a hypersurface $`D`$ in the seven-dimensional phase space $`M_m`$, composed of the boundary region $`\widehat{D}`$ in spacetime together with the entire momentum space $`P_m(x)`$ at each point of $`\widehat{D}`$. Equation (89) becomes

$$N[D]=_Df\omega ,$$
(103)
the integral we will evaluate in two different ways.

On the one hand, as derived in Eqs. (90)-(93), the integral in Eq. (103) can be written

$$N[D]=_DL_m[f]\mathrm{\Omega }.$$
(104)
With the Boltzmann equation, Eq. (96),

$$N[D]=_D[f]\mathrm{\Omega }=_{\widehat{D}}\eta \left(_{P_m(x)}[f]\pi _m\right).$$
(105)
On the other hand, the integral in Eq. (103) can be evaluated directly. As a first step, integrate over momentum space at each point of $`\widehat{D}`$. As discussed in Subsec. II.2, in such integrals over 3-dimensional regions in momentum space only the first term of $`\omega `$ in Eq. (34) contributes (see in particular the endnote in the sentence following Eq. (36)). Using Eq. (36) to express the first term of $`\omega `$, Eq. (103) can be expressed as

$`N[D]`$ $`=`$ $`{\displaystyle _{\widehat{D}}}\sigma _\mu \left({\displaystyle _{P_m(x)}}_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}f\pi _m\right)`$ (106)
$`=`$ $`{\displaystyle _{\widehat{D}}}\sigma _\mu N^\mu (x),`$ (107)
where the definition (Eq. (99)) of the particle number current has been employed in the second step. Equation (107) is a closed surface integral in spacetime; by the divergence theorem (a special case of the generalized Stokes theorem), it can be converted into a spacetime volume integral:

$`N[D]`$ $`=`$ $`{\displaystyle _{\widehat{D}}}\eta N_{}^{\mu }{}_{;\mu }{}^{}`$ (108)
$`=`$ $`{\displaystyle _{\widehat{D}}}\eta {\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{}{x^\mu }}\left(\sqrt{g}N^\mu \right).`$ (109)
Because the spacetime region $`\widehat{D}`$ is arbitrary, Eq. (100) follows from comparison of Eqs. (105) and (109).

Now we can see why this type of proof does not yield a conservative reformulation of the Boltzmann equation. In particular, we can see why no direct insight is gained into the fate of the non-conservative momentum derivative terms in the Boltzmann equation (Eq. (97)) upon integration over momentum space. The left-hand side of Eq. (97) arises from the action of the Liouville vector on $`f`$. Similarly, the phase space surface element $`\omega `$ of Eq. (34) is given by the contraction of the Liouville vector with the phase space volume element $`\mathrm{\Omega }`$. Inspection shows that the second term of Eq. (34) is closely related to the momentum derivative terms in the Boltzmann equation. But in the above derivation of the number balance equation, an expression for $`N^\mu `$ is obtained at an early stage by integration over momentum space; this integration kills these terms of $`\omega `$, so that they are no longer present by the time the spacetime divergence theorem is applied. The power and elegance of exterior calculus has hidden the algebraic details of the connection between the Boltzmann equation and the number balance equation!

This is remedied in Section IV. Instead of choosing a specialized hypersurface involving all of momentum space and evaluating the surface integral of the left-hand side of Eq. (90), the path to conservative reformulations of the Boltzmann equations involves analysis of the volume integral on the right-hand side of Eq. (90); for here all the terms of $`\omega `$ have been brought inside the exterior derivative, including the ones related to the momentum derivative terms appearing in the Boltzmann equation.

### III.2 The distribution function, Boltzmann equation, and balance equations: Lagrangian coordinates in spherical symmetry to $`O(v)`$

The first result to specialize from the previous subsection is the Boltzmann equation. An aspect of the derivation that can be demonstrated by explicit calculation is the assertion

$$df\omega =L_m[f]\mathrm{\Omega }$$
(110)
contained in Eqs. (91)-(92). In spherical symmetry $`f=f(t,m,ϵ,\mu )`$, so that its gradient is

$$df=\frac{f}{t}dt+\frac{f}{m}dm+\frac{f}{ϵ}dϵ+\frac{f}{\mu }d\mu .$$
(111)
Forming the wedge product with $`\omega `$ of Eq. (86) results in

$`df\omega `$ $`=`$ $`\{{\displaystyle \frac{\mathrm{sin}\theta ϵ^2}{4\pi \rho }}{\displaystyle \frac{f}{t}}+r^2\mathrm{sin}\theta ϵ^2\mu {\displaystyle \frac{f}{m}}+{\displaystyle \frac{\mathrm{sin}\theta ϵ^3}{4\pi \rho }}[\mu ^2({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}){\displaystyle \frac{v}{r}}]{\displaystyle \frac{f}{ϵ}}+`$ (112)
$`{\displaystyle \frac{\mathrm{sin}\theta ϵ^2(1\mu ^2)}{4\pi \rho }}[{\displaystyle \frac{1}{r}}+\mu ({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}})]{\displaystyle \frac{f}{\mu }}\}\times `$
$`dtdmd\theta d\varphi dϵd\mu d\phi ,`$
where all other terms vanish because the the wedge product of a 1-form with itself vanishes due to antisymmetry. Comparison with Eqs. (75) and (83) for $`L_m`$ and $`\mathrm{\Omega }`$ then shows that $`df\omega =L_m[f]\mathrm{\Omega }`$, as was to be demonstrated. The Boltzmann equation, Eq. (97), specializes to

$`ϵ{\displaystyle \frac{f}{t}}+4\pi r^2\rho ϵ\mu {\displaystyle \frac{f}{m}}+ϵ^2\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{f}{ϵ}}+`$
$`ϵ(1\mu ^2)\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]{\displaystyle \frac{f}{\mu }}`$ $`=`$ $`[f].`$ (113)
This agrees with Eq. (20) of Ref. Castor (1972).

Next we specialize Eq. (100) for particle number balance. Equations (99), (51)-(54), (70)-(73), and (81), together with the fact that $`f=f(t,m,ϵ,\mu )`$ in spherical symmetry, imply that

$`N^t`$ $`=`$ $`{\displaystyle f\mathrm{\hspace{0.33em}2}\pi ϵ^2𝑑ϵ𝑑\mu }\rho H^N,`$ (114)
$`N^m`$ $`=`$ $`4\pi r^2\rho {\displaystyle f\mathrm{\hspace{0.33em}2}\pi ϵ^2\mu 𝑑\mu }4\pi r^2\rho ^2G^N`$ (115)
are the only nonvanishing components of the particle number vector. With these expressions and Eqs. (48) and (81), Eq. (100) becomes

$$\frac{H^N}{t}+\frac{}{m}\left(4\pi r^2\rho G^N\right)=\frac{1}{\rho }[f]ϵ𝑑ϵ𝑑\mu 𝑑\phi .$$
(116)
This agrees with Eq. (28) (with lapse function $`\alpha `$ set to 1 for flat spacetime) of Ref. Liebendörfer et al. (2001).

Finally we specialize Eq. (102) for particle number 4-momentum balance. Equations (99), (51)-(54), (70)-(73), and (81), together with the fact that $`f=f(t,m,ϵ,\mu )`$ in spherical symmetry, imply that

$`T^{tt}`$ $`=`$ $`{\displaystyle f\mathrm{\hspace{0.33em}2}\pi ϵ^3𝑑ϵ𝑑\mu }\rho H^E,`$ (117)
$`T^{tm}=T^{mt}`$ $`=`$ $`4\pi r^2\rho {\displaystyle f\mathrm{\hspace{0.33em}2}\pi ϵ^3\mu 𝑑\mu }4\pi r^2\rho ^2G^E,`$ (118)
$`T^{mm}`$ $`=`$ $`(4\pi r^2\rho )^2{\displaystyle f\mathrm{\hspace{0.33em}2}\pi ϵ^3\mu ^2𝑑\mu }(4\pi r^2\rho )^2\rho P^E,`$ (119)
$`T^{\theta \theta }`$ $`=`$ $`{\displaystyle \frac{1}{r^2}}{\displaystyle f\pi ϵ^3(1\mu ^2)𝑑\mu }{\displaystyle \frac{\rho }{2r^2}}(H^EP^E),`$ (120)
$`T^{\varphi \varphi }`$ $`=`$ $`{\displaystyle \frac{1}{r^2\mathrm{sin}^2\theta }}{\displaystyle f\pi ϵ^3(1\mu ^2)𝑑\mu }{\displaystyle \frac{\rho }{2r^2\mathrm{sin}^2\theta }}(H^EP^E)`$ (121)
are the only nonvanishing components of the particle stress-energy tensor. With these expressions and Eqs. (48), (55)-(65), (51)-(54), (70)-(73), and (81), the $`t`$ and $`m`$ components of Eq. (102) become

$`{\displaystyle \frac{H^E}{t}}+{\displaystyle \frac{}{m}}\left(4\pi r^2\rho G^E\right)+{\displaystyle \frac{v}{r}}(H^EP^E)`$ $``$ $`\left({\displaystyle \frac{2v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)P^E`$ (122)
$`=`$ $`{\displaystyle \frac{1}{\rho }}{\displaystyle [f]ϵ^2𝑑ϵ𝑑\mu 𝑑\phi },`$
$`{\displaystyle \frac{1}{r^2\rho }}{\displaystyle \frac{}{t}}(r^2\rho G^E)+{\displaystyle \frac{1}{4\pi r^2\rho }}{\displaystyle \frac{}{m}}\left((4\pi r^2\rho )^2P^E\right){\displaystyle \frac{1}{r}}(H^EP^E)`$ $``$ $`2\left({\displaystyle \frac{2v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)G^E`$
$`(4\pi r^2\rho )\left({\displaystyle \frac{1}{2\pi r^3\rho }}+{\displaystyle \frac{\mathrm{ln}\rho }{m}}\right)P^E`$ $`=`$ $`{\displaystyle \frac{1}{\rho }}{\displaystyle [f]ϵ^2\mu 𝑑ϵ𝑑\mu 𝑑\phi }.`$ (123)
With the help of Eq. (45), the $`m`$ component equation can be expressed

$`{\displaystyle \frac{G^E}{t}}+{\displaystyle \frac{}{m}}\left(4\pi r^2\rho P^E\right){\displaystyle \frac{1}{r}}(H^EP^E)`$ $``$ $`\left({\displaystyle \frac{2v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)G^E`$ (124)
$`=`$ $`{\displaystyle \frac{1}{\rho }}{\displaystyle [f]ϵ^2\mu 𝑑ϵ𝑑\mu 𝑑\phi }.`$
Equations (122) and (124) agree with the expressions following Eq. (31) of Ref. Liebendörfer et al. (2001). As expected from the discussion in Sec. I, these 4-momentum balance equations have non-conservative terms arising from the non-vanishing connection coefficients of Eqs. (55)-(65).

It was also noted in Sec. I that in Eulerian coordinates in flat spacetime, the connection coefficients $`\mathrm{\Gamma }_{}^{\stackrel{~}{t}}{}_{\mu \nu }{}^{}`$ vanish, so that a conserved energy can be defined. Equation (102) for particle 4-momentum balance is a vector equation; the $`\stackrel{~}{t}`$ component of a vector $`V`$ is obtained from the $`t`$ and $`m`$ components by the transformation

$$V^{\stackrel{~}{t}}=\frac{\stackrel{~}{t}}{t}V^t+\frac{t}{m}V^m,$$
(125)
where the transformation coefficients are given by Eq. (45). Equation (122) is actually $`\rho ^1`$ times the $`t`$ component of 4-momentum balance, and Eq. (124) is $`(4\pi r^2\rho ^2)^1`$ times the $`m`$ component. These factors must be restored before plugging into Eq. (125); the final result for the $`\stackrel{~}{t}`$ component of 4-momentum balance is

$`{\displaystyle \frac{}{t}}\left(H^E+vG^E\right)+{\displaystyle \frac{}{m}}\left[4\pi r^2\rho \left(G^E+vP^E\right)\right]`$ $``$ $`{\displaystyle \frac{1}{r^2\rho }}{\displaystyle \frac{}{t}}\left(r^2\rho v\right)G^E`$ (126)
$`=`$ $`{\displaystyle \frac{1}{\rho }}{\displaystyle [f](1+v\mu )ϵ^2𝑑ϵ𝑑\mu 𝑑\phi }.`$
(The $`O(v)`$ baryon conservation expression

$$\frac{\mathrm{ln}\rho }{t}+\frac{2v}{r}+4\pi r^2\rho \frac{v}{m}=0,$$
(127)
obtained from Eqs. (40) and (45)-(46), has been employed.) Equation (126) agrees with Eq. (32) of Ref. Liebendörfer et al. (2001). The non-conservative third term on the left-hand side is $`O(v^2)`$ and ought to be dropped at the level of approximation we are using. However, it is retained in Ref. Liebendörfer et al. (2001) for a practical reason. In supernova simulations velocities can exceed values for which the $`O(v)`$ approximation is strictly valid; and while use of the $`O(v)`$ formalism might be questioned on physical grounds, the $`O(v^2)`$ non-conservative term can still be used to check that the numerical implementation of the $`O(v)`$ Boltzmann equation is consistent with Eq. (126).

## IV Conservative formulations of particle kinetics

We are now in a position to present conservative formulations of kinetic theory. We seek expressions closely tied to Eq. (100) for particle number balance and Eq. (102) for particle 4-momentum balance. In addition to deriving these conservative formulations in the general case, we show the relationship between them and explain the care that must be exercised in finite differencing the number balance equation in order to make it consistent with the energy balance equation. We also specialize these results to Lagrangian coordinates in spherical symmetry to $`O(v)`$, and see that the “number conservative” Boltzmann equation—arrived at in the past by guesswork Mezzacappa and Bruenn (1993b)—emerges naturally from our formalism.

### IV.1 Conservative formulations of particle kinetics: the general case

For an expression related to particle number balance, the derivation of the Boltzmann equation in the previous section had us closer to the desired result than might first be realized. In rushing headlong towards an equation for $`f`$, the key relation is easily overlooked: It is Eq. (90), the result of the generalized Stokes’ theorem. The integrand on the right-hand side of Eq. (90), $`d(f\omega )`$, is conservative: Having been obtained from Stokes’ theorem (the generalization) of the divergence theorem, it has everything inside the exterior derivative. Being an exterior derivative, it is too abstract to be directly useful; but massaging it just enough to bring it into the form $`d(f\omega )=[f]\mathrm{\Omega }`$, where $`\mathrm{\Omega }`$ is the volume element in phase space, we shall see that $`[f]`$ is in fact a conservative differential operator in the familiar, elementary sense.

Now we take a detailed look at Eq. (90), in particular the exterior derivative

$`d(f\omega )`$ $`=`$ $`d\left(\sqrt{g}p^{\widehat{\mu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}fϵ_{\mu |\nu \rho \sigma |}dx^{\nu \rho \sigma }\right)\pi _m`$ (128)
$`\left(\sqrt{g}p^{\widehat{\mu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}f\sigma _\mu \right)d\pi _m+`$
$`d\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{i}}}{p^{\widehat{j}}}}fϵ_{\widehat{0}\widehat{i}\left|\widehat{k}\widehat{n}\right|}du^{\widehat{k}\widehat{n}}\right)\eta +`$
$`\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{i}}}{p^{\widehat{j}}}}fϵ_{\widehat{0}\widehat{i}\left|\widehat{k}\widehat{n}\right|}du^{\widehat{k}\widehat{n}}\right)d\eta ,`$
where we have used Eqs. (34), (35), (26), and (27) and employed the vertical bar notation introduced in Eq. (24). The exterior derivatives will be expressed in terms of the basis $`\{dx^\mu ,du^{\widehat{i}}\}`$ on $`M_m`$. First we note that in our chosen momentum coordinates, the second term in Eq. (128) vanishes because $`d\pi _m=0`$: This is because $`\pi _m`$ as expressed in Eq. (31) has no dependence on $`\{x^\mu \}`$; and while it depends on $`\left\{u^{\widehat{i}}\right\}`$, adding another momentum $`1`$-form to the wedge product $`du^{\widehat{1}\widehat{2}\widehat{3}}`$ would cause it to vanish due to its antisymmetry. Similarly, the fourth term vanishes because $`d\eta =0`$. To understand what happens to the first term, consider one of the terms in the sum over $`\mu `$:

$$d\left(\sqrt{g}p^{\widehat{\mu }}_{}^{0}{}_{\widehat{\mu }}{}^{}fϵ_{0|123|}dx^{123}\right)\pi _m=\frac{}{x^0}\left(\sqrt{g}p^{\widehat{\mu }}_{}^{0}{}_{\widehat{\mu }}{}^{}f\right)dx^{0123}\pi _m.$$
(129)
While derivatives with respect to $`\{x^i\}`$ and $`\left\{u^{\widehat{i}}\right\}`$ are nonvanishing, only the wedge product of $`dx^0`$ with $`dx^{123}`$ is non-vanishing, and the wedge product in $`\pi _m`$ does not admit another momentum $`1`$-form. The other terms in the sum over $`\mu `$ are similar. To understand the third term in Eq. (128), consider a particular term in the sum over $`\widehat{i}`$:

$`d\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{1}}}{p^{\widehat{j}}}}fϵ_{\widehat{0}\widehat{1}\left|\widehat{2}\widehat{3}\right|}du^{\widehat{2}\widehat{3}}\right)\eta =`$
$`{\displaystyle \frac{}{u^{\widehat{1}}}}\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{1}}}{p^{\widehat{j}}}}f\right)du^{\widehat{1}\widehat{2}\widehat{3}}\eta ,`$ (130)
and similarly for the other terms in the sum over $`\widehat{i}`$. With Eq. (33), all terms in Eq. (128) can be assembled to form the expression

$$d(f\omega )=[f]\mathrm{\Omega },$$
(131)
where

$`[f]`$ $``$ $`{\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{}{x^\mu }}\left(\sqrt{g}_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}f\right)+`$ (132)
$`E(𝐩)\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|^1{\displaystyle \frac{}{u^{\widehat{i}}}}\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{i}}}{p^{\widehat{j}}}}f\right).`$
In the derivation of the Boltzmann equation, Eq. (131) can be used to replace Eq. (93) with a similar expression in which $`L_m[f]`$ is replaced by $`[f]`$, with the result that

$$[f]=[f].$$
(133)
This is our conservative formulation of particle number kinetics.

Next we seek a formulation related to particle 4-momentum balance. With arbitrary $`v_\mu `$, we evaluate the exterior derivative $`d(v_\mu p^\mu f\omega )`$ in two separate ways. (To simplify formulae, we often write $`p^\mu `$ instead of $`_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}`$, but because of our chosen momentum coordinates on $`M_m`$ the latter expression must often be used in computational steps.) First, employing Eqs. (131) and (133),

$$d(v_\mu p^\mu f\omega )=d(v_\mu p^\mu )f\omega +v_\mu p^\mu [f]\mathrm{\Omega }.$$
(134)
Computation shows that

$`d(v_\mu p^\mu )f\omega `$ $`=`$ $`\left(p^\mu p^\nu f{\displaystyle \frac{v_\mu }{x^\nu }}+v_\mu fp^\nu p^{\widehat{\mu }}{\displaystyle \frac{_{}^{\mu }{}_{\widehat{\mu }}{}^{}}{x^\nu }}\right)\mathrm{\Omega }+`$ (135)
$`v_\mu \left(_{}^{\mu }{}_{\widehat{0}}{}^{}{\displaystyle \frac{p^{\widehat{0}}}{p^{\widehat{j}}}}+_{}^{\mu }{}_{\widehat{i}}{}^{}{\displaystyle \frac{p^{\widehat{i}}}{p^{\widehat{j}}}}\right)\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right)f\mathrm{\Omega }.`$
Because of the mass shell constraint (Eq. (19)), $`p^{\widehat{0}}`$ is considered a function of the $`\left\{p^{\widehat{i}}\right\}`$. The geodesic equation (equation (16)) can be used to show that

$$\frac{p^{\widehat{0}}}{p^{\widehat{j}}}\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right)=eF_{}^{\widehat{0}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{0}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}.$$
(136)
This, together with Eq. (17) for the transformation of the connection coefficients, allows Eq. (135) to be written

$$d(v_\mu p^\mu )f\omega =\left(p^\mu p^\nu f\frac{v_\mu }{x^\nu }+v_\mu (eF_{}^{\mu }{}_{\nu }{}^{}p^\nu \mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}p^\nu p^\rho )f\right)\mathrm{\Omega }.$$
(137)
Second, we evaluate $`d(v_\mu p^\mu f\omega )`$ directly. A computation similar to that leading to Eqs. (131)-(132) yields

$$d(v_\mu p^\mu f\omega )=\left(p^\mu p^\nu f\frac{v_\mu }{x^\nu }+v_\mu 𝕋^\mu [f]\right)\mathrm{\Omega },$$
(138)
where

$`𝕋^\mu [f]`$ $``$ $`{\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{}{x^\nu }}\left(\sqrt{g}_{}^{\mu }{}_{\widehat{\mu }}{}^{}_{}^{\nu }{}_{\widehat{\nu }}{}^{}p^{\widehat{\mu }}p^{\widehat{\nu }}f\right)+`$ (139)
$`E(𝐩)\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|^1{\displaystyle \frac{}{u^{\widehat{i}}}}\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{i}}}{p^{\widehat{j}}}}_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}f\right).`$
Recalling that $`v_\mu `$ is arbitrary, Eqs. (134) and (137)-(139) can be combined into

$$𝕋^\mu [f]=\left(eF_{}^{\mu }{}_{\nu }{}^{}_{}^{\nu }{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}\mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}_{}^{\nu }{}_{\widehat{\nu }}{}^{}_{}^{\rho }{}_{\widehat{\rho }}{}^{}p^{\widehat{\nu }}p^{\widehat{\rho }}\right)f+_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}[f],$$
(140)
a “conservative” formulation of particle 4-momentum kinetics.

To illuminate the practical problem of accurately accounting for both particle number and energy, we consider the relationship between Eq. (133) for particle number kinetics and the $`t`$ component of Eq. (140) for particle 4-momentum kinetics. (This is an issue when one follows the general approach of, for example, Refs. Mezzacappa et al. (2001); Liebendörfer et al. (2001, 2002), in which a particle number distribution equation is solved and a particle energy equation is used as a consistency check.) Specifically, we need to relate the spatial and momentum divergence terms of these equations to each other and identify terms that cancel.

First we relate the spatial and momentum divergence terms of Eq. (133) and the $`t`$ component of Eq. (140). Their spatial divergence terms,

$`_x[f]`$ $``$ $`{\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{}{x^\mu }}\left(\sqrt{g}_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}f\right),`$ (141)
$`\left(𝕋^t\right)_x[f]`$ $``$ $`{\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{}{x^\mu }}\left(\sqrt{g}_{}^{t}{}_{\widehat{\nu }}{}^{}_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\nu }}p^{\widehat{\mu }}f\right),`$ (142)
are related by

$$_{}^{t}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}_x[f]=\left(𝕋^t\right)_x[f]fp^{\widehat{\mu }}p^{\widehat{\nu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}\frac{}{x^\mu }\left(_{}^{t}{}_{\widehat{\nu }}{}^{}\right).$$
(143)
Their momentum divergence terms,

$`_p[f]`$ $``$ $`E(𝐩)\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|^1{\displaystyle \frac{}{u^{\widehat{i}}}}\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\mu }\widehat{\rho }}{}^{}p^{\widehat{\mu }}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{i}}}{p^{\widehat{j}}}}f\right),`$ (144)
$`\left(𝕋^t\right)_p[f]`$ $``$ $`E(𝐩)\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|^1{\displaystyle \frac{}{u^{\widehat{i}}}}\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\mu }\widehat{\rho }}{}^{}p^{\widehat{\mu }}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{i}}}{p^{\widehat{j}}}}_{}^{t}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}f\right),`$ (145)
are related by

$$_{}^{t}{}_{\widehat{\nu }}{}^{}p^{\widehat{\nu }}_p[f]=\left(𝕋^t\right)_p[f]f\left(eF_{}^{\widehat{j}}{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\mu }\widehat{\rho }}{}^{}p^{\widehat{\mu }}p^{\widehat{\rho }}\right)\frac{u^{\widehat{i}}}{p^{\widehat{j}}}_{}^{t}{}_{\widehat{\nu }}{}^{}\frac{p^{\widehat{\nu }}}{u^{\widehat{i}}}.$$
(146)
Comparison of Eqs. (143) and (146) with Eq. (133) and the $`t`$ component of Eq. (140) shows that the following equation must be valid:

$$E_S+E_M=eF_{}^{t}{}_{\mu }{}^{}_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}\mathrm{\Gamma }_{}^{t}{}_{\mu \rho }{}^{}_{}^{\mu }{}_{\widehat{\mu }}{}^{}_{}^{\rho }{}_{\widehat{\rho }}{}^{}p^{\widehat{\mu }}p^{\widehat{\rho }},$$
(147)
where

$$E_Sfp^{\widehat{\mu }}p^{\widehat{\nu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}\frac{}{x^\mu }\left(_{}^{t}{}_{\widehat{\nu }}{}^{}\right)$$
(148)
is the “extra” term of Eq. (143) relating the spatial divergence terms of the number and energy balance equations, and

$$E_Mf\left(eF_{}^{\widehat{j}}{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\mu }\widehat{\rho }}{}^{}p^{\widehat{\mu }}p^{\widehat{\rho }}\right)\frac{u^{\widehat{i}}}{p^{\widehat{j}}}_{}^{t}{}_{\widehat{\nu }}{}^{}\frac{p^{\widehat{\nu }}}{u^{\widehat{i}}}$$
(149)
is the “extra” term in Eq. (146) relating the momentum divergence terms of the number and energy balance equations. Noting that

$$\frac{u^{\widehat{i}}}{p^{\widehat{j}}}_{}^{t}{}_{\widehat{\nu }}{}^{}\frac{p^{\widehat{\nu }}}{u^{\widehat{i}}}=\frac{p^{\widehat{0}}}{p^{\widehat{j}}}_{}^{t}{}_{\widehat{0}}{}^{}+\delta _{\widehat{j}}^{\widehat{k}}_{}^{t}{}_{\widehat{k}}{}^{},$$
(150)
Eq. (136) can be used to rewrite Eq. (149) as

$$E_M=f\left(eF_{}^{\widehat{\nu }}{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}\mathrm{\Gamma }_{}^{\widehat{\nu }}{}_{\widehat{\mu }\widehat{\rho }}{}^{}p^{\widehat{\mu }}p^{\widehat{\rho }}\right)_{}^{t}{}_{\widehat{\nu }}{}^{}.$$
(151)
Use of Eq. (17) and two applications of the identity

$$_{}^{\mu }{}_{\widehat{\rho }}{}^{}\frac{}{x^\sigma }\left(_{}^{\widehat{\rho }}{}_{\nu }{}^{}\right)=_{}^{\widehat{\rho }}{}_{\nu }{}^{}\frac{}{x^\sigma }\left(_{}^{\mu }{}_{\widehat{\rho }}{}^{}\right)$$
(152)
lead finally to

$$E_M=f\left[eF_{}^{t}{}_{\mu }{}^{}_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}\mathrm{\Gamma }_{}^{t}{}_{\mu \rho }{}^{}_{}^{\mu }{}_{\widehat{\mu }}{}^{}_{}^{\rho }{}_{\widehat{\rho }}{}^{}p^{\widehat{\mu }}p^{\widehat{\rho }}p^{\widehat{\mu }}p^{\widehat{\nu }}_{}^{\mu }{}_{\widehat{\mu }}{}^{}\frac{}{x^\mu }\left(_{}^{t}{}_{\widehat{\nu }}{}^{}\right)\right].$$
(153)
Equations (148) and (153) make it obvious that Eq. (147) is satisfied.

In a computational approach to transport in which a particle number distribution equation is solved and a particle energy equation is used as a consistency check, care must be taken that Eq. (147) is satisfied numerically. In particular, a given finite difference representation of the spacetime divergence $`_x[f]`$ implies a corresponding finite difference representation of position-dependent quantities in $`E_S`$, and the position-dependent quantities in terms of $`E_M`$ that cancel with $`E_S`$ must have a matching finite difference representation. Similarly, a given finite difference representation of the momentum divergence $`_p[f]`$ implies a corresponding finite difference representation of momentum variables in $`E_M`$, and the momentum variables in terms of $`E_S`$ that cancel with terms of $`E_M`$ must have a matching finite difference representation.

### IV.2 Conservative formulations of particle kinetics: Lagrangian coordinates in spherical symmetry to $`O(v)`$

Using Eqs. (48), (51)-(54), (66)-(69), and (70)-(73), the “number conservative” formulation of Eq. (133) specializes to

$`{\displaystyle \frac{}{t}}\left({\displaystyle \frac{f}{\rho }}\right)+{\displaystyle \frac{}{m}}\left(4\pi r^2\rho \mu {\displaystyle \frac{f}{\rho }}\right)+{\displaystyle \frac{1}{ϵ^2}}{\displaystyle \frac{}{ϵ}}\left(ϵ^3\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{f}{\rho }}\right)+`$
$`{\displaystyle \frac{}{\mu }}\left((1\mu ^2)\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]{\displaystyle \frac{f}{\rho }}\right)`$ $`=`$ $`{\displaystyle \frac{1}{\rho ϵ}}[f],`$ (154)
which agrees with Eq. (23) of Ref. Mezzacappa and Bruenn (1993b). In Ref. Mezzacappa and Bruenn (1993b), the necessity of making $`f/\rho `$ the evolved variable in order to get a conservative particle number equation is left unexplained. In the present formalism, we see that the factor $`1/\rho `$ comes from the factor $`\sqrt{g}=\mathrm{sin}\theta /4\pi \rho `$ in Eq. (100). Multiplication of Eq. (154) by $`ϵ^2dϵd\mu d\phi `$ and integrating immediately yields Eq. (116), expressing particle number “conservation”.

Similarly, the $`t`$ and $`m`$ components of the “momentum conservative” formulation of Eq. (140) are

$`{\displaystyle \frac{}{t}}\left({\displaystyle \frac{ϵf}{\rho }}\right)`$ $`+`$ $`{\displaystyle \frac{}{m}}\left(4\pi r^2\rho ϵ\mu {\displaystyle \frac{f}{\rho }}\right)+{\displaystyle \frac{1}{ϵ^2}}{\displaystyle \frac{}{ϵ}}\left(ϵ^4\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{f}{\rho }}\right)+`$ (155)
$`{\displaystyle \frac{}{\mu }}\left((1\mu ^2)\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]{\displaystyle \frac{ϵf}{\rho }}\right)`$
$`=`$ $`\left({\displaystyle \frac{2v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)ϵ\mu ^2f{\displaystyle \frac{v}{r}}ϵ(1\mu ^2)f+{\displaystyle \frac{1}{\rho }}[f]`$
and

$`{\displaystyle \frac{1}{r^2\rho }}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{r^2\rho ϵ\mu f}{\rho }}\right)`$ $`+`$ $`{\displaystyle \frac{1}{4\pi r^2\rho }}{\displaystyle \frac{}{m}}\left((4\pi r^2\rho )^2ϵ\mu ^2{\displaystyle \frac{f}{\rho }}\right)+{\displaystyle \frac{1}{ϵ^2}}{\displaystyle \frac{}{ϵ}}\left(ϵ^4\mu \left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{f}{\rho }}\right)+`$ (156)
$`{\displaystyle \frac{}{\mu }}\left((1\mu ^2)\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]{\displaystyle \frac{ϵ\mu f}{\rho }}\right)`$
$`=`$ $`2\left({\displaystyle \frac{2v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)ϵ\mu f+(4\pi r^2\rho )\left({\displaystyle \frac{1}{2\pi r^3\rho }}+{\displaystyle \frac{\mathrm{ln}\rho }{m}}\right)ϵ\mu ^2f`$
$`+{\displaystyle \frac{1}{r}}ϵ(1\mu ^2)f+{\displaystyle \frac{\mu }{\rho }}[f].`$
With the help of Eq. (45), the $`m`$ component equation can be expressed

$`{\displaystyle \frac{}{t}}\left({\displaystyle \frac{ϵ\mu f}{\rho }}\right)`$ $`+`$ $`{\displaystyle \frac{}{m}}\left(4\pi r^2\rho ϵ\mu ^2{\displaystyle \frac{f}{\rho }}\right)+{\displaystyle \frac{1}{ϵ^2}}{\displaystyle \frac{}{ϵ}}\left(ϵ^4\mu \left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{f}{\rho }}\right)+`$ (157)
$`{\displaystyle \frac{}{\mu }}\left((1\mu ^2)\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]{\displaystyle \frac{ϵ\mu f}{\rho }}\right)`$
$`=`$ $`\left({\displaystyle \frac{2v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)ϵ\mu f+{\displaystyle \frac{1}{r}}ϵ(1\mu ^2)f+{\displaystyle \frac{\mu }{\rho }}[f].`$
Multiplication of Eqs. (155)-(157) by $`ϵ^2dϵd\mu d\phi `$ and integrating immediately yields Eqs. (122)-(124), expressions for particle 4-momentum balance in the comoving frame.

Just as the comoving frame $`t`$ and $`m`$ momentum-integrated balance equations—Eqs. (122) and (124)—can be combined to form the lab frame $`\stackrel{~}{t}`$ “conservation” equation, Eq. (126), Eqs. (155) and (157) can be combined as

$`{\displaystyle \frac{}{t}}\left[{\displaystyle \frac{ϵ(1+v\mu )f}{\rho }}\right]`$ $`+`$ $`{\displaystyle \frac{}{m}}\left(4\pi r^2\rho ϵ\mu (1+v\mu ){\displaystyle \frac{f}{\rho }}\right){\displaystyle \frac{1}{r^2\rho }}{\displaystyle \frac{}{t}}\left(r^2\rho v\right)ϵ\mu {\displaystyle \frac{f}{\rho }}`$
$`{\displaystyle \frac{1}{ϵ^2}}{\displaystyle \frac{}{ϵ}}\left(ϵ^4(1+v\mu )\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{f}{\rho }}\right)+`$
$`{\displaystyle \frac{}{\mu }}\left((1\mu ^2)\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]ϵ(1+v\mu ){\displaystyle \frac{f}{\rho }}\right)={\displaystyle \frac{1}{\rho }}(1+v\mu )[f].`$
Multiplication of this equation by $`ϵ^2dϵd\mu d\phi `$ and integrating immediately yields Eq. (126) for lab frame energy “conservation”.

We now examine the relationship between Eq. (154), a conservative formulation of particle number kinetics, and Eq. (LABEL:energyConservativeBoltzmannEulerian), a conservative formulation of particle energy kinetics. We multiply Eq. (154) by $`(1+v\mu )`$ and consider what it takes to get the terms in Eq. (LABEL:energyConservativeBoltzmannEulerian). In doing so, we examine carefully only $`O(v)`$ terms. First the time derivative terms:

$$ϵ(1+v\mu )\frac{}{t}\left(\frac{f}{\rho }\right)=\frac{}{t}\left[\frac{ϵ(1+v\mu )f}{\rho }\right]+O(v^2).$$
(159)
Next we compare the mass derivative terms:

$$ϵ(1+v\mu )\frac{}{m}\left(4\pi r^2\rho \mu \frac{f}{\rho }\right)=\frac{}{m}\left(4\pi r^2\rho ϵ\mu (1+v\mu )\frac{f}{\rho }\right)4\pi r^2ϵ\mu ^2\frac{v}{m}\frac{f}{\rho }.$$
(160)
In the notation of the previous subsection, the “extra terms” from the spacetime divergence are

$$E_S=4\pi r^2ϵ\mu ^2\frac{v}{m}\frac{f}{\rho }+O(v^2).$$
(161)
Next we relate the momentum divergence terms in the number and energy equations, looking first at the energy derivatives,

$`ϵ(1+v\mu ){\displaystyle \frac{1}{ϵ^2}}{\displaystyle \frac{}{ϵ}}\left(ϵ^3\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{f}{\rho }}\right)`$ $`=`$
$`{\displaystyle \frac{1}{ϵ^2}}{\displaystyle \frac{}{ϵ}}\left(ϵ^4(1+v\mu )\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{f}{\rho }}\right)`$ $``$ $`ϵ\left[\mu ^2\left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right){\displaystyle \frac{v}{r}}\right]{\displaystyle \frac{f}{\rho }}+O(v^2),`$
and then at the angle derivatives:

$`ϵ(1+v\mu ){\displaystyle \frac{}{\mu }}\left((1\mu ^2)\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]{\displaystyle \frac{f}{\rho }}\right)`$ $`=`$
$`{\displaystyle \frac{}{\mu }}\left((1\mu ^2)\left[{\displaystyle \frac{1}{r}}+\mu \left({\displaystyle \frac{3v}{r}}+{\displaystyle \frac{\mathrm{ln}\rho }{t}}\right)\right]ϵ(1+v\mu ){\displaystyle \frac{f}{\rho }}\right)`$ $``$ $`{\displaystyle \frac{ϵv(1\mu ^2)}{r}}{\displaystyle \frac{f}{\rho }}+O(v^2).`$ (163)
In the notation of the previous subsection, the “extra terms” from the momentum divergence are

$$E_M=ϵ\left[\mu ^2\left(\frac{3v}{r}+\frac{\mathrm{ln}\rho }{t}\right)\frac{v}{r}\right]\frac{f}{\rho }+\frac{ϵv(1\mu ^2)}{r}\frac{f}{\rho }+O(v^2).$$
(164)
Having $`E_S`$ and $`E_M`$, we are ready to verify Eq. (147), which must be satisfied for consistency between the particle number and energy balance equations. Having specified electrically neutral particles, and having chosen to work with the lab frame energy expression of Eq. (LABEL:energyConservativeBoltzmannEulerian) (in which the connection coefficients vanish), the right-hand side of Eq. (147) vanishes. Employing Eq. (127) for baryon conservation, from Eqs. (161) and (164) it is easy to see that

$$E_S+E_M=0$$
(165)
is indeed satisfied analytically. But in solving the “number conservative” Eq. (154) numerically, consistency with with the “energy conservative” Eq. (LABEL:energyConservativeBoltzmannEulerian) requires that Eq. (165) be satisfied numerically as well. Ref. Liebendörfer et al. (2002) provides an example of a finite difference representation of Eq. (154) that satisfies this criterion.

## V Conclusion

In this section we summarize our conservative formulations of kinetic theory, comment on their relation to moment formalisms, and discuss their possible application in the core-collapse supernova environment.

Having in mind computational radiative transfer in astrophysical environments, we have sought formulations of relativistic kinetic theory with the following properties: (1) they are expressed in terms of global, Eulerian (or “lab-frame”) spacetime coordinates $`\{x^\mu \}`$; (2) they are expressed in terms of convenient three-momentum coordinates $`\left\{u^{\widehat{i}}\right\}`$ (e.g. spherical polar), which are taken from the orthonormal momentum components $`\left\{p^{\widehat{i}}\right\}`$ measured by an observer comoving with the medium; and (3) they are “conservative”, having transparent connections to total particle number and 4-momentum balance as expressed in Eqs. (100) and (102).

To express our formulations having these properties, we here introduce the specific particle number flux vector

$$𝒩^\mu _{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}f$$
(166)
and the specific particle stress-energy tensor

$$𝒯^{\mu \nu }_{}^{\mu }{}_{\widehat{\mu }}{}^{}_{}^{\nu }{}_{\widehat{\nu }}{}^{}p^{\widehat{\mu }}p^{\widehat{\nu }}f,$$
(167)
where the transformation to the comoving frame $`_{}^{\widehat{\mu }}{}_{\mu }{}^{}`$ is given by equation (14). (While the adjective “specific” often denotes a quantity measured per unit mass, in this context we use it to denote the particle flux and stress-energy in a given invariant momentum space volume element.) While the distribution function $`f`$ of a given particle type of mass $`m`$ and charge $`e`$ obeys the Boltzmann equation (Eq. (97)), the specific particle number flux and stress-energy satisfy the conservative equations

$`{\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{}{x^\mu }}\left(\sqrt{g}𝒩^\mu \right)+`$
$`E(𝐩)\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|^1{\displaystyle \frac{}{u^{\widehat{i}}}}\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\mu }}{}^{}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\mu }\widehat{\nu }}{}^{}p^{\widehat{\nu }}\right){\displaystyle \frac{u^{\widehat{i}}}{p^{\widehat{j}}}}_{}^{\widehat{\mu }}{}_{\mu }{}^{}𝒩^\mu \right)`$
$`=[f],`$ (168)
$`{\displaystyle \frac{1}{\sqrt{g}}}{\displaystyle \frac{}{x^\nu }}\left(\sqrt{g}𝒯^{\mu \nu }\right)+`$
$`E(𝐩)\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|^1{\displaystyle \frac{}{u^{\widehat{i}}}}\left({\displaystyle \frac{1}{E(𝐩)}}\left|det\left[{\displaystyle \frac{𝐩}{𝐮}}\right]\right|\left(eF_{}^{\widehat{j}}{}_{\widehat{\nu }}{}^{}\mathrm{\Gamma }_{}^{\widehat{j}}{}_{\widehat{\nu }\widehat{\rho }}{}^{}p^{\widehat{\rho }}\right){\displaystyle \frac{u^{\widehat{i}}}{p^{\widehat{j}}}}_{}^{\widehat{\nu }}{}_{\nu }{}^{}𝒯^{\mu \nu }\right)`$
$`=F_{}^{\mu }{}_{\nu }{}^{}𝒩^\nu \mathrm{\Gamma }_{}^{\mu }{}_{\nu \rho }{}^{}𝒯^{\nu \rho }+_{}^{\mu }{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}[f].`$ (169)
In stating that Eqs. (168) and (169) constitute conservative formulations of particle kinetics, we mean that the connection to the balance equations of Eqs. (100) and (102) is transparent, in the following sense. We can use elementary calculus to form the familiar invariant momentum space volume element

$$\frac{d^3p}{E(𝐩)}=\frac{1}{E(𝐩)}\left|det\left[\frac{𝐩}{𝐮}\right]\right|d^3u,$$
(170)
where a transformation from orthonormal momentum components $`\left\{p^{\widehat{i}}\right\}`$ to some other set of coordinates (e.g. momentum space spherical coordinates $`\left\{u^{\widehat{i}}\right\}=\{|𝐩|,\vartheta ,\phi \}`$) has been performed. Multiplying Eqs. (168) and (169) by Eq. (170) and integrating, the terms with momentum space derivatives are obviously transformed into vanishing surface terms; the results are Eqs. (100) and (102) for total particle number and 4-momentum balance.

In terms of differential forms, the procedure for obtaining the conservative formulations of kinetic theory is straightforward. First, express the volume element $`\mathrm{\Omega }`$ in the phase space for particles of definite mass in terms of the desired spacetime and 3-momentum coordinates. Next, by contraction with the Liouville vector, form the hypersurface element $`\omega =L_m\mathrm{\Omega }`$. Then bring the exterior derivative $`d(f\omega )`$ into the form $`[f]\mathrm{\Omega }`$ by direct computation; Eq. (168) results on comparison with the Boltzmann equation. This result can then be used in conjunction with an evaluation of $`d(v_\mu p^\mu f\omega )`$ for arbitrary $`v_\mu `$ to obtain Eq. (169). The reason the procedure is straightforward is that the “heavy lifting” of transforming the Boltzmann equation into conservative forms is handled by two “levers” of considerable power: the generalized Stokes theorem, and the key relation $`d\omega =0`$, which is closely related to the relativistic Liouville theorem.

A concrete example of our formalism is provided in the appendix, which contains $`O(v)`$ equations for the specific particle number density, specific particle energy density, and their angular moments—all in flat spacetime, but in coordinates sufficiently general to represent rectangular, spherical, and cylindrical coordinate systems.

We now comment on the connection of Eqs. (168) and (169) to moment formalisms. In the usual treatments, if one writes the distribution function as a function of momentum variables as measured in a given frame (lab or comoving), it is natural to form moments by multiplying the distribution function by, for example, energies and angles measured in that frame, and integrating. Lo and behold, it turns out that these moments are number densities and fluxes, and energy/momentum densities and fluxes: components of a particle number flux vector and stress-energy tensor, measured in the same frame chosen to measure angles and energies. Traditional, then, are treatments in which the components of conserved tensors as measured in a given frame are expressed as functions of momentum variables as measured in that same frame.

But this traditional approach to moments may not be the most convenient, and Eqs. (168) and (169) provide an attractive alternative. For example, Liebendörfer et al. Liebendörfer et al. (2001, 2002) form moments in the traditional way, resulting in components of conserved tensors as measured in an orthonormal comoving frame. But it is the tensor components in the lab frame that one would like to check, either because it is the coordinate basis (at least in a spatially multidimensional simulation) and therefore natural to deal with, or because it is the basis in which energy is conserved (in a simulation in comoving coordinates in spherical symmetry). The required transformation of the tensor components between these frames leads to the numerical complexity mentioned in Sec. I, and discussed further in the latter parts of Subsecs. IV.1 and IV.2. In contrast, when the coordinate basis is in the lab frame as expected in spatially multidimensional simulations, integration of Eqs. (168) and (169) over momentum variables leads directly to the tensor components in the desired coordinate basis, even though the specific particle number and stress-energy are functions of comoving frame momentum variables. The insight here is that the frame in which the tensor components are measured need not be the same as that employed to obtain the momentum variables used to parametrize the particle distributions. The appendix provides an example of a moment formalism of this kind.

In simulations of systems like core-collapse supernovae—in which careful attention to energetics is critical—a number of possible approaches, based on the conservative formulations of kinetic theory presented in this paper, might be suggested.

First, the general approach of, for example, Refs. Mezzacappa et al. (2001); Liebendörfer et al. (2001, 2002) could be followed, in which the conservative particle number distribution equation (Eq. (168)) is solved, and the conservative particle energy equation (the time component of Eq. (169)) is used as a consistency check. The quantitiy $`\sqrt{g}𝒩^0`$, which might be called the specific particle number density, would be the primary neutrino distribution variable: It is the contribution of each comoving frame momentum bin to the lab frame particle number density. This approach makes number conservation a somewhat natural outcome, but energy conservation would require finite-differenced representations of various quantities to be “matched” in order that Eq. (147) be satisified numerically. This might be considered the most rigorous and self-consistent method. (Note that if one solves the “plain”, non-conservative Boltzmann equation—Eq. (97)—for the scalar distribution function $`f`$ as a function of comoving frame momentum variables, conservation of neither lab frame particle number or energy is straightforward. The same is true of methods (e.g., Ref. Burrows et al. (2000)) based on a non-conservative form of the transport equation for the comoving frame specific intensity.)

In order to avoid the intricate finite differencing of numerous terms demanded by this method, a second option would be the use of $`\sqrt{g}𝒯^{00}`$, which might be called the specific particle energy density, as the primary neutrino distribution variable. Designing a code around a differenced version of the $`\mu =0`$ component of Eq. (169) would make accurate accounting of total neutrino energy (as represented by the $`\mu =0`$ component of Eq. (102)) relatively straightfoward. Of course, the neutrino number balance equation expressed in terms of $`\sqrt{g}𝒯^{00}`$ would have numerical errors unless certain finite differencings were carefully designed. But with respect to the crucial energetics of the physical system, it is worth noting that there are a couple of factors mitigating the impact of errors in number conservation in comparison with errors in energy conservation. Errors in number conservation translate into errors in the electron fraction $`Y_e`$, which affect energy conservation through the equation of state, but: (1) Only $`\nu _e`$ and $`\left(\overline{\nu }\right)_e`$ affect $`Y_e`$, while all species impact the energy budget. Better to have two species contributing to error rather than six! (2) The effects of $`\nu _e`$ and $`\left(\overline{\nu }\right)_e`$ on $`Y_e`$ are opposite in sign (unlike their contributions to energy), so that to the extent that their distributions are similar, the impact of their errors on $`Y_e`$ may approximately cancel.

A third possibility would be to solve for both the specific particle energy density $`\sqrt{g}𝒯^{00}`$ and the specific particle number density $`\sqrt{g}𝒩^0`$. (Rampp and Janka Rampp and Janka (2002) solve for both number and energy distributions, but in the comoving frame; this limits the utility of their approach with respect to accurate tracking of lab frame quantities.) Instead of pre-defining both the boundaries and center values of bins in energy space, one could define the boundaries only and use the values of $`\sqrt{g}𝒯^{00}`$ and $`\sqrt{g}𝒩^0`$ (along with the transformation $`_{}^{\widehat{\mu }}{}_{\mu }{}^{}`$) to obtain center values of the energy bins in each spatial zone and each time step. The consistency of the solutions would be arguably reasonable as long as the derived center values of the energy bins do not wander outside the pre-defined bin boundaries.

*

## Appendix A Flat spacetime, $`O(v)`$ equations in a general coordinate system

In this appendix we present an application of the conservative formulations of kinetic theory derived in this paper: $`O(v)`$ equations for the specific particle number density and specific particle energy density, as well as angular integrals of these equations—all in flat spacetime, but in coordinates sufficiently general to represent rectangular, spherical, and cylindrical coordinate systems. The angle-integrated equations constitute “monochromatic” moment formalisms that provide an alternative to traditional variable Eddington factor methods of handling radiation transport. We specialize to massless, electrically neutral particles.

We begin by describing our spacetime and momentum space coordinate systems. While we assume flat spacetime, in order to accomodate curvilinear coordinate systems we employ a general spacetime coordinate labeling $`(x^\mu )=(x^1,x^2,x^3,t)^T`$. The line element is

$$ds^2=g_{\mu \nu }dx^\mu dx^\nu ,$$
(171)
with

$$(g_{\mu \nu })=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& a^2(x^1)& 0& 0\\ 0& 0& b^2(x^1)c^2(x^2)& 0\\ 0& 0& 0& 1\end{array}\right).$$
(172)
In this matrix expression for the metric components, rows and columns are ordered $`1,2,3,0`$. In Cartesian coordinates, $`(x^1,x^2,x^3)=(x,y,z)`$, and $`(a,b,c)=(1,1,1)`$. In spherical coordinates, $`(x^1,x^2,x^3)=(r,\theta ,\varphi )`$, and $`(a,b,c)=(r,r,\mathrm{sin}\theta )`$. In cylindrical coordinates, $`(x^1,x^2,x^3)=(r,z,\varphi )`$, and $`(a,b,c)=(1,r,1)`$. Our spacetime coordinate systems are “lab frames,” so that the equations we derive are Eulerian. Orthonormal “lab frame” coordinates, indicated by barred indices, are obtained by the transformation

$$dx^{\overline{\mu }}=e_{}^{\overline{\mu }}{}_{\mu }{}^{}dx^\mu ,$$
(173)
with

$$\left(e_{}^{\overline{\mu }}{}_{\mu }{}^{}\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& a(x^1)& 0& 0\\ 0& 0& b(x^1)c(x^2)& 0\\ 0& 0& 0& 1\end{array}\right).$$
(174)
The transformation to an orthonormal frame comoving with the fluid, indicated by indices adorned with a hat ($`\widehat{}`$), is

$$dx^{\widehat{\mu }}=\mathrm{\Lambda }_{}^{\widehat{\mu }}{}_{\overline{\mu }}{}^{}dx^{\overline{\mu }},$$
(175)
where to $`O(v)`$

$$\left(\mathrm{\Lambda }_{}^{\widehat{\mu }}{}_{\overline{\mu }}{}^{}\right)=\left(\begin{array}{cccc}1& 0& 0& v_{\overline{1}}\\ 0& 1& 0& v_{\overline{2}}\\ 0& 0& 1& v_{\overline{3}}\\ v_{\overline{1}}& v_{\overline{2}}& v_{\overline{3}}& 1\end{array}\right).$$
(176)
The bars on the velocity indices are a reminder that the fluid velocity is expressed in the orthonormal lab frame coordinate system. The combined transformation from the lab coordinate frame to the orthonormal comoving frame is

$$dx^{\widehat{\mu }}=_{}^{\widehat{\mu }}{}_{\mu }{}^{}dx^\mu ,$$
(177)
where

$$_{}^{\widehat{\mu }}{}_{\mu }{}^{}=\mathrm{\Lambda }_{}^{\widehat{\mu }}{}_{\overline{\mu }}{}^{}e_{}^{\overline{\mu }}{}_{\mu }{}^{}.$$
(178)
The neutrino 4-momentum is described in terms of its comoving frame components, $`\left(p^{\widehat{\mu }}\right)=(p^{\widehat{1}},p^{\widehat{2}},p^{\widehat{3}},p^{\widehat{0}})^T`$. Only three momentum variables are independent; we choose polar coordinates in momentum space, defined by

$`p^{\widehat{1}}`$ $`=`$ $`ϵ\mathrm{cos}\vartheta ,`$ (179)
$`p^{\widehat{2}}`$ $`=`$ $`ϵ\mathrm{sin}\vartheta \mathrm{cos}\phi ,`$ (180)
$`p^{\widehat{3}}`$ $`=`$ $`ϵ\mathrm{sin}\vartheta \mathrm{sin}\phi .`$ (181)
The radiation field is a function of the variables $`t,x^1,x^2,x^3,ϵ,\vartheta ,\phi `$. The invariant spacetime volume element in Eq. (28) is $`abcdx^1dx^2dx^3dt`$, and the invariant volume element on the mass shell in momentum space in Eq. (31) is $`ϵ\mathrm{sin}\vartheta dϵd\vartheta d\phi `$.

We first present an equation for the specific particle density $`𝒩`$, defined by

$$𝒩=p^tf=_{}^{t}{}_{\widehat{\mu }}{}^{}p^{\widehat{\mu }}f,$$
(182)
where $`f`$ is the invariant particle distribution function defined in Sec. III. The specific particle density is related to the lab frame particle density $`𝔫N^t`$ by

$$𝔫=𝒩ϵ\mathrm{sin}\vartheta dϵd\vartheta d\phi .$$
(183)
From Eq. (168), we find that the specific particle density satisfies

$$𝔻_t[𝒩]+𝔻_{x^1}[𝒩]+𝔻_{x^2}[𝒩]+𝔻_{x^3}[𝒩]+𝔻_\vartheta [𝒩]+𝔻_\phi [𝒩]+𝕆_ϵ[𝒩]+𝕆_\vartheta [𝒩]+𝕆_\phi [𝒩]=,$$
(184)
where $``$ is the invariant collision integral appearing in the Boltzmann equation (Eq. (97)), and

$`𝔻_t[𝒩]`$ $`=`$ $`{\displaystyle \frac{𝒩}{t}},`$ (185)
$`𝔻_{x^1}[𝒩]`$ $`=`$ $`{\displaystyle \frac{1}{ab}}{\displaystyle \frac{}{x^1}}(ab((1c_\vartheta ^2)v_{\overline{1}}+c_\vartheta (1s_\vartheta c_\phi v_{\overline{2}}s_\vartheta s_\phi v_{\overline{3}}))𝒩),`$ (186)
$`𝔻_{x^2}[𝒩]`$ $`=`$ $`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{}{x^2}}(c((1s_\vartheta ^2c_\phi ^2)v_{\overline{2}}+c_\phi s_\vartheta (1c_\vartheta v_{\overline{1}}s_\vartheta s_\phi v_{\overline{3}}))𝒩),`$ (187)
$`𝔻_{x^3}[𝒩]`$ $`=`$ $`{\displaystyle \frac{1}{bc}}{\displaystyle \frac{}{x^3}}(((1s_\vartheta ^2s_\phi ^2)v_{\overline{3}}+s_\vartheta s_\phi (1c_\vartheta v_{\overline{1}}s_\vartheta c_\phi v_{\overline{2}}))𝒩),`$ (188)
$`𝔻_\vartheta [𝒩]`$ $`=`$ $`{\displaystyle \frac{1}{s_\vartheta }}{\displaystyle \frac{}{\vartheta }}\left(\left({\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}s_\vartheta ^2c_\phi ^2+{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}s_\vartheta ^2s_\phi ^2\right)𝒩\right),`$ (189)
$`𝔻_\phi [𝒩]`$ $`=`$ $`{\displaystyle \frac{1}{s_\vartheta }}{\displaystyle \frac{}{\phi }}\left(\left({\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}c_\vartheta s_\vartheta c_\phi s_\phi +{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}c_\vartheta s_\vartheta c_\phi s_\phi +{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}s_\vartheta ^2s_\phi \right)𝒩\right),`$ (190)
$`𝕆_ϵ[𝒩]`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{ϵ}}(ϵ^2({\displaystyle \frac{v_{\overline{1}}}{x^1}}c_\vartheta ^2+{\displaystyle \frac{v_{\overline{2}}}{x^1}}c_\vartheta s_\vartheta c_\phi +{\displaystyle \frac{v_{\overline{3}}}{x^1}}c_\vartheta s_\vartheta s_\phi +`$ (191)
$`{\displaystyle \frac{1}{a}}\left({\displaystyle \frac{v_{\overline{1}}}{x^2}}c_\vartheta s_\vartheta c_\phi +{\displaystyle \frac{v_{\overline{2}}}{x^2}}s_\vartheta ^2c_\phi ^2+{\displaystyle \frac{v_{\overline{3}}}{x^2}}s_\vartheta ^2c_\phi s_\phi \right)+`$
$`{\displaystyle \frac{1}{bc}}\left({\displaystyle \frac{v_{\overline{1}}}{x^3}}c_\vartheta s_\vartheta s_\phi +{\displaystyle \frac{v_{\overline{2}}}{x^3}}s_\vartheta ^2c_\phi s_\phi +{\displaystyle \frac{v_{\overline{3}}}{x^3}}s_\vartheta ^2s_\phi ^2\right)+`$
$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}(s_\vartheta ^2c_\phi ^2v_{\overline{1}}c_\vartheta s_\vartheta c_\phi v_{\overline{2}})+`$
$`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}(s_\vartheta ^2s_\phi ^2v_{\overline{1}}c_\vartheta s_\vartheta s_\phi v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}(s_\vartheta ^2s_\phi ^2v_{\overline{2}}s_\vartheta ^2c_\phi s_\phi v_{\overline{3}}))𝒩),`$
$`𝕆_\vartheta [𝒩]`$ $`=`$ $`{\displaystyle \frac{1}{s_\vartheta }}{\displaystyle \frac{}{\vartheta }}(({\displaystyle \frac{v_{\overline{1}}}{x^1}}c_\vartheta s_\vartheta ^2+{\displaystyle \frac{v_{\overline{2}}}{x^1}}c_\vartheta ^2s_\vartheta c_\phi +{\displaystyle \frac{v_{\overline{3}}}{x^1}}c_\vartheta ^2s_\vartheta s_\phi +`$ (192)
$`{\displaystyle \frac{1}{a}}\left({\displaystyle \frac{v_{\overline{1}}}{x^2}}s_\vartheta ^3c_\phi +{\displaystyle \frac{v_{\overline{2}}}{x^2}}c_\vartheta s_\vartheta ^2c_\phi ^2+{\displaystyle \frac{v_{\overline{3}}}{x^2}}c_\vartheta s_\vartheta ^2c_\phi s_\phi \right)+`$
$`{\displaystyle \frac{1}{bc}}\left({\displaystyle \frac{v_{\overline{1}}}{x^3}}s_\vartheta ^3s_\phi +{\displaystyle \frac{v_{\overline{2}}}{x^3}}c_\vartheta s_\vartheta ^2c_\phi s_\phi +{\displaystyle \frac{v_{\overline{3}}}{x^3}}c_\vartheta s_\vartheta ^2s_\phi ^2\right)+`$
$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}(c_\phi s_\vartheta v_{\overline{2}}+c_\phi s_\vartheta ^3s_\phi ^2v_{\overline{2}}c_\phi ^2s_\vartheta ^3s_\phi v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}(c_\phi s_\vartheta ^3s_\phi ^2v_{\overline{2}}+s_\vartheta s_\phi v_{\overline{3}}+c_\phi ^2s_\vartheta ^3s_\phi v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}(c_\vartheta s_\vartheta ^2s_\phi ^2v_{\overline{2}}c_\vartheta s_\vartheta ^2c_\phi s_\phi v_{\overline{3}}))𝒩),`$
$`𝕆_\phi [𝒩]`$ $`=`$ $`{\displaystyle \frac{1}{s_\vartheta }}{\displaystyle \frac{}{\phi }}(({\displaystyle \frac{v_{\overline{2}}}{x^1}}c_\vartheta s_\phi +{\displaystyle \frac{v_{\overline{3}}}{x^1}}c_\vartheta c_\phi +`$ (193)
$`{\displaystyle \frac{1}{a}}\left({\displaystyle \frac{v_{\overline{2}}}{x^2}}s_\vartheta c_\phi s_\phi +{\displaystyle \frac{v_{\overline{3}}}{x^2}}s_\vartheta c_\phi ^2\right)+`$
$`{\displaystyle \frac{1}{bc}}\left({\displaystyle \frac{v_{\overline{2}}}{x^3}}s_\vartheta s_\phi ^2+{\displaystyle \frac{v_{\overline{3}}}{x^3}}s_\vartheta c_\phi s_\phi \right)+`$
$`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}(s_\vartheta ^3c_\phi s_\phi v_{\overline{1}}c_\vartheta s_\vartheta ^2c_\phi ^2s_\phi v_{\overline{2}}+c_\vartheta ^3c_\phi v_{\overline{3}}+c_\vartheta s_\vartheta ^2c_\phi ^3v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}(c_\vartheta s_\vartheta ^2s_\phi v_{\overline{1}}+c_\vartheta ^2s_\vartheta c_\phi s_\phi v_{\overline{2}}+s_\vartheta v_{\overline{3}}+c_\vartheta ^2s_\vartheta s_\phi ^2v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}(s_\vartheta ^3c_\phi s_\phi v_{\overline{1}}c_\vartheta ^3s_\phi v_{\overline{2}}c_\vartheta s_\vartheta ^2s_\phi ^3v_{\overline{2}}+c_\vartheta s_\vartheta ^2c_\phi s_\phi ^2v_{\overline{3}}))𝒩).`$
In these expressions we have employed the notation $`s_\vartheta =\mathrm{sin}\vartheta `$, $`c_\vartheta =\mathrm{cos}\vartheta `$. In accordance with the $`O(v)`$ limit, we have ignored acceleration terms and terms nonlinear in velocities and/or velocity spatial derivatives.

Next we consider an equation for the specific particle energy density $``$, defined by

$$=p^tp^tf=_{}^{t}{}_{\widehat{\mu }}{}^{}_{}^{t}{}_{\widehat{\nu }}{}^{}p^{\widehat{\mu }}p^{\widehat{\nu }}f.$$
(194)
The specific particle energy density is related to the lab frame particle energy density $`𝔢T^{tt}`$ by

$$𝔢=ϵ\mathrm{sin}\vartheta dϵd\vartheta d\phi .$$
(195)
From Eq. (169), we find that the specific particle density satisfies

$$𝔻_t[]+𝔻_{x^1}[]+𝔻_{x^2}[]+𝔻_{x^3}[]+𝔻_\vartheta []+𝔻_\phi []+𝕆_ϵ[]+𝕆_\vartheta []+𝕆_\phi []=ϵ(1+v_{\overline{1}}c_\vartheta +v_{\overline{2}}s_\vartheta c_\phi +v_{\overline{3}}s_\vartheta s_\phi ).$$
(196)
The terms on the left-hand side are given by Eqs. (185)-(193) with $`𝒩`$ replaced by $``$.555Because the relations between $`𝒩`$, $``$ and $`f`$ have different dependencies on velocity, expressions of Eqs. (184) and (196) in terms of $`f`$ would differ by terms of $`O(v^2)`$ and higher. While in principle these discrepancies vanish to $`O(v)`$, a simulation solving for both $`𝒩`$ and $``$ in which $`v^2`$ turns out to be non-negligible compared to unity can be expected to exhibit some inconsistency between values of $`f`$ derived from Eq. (182) and Eq. (194).  The $`t`$ component of the second term of Eq. (169), which we refer to as the “acceleration term”, vanishes for the flat spacetime metric of Eq. (172).

We now consider angle-integrated versions of Eqs. (184) and (196), which constitute equations for “monochromatic moments” of the radiation field. These moments are functions of the variables $`t,x^1,x^2,x^3,ϵ`$. In deriving and presenting these expressions, it is convenient to define

$`n_{\widehat{1}}`$ $`=`$ $`\mathrm{cos}\vartheta ,`$ (197)
$`n_{\widehat{2}}`$ $`=`$ $`\mathrm{sin}\vartheta \mathrm{cos}\phi ,`$ (198)
$`n_{\widehat{3}}`$ $`=`$ $`\mathrm{sin}\vartheta \mathrm{sin}\phi ,`$ (199)
the orthonormal comoving frame components of the unit 3-vector specifying the particle direction. In accordance with the $`O(v)`$ limit, we ignore acceleration terms and terms nonlinear in velocities and/or velocity spatial derivatives.

Our first equation of this type dictates the evolution of the lab frame monochromatic particle density $`𝔑`$, defined by

$$𝔑=𝒩\mathrm{sin}\vartheta d\vartheta d\phi .$$
(200)
It is related to the lab frame particle density $`𝔫N^t`$ by

$$𝔫=𝔑ϵ𝑑ϵ,$$
(201)
and its evolution is governed by

$$𝔻_t[𝔑]+𝔻_{x^1}[𝔑]+𝔻_{x^2}[𝔑]+𝔻_{x^3}[𝔑]+𝕆_ϵ[𝔑]=\mathrm{sin}\vartheta d\vartheta d\phi ,$$
(202)
where

$`𝔻_t[𝔑]`$ $`=`$ $`{\displaystyle \frac{𝔑}{t}},`$ (203)
$`𝔻_{x^1}[𝔑]`$ $`=`$ $`{\displaystyle \frac{1}{ab}}{\displaystyle \frac{}{x^1}}(ab(f_{\widehat{1}}+v_{\overline{1}}f_{\widehat{1}\widehat{1}}v_{\overline{1}}f_{\widehat{1}\widehat{2}}v_{\overline{2}}f_{\widehat{1}\widehat{3}}v_{\overline{3}})𝔑),`$ (204)
$`𝔻_{x^2}[𝔑]`$ $`=`$ $`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{}{x^2}}(c(f_{\widehat{2}}+v_{\overline{2}}f_{\widehat{1}\widehat{2}}v_{\overline{1}}f_{\widehat{2}\widehat{2}}v_{\overline{2}}f_{\widehat{2}\widehat{3}}v_{\overline{3}})𝔑),`$ (205)
$`𝔻_{x^3}[𝔑]`$ $`=`$ $`{\displaystyle \frac{1}{bc}}{\displaystyle \frac{}{x^3}}((f_{\widehat{3}}+v_{\overline{3}}f_{\widehat{1}\widehat{3}}v_{\overline{1}}f_{\widehat{2}\widehat{3}}v_{\overline{2}}f_{\widehat{3}\widehat{3}}v_{\overline{3}})𝔑),`$ (206)
$`𝕆_ϵ[𝔑]`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{ϵ}}(ϵ^2({\displaystyle \frac{v_{\overline{1}}}{x^1}}f_{\widehat{1}\widehat{1}}+{\displaystyle \frac{v_{\overline{2}}}{x^1}}f_{\widehat{1}\widehat{2}}+{\displaystyle \frac{v_{\overline{3}}}{x^1}}f_{\widehat{1}\widehat{3}}+`$ (207)
$`{\displaystyle \frac{1}{a}}\left({\displaystyle \frac{v_{\overline{1}}}{x^2}}f_{\widehat{1}\widehat{2}}+{\displaystyle \frac{v_{\overline{2}}}{x^2}}f_{\widehat{2}\widehat{2}}+{\displaystyle \frac{v_{\overline{3}}}{x^2}}f_{\widehat{2}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{bc}}\left({\displaystyle \frac{v_{\overline{1}}}{x^3}}f_{\widehat{1}\widehat{3}}+{\displaystyle \frac{v_{\overline{2}}}{x^3}}f_{\widehat{2}\widehat{3}}+{\displaystyle \frac{v_{\overline{3}}}{x^3}}f_{\widehat{3}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}(f_{\widehat{2}\widehat{2}}v_{\overline{1}}f_{\widehat{1}\widehat{2}}v_{\overline{2}})+`$
$`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}(f_{\widehat{3}\widehat{3}}v_{\overline{1}}f_{\widehat{1}\widehat{3}}v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}(f_{\widehat{3}\widehat{3}}v_{\overline{2}}f_{\widehat{2}\widehat{3}}v_{\overline{3}}))𝔑).`$
In these expressions, we have defined

$`f_{\widehat{i}}`$ $``$ $`{\displaystyle \frac{1}{𝔑}}{\displaystyle 𝒩n_{\widehat{i}}\mathrm{sin}\vartheta d\vartheta d\phi },`$
$`f_{\widehat{i}\widehat{j}}`$ $``$ $`{\displaystyle \frac{1}{𝔑}}{\displaystyle 𝒩n_{\widehat{i}}n_{\widehat{j}}\mathrm{sin}\vartheta d\vartheta d\phi },`$ (208)
which are akin to Eddington factors in traditional moment approaches to radiation transport.

Finally, we present a set of coupled equations for the lab frame monochromatic particle energy density $`𝔈`$ and orthonormal lab frame monochromatic particle momentum density $`𝔓^{\overline{i}}`$, defined by

$`𝔈`$ $`=`$ $`{\displaystyle \mathrm{sin}\vartheta d\vartheta d\phi },`$ (209)
$`𝔓^{\overline{i}}`$ $`=`$ $`e_{}^{\overline{i}}{}_{\mu }{}^{}{\displaystyle fp^tp^\mu \mathrm{sin}\vartheta d\vartheta d\phi }.`$ (210)
They are related to the lab frame particle energy density $`𝔢T^{tt}`$ and orthonormal lab frame particle momentum density $`𝔭^{\overline{i}}e_{}^{\overline{i}}{}_{\mu }{}^{}T^{t\mu }`$ by

$`𝔢`$ $`=`$ $`{\displaystyle 𝔈ϵ𝑑ϵ},`$ (211)
$`𝔭^{\overline{i}}`$ $`=`$ $`{\displaystyle 𝔓^{\overline{i}}ϵ𝑑ϵ}.`$ (212)
The evolution of the lab frame monochromatic particle energy density is governed by

$$𝔻_t[𝔈]+𝔻_{x^1}[𝔈]+𝔻_{x^2}[𝔈]+𝔻_{x^3}[𝔈]+𝕆_ϵ[𝔈]=ϵ(1+v_{\overline{1}}n_{\widehat{1}}+v_{\overline{2}}n_{\widehat{2}}+v_{\overline{3}}n_{\widehat{3}})\mathrm{sin}\vartheta d\vartheta d\phi ,$$
(213)
where

$`𝔻_t[𝔈]`$ $`=`$ $`{\displaystyle \frac{𝔈}{t}},`$ (214)
$`𝔻_{x^1}[𝔈]`$ $`=`$ $`{\displaystyle \frac{1}{ab}}{\displaystyle \frac{}{x^1}}\left(ab𝔓^{\overline{1}}\right),`$ (215)
$`𝔻_{x^2}[𝔈]`$ $`=`$ $`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{}{x^2}}\left(c𝔓^{\overline{2}}\right),`$ (216)
$`𝔻_{x^3}[𝔈]`$ $`=`$ $`{\displaystyle \frac{1}{bc}}{\displaystyle \frac{𝔓^{\overline{3}}}{x^3}},`$ (217)
$`𝕆_ϵ[𝔈]`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{ϵ}}(ϵ^2({\displaystyle \frac{v_{\overline{1}}}{x^1}}h_{\widehat{1}\widehat{1}}+{\displaystyle \frac{v_{\overline{2}}}{x^1}}h_{\widehat{1}\widehat{2}}+{\displaystyle \frac{v_{\overline{3}}}{x^1}}h_{\widehat{1}\widehat{3}}+`$ (218)
$`{\displaystyle \frac{1}{a}}\left({\displaystyle \frac{v_{\overline{1}}}{x^2}}h_{\widehat{1}\widehat{2}}+{\displaystyle \frac{v_{\overline{2}}}{x^2}}h_{\widehat{2}\widehat{2}}+{\displaystyle \frac{v_{\overline{3}}}{x^2}}h_{\widehat{2}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{bc}}\left({\displaystyle \frac{v_{\overline{1}}}{x^3}}h_{\widehat{1}\widehat{3}}+{\displaystyle \frac{v_{\overline{2}}}{x^3}}h_{\widehat{2}\widehat{3}}+{\displaystyle \frac{v_{\overline{3}}}{x^3}}h_{\widehat{3}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}(h_{\widehat{2}\widehat{2}}v_{\overline{1}}h_{\widehat{1}\widehat{2}}v_{\overline{2}})+`$
$`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}(h_{\widehat{3}\widehat{3}}v_{\overline{1}}h_{\widehat{1}\widehat{3}}v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}(h_{\widehat{3}\widehat{3}}v_{\overline{2}}h_{\widehat{2}\widehat{3}}v_{\overline{3}}))𝔈).`$
The evolution of the first component of the orthonormal lab frame monochromatic particle momentum density is governed by

$$𝔻_t[𝔓^{\overline{1}}]+𝔻_{x^1}[𝔓^{\overline{1}}]+𝔻_{x^2}[𝔓^{\overline{1}}]+𝔻_{x^3}[𝔓^{\overline{1}}]+𝔽_{22}[𝔓^{\overline{1}}]+𝔽_{33}[𝔓^{\overline{1}}]+𝕆_ϵ[𝔓^{\overline{1}}]=ϵ(n_{\widehat{1}}+v_{\overline{1}})\mathrm{sin}\vartheta d\vartheta d\phi ,$$
(219)
where

$`𝔻_t[𝔓^{\overline{1}}]`$ $`=`$ $`{\displaystyle \frac{𝔓^{\overline{1}}}{t}},`$ (220)
$`𝔻_{x^1}[𝔓^{\overline{1}}]`$ $`=`$ $`{\displaystyle \frac{1}{ab}}{\displaystyle \frac{}{x^1}}\left(ab\left(2v_{\overline{1}}𝔓^{\overline{1}}+(h_{\widehat{1}\widehat{1}}2(h_{\widehat{1}\widehat{1}\widehat{1}}v_{\overline{1}}+h_{\widehat{1}\widehat{1}\widehat{2}}v_{\overline{2}}+h_{\widehat{1}\widehat{1}\widehat{3}}v_{\overline{3}}))𝔈\right)\right),`$ (221)
$`𝔻_{x^2}[𝔓^{\overline{1}}]`$ $`=`$ $`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{}{x^2}}\left(c\left(𝔓^{\overline{2}}v_{\overline{1}}+𝔓^{\overline{1}}v_{\overline{2}}+(h_{\widehat{1}\widehat{2}}2(h_{\widehat{1}\widehat{1}\widehat{2}}v_{\overline{1}}+h_{\widehat{1}\widehat{2}\widehat{2}}v_{\overline{2}}+h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{3}}))𝔈\right)\right),`$ (222)
$`𝔻_{x^3}[𝔓^{\overline{1}}]`$ $`=`$ $`{\displaystyle \frac{1}{bc}}{\displaystyle \frac{}{x^3}}\left(𝔓^{\overline{3}}v_{\overline{1}}+𝔓^{\overline{1}}v_{\overline{3}}+(h_{\widehat{1}\widehat{3}}2(h_{\widehat{1}\widehat{1}\widehat{3}}v_{\overline{1}}+h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{2}}+h_{\widehat{1}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈\right),`$ (223)
$`𝔽_{22}[𝔓^{\overline{1}}]`$ $`=`$ $`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}\left(2𝔓^{\overline{2}}v_{\overline{2}}+(h_{\widehat{2}\widehat{2}}2(h_{\widehat{1}\widehat{2}\widehat{2}}v_{\overline{1}}+h_{\widehat{2}\widehat{2}\widehat{2}}v_{\overline{2}}+h_{\widehat{2}\widehat{2}\widehat{3}}v_{\overline{3}}))𝔈\right),`$ (224)
$`𝔽_{33}[𝔓^{\overline{1}}]`$ $`=`$ $`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}\left(2𝔓^{\overline{3}}v_{\overline{3}}+(h_{\widehat{3}\widehat{3}}2(h_{\widehat{1}\widehat{3}\widehat{3}}v_{\overline{1}}+h_{\widehat{2}\widehat{3}\widehat{3}}v_{\overline{2}}+h_{\widehat{3}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈\right),`$ (225)
$`𝕆_ϵ[𝔓^{\overline{1}}]`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{ϵ}}(ϵ^2({\displaystyle \frac{v_{\overline{1}}}{x^1}}h_{\widehat{1}\widehat{1}\widehat{1}}+{\displaystyle \frac{v_{\overline{2}}}{x^1}}h_{\widehat{1}\widehat{1}\widehat{2}}+{\displaystyle \frac{v_{\overline{3}}}{x^1}}h_{\widehat{1}\widehat{1}\widehat{3}}+`$ (226)
$`{\displaystyle \frac{1}{a}}\left({\displaystyle \frac{v_{\overline{1}}}{x^2}}h_{\widehat{1}\widehat{1}\widehat{2}}+{\displaystyle \frac{v_{\overline{2}}}{x^2}}h_{\widehat{1}\widehat{2}\widehat{2}}+{\displaystyle \frac{v_{\overline{3}}}{x^2}}h_{\widehat{1}\widehat{2}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{bc}}\left({\displaystyle \frac{v_{\overline{1}}}{x^3}}h_{\widehat{1}\widehat{1}\widehat{3}}+{\displaystyle \frac{v_{\overline{2}}}{x^3}}h_{\widehat{1}\widehat{2}\widehat{3}}+{\displaystyle \frac{v_{\overline{3}}}{x^3}}h_{\widehat{1}\widehat{3}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}(h_{\widehat{1}\widehat{2}\widehat{2}}v_{\overline{1}}h_{\widehat{1}\widehat{1}\widehat{2}}v_{\overline{2}})+`$
$`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}(h_{\widehat{1}\widehat{3}\widehat{3}}v_{\overline{1}}h_{\widehat{1}\widehat{1}\widehat{3}}v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}(h_{\widehat{1}\widehat{3}\widehat{3}}v_{\overline{2}}h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{3}}))𝔈).`$
The evolution of the second component of the orthonormal lab frame monochromatic particle momentum density is governed by

$$𝔻_t[𝔓^{\overline{2}}]+𝔻_{x^1}[𝔓^{\overline{2}}]+𝔻_{x^2}[𝔓^{\overline{2}}]+𝔻_{x^3}[𝔓^{\overline{2}}]+𝔽_{12}[𝔓^{\overline{2}}]+𝔽_{33}[𝔓^{\overline{2}}]+𝕆_ϵ[𝔓^{\overline{2}}]=ϵ(n_{\widehat{2}}+v_{\overline{2}})\mathrm{sin}\vartheta d\vartheta d\phi ,$$
(227)
where

$`𝔻_t[𝔓^{\overline{2}}]`$ $`=`$ $`{\displaystyle \frac{𝔓^{\overline{2}}}{t}},`$ (228)
$`𝔻_{x^1}[𝔓^{\overline{2}}]`$ $`=`$ $`{\displaystyle \frac{1}{ab}}{\displaystyle \frac{}{x^1}}\left(ab\left(𝔓^{\overline{2}}v_{\overline{1}}+𝔓^{\overline{1}}v_{\overline{2}}+(h_{\widehat{1}\widehat{2}}2(h_{\widehat{1}\widehat{1}\widehat{2}}v_{\overline{1}}+h_{\widehat{1}\widehat{2}\widehat{2}}v_{\overline{2}}+h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{3}}))𝔈\right)\right),`$ (229)
$`𝔻_{x^2}[𝔓^{\overline{2}}]`$ $`=`$ $`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{}{x^2}}\left(c\left(2𝔓^{\overline{2}}v_{\overline{2}}+(h_{\widehat{2}\widehat{2}}2(h_{\widehat{1}\widehat{2}\widehat{2}}v_{\overline{1}}+h_{\widehat{2}\widehat{2}\widehat{2}}v_{\overline{2}}+h_{\widehat{2}\widehat{2}\widehat{3}}v_{\overline{3}}))𝔈\right)\right),`$ (230)
$`𝔻_{x^3}[𝔓^{\overline{2}}]`$ $`=`$ $`{\displaystyle \frac{1}{bc}}{\displaystyle \frac{}{x^3}}\left(𝔓^{\overline{3}}v_{\overline{2}}+𝔓^{\overline{2}}v_{\overline{3}}+(h_{\widehat{2}\widehat{3}}2(h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{1}}+h_{\widehat{2}\widehat{2}\widehat{3}}v_{\overline{2}}+h_{\widehat{2}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈\right),`$ (231)
$`𝔽_{12}[𝔓^{\overline{2}}]`$ $`=`$ $`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}\left(𝔓^{\overline{2}}v_{\overline{1}}+𝔓^{\overline{1}}v_{\overline{2}}+(h_{\widehat{1}\widehat{2}}2(h_{\widehat{1}\widehat{1}\widehat{2}}v_{\overline{1}}+h_{\widehat{1}\widehat{2}\widehat{2}}v_{\overline{2}}+h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{3}}))𝔈\right),`$ (232)
$`𝔽_{33}[𝔓^{\overline{2}}]`$ $`=`$ $`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}\left(2𝔓^{\overline{3}}v_{\overline{3}}+(h_{\widehat{3}\widehat{3}}2(h_{\widehat{1}\widehat{3}\widehat{3}}v_{\overline{1}}+h_{\widehat{2}\widehat{3}\widehat{3}}v_{\overline{2}}+h_{\widehat{3}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈\right),`$ (233)
$`𝕆_ϵ[𝔓^{\overline{2}}]`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{ϵ}}(ϵ^2({\displaystyle \frac{v_{\overline{1}}}{x^1}}h_{\widehat{1}\widehat{1}\widehat{2}}+{\displaystyle \frac{v_{\overline{2}}}{x^1}}h_{\widehat{1}\widehat{2}\widehat{2}}+{\displaystyle \frac{v_{\overline{3}}}{x^1}}h_{\widehat{1}\widehat{2}\widehat{3}}+`$ (234)
$`{\displaystyle \frac{1}{bc}}\left({\displaystyle \frac{v_{\overline{1}}}{x^3}}h_{\widehat{1}\widehat{2}\widehat{3}}+{\displaystyle \frac{v_{\overline{2}}}{x^3}}h_{\widehat{2}\widehat{2}\widehat{3}}+{\displaystyle \frac{v_{\overline{3}}}{x^3}}h_{\widehat{2}\widehat{3}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{a}}\left({\displaystyle \frac{v_{\overline{1}}}{x^2}}h_{\widehat{1}\widehat{2}\widehat{2}}+{\displaystyle \frac{v_{\overline{2}}}{x^2}}h_{\widehat{2}\widehat{2}\widehat{2}}+{\displaystyle \frac{v_{\overline{3}}}{x^2}}h_{\widehat{2}\widehat{2}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}(h_{\widehat{2}\widehat{2}\widehat{2}}v_{\overline{1}}h_{\widehat{1}\widehat{2}\widehat{2}}v_{\overline{2}})+`$
$`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}(h_{\widehat{2}\widehat{3}\widehat{3}}v_{\overline{1}}h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}(h_{\widehat{2}\widehat{3}\widehat{3}}v_{\overline{2}}h_{\widehat{2}\widehat{2}\widehat{3}}v_{\overline{3}}))𝔈).`$
The evolution of the third component of the orthonormal lab frame monochromatic particle momentum density is governed by

$$𝔻_t𝔓^{\overline{3}}+𝔻_{x^1}[𝔓^{\overline{3}}]+𝔻_{x^2}[𝔓^{\overline{3}}]+𝔻_{x^3}[𝔓^{\overline{3}}]+𝔽_{13}[𝔓^{\overline{3}}]+𝔽_{23}[𝔓^{\overline{3}}]+𝕆_ϵ[𝔓^{\overline{3}}]=ϵ(n_{\widehat{3}}+v_{\overline{3}})\mathrm{sin}\vartheta d\vartheta d\phi ,$$
(235)
where

$`𝔻_t[𝔓^{\overline{3}}]`$ $`=`$ $`{\displaystyle \frac{𝔓^{\overline{3}}}{t}},`$ (236)
$`𝔻_{x^1}[𝔓^{\overline{3}}]`$ $`=`$ $`{\displaystyle \frac{1}{ab}}{\displaystyle \frac{}{x^1}}\left(ab\left(𝔓^{\overline{3}}v_{\overline{1}}+𝔓^{\overline{1}}v_{\overline{3}}+(h_{\widehat{1}\widehat{3}}2(h_{\widehat{1}\widehat{1}\widehat{3}}v_{\overline{1}}+h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{2}}+h_{\widehat{1}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈\right)\right),`$ (237)
$`𝔻_{x^2}[𝔓^{\overline{3}}]`$ $`=`$ $`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{}{x^2}}\left(c\left(𝔓^{\overline{3}}v_{\overline{2}}+𝔓^{\overline{2}}v_{\overline{3}}+(h_{\widehat{2}\widehat{3}}2(h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{1}}+h_{\widehat{2}\widehat{2}\widehat{3}}v_{\overline{2}}+h_{\widehat{2}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈\right)\right),`$ (238)
$`𝔻_{x^3}[𝔓^{\overline{3}}]`$ $`=`$ $`{\displaystyle \frac{1}{bc}}{\displaystyle \frac{}{x^3}}\left(2𝔓^{\overline{3}}v_{\overline{3}}+(h_{\widehat{3}\widehat{3}}2(h_{\widehat{1}\widehat{3}\widehat{3}}v_{\overline{1}}+h_{\widehat{2}\widehat{3}\widehat{3}}v_{\overline{2}}+h_{\widehat{3}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈\right),`$ (239)
$`𝔽_{13}[𝔓^{\overline{3}}]`$ $`=`$ $`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}\left(𝔓^{\overline{3}}v_{\overline{1}}+𝔓^{\overline{1}}v_{\overline{3}}+(h_{\widehat{1}\widehat{3}}2(h_{\widehat{1}\widehat{1}\widehat{3}}v_{\overline{1}}+h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{2}}+h_{\widehat{1}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈\right),`$ (240)
$`𝔽_{23}[𝔓^{\overline{3}}]`$ $`=`$ $`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}\left(𝔓^{\overline{3}}v_{\overline{2}}+𝔓^{\overline{2}}v_{\overline{3}}+(h_{\widehat{2}\widehat{3}}2(h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{1}}+h_{\widehat{2}\widehat{2}\widehat{3}}v_{\overline{2}}+h_{\widehat{2}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈\right),`$ (241)
$`𝕆_ϵ[𝔓^{\overline{3}}]`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{ϵ}}(ϵ^2({\displaystyle \frac{v_{\overline{1}}}{x^1}}h_{\widehat{1}\widehat{1}\widehat{3}}+{\displaystyle \frac{v_{\overline{2}}}{x^1}}h_{\widehat{1}\widehat{2}\widehat{3}}+{\displaystyle \frac{v_{\overline{3}}}{x^1}}h_{\widehat{1}\widehat{3}\widehat{3}}+`$ (242)
$`{\displaystyle \frac{1}{a}}\left({\displaystyle \frac{v_{\overline{1}}}{x^2}}h_{\widehat{1}\widehat{2}\widehat{3}}+{\displaystyle \frac{v_{\overline{2}}}{x^2}}h_{\widehat{2}\widehat{2}\widehat{3}}+{\displaystyle \frac{v_{\overline{3}}}{x^2}}h_{\widehat{2}\widehat{3}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{bc}}\left({\displaystyle \frac{v_{\overline{1}}}{x^3}}h_{\widehat{1}\widehat{3}\widehat{3}}+{\displaystyle \frac{v_{\overline{2}}}{x^3}}h_{\widehat{2}\widehat{3}\widehat{3}}+{\displaystyle \frac{v_{\overline{3}}}{x^3}}h_{\widehat{3}\widehat{3}\widehat{3}}\right)+`$
$`{\displaystyle \frac{1}{a}}{\displaystyle \frac{a}{x^1}}(h_{\widehat{2}\widehat{2}\widehat{3}}v_{\overline{1}}h_{\widehat{1}\widehat{2}\widehat{3}}v_{\overline{2}})+`$
$`{\displaystyle \frac{1}{b}}{\displaystyle \frac{b}{x^1}}(h_{\widehat{3}\widehat{3}\widehat{3}}v_{\overline{1}}h_{\widehat{1}\widehat{3}\widehat{3}}v_{\overline{3}})+`$
$`{\displaystyle \frac{1}{ac}}{\displaystyle \frac{c}{x^2}}(h_{\widehat{3}\widehat{3}\widehat{3}}v_{\overline{2}}h_{\widehat{2}\widehat{3}\widehat{3}}v_{\overline{3}}))𝔈).`$
In these energy and momentum equations, we have defined

$`h_{\widehat{i}\widehat{j}}`$ $``$ $`{\displaystyle \frac{1}{𝔈}}{\displaystyle n_{\widehat{i}}n_{\widehat{j}}\mathrm{sin}\vartheta d\vartheta d\phi },`$ (243)
$`h_{\widehat{i}\widehat{j}\widehat{k}}`$ $``$ $`{\displaystyle \frac{1}{𝔈}}{\displaystyle n_{\widehat{i}}n_{\widehat{j}}n_{\widehat{k}}\mathrm{sin}\vartheta d\vartheta d\phi },`$ (244)
which are akin to Eddington factors in traditional moment approaches to radiation transport. We point out that $`f_{\widehat{i}}`$ is not equal to $`𝔓^{\overline{i}}/𝔈`$, and that the second angular moment factors $`f_{\widehat{i}\widehat{j}}`$ and $`h_{\widehat{i}\widehat{j}}`$ are not equal to each other. This is because $`𝒩`$ and $``$ have different angular dependencies, as can be seen from Eqs. (182) and (194). This complication is a result of taking moments of lab frame number and energy distributions with respect to comoving frame angular variables.

###### Acknowledgements.

We are grateful to Matthias Liebendörfer, Bronson Messer, and Raph Hix for insightful discussions about various aspects of this work, and helpful suggestions on the manuscript. We thank Dimitri Mihalas and Steve Bruenn for discussions about momentum space coordinate systems. This work was supported by the U.S. Department of Energy (DoE) HENP Scientific Discovery Through Advanced Computing Program; Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the DoE under contract DE-AC05-00OR22725; the Joint Institute for Heavy Ion Research; and a DoE PECASE grant.