Conservation Laws and Expectation Value Matching

In the S approach, the complete dynamics in the phase space of the transported particles is determined by one kinetic equation (15). Its integration over the momentum phase space is equivalent to a continuity equation (24), its integration weighed with particle energies, , is equivalent to an energy equation (27), and its integration weighed with the particle direction cosine and energy would lead to a momentum equation. These derived equations certainly fulfill the macroscopic conservation laws as given in Eq. (7). But it is a challenge to obtain the same level of consistency in a finite difference representation of the Boltzmann equation. We illustrate this with the following little example:

For simplicity we assume that a distribution function, , of particles propagating with speed only depends on time, , and on one dimension in the phase space, . The simplest form of a transport equation would then relate the time derivative of the distribution function to an advection term,

$$\frac{\partial f}{\partial t}+c\frac{\partial f}{\partial x}=0.$$ (31)

If we now ask for the time evolution of the expectation value of the distribution function with respect to an operator, , we apply integrations by parts to find

$$\frac{\partial }{c\partial t}\int gfdx=\int \left( \frac{\pa... ...tial g}{\partial x}\right) fdx-\left[ gf\right] _{\partial x}.$$ (32)

Here, the time evolution of the expectation value is described by other expectation values and the value of at the boundary of the domain of integration. This equation has an immediate physical interpretation: The change of the expectation value along the characteristics of the particle flow is simply given by the change of the weight along the characteristics of the particle flow because the phase space volume stays constant,

$$\left( \frac{\partial }{c\partial t}+\frac{\partial }{\parti... ...tial g}{c\partial t}+\frac{\partial g}{\partial x}\right) fdx.$$ (33)

The second term on the left hand side is equivalent to the boundary term in Eq. (). Let us now assume that numerical stability requires upwind differencing of the advection term in Eq. (), and, for simplicity, that the wind is unidirectional. A finite difference representation would then read like

$$\frac{\partial f_{i'}}{\partial t}+c\frac{f_{i'}-f_{i'-1}}{dx_{i'}}=0.$$ (34)

We use the convention that a prime at the index points to the zone center value while an integer index points to the zone edge. We simplify the example once more and assume that the operator is not time-dependent. With this, we track above integration by parts in the finite difference representation and calculate the evolution of the expectation value of operator ,

$$\frac{\partial }{c\partial t}\sum ^{n}_{i=1}g_{i'}f_{i'}dx_{i'} = \sum ^{n}_{i=1}g_{i'}\frac{\partial f_{i'}}{c\partial t}dx_{i'}=-\sum ^{n}_{i=1}g_{i'}\left( f_{i'}-f_{i'-1}\right)$$

$$= -\sum ^{n}_{i=1}g_{i'}f_{i'}+\sum ^{n-1}_{i=0}g_{i'+1}f_{i'}$$

$$= \sum ^{n}_{i=1}\frac{g_{i'+1}-g_{i'}}{dx_{i'}}f_{i'}dx_{i'}-\left[ g_{n'+1}f_{n'}-g_{0'+1'}f_{0'}\right] .$$ (35)

Firstly, we remark that there is a discrete analogue to the integration by parts and that Eq. () has exactly the same structure as Eq. (). Secondly, we note that the choice of finite differencing for the advection term in Eq. () determines the finite differencing of the weight function in the expectation value on the right hand side of Eq. (). In the following, we will frequently use this type of investigation to optimize for requirements (3) and (4): First we choose a stable finite difference representation of a term in the transport equation. Then we evaluate relevant expectation values by performing discrete ``integrations by parts'', analogously to this simple example. The emerging finite difference representations of the expectation values along one phase space dimension are then used as prefactors for the finite differencing of terms related to other phase space dimensions in order to build a consistent finite difference representation of the full transport equation.

In order to construct a conservative finite difference representation of Eq. (15), we analyze the structure of the conservation equations in the continuum world, where insight is not buried in discretization indices. The emergence of lepton number conservation in Eq. (24) is straightforward, because the observer corrections are already written in a form that allows an immediate integration over energy or the angle cosine. Number conservation will also emerge naturally in the finite difference representation when we finite difference the same basic structure of the equation. The energy conservation equation in the frame of a distant observer, however, introduces the weight in Eq. (26). It depends on all four phase space dimensions and leads to the contribution of several expectation values of the distribution function in the evaluation of the energy conservation equation,

$$\int \left( \Gamma +u\mu \right) E\left[ C_{t}+D_{a}+D_{\mu }+D_{E}+O_{E}+O_{\mu }-C_{c}\right] d\mu E^{2}dE=0.$$ (36)

All terms in the Boltzmann equation (15), , , , , , and , are written as the derivative of an expression with respect to time, rest mass, angle cosine, or particle energy. We now perform an integration by parts with respect to these integration variables along the lines of Eq. (). A multitude of correction terms proportional to arise. As in Eq. (), they describe the evolution of the weight along the characteristic of the particle flow. However, we showed in Liebendoerfer_Mezzacappa_Thielemann_01 that is nearly a constant of motion along the characteristics. Therefore, it is natural that most of the partial terms actually cancel in the total evolution of the radiation energy in the frame of a distant observer. If, like in Eq. (), the integrations by parts are also carefully followed in the finite difference representation, mutual cancellations of important expectation values can be forced to be exact--independently of the resolution. Because these canceling terms individually reach large values around and after bounce, this expectation value matching is an essential step for the conservation of energy in the finite difference representation. We perform the integrations by parts in Eq. () and start with the identification of canceling contributions in the following overview of contributions from , , , , , , in the transport equation (15):

$$C_{t}: \frac{\partial }{\alpha \partial t}\left( \Gamma J\right) -\frac{... ...ial }{\alpha \partial t}\left( uH\right) -\frac{\partial u}{\alpha \partial t}H$$

$$D_{a}: \frac{\partial }{\alpha \partial a}\left[ 4\pi r^{2}\alpha \rho \Gamma H\right] -4\pi r^{2}\rho \frac{\partial \Gamma }{\partial a}H$$

$$+\frac{\partial }{\alpha \partial a}\left[ 4\pi r^{2}\alpha \rho uK\right] -4\pi r^{2}\rho \frac{\partial u}{\partial a}K$$

$$D_{\mu }: -\Gamma \frac{u}{r}\left( J-K\right) +\Gamma u\frac{\partial \Phi }{\partial r}\left( J-K\right)$$

$$D_{E}: \Gamma \frac{\partial \Phi }{\partial r}\left( \Gamma H+uK\right)$$

$$O_{E}: -\left[ \frac{\partial \ln \rho }{\alpha \partial t}+\frac{2u}{r}... ...Q\right) +\Gamma \frac{u}{r}\left( J-K\right) +\frac{u^{2}}{r}\left( H-Q\right)$$

$$O_{\mu }: -\left[ \frac{\partial \ln \rho }{\alpha \partial t}+\frac{2u}{r}\right] u\left( H-Q\right) -\frac{u^{2}}{r}\left( H-Q\right)$$

$$C_{c}: \Gamma \int \left( \frac{j}{\rho }-\chi F\right) E^{3}dEd\mu -u\int \chi FE^{3}dE\mu d\mu .$$ (37)

We will neglect terms that are nonlinear in the radiation field, i.e. terms that describe a gravitational interaction of the radiation field with itself. From the hydrodynamics equations in Misner_Sharp_64,May_White_66 we then derive the following useful relationships for the coefficients in Eq. ().

$$\frac{1}{\Gamma }\frac{\partial \Gamma }{\alpha \partial t} = u\frac{\partial \Phi }{\partial r}$$ (38)

$$-\left( \frac{\partial \ln \rho }{\alpha \partial t}+\frac{2u}{r}\right) = \frac{4\pi r^{2}\rho }{\Gamma }\frac{\partial u}{\partial a}$$ (39)

$$\frac{\partial u}{\alpha \partial t} = \Gamma ^{2}\frac{\partial \Phi }{\partial r}-\frac{m}{r^{2}}-4\pi rp.$$ (40)

We may now identify cancellations in Eq. () and label them for later reference: The fourth term of cancels with the first term of ( ) by Eq. (). The first and second term of cancel with the third and fourth term in ( ). These are the only O cancellations and therefore the most important ones. In an analogous notation, we find the following higher order cancellations: ( ) by Eq. (), ( ), ( ), ( ). The remaining terms, , reduce by Eq. () and the definition of to the general relativistic term,

$$\left( -\frac{\partial u}{\alpha \partial t}-\Gamma \frac{\p... ...}{\partial r}\right) H=4\pi r\rho \left( 1+e+p/\rho \right) H.$$ (41)

They enter in this form the radiation energy conservation equation (27).

We will adopt the following strategy in the finite differencing of Eq. (15) (the reader is invited to draw lines in Eq. () to visualize the relationship between the different terms): (i) The finite differencing of and are straightforward. (ii) Appreciable experience has been gathered in previous work (see e.g. Lewis_Miller_84) with the finite differencing of the O terms in and . They are assumed to be well-chosen in Mezzacappa_Bruenn_93a and therefore not subject to changes. (iii) Based on this, the cancellation ( ) dictates the finite difference representation of in and . The cancellation ( ) sets the finite difference representation of in and . (iv) The cancellation ( ) dictates the finite difference representation of in and . The cancellation ( ) propagates the finite difference representation to in . (v) The evaluation of according to Eq. () constrains the finite difference representation of in and therewith defines . (vi) And finally we choose a different finite difference representation for in and , as suggested by the cancellation ( ). By these six chains, all terms in Eq. (15) become constrained.