Conservative Moment Equations for Neutrino Radiation Transport with Limited Relativity
aa r X i v : . [ a s t r o - ph . S R ] D ec D RAFT VERSION S EPTEMBER
18, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
CONSERVATIVE MOMENT EQUATIONS FOR NEUTRINO RADIATION TRANSPORT WITH LIMITED RELATIVITY E IRIK E NDEVE , C HRISTIAN
Y. C
ARDALL , AND A NTHONY M EZZACAPPA
Draft version September 18, 2018
ABSTRACTWe derive conservative, multidimensional, energy-dependent moment equations for neutrino transport incore-collapse supernovae and related astrophysical systems, with particular attention to the consistency of con-servative four-momentum and lepton number transport equations. After taking angular moments of conserva-tive formulations of the general relativistic Boltzmann equation, we specialize to a conformally flat spacetime,which also serves as the basis for four further limits. Two of these—the multidimensional special relativisticcase, and a conformally flat formulation of the spherically symmetric general relativistic case—are given in ap-pendices for the sake of comparison with extant literature. The third limit is a weak-field, ‘pseudo-Newtonian’approach (Kim et al. 2009, 2012) in which the source of the gravitational potential includes the trace of thestress-energy tensor (rather than just the mass density), and all orders in fluid velocity v are retained. Our pri-mary interest here is in the fourth limit: ‘ O ( v ) ’ moment equations for use in conjunction with Newtonian self-gravitating hydrodynamics. We show that the concept of ‘ O ( v ) ’ transport requires care when dealing with bothconservative four-momentum and conservative lepton number transport, and present two self-consistent op-tions: ‘ O ( v ) -plus’ transport, in which an O ( v ) energy equation combines with an O ( v ) momentum equationto give an O ( v ) number equation; and ‘ O ( v ) -minus’ transport, in which an O ( v ) energy equation combineswith an O (1) momentum equation to give an O ( v ) number equation. Subject headings: neutrinos — radiative transfer — supernovae: general INTRODUCTION
Detailed simulations of neutrino transport are central tounderstanding the core-collapse supernova explosion mech-anism (e.g., Mezzacappa 2005; Woosley & Janka 2005;Kotake et al 2006; Janka 2012). Neutrinos of all flavors (i.e.,electron, mu, and tau neutrinos and antineutrinos) carry away ∼ of the gravitational energy released ( ∼ erg) dur-ing core collapse of a massive star ( M & M ⊙ ). As thedominant energy source (in the case of slowly-rotating pro-genitors), it is likely that the neutrinos streaming from theneutrinospheres of a proto-neutron star (PNS; the collapsedcore of a massive star) deposit enough energy into the stellarfluid to power the supernova: about 5-10% ( ∼ erg) of theneutrino energy emitted in the first second after core bounce( & erg) must be transferred to give rise to a neutrino-driven supernova.Certainly, robust and accurate numerical methods must bedeployed. Numerical simulations involving multi-physicscodes (coupling neutrino radiation transport, magnetohydro-dynamics, nuclear reaction kinetics and equations of state,and gravity, preferably—and ultimately—in general relativ-ity) must be evolved for ∼ O (10 - ) timesteps to study indetail the energy transfer from the neutrino radiation field tothe supernova matter, and the onset and development of theexplosion. In particular, conserved quantities (i.e., total en-ergy and lepton number) must be preserved within tolerablemargins, before firm conclusions can be drawn from the sim-ulation outcomes.So far, only spherically symmetric simulations havereached the desired level of realism in neutrino transport, Computer Science and Mathematics Division, Oak Ridge NationalLaboratory, Oak Ridge, TN 37831-6354, USA; [email protected] Physics Division, Oak Ridge National Laboratory, Oak Ridge, TN37831-6354, USA Department of Physics and Astronomy, University of Tennessee,Knoxville, TN 37996-1200, USA where general relativistic, angle- and energy-dependent neu-trino radiation hydrodynamics simulations can be performed(e.g., Liebend¨orfer et al. 2001, 2004; Sumiyoshi et al.2005). However, spherically symmetric simulations do notresult in neutrino-driven explosions (Rampp & Janka 2000;Liebend¨orfer et al. 2001; Thompson et al. 2003, except inthe case of the lightest progenitors with O-Ne-Mg cores;Kitaura et al. 2006), and, based on results from 2D (ax-isymmetric) simulations, there is now an emerging consen-sus that multidimensional effects due to hydrodynamics insta-bilities are important in shaping the core-collapse supernovaexplosion (e.g., Herant et al. 1994; Burrows et al. 1995;Janka & Muller 1996; Fryer & Warren 2002; Burrows et al.2006, 2007; Bruenn et al. 2009; Marek & Janka 2009;Suwa et al. 2010; M¨uller et al. 2012a,b; Bruenn et al.2012).Neutrino transport in multidimensional simulations is farless mature. Discrepant results from 2D simulations by dif-ferent research groups, in part due to the physics includedin the models and the approximations made, remain to befully investigated. It appears clear, however, that energy-dependent neutrino transport (including “observer correc-tions” due to fluid motions and inelastic scattering; e.g.,Lentz et al. 2012a) as well as a comprehensive set of weak-interactions are required (e.g., Lentz et al. 2012b). 3D sim-ulations are still in their infancy, although recent progresshas been made (e.g., Fryer & Young 2007; Takiwaki et al.2012).Current approaches to neutrino transport in multidimen-sional simulations include the so-called ray-by-ray ap-proximation with the multigroup (energy-dependent) flux-limited diffusion approximation (e.g., Bruenn et al. 2009,2012), the two-moment model with approximate Boltz-mann closure (e.g., Marek & Janka 2009; M¨uller et al.2012a), and the so-called isotropic diffusion source approx-imation (e.g., Liebend¨orfer et al. 2009; Suwa et al. 2010). Endeve et al.Other recent approaches to neutrino radiation transport in-clude multidimensional, multigroup flux-limited diffusion—valid to O ( v ) (e.g., Swesty & Myra 2009; Zhang et al.2012) or to O (1) (e.g., Burrows et al. 2006)—or multi-dimensional, energy-integrated (grey) or multigroup two-moment models with analytic closure (e.g., Kuroda et al.2012; Obergaulinger & Janka 2011; Shibata et al. 2011).Ott et al. (2008) simulated the post-bounce phase of core-collapse supernovae in axial symmetry using O (1) multi-group and multiangle neutrino transport (see also the recentdevelopments by Sumiyoshi & Yamada 2012).Lack of sufficient computational resources is the primaryreason for the various approximations made in multidimen-sional simulations of neutrino transport in core-collapse su-pernovae. Gravitational collapse, a steep density cliff atthe contracting surface of the PNS, and regions of post-shock turbulence recommend regions of high spatial res-olution. Moreover, differences in the evolution of im-portant hydrodynamic phenomena in two and three spatialdimensions—for example, the Standing Accretion Shock In-stability (e.g., Blondin et al. 2003; Blondin & Mezzacappa2007), or turbulence (Ishihara et al. 2009; Boffetta & Ecke2012)—suggest that simulations should ideally (and eventu-ally) be performed in 3D. Recent explorations using so-calledneutrino light-bulb models emphasize differences between 2Dand 3D simulations (e.g., Nordhaus et al. 2010; Hanke et al.2012; Burrows et al. 2012; Couch 2012). Models based onevolving angular moments of the neutrino distribution func-tion are therefore attractive for 3D simulations, where the so-lution space must be reduced.With an eye toward future multidimensional simula-tions of angle- and energy-dependent neutrino transport,Cardall & Mezzacappa (2003, hereafter CM03) derived con-servative general relativistic kinetic equations intended forimplementation in numerical codes for simulations of neu-trino transport in core-collapse supernovae. They adopted dis-tinct position and momentum space coordinates, by employ-ing spacetime position coordinates associated with a globalcoordinate basis and momentum space coordinates associatedwith an orthonormal comoving basis. The kinetic equationsderived by CM03 describe the evolution of the specific par-ticle number density and four-momentum (six-dimensionalphase space densities). Such conservative equations allowfor convenient treatment of neutrino-matter interactions in thecomoving-frame, and may be particularly useful when con-structing numerical methods for the multidimensional Boltz-mann equation with desirable conservation properties (i.e.,conserves total lepton number and energy). However, Boltz-mann neutrino transport in 3D core-collapse supernova simu-lations will not be possible for some time to come (exascalecomputational resources will be required, even with moderatephase space resolution), and approximate methods will con-tinue to play a dominant role in the foreseeable future.In preparation for development of numerical methods forneutrino radiation transport in our astrophysical simulationcode GenASiS (Cardall et al. 2012a,b), we derive conserva-tive , multidimensional, monochromatic moment equations forneutrino radiation transport from first principles, employingthe general relativistic framework of CM03. We derive a two-moment model, which evolves angular moments of the neu-trino distribution function; i.e., the neutrino energy and mo-mentum (as opposed to only the energy in the flux-limiteddiffusion approach). The radiation moments are functionsof spacetime position components associated with a global coordinate basis x µ and the radiation energy measured by acomoving observer ǫ . First we specialize the moment equa-tions to the case of a conformally flat spacetime. Then wespecialize the equations further to the pseudo-Newtonian (cf.Kim et al. 2012) and the Newtonian gravity, O ( v ) limits.Special relativistic equations are presented in Appendix A,while we present general relativistic moment equations forspherically symmetric spacetimes in Appendix B. (Generalrelativistic moment equations employing the 3+1 formula-tion have been presented separately in Cardall et al. 2012;see also Shibata et al. 2011.) We elucidate the relationshipbetween the moment equations for neutrino four-momentumand neutrino number in general (see also Cardall et al. 2012),and we pay special attention to consistency between the con-servative moment equations for neutrino four-momentum andthe conservative equation for the neutrino number density inthe non-relativistic limits. Detailed knowledge of this rela-tionship may be useful when constructing consistent numer-ical methods for neutrino radiation transport, where leptonnumber conservation (in addition to total energy conserva-tion) must also be considered. Our conservative momentequations correspond to similar non-conservative momentequations presented by other authors (e.g., Buchler 1979,1983, 1986; Kaneko et al. 1984; Munier & Weaver 1986b;M¨uller et al. 2010).The O ( v ) limit of the radiation moment equations is notuniquely defined, and care should be taken to ensure con-sistency between the conservative equations for the radiationfour-momentum and the neutrino number equation (see alsoCardall et al. 2005, for a detailed discussion). In the non-relativistic limit, we find that exact consistency between theconservative moment equations for neutrino four-momentumand the conservative neutrino number equation can be ob-tained by adopting different orders of v for the radiation en-ergy and momentum equations, respectively. As an illustra-tive example, in the laboratory frame the neutrino numberdensity E N and the neutrino energy density E and momentumdensity F i are to O ( v ) related by ǫ E N = E − v i F i , (1)where v i is the fluid three-velocity. The evolution equationsfor the respective quantities are similarly related. By approxi-mating the conservative (lab-frame) radiation energy and mo-mentum equations to O ( v ) and O (1) , respectively, we ob-tain exact consistency with the conservative O ( v ) neutrinonumber equation, and we avoid the introduction of ‘extra’ O ( v ) terms. However, because we use an O (1) momentumequation, this description is formally only accurate to O (1) (we call this the O ( v ) -minus approximation to indicate thatvelocity-dependent terms are included, and to distinguish itfrom ‘traditional’ O (1) descriptions, where no moving fluideffects are included). We obtain a description that is formallyaccurate to O ( v ) by adopting O ( v ) and O ( v ) radiation en-ergy and momentum equations, respectively. This descriptionis exactly consistent with the conservative O ( v ) neuterinonumber equation. We refer to this system of moment equa-tions as the O ( v ) -plus radiation moment equations.We have organized the paper as follows: we briefly summa-rize general relativistic kinetic equations in Section 2. Angu-lar moments of the distribution are defined in Section 3, andwe derive conservative evolution equations for the so-calledzeroth- and first-order angular moments in Section 4. (We de-rive conservative monochromatic equations for the neutrinoONSERVATIVE RADIATION MOMENT EQUATIONS 3number density and four-momentum.) We also elucidate therelationship between the stress-energy equation and the num-ber equation in Section 4. In Section 5 we specialize themoment equations for a sufficiently general (conformally flat)spacetime metric to accommodate the pseudo-Newtonian ra-diation moment equations; the Newtonian gravity, O ( v ) equa-tions; the special relativistic moment equations; and the gen-eral relativistic moment equations for spherically symmetricspacetimes. The pseudo-Newtonian moment equations arepresented in Section 6. We introduce the Newtonian grav-ity, O ( v ) -plus and O ( v ) -minus approximations in Section 7.A summary and discussion is given in Section 8. Special rela-tivistic moment equations are presented in Appendix A. Gen-eral relativistic moment equations for spherically symmetricspacetimes are presented in Appendix B, where we also com-pare the equations with those solved by M¨uller et al. (2010). GENERAL RELATIVISTIC KINETIC EQUATIONS
In this section we briefly summarize general relativistic ki-netic equations, and define the variables associated with theequations (see for example Lindquist 1966; Ehlers 1971;Israel 1972, for detailed treatments of relativistic kinetic the-ory). The kinetic equations form the basis for our derivationof angular moment equations for neutrino radiation transport.To this end, we consider a single species of massless, elec-trically neutral particles (i.e., classical neutrinos). We adopt a‘geometrized’ unit system in which the vacuum speed of light,the gravitational constant, and the Planck constant are unity,and we adopt the usual Einstein summation convention andlet repeated Greek indices run from to , and repeated Latinindices run from to .The classical neutrino distribution function f ( x µ , p ˆ µ ) —thephase space particle density—is governed by the Boltzmannequation. In its geometric form, it states that the change in f along a phase space trajectory, parametrized by λ , equalsthe density C [ f ] of point-like collisions that add or removeparticles from the trajectory; i.e., dfdλ = C [ f ] . (2)For practical computations it is necessary to introduce phasespace coordinates. By introducing spacetime coordinates x µ ,and momentum space coordinates p ˆ ı , Equation (2) becomes dx µ dλ ∂f∂x µ + dp ˆ ı dλ ∂f∂p ˆ ı = C [ f ] . (3)(By the mass shell constraint p ˆ µ p ˆ µ = 0 , only three of thefour-momentum components are independent, and the con-traction over momentum coordinates in Equation (3) onlyruns over three indices.) The geodesic equations describingthe trajectory can be written as dx µ dλ = L µ ˆ µ p ˆ µ , (4) dp ˆ µ dλ = − Γ ˆ µ ˆ ν ˆ ρ p ˆ ν p ˆ ρ . (5)Then, the general relativistic Boltzmann equation reads (e.g.,Lindquist (1966); Mezzacappa & Matzner (1989); CM03) L µ ˆ µ p ˆ µ ∂f∂x µ − Γ ˆ ı ˆ ν ˆ ρ p ˆ ν p ˆ ρ ∂f∂p ˆ ı = C [ f ] . (6)In the above we have used our freedom in choosing dis-tinct spacetime and momentum space coordinates: { x µ } = { t, x i } are spacetime position components in a global coor-dinate (holonomic) basis , while { p ˆ µ } = ǫ { , n ˆ ı } are four-momentum components in a comoving orthonormal basis ,where n ˆ ı are components of the spatial unit vector parallelto p ˆ ı , with basis tangent to the spatial coordinate lines in acomoving orthonormal basis.A ‘comoving-frame’ at a spacetime point is an inertial ref-erence frame whose four-velocity instantaneously coincideswith that of a moving fluid element at that point, and is theframe where interactions between the radiation field and thefluid are most easily described (Mihalas & Mihalas 1999).In the comoving frame, the four-velocity of a comoving ob-server is { u ˆ µ } = { , , , } . Quantities with unaccentedindices are defined with respect to the global coordinate ba-sis, while quantities whose indices are accented with a hat(comoving-frame quantities) are defined with respect to theorthonormal comoving basis. The composite transformation L µ ˆ µ ≡ e µ ¯ µ Λ ¯ µ ˆ µ transforms four-vectors defined with respectto the orthonormal comoving basis into four-vectors definedwith respect to the global coordinate basis. Thus, the four-velocity of a comoving observer with respect to the globalcoordinate basis is u µ = L µ ˆ µ u ˆ µ . The composite transfor-mation consists of a Lorentz transformation Λ ¯ µ ˆ µ from the or-thonormal comoving basis (comoving-frame) to an in generalnoncomoving orthonormal tetrad basis (the orthonormal lab-frame; e.g., u ¯ µ = Λ ¯ µ ˆ µ u ˆ µ ), followed by a transformation e µ ¯ µ from the orthonormal tetrad basis to the global coordinatebasis (e.g., u µ = e µ ¯ µ u ¯ µ ). Quantities with indices accentedwith a bar (lab-frame quantities) are defined with respect tothe orthonormal tetrad basis. The “tetrad transformation” lo-cally transforms the metric tensor g µν into the Minkowskimetric tensor; e µ ¯ µ e ν ¯ ν g µν = η ¯ µ ¯ ν ≡ diag [ − , , , . (TheMinkowski metric is preserved by the Lorentz transforma-tion, and the composite transformation is also a tetrad trans-formation.) The inverse composite transformation is denoted L ˆ µµ ≡ Λ ˆ µ ¯ µ e ¯ µµ ; i.e., L µ ˆ µ L ˆ µν = δ µν , where δ µν is the Kroneckerdelta.The ‘connection coefficients’ associated with the orthonor-mal comoving basis Γ ˆ µ ˆ ν ˆ ρ (cf. Equation (6)) are expressed interms of the connection coefficients associated with the globalspacetime coordinate basis as Γ ˆ µ ˆ ν ˆ ρ = L ˆ µµ L ν ˆ ν L ρ ˆ ρ Γ µνρ + L ˆ µµ L ρ ˆ ρ ∂L µ ˆ ν ∂x ρ (7)where the coordinate basis connection coefficients are givenin terms of the spacetime metric g µν as Γ µνρ = 12 g µσ (cid:16) ∂g νσ ∂x ρ + ∂g σρ ∂x ν − ∂g νρ ∂x σ (cid:17) , (8)and g µν is the contravariant metric tensor; g µν g νρ = δ µρ .Note now that the factors outside the momentum space deriva-tive of f in the second term on the left-hand side of Equation(6) are proportional to space and time derivatives of the fluidthree-velocity and components of the metric tensor. Theseterms account for Doppler and Einstein shifts. See for ex-ample Lentz et al. (2012a) for a general discussion of theseterms and a demonstration of their importance in simulationsof neutrino transport in core-collapse supernovae.The collision operator C [ f ] describes the momentum spaceevolution of the distribution function f due to neutrino-matterinteractions. Neutrino-matter interactions crucial to realistic Endeve et al.simulations of CCSNe include particle creation and destruc-tion via emission and absorption, scattering (both elastic andinelastic), and neutrino pair creation and annihilation. We willnot discuss neutrino-matter interactions in any detail in thispaper. See, however, Bruenn (1985); Burrows & Thompson(2004), and the references therein for general details. See alsoLentz et al. (2012b) for a recent discussion of neutrino opac-ities in core-collapse supernova simulations.Our choice of phase space coordinates is well suited fornumerical solution of the Boltzmann equation (or equationsbased on angular moments of the Boltzmann equation). InEquation (6), the particle distribution function is parametrizedin terms of spacetime position components in a global coordi-nate basis { x µ } and momentum components in an orthonor-mal basis comoving with the fluid { p ˆ ı } . This specific choiceis motivated by our intent to develop numerical methods forcomputer simulations of neutrino transport. Neutrino-matterinteractions are on the one hand most easily handled computa-tionally in the frame comoving with the fluid, where materialproperties are isotropic (e.g., Mihalas & Mihalas 1999). Onthe other hand, when integrated over momentum space, theBoltzmann equation expresses conservation of particle num-ber and energy in the laboratory frame. Accurate account-ing of total lepton number and energy in numerical simula-tions of core-collapse supernovae is important and extremelychallenging. Convenient treatment of neutrino-matter interac-tions and global conservation are naturally expressed in thisphase space coordinate basis. However, Equation (6) is notin conservative form, and, from a practical standpoint, maynot be the best starting point for designing numerical meth-ods to be implemented in simulation codes. Numerical meth-ods based on the solution of the Boltzmann equation as it isexpressed in Equation (6) will in general not conserve leptonnumber or energy. Even for numerical methods based on anumber-conservative formulation of the Boltzmann equation,care must be taken in the discretization process in order toensure total lepton number and energy conservation withinacceptable limits, so that firm conclusions can be drawn fromthe simulations (Liebend¨orfer et al. 2004).To help facilitate the realization of lepton number and en-ergy conservation in multidimensional simulations of neu-trino transport, CM03 derived conservative general relativis-tic formulations of kinetic theory. In particular, their number-conservative reformulation of Equation (6) (cf. Eq. (168) inCM03) reads S N [ f ] + M N [ f ] = C [ f ] , (9)where the spacetime divergence is S N [ f ] = 1 √− g ∂∂x µ (cid:16) √− g L µ ˆ µ p ˆ µ f (cid:17) , (10)and the momentum space divergence is M N [ f ] = −| p | (cid:12)(cid:12)(cid:12) det h ∂ p ∂ u i(cid:12)(cid:12)(cid:12) − ×× ∂∂u ˆ ı (cid:16) | p | (cid:12)(cid:12)(cid:12) det h ∂ p ∂ u i(cid:12)(cid:12)(cid:12) Γ ˆ ˆ µ ˆ ν ∂u ˆ ı ∂p ˆ p ˆ µ p ˆ ν f (cid:17) . (11)The determinant of the metric tensor g µν is denoted g .Similarly, the four-momentum-conservative reformulation ofEquation (6) (cf. Eq. (169) in CM03) reads S µ T [ f ] + M µ T [ f ] = L µ ˆ µ p ˆ µ C [ f ] , (12) where the spacetime divergence is S µ T [ f ] = 1 √− g ∂∂x ν (cid:16) √− g L µ ˆ µ L ν ˆ ν p ˆ µ p ˆ ν f (cid:17) +Γ µρν L ρ ˆ ρ L ν ˆ ν p ˆ ρ p ˆ ν f, (13)and the momentum space divergence is M µ T [ f ] = −| p | (cid:12)(cid:12)(cid:12) det h ∂ p ∂ u i(cid:12)(cid:12)(cid:12) − ×× ∂∂u ˆ ı (cid:16) | p | (cid:12)(cid:12)(cid:12) det h ∂ p ∂ u i(cid:12)(cid:12)(cid:12) Γ ˆ ˆ ν ˆ ρ ∂u ˆ ı ∂p ˆ L µ ˆ µ p ˆ µ p ˆ ν p ˆ ρ f (cid:17) . (14)Equations (9) and (12) are conservative in the sense that,when integrated over momentum space, the momentum spacederivatives vanish, and the resulting equations—expressingnumber and four-momentum balance, respectively—are fa-miliar from position space conservation law theory (seeCM03, for a detailed discussion). Note that Equations (9)and (12) are merely conservative (nontrivial) reformulationsof Equation (6). Thus, Equations (9) and (12) are, of course,not independent. While a numerical method is based on solv-ing one of the equations, the numerical solution should ideallyremain consistent with both. A formidable challenge for fu-ture efforts will be to construct a discrete representation of themultidimensional Boltzmann equation (Equation (9) or Equa-tion (12)) that is both number and energy conservative (i.e., si-multaneously consistent with Eqs. (9) and (12))—as was doneby Liebend¨orfer et al. (2004) in the spherically symmetriccase. The conservative formulations provided by Equations(9) and (12) may be helpful for this task. Numerical meth-ods based on solving for angular moments of the distributionfunction should also retain this “dual consistency.” In this pa-per we will discuss requirements and prospects for construct-ing numerical methods for neutrino transport in CCSNe basedon moment models that simultaneously conserve total leptonnumber and total energy.In Equations (9) and (12), a change to spherical momen-tum space coordinates { u ˆ ı } = { ǫ, ϑ, ϕ } has been facili-tated. The Cartesian momentum components are expressedin terms of the spherical momentum space coordinates by { p ˆ ı } = ǫ { n ˆ1 , n ˆ2 , n ˆ3 } = ǫ { cos ϑ, sin ϑ cos ϕ, sin ϑ sin ϕ } .The Jacobian matrix associated with the transformation is ∂p ˆ ı ∂u ˆ = cos ϑ − ǫ sin ϑ ϑ cos ϕ ǫ cos ϑ cos ϕ − ǫ sin ϑ sin ϕ sin ϑ sin ϕ ǫ cos ϑ sin ϕ ǫ sin ϑ cos ϕ ! (15)whose inverse is ∂u ˆ ı ∂p ˆ = 1 ǫ ǫ cos ϑ ǫ sin ϑ cos ϕ ǫ sin ϑ sin ϕ − sin ϑ cos ϑ cos ϕ cos ϑ sin ϕ − sin ϕ/ sin ϑ cos ϕ/ sin ϑ ! . (16)The particle energy measured by a comoving observer (as-suming massless neutrinos) is | p | = ǫ , and the determinant ofthe Jacobian matrix ( ∂p ˆ ı /∂u ˆ ) is (cid:12)(cid:12)(cid:12) det h ∂ p ∂ u i(cid:12)(cid:12)(cid:12) = ǫ sin ϑ. (17)This concludes our summary of general relativistic kineticequations. ANGULAR MOMENTS OF THE DISTRIBUTION FUNCTION
ONSERVATIVE RADIATION MOMENT EQUATIONS 5The distribution function f ( x µ , p ˆ ı ) provides a detailed sta-tistical description of the particle momentum at a given space-time location. However, as a seven-dimensional object itis in general computationally prohibitive to compute withsufficient phase space resolution in multidimensional super-nova simulations using currently available computer hard-ware. The solution space must be reduced. Instead of solv-ing the Boltzmann equation directly for the particle distri-bution function, solving equations for moments of the dis-tribution function has become a popular method of reducingthe dimensionality of the transport problem to a computation-ally tractable one (e.g., Lindquist 1966; Anderson & Spiegel1972; Castor 1972; Thorne 1981; Munier & Weaver 1986b;Mihalas & Mihalas 1999). Variants of so-called momentmodels (with various degrees of sophistication) have beenused extensively for simulation of neutrino transport in core-collapse supernovae (e.g., Bruenn 1985; Herant et al. 1994;Rampp & Janka 2002; Bruenn et al. 2006; Buras et al.2006; Burrows et al. 2006, 2007; Bruenn et al. 2009;Marek & Janka 2009; Swesty & Myra 2009; M¨uller et al.2010; Obergaulinger & Janka 2011; Bruenn et al. 2012).Energy-dependent, general relativistic angular moment equa-tions are presented in this paper (see also M¨uller et al. 2010;Shibata et al. 2011; Cardall et al. 2012).The neutrino heating rate in the gain region (i.e., the re-gion behind the shock where the neutrino heating rate exceedsthe neutrino cooling rate) depends quadratically on the neu-trino energy spectrum (e.g., Mezzacappa 2005), and mainlyfor this reason do we retain the energy dependence, while weintegrate over the momentum directions. The dimensionalityof the transport problem is then reduced to five (four space-time dimensions, and one momentum space dimension). Inparticular, we wish to derive equations for the lab-frame neu-trino number density and four-momentum density expressedin terms of angular moments of the distribution function de-fined in the comoving-frame. The invariant momentum spacevolume element is d p | p | = (cid:12)(cid:12)(cid:12) det (cid:16) ∂ p ∂ u (cid:17)(cid:12)(cid:12)(cid:12) d u | p | = ǫ sin ϑ dǫ dϑ dϕ. (18)It is expressed in terms of spherical comoving-frame momen-tum space coordinates on the far right side.The number-flux four-vector (integrated over the entire mo-mentum space V p ) is defined in terms of moments of the dis-tribution function as (e.g., Lindquist 1966) N µ = Z V p p µ f d p | p | = L µ ˆ µ Z ∞ N ˆ µ ǫ dǫ, (19)where the monochromatic number-flux four-vector in thecomoving-frame is defined in terms of angular moments ofthe distribution function as N ˆ µ = 1 ǫ Z Ω p ˆ µ f d Ω . (20)The solid angle element is d Ω = sin ϑ dϑ dϕ , and the an-gular integration is carried out over the unit sphere; i.e., R Ω . . . d Ω = R π R π . . . sin ϑ dϑ dϕ . The monochromaticlab-frame number-flux four-vector is related to the corre-sponding comoving-frame four-vector by N µ = L µ ˆ µ N ˆ µ = L µ ˆ0 J N + L µ ˆ ı H ˆ ı N , (21) where the comoving-frame neutrino number density and num-ber flux density are denoted J N and H ˆ ı N , respectively. We usethe calligraphic font to distinguish energy-dependent radia-tion quantities from energy-integrated (grey) radiation quan-tities (e.g., N µ vs. N µ ). Note that what we here call thelab-frame number-flux four-vector is a function of the co-ordinate basis spacetime position components x µ and theneutrino energy measured by a comoving observer ǫ ; i.e., N µ = N µ ( x µ , ǫ ) . This parametrization of radiation trans-port variables was suggested by CM03 (see also Riffert 1986;Mezzacappa & Matzner 1989). This is also the parametriza-tion used by Shibata et al. (2011) and Cardall et al. (2012).In a similar manner, the momentum space integrated sym-metric stress-energy tensor is defined in terms of moments ofthe distribution function as (e.g., Lindquist 1966) T µν = Z V p p µ p ν f d p | p | = L µ ˆ µ L ν ˆ ν Z ∞ T ˆ µ ˆ ν ǫ dǫ, (22)where the monochromatic stress-energy tensor in thecomoving-frame is defined in terms of angular moments ofthe distribution function as T ˆ µ ˆ ν = 1 ǫ Z Ω p ˆ µ p ˆ ν f d Ω . (23)The monochromatic lab-frame stress-energy tensor T µν = L µ ˆ µ L ν ˆ ν T ˆ µ ˆ ν is also a function of the neutrino energy mea-sured by a comoving observer; i.e., T µν = T µν ( x µ , ǫ ) . Fur-thermore, we define the components of the monochromaticcomoving-frame stress-energy tensor as (cid:18) T ˆ0ˆ0 T ˆ0ˆ T ˆ ı ˆ0 T ˆ ı ˆ (cid:19) = (cid:18) J H ˆ H ˆ ı K ˆ ı ˆ (cid:19) (24)(i.e., in terms of the comoving-frame radiation energy density J , momentum density H ˆ ı , and stress K ˆ ı ˆ , respectively). Interms of the monochromatic energy density, momentum den-sity, and stress defined in Equation (24), the monochromaticlab-frame stress-energy tensor can now be written as T µν = L µ ˆ0 L ν ˆ0 J + (cid:16) L µ ˆ ı L ν ˆ0 + L ν ˆ ı L µ ˆ0 (cid:17) H ˆ ı + L µ ˆ ı L ν ˆ K ˆ ı ˆ . (25)We can also write K ˆ ı ˆ = k ˆ ı ˆ J , and define the monochromaticsymmetric rank-two “variable Eddington tensor” as k ˆ ı ˆ = R Ω n ˆ ı n ˆ f d Ω R Ω f d Ω . (26)For the moment equations we derive in the next section itis useful to define the completely symmetric third-order mo-ment U µνρ = Z V p p µ p ν p ρ f d p | p | (27) = L µ ˆ µ L ν ˆ ν L ρ ˆ ρ Z ∞ U ˆ µ ˆ ν ˆ ρ ǫ dǫ, (28)where the monochromatic third-order comoving-frame mo-ments are given by U ˆ µ ˆ ν ˆ ρ = 1 ǫ Z Ω p ˆ µ p ˆ ν p ˆ ρ f d Ω , (29) Endeve et al.with U µνρ = L µ ˆ µ L ν ˆ ν L ρ ˆ ρ U ˆ µ ˆ ν ˆ ρ . Note in particular that U ˆ0ˆ µ ˆ ν = U ˆ µ ˆ0ˆ ν = U ˆ µ ˆ ν ˆ0 = ǫ T ˆ µ ˆ ν , (30)and that we can expand the monochromatic third-order lab-frame moments in terms of the comoving-frame moments as ǫ U µνρ = J L µ ˆ0 L ν ˆ0 L ρ ˆ0 + L µ ˆ ı H ˆ ı L ν ˆ0 L ρ ˆ0 + L ν ˆ ı H ˆ ı L µ ˆ0 L ρ ˆ0 + L ρ ˆ ı H ˆ ı L µ ˆ0 L ν ˆ0 + L µ ˆ ı L ν ˆ K ˆ ı ˆ L ρ ˆ0 + L µ ˆ ı L ρ ˆ K ˆ ı ˆ L ν ˆ0 + L ν ˆ ı L ρ ˆ K ˆ ı ˆ L µ ˆ0 + L µ ˆ ı L ν ˆ L ρ ˆ k L ˆ ı ˆ ˆ k , (31)where L ˆ ı ˆ ˆ k = l ˆ ı ˆ ˆ k J , and we have introduced the completelysymmetric rank-three variable Eddington tensor l ˆ ı ˆ ˆ k = R Ω n ˆ ı n ˆ n ˆ k f d Ω R Ω f d Ω . (32)We have defined the lab-frame number-flux four-vector, thelab-frame stress-energy tensor, and the third-order lab-framemoments (a rank-three tensor). They are expressed in terms ofangular moments of the distribution function in the comoving-frame (cf. Equations (21), (25), and (31)). From the defini-tions in Equations (20) and (23), we clearly have ǫ N ˆ µ = T ˆ µ ˆ0 . This expression can be generalized to a useful covariantexpression relating the number-flux four-vector and the stress-energy tensor by noting that T ˆ µ ˆ0 = − u ˆ0 T ˆ µ ˆ0 = − u ˆ ν T ˆ µ ˆ ν (Cardall et al. 2012). Thus, we have N µ = − ǫ u ν T µν , (33)which is a covariant expression (valid in any frame; although ǫ is always the neutrino energy measured by a comoving ob-server). Similarly, we have T µν = − ǫ u ρ U µνρ , (34)which is also a covariant expression relating the stress-energytensor and the third-order tensor. GENERAL RELATIVISTIC ANGULAR MOMENT EQUATIONS
We derive general relativistic angular moment equations inthis section. The main results are the number-conservativemonochromatic lab-frame number equation (Equation (41)),which is based on Equation (9), and the four-momentum-conservative monochromatic lab-frame stress-energy equa-tion (Equation (47)), which is based on Equation (12). Thesteps are straightforward.
Number Equation
From the definition of the monochromatic number-fluxfour-vector in Equations (20) and (21), and the expression in-side the spacetime derivative in Equation (9), it is clear that weobtain the evolution equation for the lab-frame neutrino num-ber density by integrating Equation (9) over the unit sphere Ω and dividing by ǫ . For the spacetime divergence we find √− g ∂∂x µ (cid:16) √− g N µ (cid:17) . (35)With spherical momentum space coordinates, the momen- tum space divergence in Equation (9) becomes − ǫ sin ϑ ∂∂u ˆ ı (cid:16) ǫ sin ϑ Γ ˆ ˆ µ ˆ ν ∂u ˆ ı ∂p ˆ p ˆ µ p ˆ ν f (cid:17) = − ǫ ∂∂ǫ (cid:16) ǫ Γ ˆ ˆ µ ˆ ν ∂u ˆ1 ∂p ˆ p ˆ µ p ˆ ν f (cid:17) − ϑ ∂∂ϑ (cid:16) sin ϑ Γ ˆ ˆ µ ˆ ν ∂u ˆ2 ∂p ˆ p ˆ µ p ˆ ν f (cid:17) − ∂∂ϕ (cid:16) Γ ˆ ˆ µ ˆ ν ∂u ˆ3 ∂p ˆ p ˆ µ p ˆ ν f (cid:17) . (36)We only focus on the first term (containing the derivative withrespect to energy ǫ ) on the right-hand side of Equation (36).The second and third terms on the right-hand side of Equa-tion (36) vanish upon integration over momentum space an-gles ( ϑ , ϕ ). In particular, integrating over angles and dividingby ǫ , the momentum space divergence becomes − ǫ ∂∂ǫ (cid:16) ǫ Γ ˆ0ˆ µ ˆ ν Z Ω p ˆ µ p ˆ ν f d Ω (cid:17) , (37)where we have used the fact that ǫ Γ ˆ ı ˆ µ ˆ ν ∂u ˆ1 ∂p ˆ ı p ˆ µ p ˆ ν = Γ ˆ ı ˆ µ ˆ ν p ˆ ı p ˆ µ p ˆ ν = ǫ Γ ˆ0ˆ µ ˆ ν p ˆ µ p ˆ ν (38)(cf. Equation (16)). The rightmost expression in Equation(38) is derived from the mass-shell constraint ( p ˆ µ p ˆ µ = 0 ;i.e., p ˆ0 is considered a function of the independent momen-tum space coordinates { p ˆ ı } ) and Equation (5); cf. Equation(136) in CM03. Then, writing the expression inside the en-ergy derivative in Equation (37) in terms of comoving-frameangular moments of the distribution function, we have ǫ Γ ˆ0ˆ µ ˆ ν Z Ω p ˆ µ p ˆ ν f d Ω = ǫ Γ ˆ0ˆ µ ˆ ν T ˆ µ ˆ ν = − ǫ T µν ∇ µ L ˆ0 ν = ǫ T µν ∇ µ u ν , (39)where we have eliminated the connection coefficients associ-ated with the orthonormal comoving basis with Equation (7),and introduced the (coordinate basis) four-velocity of a co-moming observer u ν = L ˆ νν u ˆ ν = − L ˆ0 ν (cf. Section 2). Thecovariant derivative of a vector A ν is (cf. Landau & Lifshitz1975) ∇ µ A ν = ∂A ν ∂x µ − Γ ρνµ A ρ . (40)Thus, the number-conservative monochromatic lab-framenumber equation, expressed explicitly in terms of angular mo-ments of the particle distribution function in the comoving-frame (via Equations (21) and (25)), becomes √− g ∂∂x µ (cid:16) √− g N µ (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ T µν ∇ µ u ν (cid:17) = 1 ǫ Z Ω C [ f ] d Ω . (41)Indeed, an integration over energy shells results in √− g ∂∂x µ (cid:16) √− g N µ (cid:17) = Z V p C [ f ] d p | p | (42)(cf. Lindquist 1966).ONSERVATIVE RADIATION MOMENT EQUATIONS 7 Stress-Energy Equation
Similarly, from the definition of the stress-energy tensor inEquation (23) and the expression inside the spacetime deriva-tive in Equation (12), it is clear that we obtain the evolu-tion equation for the monochromatic lab-frame neutrino four-momentum by integrating Equation (12) over the full solid an-gle Ω and dividing by ǫ . The spacetime divergence becomes √− g ∂∂x ν (cid:16) √− g T µν (cid:17) + Γ µρν T ρν . (43)The momentum space divergence in Equation (12) is − ǫ sin ϑ ∂∂u ˆ ı (cid:16) ǫ sin ϑ Γ ˆ ˆ ν ˆ ρ ∂u ˆ ı ∂p ˆ L µ ˆ µ p ˆ µ p ˆ ν p ˆ ρ f (cid:17) = − ǫ ∂∂ǫ (cid:16) ǫ Γ ˆ ˆ ν ˆ ρ ∂u ˆ1 ∂p ˆ L µ ˆ µ p ˆ µ p ˆ ν p ˆ ρ f (cid:17) − ϑ ∂∂ϑ (cid:16) sin ϑ Γ ˆ ˆ ν ˆ ρ ∂u ˆ2 ∂p ˆ L µ ˆ µ p ˆ µ p ˆ ν p ˆ ρ f (cid:17) − ∂∂ϕ (cid:16) Γ ˆ ˆ ν ˆ ρ ∂u ˆ3 ∂p ˆ L µ ˆ µ p ˆ µ p ˆ ν p ˆ ρ f (cid:17) . (44)Again, we focus only on the term with the energy derivativeon the right-hand side of Equation (44). After integrating overthe unit sphere and dividing by ǫ , the momentum space diver-gence reduces to − ǫ ∂∂ǫ (cid:16) ǫ L µ ˆ µ Γ ˆ0ˆ ν ˆ ρ Z Ω p ˆ µ p ˆ ν p ˆ ρ f d Ω (cid:17) , (45)where the relations in Equation (38) have been used. Weseek to express the momentum space divergence in terms ofcomoving-frame angular moments of the distribution func-tion, and write ǫ L µ ˆ µ Γ ˆ0ˆ ν ˆ ρ Z Ω p ˆ µ p ˆ ν p ˆ ρ f d Ω = ǫ L µ ˆ µ Γ ˆ0ˆ ν ˆ ρ U ˆ µ ˆ ν ˆ ρ = − ǫ U µνρ ∇ ν L ˆ0 ρ = ǫ U µνρ ∇ ν u ρ , (46)where we have taken steps similar to those used to arrive atEquation (39).Then, combining all the terms, the monochromaticlab-frame stress-energy equation, expressed in terms ofcomoving-frame angular moments of the distribution function(via Equations (25) and (31)), becomes √− g ∂∂x ν (cid:16) √− g T µν (cid:17) + Γ µρν T ρν − ǫ ∂∂ǫ (cid:16) ǫ U µνρ ∇ ν u ρ (cid:17) = 1 ǫ Z Ω p µ C [ f ] d Ω . (47)Equation (47) corresponds to Equation (3.18) in Shibata et al.(2011), which was derived from the moment formalism ofThorne (1981). Integrating Equation (47) over comoving-frame energy shells gives the familiar result √− g ∂∂x ν (cid:16) √− g T µν (cid:17) +Γ µρν T ρν = Z V p p µ C [ f ] d p | p | (48)(cf. Lindquist 1966).By invoking the transformation from the orthonormal tetradbasis to the coordinate basis e µ ¯ µ we can rewrite Equation (47) in an equivalent form in terms of the stress-energy tensor as-sociated with the orthonormal tetrad basis √− g ∂∂x ν (cid:16) √− g e ν ¯ ν T ¯ µ ¯ ν (cid:17) + Γ ¯ µ ¯ ρ ¯ ν T ¯ ρ ¯ ν − ǫ ∂∂ǫ (cid:16) ǫ U ¯ µ ¯ ν ¯ ρ ∇ ¯ ν u ¯ ρ (cid:17) = 1 ǫ Z Ω p ¯ µ C [ f ] d Ω , (49)where the stress-energy tensor T ¯ µ ¯ ν is related to thecomoving-frame angular moments through the Lorentz trans-formation (Appendix A), the covariant derivative of a four-vector A ¯ ρ is ∇ ¯ ν A ¯ ρ = e ν ¯ ν ∂A ¯ ρ ∂x ν − Γ ¯ µ ¯ ρ ¯ ν A ¯ µ , (50)and the connection coefficients associated with the orthonor-mal tetrad basis are Γ ¯ µ ¯ ρ ¯ ν = e ¯ µµ e ρ ¯ ρ e ν ¯ ν Γ µρν + e ¯ µµ e ν ¯ ν ∂e µ ¯ ρ ∂x ν . (51)Going even one step further, by invoking the Lorentz trans-formation, we can rewrite Equation (49) in an equivalent formin terms of the stress-energy tensor associated with the or-thonormal comoving basis √− g ∂∂x ν (cid:16) √− g L ν ˆ ν T ˆ µ ˆ ν (cid:17) + Γ ˆ µ ˆ ρ ˆ ν T ˆ ρ ˆ ν − ǫ ∂∂ǫ (cid:16) ǫ U ˆ µ ˆ ν ˆ ρ ∇ ˆ ν u ˆ ρ (cid:17) = 1 ǫ Z Ω p ˆ µ C [ f ] d Ω , (52)where the covariant derivative of a four-vector A ˆ ρ is ∇ ˆ ν A ˆ ρ = L ν ˆ ν ∂A ˆ ρ ∂x ν − Γ ˆ µ ˆ ρ ˆ ν A ˆ µ , (53)and the connection coefficients associated with the orthonor-mal comoving basis are given in terms of connection coeffi-cients associated with the orthonormal tetrad basis by Γ ˆ µ ˆ ρ ˆ ν = Λ ˆ µ ¯ µ Λ ¯ ρ ˆ ρ Λ ¯ ν ˆ ν Γ ¯ µ ¯ ρ ¯ ν + Λ ˆ µ ¯ µ Λ ¯ ν ˆ ν e ν ¯ ν ∂ Λ ¯ µ ˆ ρ ∂x ν , (54)or in terms of the connection coefficients associated with thecoordinate basis in Equation (7).The covariant nature of the general relativistic momentequations is illustrated by Equations (47), (49), and (52). Al-though they are analytically equivalent, Equation (47) maybe preferred over Equations (49) and (52) as the basis fordeveloping a numerical method to evolve the neutrino radia-tion field in full general relativity employing the so-called 3+1decomposition (Shibata et al. 2011; Cardall et al. 2012). InSection 5 we use Equation (47) to obtain moment equa-tions valid for conformally flat spacetimes. Equations (47)and (49) are conservative equations for the lab-frame four-momentum. In the absence of gravity and neutrino-matterinteractions, Equations (47) and (49) express exact conser-vation of neutrino energy and —if Cartesian coordinates areused—momentum (Appendix A). Equation (52) is expressedin terms of the comoving-frame stress-energy tensor. Re-cently, M¨uller et al. (2010) presented numerical methodsfor neutrino radiation transport based on equations that areclosely related to Equation (52) (see Appendix B). Equa-tion (52) is nonconservative—even in the absence of gravityand neutrino-matter interactions—since the connection coef-ficients Γ ˆ µ ˆ ν ˆ ρ depend on time and space derivatives of the fluid Endeve et al.three-velocity (cf. Equation (54)). In Appendix B we useEquation (47) to obtain conservative general relativistic mo-ment equations valid for spherically symmetric spacetimes(i.e., similar to M¨uller et al. 2010, but in conservative form). The Relationship Between the Number Equation and theStress-Energy Equation
In this subsection we illustrate the relationship between themoment equations for the neutrino four-momentum densityand the moment equation for the neutrino number density;Equations (47) and (41), respectively. From Equation (33),it is evident how the neutrino number equation is related tothe stress-energy equation. Thus, contracting − ǫ − u µ withEquation (47) results in √− g ∂∂x µ (cid:16) √− g N µ (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ T µν ∇ µ u ν (cid:17) + 1 ǫ (cid:16) T νρ + 1 ǫ u µ U µνρ (cid:17) ∇ ν u ρ = 1 ǫ Z Ω C [ f ] d Ω . (55)To arrive at Equation (55), we have used Equations (33) and(34) in the first two terms on the left-hand side, and the co-variant relation u µ p µ = u ˆ µ p ˆ µ = − ǫ to obtain the colli-sion term on the right-hand side. The third term on the left-hand side of Equation (55) consists of the ‘leftover’ terms af-ter bringing − ǫ − u µ inside the spacetime and energy deriva-tives, and vanishes exactly by virtue of Equation (34). Thus,as expected, Equation (55)—the contraction of − ǫ − u µ withEquation (47)—reduces to Equation (41); i.e., the number-conservative monochromatic number equation. The solutionto Equation (47) can thus be used to construct the solution toEquation (41) via Equation (33).When carrying out the steps to obtain the number-conservative monochromatic neutrino number equation fromthe four-momentum-conservative monochromatic stress-energy equation, we observe that terms emanating from thespacetime divergences and the geometry sources cancel withterms emanating from the momentum space divergences. Theremaining terms constitute the left-hand side of the numberequation. Such detailed knowledge about how the conserva-tive number equation is ‘built into’ the four-momentum con-servative stress-energy equation may be useful when develop-ing numerical methods. In particular, when discretizing thestress-energy equation to develop a two-moment model forneutrino transport, care should be taken to mimic the cancel-lations that occur in the continuum limit in order to ensurethat the resulting numerical solution is also consistent witha discrete version of the conservative number equation, andthereby facilitating lepton number conservation as a ‘built-in’property of the numerical method. We will elaborate furtheron this issue with specific examples in Sections 5 and 7, andin Appendix A.We have derived conservative general relativistic momentequations for the number density N (Equation (41)) andthe four-momentum density T µ (e.g., Equation (47)). Theequations are in conservative form in the sense that (moduloneutrino-matter interactions and the geometry source terms inthe stress-energy equation) the time rates of change of N and T µ are governed by space and momentum space diver-gences. When integrated over comoving-frame energy bins(with ǫ dǫ as the integration measure), the equations reduceto familiar position space conservation laws (Equations (42)and (48), respectively). The equations govern the evolutionof quantities defined with respect to the global coordinate basis, which are explicitly expressed in terms of comoving-frame angular moments, N ˆ µ and T ˆ µ ˆ ν via Equations (21) and(25). Such equations may be suitable for implementation ofneutrino transport capabilities in numerical codes intendedfor simulations of core-collapse supernovae, where accept-able conservation of total energy and lepton number is critical.We have also illustrated the relationship between the numberequation and the stress-energy equation.In terms of a hierarchy of moment equations, we havetruncated the hierarchy at the level of the so-called first-order moment (the comoving-frame energy flux density H ˆ ı ;the comoving-frame energy density J is proportional to thezeroth-order moment), which results in a two-moment modelfor neutrino transport. When truncating at this particularlevel, the resulting moment equations involve higher-ordermoments (or ratios of higher-order moments to the zerothorder moment; i.e., the variable Eddington tensors) k ˆ ı ˆ and l ˆ ı ˆ ˆ k , which must be determined in order to obtain a closedsystem of equations. Due to symmetry, six unique com-ponents of k ˆ ı ˆ and ten unique components of l ˆ ı ˆ ˆ k must bedetermined (sixteen in total) per energy group. The pro-cedure to determine these higher-order moments in a clo-sure prescription that relates the higher moments to the firsttwo is referred to as the moment closure problem. (Wedo not address the moment closure problem in this paper.)Common approaches to the moment closure problem includeanalytic closure (e.g., Minerbo 1978; Levermore 1984),entropy-based closure (e.g., Cernohorsky & Bludman 1994;Levermore 1996; Smit et al. 2000; Brunner & Holloway2001; Hauck & McClarren 2010), and closure based on thesolution of an approximate (simplified) Boltzmann equation(e.g., Rampp & Janka 2002; M¨uller et al. 2010). MOMENT EQUATIONS FOR CONFORMALLY FLAT SPACETIMES
The moment equations derived in Section 4 are valid forgeneral spacetime coordinate systems. In this section wepresent moment equations suitable for simulations assum-ing a simplified spacetime metric. In particular, we imposethe conformal flatness condition on the spatial metric (cf.Wilson et al. 1996; Isenberg 2008). The equations derivedin this section are sufficiently general to accommodate thepseudo-Newtonian moment equations (Section 6), the New-tonian gravity, O ( v ) approximation of the radiation momentequations (Section 7), the special relativistic limit (AppendixA), and the fully general relativistic case for spherically sym-metric spacetimes (Appendix B). (See Shibata et al. 2011;Cardall et al. 2012, for moment equations for radiation trans-port employing the full 3+1 formulation of general relativity.)To this end, we adopt the 3+1 form of the spacetime metric g µν = (cid:18) − α + β k β k β j β i γ ij (cid:19) , (56)where α , β i , and γ ij are the lapse function, the shift vec-tor, and the spatial three-metric, respectively. The confor-mal flatness condition is imposed by performing a confor-mal transformation of the three-metric, γ ij = ψ ¯ γ ij , where ψ ( x µ ) is the conformal factor, and insisting that the con-formally related metric ¯ γ ij is diagonal. In particular, weset ¯ γ ij = diag (cid:2) , a ( x ) , b ( x ) c ( x ) (cid:3) , where the metricfunctions (or scale factors) a , b , and c are sufficiently generalto accommodate Cartesian, spherical, or cylindrical coordi-nates (see for example Table I in Cardall et al. 2005). TheONSERVATIVE RADIATION MOMENT EQUATIONS 9scale factors, and hence the conformally related metric, maydepend on the spatial position coordinate, but are taken hereto be independent of time . The lapse function, the shift vector,and the conformal factor completely determine the spacetimemetric, and can be obtained by solving a system of nonlinearelliptic equations (Wilson et al. 1996). The four-metric andits inverse lower and raise indices on four-vectors and tensors,while the spatial metric γ ij and its inverse γ ij can be used tolower and raise spatial indices of purely spatial vectors andtensors. The determinant of the spacetime metric is denoted g = − α γ , where γ = ψ ( abc ) is the determinant of thespatial metric.In the 3+1 formulation of general relativity, the four-dimensional spacetime is sliced into a “stack” of spatial(spacelike) hypersurfaces of constant (global) time coordinate(see for example the recent book by Baumgarte & Shapiro2010). In particular, the unit normal to the spacelike hyper-surfaces is denoted n µ = (cid:0) α − , − α − β i (cid:1) . With the metricspecified in Equation (56) we have n µ = g µν n ν = ( − α, .Observers whose four-velocity coincides with n µ are at restwith respect to the spacelike slice, and are referred to as Eule-rian observers (e.g., Banyuls et al. 1997). For the transforma-tion from the orthonormal tetrad basis to the global coordinatebasis we can use e µ ¯ µ = (cid:18) α − − α − β i e i ¯ ı (cid:19) , (57)where e i ¯ ı = ψ − diag (cid:2) , a − , b − c − (cid:3) (i.e., e µ ¯ µ e ν ¯ ν g µν = η ¯ µ ¯ ν ). With e µ ¯ µ given by Equation (57), we have e µ ¯0 n µ = − .That is, the spacetime orientation of the orthonormal tetradbasis has been chosen so that e µ ¯0 is aligned with the unit nor-mal of the spacelike hypersurfaces (i.e., e µ ¯0 = n µ ). The four-velocity of a comoving observer with respect to the lab-framecoodinate basis can now be computed. The result is u µ = L µ ˆ µ u ˆ µ = W ( n µ + v µ ) , (58)where v µ = (cid:0) , v i (cid:1) and v i = e i ¯ ı ¯ v ¯ ı is the coordinate ba-sis three-velocity of the observer comoving with the fluid(cf. Appendix A). Note that u µ u µ = − W (cid:0) − v i v i (cid:1) → W = (cid:0) − v i v i (cid:1) − .Following Cardall et al. (2012), we form ‘Eulerian’ and‘Lagrangian’ decompositions of four-vectors and tensors. Inparticular, the Eulerian and Lagrangian decompositions of thecoordinate basis number-flux four-vector are N µ = E N n µ + F µ N , (59) = J N u µ + H µ N , (60)where the lab-frame (or Eulerian-frame) number density andflux (the ‘Eulerian projections’ of the number-flux four-vector) are E N = N ¯0 and F µ N = e µ ¯ ı N ¯ ı , respectively. Thecomoving-frame number density and flux (the ‘Lagrangianprojections’ of the number-flux four-vector) are J N = N ˆ0 and H µ N = L µ ˆ ı N ˆ ı (cf. Equation (21)). Note that F µ N is or-thogonal to the Eulerian observer’s four-velocity ( n µ F µ N =0 ), while H µ N is orthogonal to the four-velocity of the comov-ing observer ( u µ H µ N = 0 ).Similarly, the Eulerian and Lagrangian decompositions ofthe coordinate basis stress-energy tensor are T µν = E n µ n ν + F µ n ν + F ν n µ + S µν , (61) = J u µ u ν + H µ u ν + H ν u µ + K µν , (62) where the monochromatic lab-frame radiation energy density,flux, and stress (the Eulerian projections of the stress-energytensor) are related to the components of the stress-energy ten-sor in the orthonormal tetrad basis by E = T ¯0¯0 , F µ = e µ ¯ ı T ¯ ı ¯0 , and S µν = e µ ¯ ı e ν ¯ T ¯ ı ¯ , (63)respectively. Note that n µ F µ = n µ S µν = n ν S µν = 0 .The monochromatic comoving-frame radiation energy den-sity, flux, and stress (the Lagrangian projections of the stress-energy tensor) are related to the components of the stress-energy tensor in the orthonormal comoving basis by (cf.Equation (25)) J = T ˆ0ˆ0 , H µ = L µ ˆ ı T ˆ ı ˆ0 , and K µν = L µ ˆ ı L ν ˆ T ˆ ı ˆ , (64)where u µ H µ = u µ K µν = u ν K µν = 0 .Finally, the Eulerian and Lagrangian decompositions of thethird-order moment (rank-three tensor) are U µνρ = G n µ n ν n ρ + I µ n ν n ρ + I ν n µ n ρ + I ρ n µ n ν + P µν n ρ + P µρ n ν + P νρ n µ + Q µνρ (65) = ǫ ( J u µ u ν u ρ + H µ u ν u ρ + H ν u µ u ρ + H ρ u µ u ν + K µν u ρ + K µρ u ν + K νρ u µ + L µνρ ) , (66)where G = U ¯0¯0¯0 , I µ = e µ ¯ ı U ¯ ı ¯0¯0 , P µν = e µ ¯ ı e ν ¯ U ¯ ı ¯ ¯0 , Q µνρ = e µ ¯ ı e ν ¯ e ρ ¯ k U ¯ ı ¯ ¯ k , and L µνρ = L µ ˆ ı L ν ˆ L ρ ˆ k L ˆ ı ˆ ˆ k . Again, I µ , P µν , and Q µνρ are purely spatial and orthogonal to thefour-velocity of the Eulerian observer, while L µνρ is orthog-onal to the four-velocity of the comoving observer. TheEulerian and Lagrangian decompositions of the number-fluxfour-vector, the radiation stress-energy tensor, and the third-order moments are useful when deriving conservative momentequations for the radiation field—in particular, when seek-ing conservative evolution equations for lab-frame radiationquantities, expressed in terms of comoving-frame moments.The radiation moment equations for conformally flat space-times can be obtained by specializing the 3+1 moment equa-tions in Cardall et al. (2012), or by adopting the spacetimemetric in Equation (56), with γ ij = ψ ¯ γ ij , and computingconnection coefficients directly. Here we have adopted thelatter approach, but see also Cardall et al. (2012). We obtainthe evolution equation for the monochromatic lab-frame radi-ation energy density by contracting Equation (47) with − n µ .The result is √− g ∂∂t (cid:16) √ γ E (cid:17) + 1 √− g ∂∂x i (cid:16) √ γ h α F i − β i E i(cid:17) + F i ∂ ln α∂x i + S ii ∂ ln ψ ∂τ − S ij α (cid:16) ¯ ∇ i β j − ∂ ln ψ ∂x i β j (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ F ǫ (cid:17) = − n µ ǫ Z Ω p µ C [ f ] d Ω , (67)0 Endeve et al.where we have defined the energy space energy flux F ǫ = − n µ U µνρ ∇ ν u ρ = W n I i ∂v i ∂τ + P ji ∂v i ∂x j + 12 ψ P ij v k ∂ ¯ γ ij ∂x k + (cid:16) I i − G v i (cid:17) ∂ ln α∂x i + P ii (cid:16) ∂ ln ψ ∂τ + v j ∂ ln ψ ∂x j (cid:17) − P ij α (cid:16) ¯ ∇ i β j − ∂ ln ψ ∂x i β j (cid:17) + v i I j α ∂β j ∂x i o + (cid:16) I i v i − G (cid:17) ∂W∂τ + (cid:16) P ji v i − I j (cid:17) ∂W∂x j . (68)We have also defined the covariant derivative with respect tothe conformally related three-metric ¯ ∇ i β j = ∂β j ∂x i − ¯Γ kji β k , (69)where ¯Γ kij = 12 ¯ γ kl (cid:18) ∂ ¯ γ li ∂x j + ∂ ¯ γ lj ∂x i − ∂ ¯ γ ij ∂x l (cid:19) , (70)and the proper time derivative along constant coordinate lines ∂∂τ = 1 α ∂∂t − β i α ∂∂x i . (71)We obtain the monochromatic radiation momentum equa-tion by contracting Equation (47) with the projection operator γ iµ = g iµ + n i n µ . The result is √− g ∂∂t (cid:16) √ γ F i (cid:17) + 1 √− g ∂∂x j (cid:16) √ γ h α S ji − β j F i i(cid:17) − ψ S jk ∂ ¯ γ jk ∂x i + E ∂ ln α∂x i − S jj ∂ ln ψ ∂x i − F j α ∂β j ∂x i − ǫ ∂∂ǫ (cid:16) ǫ S ǫi (cid:17) = γ iµ ǫ Z Ω p µ C [ f ] d Ω , (72)where we have defined the energy space momentum flux S ǫi = γ iµ U µνρ ∇ ν u ρ = W n P ij ∂v j ∂τ + Q kij ∂v j ∂x k + 12 ψ Q jki v l ∂ ¯ γ jk ∂x l + (cid:16) P ji − I i v j (cid:17) ∂ ln α∂x j + Q ji j (cid:16) ∂ ln ψ ∂τ + v k ∂ ln ψ ∂x k (cid:17) − Q jki α (cid:16) ¯ ∇ j β k − ∂ ln ψ ∂x j β k (cid:17) + v j P ik α ∂β k ∂x j o + (cid:16) P ij v j − I i (cid:17) ∂W∂τ + (cid:16) Q kij v j − P ki (cid:17) ∂W∂x k . (73)Equations (67) and (72) are conservative equations for the lab-frame radiation energy density and momentum density. Theindependent variables are the coordinate basis spacetime po-sition components x µ and the neutrino energy measured bya comoving observer ǫ . When integrated over the comoving-frame energy, the equations reduce to familiar position spaceconservation laws. They express exact energy and momentumconservation in the absence of neutrino-matter interactionsand geometry sources (due to gravity and curvilinear coordi-nates). The terms in the energy derivatives (cf. Equations (68)and (73)) result in changes in the radiation energy spectrummeasured in the comoving-frame from gravitational redshifts and acceleration of the comoving observer. Note that Equa-tions (67) and (72) differ from the 3+1 general relativistic mo-ment equations given by Shibata et al. (2011); Cardall et al.(2012), in that they are not expressed in terms of the extrinsiccurvature K ij (an evolved quantity in 3+1 general relativity;e.g., Baumgarte & Shapiro 2010).For the terms inside the energy derivatives in Equations(67) and (72), we can use Equation (34), and write u µ = W ( n µ + v µ ) , to relate Eulerian projections of the third-order moments (Equation (65)) in terms of Eulerian projec-tions of the stress-energy tensor (Equation (61); Cardall et al.2012). In particular, we have W I i = ǫ n F i + S ij v j + Wǫ Q ijk v j v k o , (74) W P ij = ǫ n S ij + Wǫ Q ijk v k o , (75) W (cid:16) I i − G v i (cid:17) = ǫ n(cid:16) F i − E v i (cid:17) + (cid:16) S ij − v i F j (cid:17) v j + (cid:16) Wǫ Q ijk − v i S jk (cid:17) v j v k − v i Wǫ Q jkl v j v k v l o , (76) W (cid:16) P ij − I i v j (cid:17) = ǫ n(cid:16) S ij − F i v j (cid:17) + (cid:16) Wǫ Q ijk − S ik v j (cid:17) v k − v j Wǫ Q ikl v k v l o , (77) W (cid:16) v i I i − G (cid:17) = − ǫ E , (78) W (cid:16) v j P ij − I i (cid:17) = − ǫ F i , (79) W (cid:16) v k Q ijk − P ij (cid:17) = − ǫ S ij . (80)The neutrino number equation can be obtained from Equa-tions (67) and (72). From Equation (33) we have E N = Wǫ (cid:16) E − v i F i (cid:17) , (81) F i N = Wǫ (cid:16) F i − v j S ij (cid:17) . (82)As was done in Section 4.3, we obtain a conservative equationfor the monochromatic lab-frame neutrino number density byadding ǫ − W times Equation (67) and − ǫ W v i contractedwith Equation (72)—and subsequently bringing the necessaryterms inside the time, space, and energy derivatives, respec-tively. The result is √− g ∂∂t (cid:16) √ γ E N (cid:17) + 1 √− g ∂∂x i (cid:16) √ γ h α F i N − β i E N i(cid:17) − ǫ ∂∂ǫ (cid:16) ǫ F ǫ N (cid:17) = 1 ǫ Z Ω C [ f ] d Ω , (83)ONSERVATIVE RADIATION MOMENT EQUATIONS 11where the energy space number flux is F ǫ N = T µν ∇ µ u ν = Wǫ (cid:16) F ǫ − v i S ǫi (cid:17) = W n F i ∂v i ∂τ + S ji ∂v i ∂x j + 12 ψ S ij v k ∂ ¯ γ ij ∂x k + (cid:16) F i − E v i (cid:17) ∂ ln α∂x i + S ii (cid:16) ∂ ln ψ ∂τ + v j ∂ ln ψ ∂x j (cid:17) − S ij α (cid:16) ¯ ∇ i β j − ∂ ln ψ ∂x i β j (cid:17) + v i F j α ∂β j ∂x i o − (cid:16) E − v i F i (cid:17) ∂W∂τ − (cid:16) F j − v i S ji (cid:17) ∂W∂x j . (84)When deriving Equation (83) from Equations (67) and (72),terms emanating from the time derivatives, space derivatives,and the geometry sources cancel exactly with terms emanat-ing from the energy derivatives (see also Cardall et al. 2012).Ideally, a discrete representation of the neutrino energy andmomentum equations can be constructed so that a conserva-tive neutrino number equation can be analogously obtained inthe discrete limit. The discretization of the left-hand sides ofthe energy and momentum equations is then consistent withneutrino number conservation. (A similar consistency mustalso be considered for the right-hand side ot the equations.)Clearly, in order to construct such a discretization, the indi-vidual terms in the neutrino energy and momentum equationscannot be discretized independently . PSEUDO-NEWTONIAN MOMENT EQUATIONS FOR NEUTRINORADIATION TRANSPORT
In this section, we specialize the moment equations pre-sented in Section 5 to the ‘pseudo-Newtonian’ limit by adopt-ing a spacetime metric consistent with the line element (cf.Kim et al. 2009, 2012) ds = − ( 1 + 2 Φ ) dt + ( 1 + 2 Φ ) − ¯ γ ij dx i dx j , (85)where Φ ≪ is a pseudo-Newtonian gravitational po-tential. In this approximation, all orders in the fluid ve-locity v are retained, and the pseudo-Newtonian gravita-tional potential is obtained by solving a Poisson equationin which the trace of the stress-energy tensor—not just therest-mass density—contributes to the source on the right-hand side. We allow for common curvilinear spatial coor-dinates (i.e., Cartesian, spherical, and cylidrical) by setting ¯ γ ij = diag [ 1 , a ( x ) , b ( x ) c ( x ) ] (Section 5). The spa-tial metric γ ij = (1 + 2 Φ) − ¯ γ ij and its inverse are used tolower and raise indices on (spatial) vectors and tensors. Thedeterminant of the spacetime metric becomes √− g = ( 1 + 2 Φ ) − √ ¯ γ. (86)Similarly, we have √ γ = (1 + 2 Φ) − / √ ¯ γ , where √ ¯ γ = abc . Note that the spacetime metric in Equation (85) dif-fers slightly from the Newtonian one given by Misner et al.(1973); Schutz (1985), but agrees to O (Φ) ; i.e., (1+2Φ) − =(1 − O (Φ ) . Using the spacetime metric in Equation(85) turns out to be algebraically advantageous. Note in par-ticular that √− g = ˜ α √ γ , where ˜ α = (1 + 2 Φ) / . Wenote in passing that Kim et al. (2009, 2012) recently com-bined relativistic hydrodynamics with the weak-field (Newto-nian) limit of general relativity (cf. Equation (85)) to studyequilibrium solutions of rotating relativistic stars and foundgood agreement with general relativistic computations in the reported cases. (See also Takiwaki et al. 2009, who adopteda similar approach in magneto-rotational core-collapse simu-lations.)We obtain the pseudo-Newtonian radiation moment equa-tions from Equations (67) and (72) by setting α = ˜ α ≡ (1 + 2 Φ) / , β i = 0 , and ψ = ˜ ψ ≡ (1 + 2 Φ) − / , andretaining all terms of O (Φ) . We state the results for easy ref-erence. The radiation energy equation becomes √− g ∂∂t (cid:16) √ γ E (cid:17) + 1 √− g ∂∂x i (cid:16) √− g F i (cid:17) + F i ∂ Φ ∂x i −S ii ∂ Φ ∂t − ǫ ∂∂ǫ (cid:16) ǫ F ǫ (cid:17) = ˜ αǫ Z Ω p C [ f ] d Ω , (87)where the lab-frame radiation energy density, momentumdensity, and stress ( E , F i , and S ij , respectively) are related tothe comoving-frame moments in Equations (A9)-(A11). Theenergy space energy flux (cf. Equation (68)) reduces to F ǫ = W n I i α ∂v i ∂t + P ji ∂v i ∂x j + 12 ˜ ψ P ij v k ∂ ¯ γ ij ∂x k + (cid:16) I i − G v i (cid:17) ∂ Φ ∂x i − P ii (cid:16) ∂ Φ ∂t + v j ∂ Φ ∂x j (cid:17)o + (cid:16) I i v i − G (cid:17) α ∂W∂t + (cid:16) P ji v i − I j (cid:17) ∂W∂x j . (88)Similarly, the radiation momentum equation becomes √− g ∂∂t (cid:16) √ γ F i (cid:17) + 1 √− g ∂∂x j (cid:16) √− g S ji (cid:17) − ˜ ψ S jk ∂ ¯ γ jk ∂x i + E ∂ Φ ∂x i + S jj ∂ Φ ∂x i − ǫ ∂∂ǫ (cid:16) ǫ S ǫi (cid:17) = γ ij ǫ Z Ω p j C [ f ] d Ω , (89)where the energy space momentum flux is given by (cf. Equa-tion (73)) S ǫi = W n P ij α ∂v j ∂t + Q kij ∂v j ∂x k + 12 ˜ ψ Q jki v l ∂ ¯ γ jk ∂x l + (cid:16) P ji − I i v j (cid:17) ∂ Φ ∂x j − Q ji j (cid:16) ∂ Φ ∂t + v k ∂ Φ ∂x k (cid:17)o + (cid:16) P ij v j − I i (cid:17) α ∂W∂t + (cid:16) Q kij v j − P ki (cid:17) ∂W∂x k . (90)Equations (87) and (89) are valid to all orders in v , but lim-ited to weak gravitational fields ( Φ ≪ ). In the case of nogravitational fields ( Φ ≡ ) they reduce to the special rela-tivistic moment equations (cf. Equations (A22) and (A24) inAppendix A). Moreover, Equations (87) and (89) are consis-tent with the conservative neutrino number equation given by √− g ∂∂t (cid:16) √ γ E N (cid:17) + 1 √− g ∂∂x i (cid:16) √− g F i N (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ F ǫ N (cid:17) = 1 ǫ Z Ω C [ f ] d Ω , (91)2 Endeve et al.where the energy space number flux is (cf. Equation (84)) F ǫ N = W n F i α ∂v i ∂t + S ji ∂v i ∂x j + 12 ˜ ψ S ij v k ∂ ¯ γ ij ∂x k + (cid:16) F i − E v i (cid:17) ∂ Φ ∂x i − S ii (cid:16) ∂ Φ ∂t + v j ∂ Φ ∂x j (cid:17)o − (cid:16) E − v i F i (cid:17) α ∂W∂t − (cid:16) F j − v i S ji (cid:17) ∂W∂x j . (92)The lab-frame number density and number flux are relatedto the corresponding comoving-frame moments in Equations(A6) and (A7), respectively.Equations (87), (89), and (91) simplify even further whenslowly varying gravitational fields are considered ( ∂ Φ /∂t =0 ; cf. Kim et al. 2012). NON-RELATIVISTIC SELF-GRAVITATING NEUTRINORADIATION HYDRODYNAMICS
In this section we detail the full set of neutrino radia-tion hydrodynamics equations intended for our planned non-relativistic simulations of core-collapse supernovae and re-lated systems. The equations are deduced from the relativisticequations derived in previous sections. Per our initial dis-cussion in Section 1, in order to formulate a system of mo-ment equations that are consistent with a conservative equa-tion for the neutrino number density, we adopt different ordersof v for the radiation energy and momentum equations. Weconsider two possible cases: the ‘ O ( v ) -plus’ and the ‘ O ( v ) -minus’ moment equations. The O ( v ) -plus moment equationsconsist of an O ( v ) energy equation and an O ( v ) momentumequation, and are consistent with the conservative O ( v ) neu-trino number equation. Similarly, the O ( v ) -minus momentequations consist of an O ( v ) energy equation and an O (1) momentum equation, and are consistent with the conservative O ( v ) neutrino number equation. Dimensional analysis sug-gests that gravitational effects on the radiation field are for-mally O ( v ) . Therefore, in the O ( v ) -plus moment equationsfor neutrino radiation transport, we also retain terms due toa Newtonian gravitational potential. These ‘gravitational red-shift’ terms may play a non-negligible role in the post-bouncesupernova environment (e.g., Bruenn et al. 2001). More-over, by including gravitational redshift effects, the O ( v ) -plus moment equations are conceptually more similar to thefully relativistic moment equations (e.g., Shibata et al. 2011;Cardall et al. 2012). O ( v ) -Plus Moment Equations for Neutrino Transport In this subsection we detail radiation moment equations for-mally valid to O ( v ) , including effects due to a Newtoniangravitational field. We obtain these Newtonian, O ( v ) -plus ra-diation moment equations from the pseudo-Newtonian equa-tions in Section 6—Equations (87) and (89)—by consideringnon-relativistic fluid velocities. In the radiation moment equa-tions presented in this section, we also omit terms contain-ing the time derivative of the gravitational potential. Further-more, we retain at least all terms that are linear in the fluidvelocity. We retain O ( v ) terms in the radiation energy equa-tion, while we retain O ( v ) terms in the radiation momentumequation. As discussed in the introduction, and in detail be-low, we adopt this ‘mixed’ order of the radiation energy andmomentum equations in order to achieve satisfactory consis-tency with the conservative neutrino number equation. Sincewe assume v ≪ , we set W = ˜ W ≡ v , where v = v i v i . Then, the monochromatic O ( v ) radiation energyequation becomes (cf. Equation (87)) √− g ∂∂t (cid:16) √ γ ˜ E (cid:17) + 1 √− g ∂∂x i (cid:16) √− g ˜ F i (cid:17) + F i ∂ Φ ∂x i − ǫ ∂∂ǫ (cid:16) ǫ ˜ F ǫ (cid:17) = ˜ C E , (93)where the O ( v ) lab-frame radiation energy density and mo-mentum density are related to the comoving-frame moments( J , H i , and K ij ) by ˜ E = E + v i (cid:16) v i J + v j K ij (cid:17) , (94) ˜ F i = F i + v j (cid:16) v i H j + 12 H i v j (cid:17) , (95)and the O ( v ) radiation energy density, momentum density,and stress are related to the comoving-frame moments by E = J + 2 v i H i , (96) F i = H i + v i J + v j K ij , (97) S ij = K ij + v i H j + H i v j (98)(cf. Equations (A9)-(A11)). Radiation quantities adornedwith a tilde (e.g., ˜ E ) are accurate to O ( v ) . To O ( v ) we alsohave H i = e i ¯ ı δ ¯ ı ˆ ı H ˆ ı , (99) K ij = e i ¯ ı e j ¯ δ ¯ ı ˆ ı δ ¯ ˆ K ˆ ı ˆ , (100) L ijk = e i ¯ ı e j ¯ e k ¯ k δ ¯ ı ˆ ı δ ¯ ˆ δ ¯ k ˆ k L ˆ ı ˆ ˆ k , (101)(cf. Equations (A12), (A13), and (A15)). The comoving-frame moments J , H ˆ ı , K ˆ ı ˆ , and L ˆ ı ˆ ˆ k are defined in Section 3.The energy space energy flux, written in terms of comoving-frame moments, is (cf. Equation (88)) ˜ F ǫ = F ǫ + ǫ n(cid:16) J v j + 2 v i K ij (cid:17) α ∂v j ∂t + (cid:16) v j H k + H j v k + v i L kij (cid:17) ∂v j ∂x k + 12 ˜ ψ (cid:16) v j H k + H j v k + v i L jki (cid:17) v l ∂ ¯ γ jk ∂x l + (cid:16) v i K ji − K ii v j (cid:17) ∂ Φ ∂x j − ˜ E α ∂ ln ˜ W∂t − ˜ F i ∂ ln ˜ W∂x i o , (102)where the O ( v ) and O (Φ) terms are defined separately in F ǫ = ǫ n H j α ∂v j ∂t + K kj ∂v j ∂x k + 12 ˜ ψ K jk v l ∂ ¯ γ jk ∂x l + H j ∂ Φ ∂x j o . (103)In Equation (102), except for the last two terms, we havedropped terms of order equal to or higher than O ( v ) and O ( v Φ) . Similarly, the monochromatic O ( v ) radiation mo-mentum equation becomes √− g ∂∂t (cid:16) √ γ F i (cid:17) + 1 √− g ∂∂x j (cid:16) √− g S ji (cid:17) − ˜ ψ S jk ∂ ¯ γ jk ∂x i + (cid:16) J + K jj (cid:17) ∂ Φ ∂x i − ǫ ∂∂ǫ (cid:16) ǫ ˜ S ǫi (cid:17) = C F i , (104)ONSERVATIVE RADIATION MOMENT EQUATIONS 13where energy space momentum flux is ˜ S ǫi = S ǫi − ǫ n F i α ∂ ln ˜ W∂t + S ji ∂ ln ˜ W∂x j o , (105)and the O ( v ) and O (Φ) terms are contained in S ǫi = ǫ n K ij α ∂v j ∂t + L kij ∂v j ∂x k + 12 ˜ ψ L jki v l ∂ ¯ γ jk ∂x l + K ji ∂ Φ ∂x j o . (106)In Equation (105), except for the last two terms, we havedropped terms of order equal to or higher than O ( v ) and O ( v Φ) . Note that ˜ S ǫi is not accurate to O ( v ) , but is adornedwith a tilde to signify that some higher-order terms have beenretained. Also note that we have retained the Φ -dependencein the determinant of the spacetime metric (cf. Equation (86))appearing in the moment equations, and retained factors ˜ α and ˜ ψ , which give rise to higher-order terms in some cases (e.g.,Equation (106)). These should all be kept for consistency withthe conservative O ( v ) number equation (see below).In Equations (93) and (104), the collision terms on the right-hand side are, to O ( v ) and O ( v ) , respectively, ˜ C E = ˜ W Z Ω C [ f ] d Ω + ¯ v ˆ ı Z Ω n ˆ ı C [ f ] d Ω , (107) C F i = v i Z Ω C [ f ] d Ω + γ ij e j ¯ δ ¯ ˆ Z Ω n ˆ C [ f ] d Ω . (108)( ¯ v ˆ ı is the fluid three-velocity with respect to the orthonormaltetrad basis; Appendix A.)Equations (93) and (104) are conservative equations for themonochromatic lab-frame radiation energy density and mo-mentum density, respectively. Equation (93) expresses ex-act conservation of radiation energy in the absence of gravityand neutrino-matter interactions, and—if in addition, Carte-sian coordinates are used—Equation (104) expresses exactconservation of radiation momentum. In the strict O ( v ) limit , Equations (93) and (104) can be compared with corre-sponding equations derived by Buchler (1979, his Equations(9) and (10), respectively). See also Kaneko et al. (1984);Munier & Weaver (1986b); Buras et al. (2006). We obtainthe strict O ( v ) limit of Equations (93) and (104) by set-ting Φ = 0 , letting ˜ E , ˜ F i , ˜ F ǫ , ˜ S ǫi → E , F i , F ǫ , S ǫi , and ˜ C E → C E by setting ˜ W = 1 in Equation (107). Buchler’sequations evolve the comoving-frame moments and are non-conservative . However, to O ( v ) , the terms inside the energyderivatives in Equations (93) and (104) agree with the energyderivative terms in Buchler’s energy and momentum equa-tions. Shibata et al. (2011) also list conservative momentequations for neutrino radiation transport in the ‘slow-motion’limit (cf. their Equations (8.5) and (8.6)). They do not includegravitational effects, or any higher-order terms in the fluidvelocity. They also do not include the terms involving thederivative of the fluid velocity with respect to time in the en-ergy derivative terms, which can be important when λ/τ . (e.g., in optically thin regions; Buchler 1979), where λ and τ are typical length and time scales, respectively.By integrating Equations (93) and (104) over the comoving-frame energy (with ǫ dǫ as the measure of integration) we By the “strict” O ( v ) limit we mean that all (or most) terms that are atmost linear in the fluid velocity have been retained in the radiation energy and momentum equations. As we emphasize in this section, in this strict O ( v ) limit, the radiation moment equations are not fully consistent with a conservative equation for the neutrino number density. obtain the energy-integrated (grey) Eulerian-frame radiationenergy equation, √− g ∂ (cid:16) √ γ ˜ E (cid:17) ∂t + 1 √− g ∂ (cid:16) √− g ˜ F i (cid:17) ∂x i = − F i ∂ Φ ∂x i + ˜ C E , (109)and the grey Eulerian-frame radiation momentum equation, √− g ∂ (cid:16) √ γ F i (cid:17) ∂t + 1 √− g ∂ (cid:16) √− g S ji (cid:17) ∂x j − ˜ ψ S jk ∂ ¯ γ jk ∂x i = − (cid:16) J + K jj (cid:17) ∂ Φ ∂x i + C F i , (110)where the grey radiation variables are (cid:8) ˜ E, ˜ F i , F i , S ij , J, K ij , ˜ C E , C F i (cid:9) = Z ∞ (cid:8) ˜ E , ˜ F i , F i , S ij , J , K ij , ˜ C E , C F i (cid:9) ǫ dǫ. (111)Note that we use the font type to distinguish the energy-dependent radiation variables (denoted with calligraphic font)from the grey radiation variables. The Eulerian-frame greyradiation energy density, momentum density, and stress ( E , F i , and S ij , respectively) are related to the correspondingcomoving-frame quantities ( J , H i , and K ij ) via relationsanalogous to those listed in Equations (96)-(98). In the strict O ( v ) limit, for Cartesian coordinates, and with Φ = 0 , Equa-tions (109) and (110) correspond to Equations (32a) and (32b)given by Lowrie et al. (2001), which they refer to as the ‘cor-rect’ O ( v ) Eulerian-frame radiation energy and momentumequations, respectively. One requirement for the correct O ( v ) limit is that the hyperbolic wave speeds associated with thesystem of equations is bounded by the speed of light. Thepresence of geometry sources due to the gravitational fieldor the use of curvilinear coordinates does not destroy the hy-perbolic character of the equations, since they do not con-tain any differential operators acting on the radiation variables(Banyuls et al. 1997).Equations (93) and (104) are correct to O ( v ) and O ( v ) ,respectively. We have explicitly retained higher-order (non-linear) terms in the fluid velocity inside the energy deriva-tive of both equations. In the energy equation, we have alsoretained terms containing the dot product of the fluid three-velocity with the gradient of the gravitational potential. How-ever, when combined, the system is only accurate to O ( v ) .Hence, we refer to these equations as the ‘ O ( v ) -plus’ approx-imation of the radiation moment equations. These higher-order terms must (by definition) be small in the O ( v ) limitwe are interested in here. However, by retaining these terms,the solution to Equations (93) and (104) becomes consistentwith the conservative O ( v ) neutrino number equation. Thisconsistency may help enable exact lepton number conserva-tion in the two-moment model for neutrino radiation transportin the O ( v ) limit. We elaborate further on the details here.Equation (81) relates the Eulerian-frame neutrino numberdensity to the corresponding neutrino energy density and mo-mentum density. The equations governing their time evolu-tion are similarly related. In particular, by adding ǫ − ˜ W timesEquation (93) and − ǫ − ˜ W v i contracted with Equation (104)4 Endeve et al.we obtain √− g ∂∂t (cid:16) √ γ ˜ Wǫ h ˜ E − v i F i i(cid:17) + 1 ǫ h ˜ W F i α ∂v i ∂t − (cid:16) ˜ E − v i F i (cid:17) α ∂ ˜ W∂t i (112)from the time derivative terms. Remember ˜ W = 1 + v .Similarly we obtain √− g ∂∂x j (cid:16) √− g ˜ Wǫ h ˜ F j − v i S ji i(cid:17) + 1 ǫ h ˜ W S ji ∂v i ∂x j − (cid:16) ˜ F j − v i S ji (cid:17) ∂ ˜ W∂x j i (113)from the space derivative terms, while the combination of thegeometry sources results in ˜ Wǫ n
12 ˜ ψ S ij v k ∂ ¯ γ ij ∂x k + (cid:16) H j + v i K ji + K ii v j (cid:17) ∂ Φ ∂x j o . (114)From the energy derivative terms we obtain − ǫ ∂∂ǫ (cid:16) ǫ ˜ F ǫ N (cid:17) − ǫ ˜ F ǫ N , (115)where the number-flux in energy space is obtained from ˜ F ǫ N = ˜ Wǫ (cid:16) ˜ F ǫ − v i ˜ S ǫi (cid:17) = ˜ W n(cid:16) H i + J v i + v j K ij (cid:17) α ∂v i ∂t + (cid:16) K ji + v i H j + H i v j (cid:17) ∂v i ∂x j + 12 ˜ ψ (cid:16) K ij + v i H j + H i v j (cid:17) v k ∂ ¯ γ ij ∂x k + (cid:16) H j + v i K ji + K ii v j (cid:17) ∂ Φ ∂x j o − (cid:16) ˜ E − v i F i (cid:17) α ∂ ˜ W∂t − (cid:16) ˜ F j − v i S ji (cid:17) ∂ ˜ W∂x j . (116)Finally, the collision terms combine to give ˜ C N = ˜ Wǫ (cid:16) ˜ C E − v i C F i (cid:17) = (cid:16) − v (cid:17) ǫ Z Ω C [ f ] d Ω . (117)(Note that v i e i ¯ ı δ ¯ ı ˆ ı = ¯ v ˆ ı ; Appendix A.)When adding Equations (112)-(115), which result in theleft-hand side of the conservative neutrino number equation(Equation (118) below), we note that the second term in Equa-tion (115) (cf. Equation (116)) cancels exactly with the ‘left-over’ terms in Equations (112) and (113), and Equation (114).In particular, the first term in the second line of Equation (112)cancels with terms emanating from the second line on theright-hand side of Equation (116). Similarly, the first termin the second line of Equation (113) cancels with terms em-anating from the third line on the right-hand side of Equa-tion (116). The terms emanating from the geometry sources(Equation (114)) cancel with terms emanating from the fourthand fifth line on the right-hand side of Equation (116). Finally,the terms involving time and space derivatives of the approxi-mate Lorentz factor ˜ W cancel with the terms emanating from the last line on the right-hand side of Equation (116). More-over, when contracted with − ˜ W ǫ − v i , several terms in theenergy derivative of the radiation momentum equation cancelexactly with some of the O ( v ) terms in the energy derivativeof the ˜ W ǫ − -multiplied radiation energy equation (these can-cellations are only obtained by retaining higher-order termsin Equation (93)). We also find that the O ( v ) terms involv-ing the time and space derivatives of ˜ W in Equations (102)and (105) are needed to cancel with the corresponding termsin Equations (112) and (113), obtained after pulling ˜ W insidethe time and space derivatives.By combining all the terms (Equations (112)-(117)) we ob-tain √− g ∂∂t (cid:16) √ γ ˜ E N (cid:17) + 1 √− g ∂∂x i (cid:16) √− g ˜ F i N (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ ˜ F ǫ N (cid:17) = ˜ C N , (118)where the Eulerian-frame O ( v ) number density and numberflux density are (cf. Equations (A6) and (A7)) ˜ E N = ˜ Wǫ (cid:16) ˜ E − v i F i (cid:17) = ˜ W J N + v i H i N + O ( v ) , (119) ˜ F i N = ˜ Wǫ (cid:16) ˜ F i − v j S ij (cid:17) = H i N + v i J N + O ( v ) , (120)where we have used Equations (94)-(95) and (96)-(98). Wehave H i N = e i ¯ ı δ ¯ ı ˆ ı H ˆ ı N , and the comoving-frame moments J N and H ˆ ı N are defined in Section 3. We have retained the O ( v ) term appearing in the collision term on the right-hand side ofEquation (118) (defined in Equation (117)), but this term cansafely be dropped in practical computations.Equation (118) is a conservative equation for the Eulerian-frame neutrino number density. It is valid to O ( v ) . Togetherwith an equation for the electron density, it states that the totallepton number is conserved during lepton number exchangewith the fluid (see discussion in Section 7.4). In the absenceof neutrino-matter interactions, Equation (118) states conser-vation of particle number. It is obtained analytically from themonochromatic radiation energy and momentum equations.Modulo terms of O ( v ) or higher (cf. Equations (117), (119),and (120)), it can also be obtained directly from Equation(83).A numerical solution to Equations (93) and (104) shouldalso be consistent with Equation (118) in order to ensurelepton number conservation in simulations of neutrino radi-ation transport in core-collapse supernovae and related sys-tems (Section 7.4). Ideally, the discretized neutrino energyand momentum equations can be similarly combined to obtaina discrete representation of the conservative neutrino numberequation. Note that we arrived at Equation (118) in an an-alytically exact manner (due to exact cancellations; we didnot throw away any O ( v ) terms outside the time, space, orenergy derivatives). By retaining the higher-order terms inEquation (93) and Equation (104), the two-moment modelof neutrino radiation transport is consistent with the conser-vative O ( v ) neutrino number equation. If the higher-orderterms in Equations (93) and (104) are not retained (i.e., in thestrict O ( v ) limit), there will be additional O ( v ) terms in thenumber equation derived from the O ( v ) energy and momen-tum equations, and the consistency of the two-moment modelwith the number equation is only O ( v ) —which may be ac-ONSERVATIVE RADIATION MOMENT EQUATIONS 15ceptable in practical numerical computations. However, notethat by omitting the ∂v i /∂t -term inside the energy derivativein the energy equation, the consistency with the conservativeneutrino number equation can reduce to O (1) . Thus, omit-ting this term in numerical simulations of core-collapse super-novae can potentially result in severe consequences for leptonnumber conservation.By integrating Equation (118) over the comoving-frameenergy we obtain the Eulerian-frame grey neutrino numberequation √− g ∂∂t (cid:16) √ γ ˜ E N (cid:17) + 1 √− g ∂∂x i (cid:16) √− g ˜ F i N (cid:17) = ˜ C N , (121)where (cid:8) ˜ E N , ˜ F i N , ˜ C N (cid:9) = Z ∞ (cid:8) ˜ E N , ˜ F i N , ˜ C N (cid:9) ǫ dǫ. (122)We have presented conservative equations for themonochromatic lab-frame radiation energy density and mo-mentum density—valid to O ( v ) and O ( v ) , respectively—forcommon curvilinear coordinates (Cartesian, spherical polar,and cylindrical); Equations (93) and (104), respectively. Wehave also demonstrated how these equations (and thereforealso their solutions) are fully consistent with the conservativeequation for the monochromatic O ( v ) lab-frame radiationparticle density (Equation (118)). Note that for a consistentdescription of neutrino radiation hydrodynamics with the O ( v ) -plus moment equations, the hydrodynamics equationsmay also have to be promoted to include higher-order termsin the fluid velocity. Further Simplifications of the Neutrino RadiationTransport Equations: O ( v ) -Minus Moment Equations In this subsection we further specialize the radiation mo-ment equations presented in Section 7.1 by introducingthe O ( v ) -minus moment equations. Although apparentlyless complex than the fully relativistic moment equations(cf. Shibata et al. 2011; Cardall et al. 2012, Section 5),the O ( v ) -plus moment equations presented in Section 7.1are still nontrivial to discretize for numerical computations.(Their complexity rivals that of the pseudo-Newtonian mo-ment equations presented in Section 6.) Therefore, we pro-pose a further simplification, which involves solving the O ( v ) energy equation and the O (1) momentum equation, as a use-ful intermediate (first) step beyond the O (1) moment for-malism previously used by some to model neutrino radiationtransport in core-collapse supernovae (e.g., Burrows et al.2006). The resulting system is formally valid to O (1) . How-ever, the moment equations are consistent with the conser-vative O ( v ) neutrino number equation. Moreover, the equa-tions presented here evolve the radiation energy density and momentum density, and, with a proper closure prescrip-tion, may be an improvement over the equations solved incommon O ( v ) multigroup flux-limited diffusion approaches(e.g., Swesty & Myra 2009; Zhang et al. 2012). In the non-relativistic, Newtonian gravity limit, we have Φ , c s , v ≪ ,where c s is the sound speed. Thus, we ignore gravitationaleffects on the radiation field in this subsection (i.e., Φ = 0 ).Then, the monochromatic O ( v ) lab-frame radiation energy equation becomes (cf. Equation (93)) ∂ E ∂t + 1 √ ¯ γ ∂∂x i (cid:16) √ ¯ γ F i (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ F ǫ (cid:17) = Z Ω C [ f ] d Ω + ¯ v ˆ ı Z Ω n ˆ ı C [ f ] d Ω = C E , (123)where the O ( v ) lab-frame radiation energy density and mo-mentum density, E and F i , are given by Equations (96) and(97), respectively. The energy space energy flux F ǫ is givenby Equation (103), with Φ = 0 .Similarly, the monochromatic O (1) radiation momentumequation becomes (cf. Equation (104)) ∂ H i ∂t + 1 √ ¯ γ ∂∂x j (cid:16) √ ¯ γ K ji (cid:17) − K jk ∂ ¯ γ jk ∂x i = ¯ γ ij e j ¯ δ ¯ ˆ Z Ω n ˆ C [ f ] d Ω = C H i . (124)In Equation (124), the energy derivative terms vanish sincewe have dropped all the velocity-dependent terms in additionto the gravitational terms (cf. Equation (105)).Equations (123) and (124) are consistent with the conser-vative O ( v ) neutrino number equation: by adding ǫ − timesEquation (123) and − ǫ − v i contracted with Equation (124)we obtain ∂ E N ∂t + 1 √ ¯ γ ∂∂x i (cid:16) √ ¯ γ F i N (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ F ǫ N (cid:17) = 1 ǫ Z Ω C [ f ] d Ω , (125)where the O ( v ) lab-frame number density and number fluxdensity, and the number flux in energy space are E N = J N + v i H i N , (126) F i N = H i N + v i J N , (127) F ǫ N = H i ∂v i ∂t + K ji ∂v i ∂x j + 12 K ij v k ∂ ¯ γ ij ∂x k . (128)Equation (125) is the conservative O ( v ) number equation,and is obtained in a manner similar to the O ( v ) number equa-tion detailed in Section 7.1. On the left-hand side, the extraterms emanating from pulling the fluid three-velocity insidethe time and space derivatives of the momentum equation can-cel with the extra terms emanating from pulling ǫ − inside theenergy derivative of the energy equation. On the right-handside, the velocity-dependent term in the collision term of theenergy equation cancels with the contraction of − ǫ − v i withthe collision term of the momentum equation. However, byreducing the order of the energy and momentum equations(to O ( v ) and O (1) , respectively) the number of cancellationsthat occur is dramatically reduced. Because of these simplifi-cations, Equations (123) and (124), which are consistent withthe conservative O ( v ) neutrino number equation, may be asuitable starting point for developing lepton number conserva-tive numerical methods for neutrino radiation hydrodynamicsbased on the two-moment model. Non-Relativistic, Self-Gravitating Hydrodynamics
For self-consistent Newtonian, non-relativistic simulationsof self-gravitating neutrino radiation hydrodynamics, the mo-ment equations for the radiation field in Section 7.2 (one setfor each of the neutrino species) must be coupled to the equa-tions of self-gravitating hydrodynamics. For completeness6 Endeve et al.we list the non-relativistic hydrodynamics equations includ-ing self-gravity in this subsection. (For brevity do we notconsider nuclear reactions here. However, see for examplePlewa & M¨uller 1999.) We consider a perfect fluid; i.e., weignore fluid viscosity and thermal conduction.The equations describing a self-gravitating perfect fluid(e.g., Landau & Lifshitz 1959) include the mass conservationequation ∂ρ∂t + 1 √ ¯ γ ∂∂x i (cid:16) √ ¯ γρ v i (cid:17) = 0 , (129)the fluid momentum equation ∂F f i ∂t + 1 √ ¯ γ ∂∂x j (cid:16) √ ¯ γS j f i (cid:17) − S jk f ∂ ¯ γ jk ∂x i = − ρ ∂ Φ ∂x i − C H i , (130)and the fluid energy equation ∂E f ∂t + 1 √ ¯ γ ∂∂x i (cid:16) √ ¯ γ h E f + p i v i (cid:17) = − ρv i ∂ Φ ∂x i − C E . (131)For simulations involving a nuclear equation of state, Equa-tions (129)-(131) must be supplied with a balance equationfor the electron number density ∂n e ∂t + 1 √ ¯ γ ∂∂x i (cid:16) √ ¯ γ n e v i (cid:17) = − (cid:16) C N ν e − C N ¯ ν e (cid:17) . (132)In Equations (129)-(132), the fluid energy density, momen-tum density and stress are E f = e + ρ v i v i , F f i = ρ v i ,and S i f j = ρ v i v j + δ ij p , where ρ = ¯ m b n b , v i , p , e , and n e denote the mass density, the i th component of the fluid three-velocity, the fluid pressure and internal energy density, andthe electron density (electrons minus positrons), respectively.The average baryon mass and the baryon density are denoted ¯ m b and n b . Equations (129)-(132) are closed with the spec-ification of an equation of state p = p ( ρ, T, Y e ) , where T and Y e = n e /n b are the fluid temperature and electron frac-tion, respectively. In Equations (130) and (131), the collisionterms include energy and momentum exchange with all neu-trino species s ; i.e., (cid:8) C E , C H i (cid:9) = X s (cid:8) C Es , C Hs i (cid:9) . (133)Here we consider electron lepton number exchange betweenthe fluid and the neutrino radiation field. Only interactionsinvolving electron neutrinos ( ν e ) and electron antineutrinos( ¯ ν e ) contribute to the right-hand side of Equation (132).The Newtonian gravitational potential is obtained by solv-ing Poisson’s equation √ ¯ γ ∂∂x i (cid:16) √ ¯ γ ¯ γ ij ∂ Φ ∂x j (cid:17) = 4 π ρ. (134)Note that in the non-relativistic, Newtonian gravity limit onlythe mass density contributes to the gravitational field (as op-posed to all types of stress-energy in general relativity); i.e.,we assume v ≪ and e, p, J, K ii ≪ ρ .Using Equation (134), we can rewrite the gravitational forcein Equation (130) as ρ ∂ Φ ∂x i = 1 √ ¯ γ ∂∂x j (cid:16) √ ¯ γ S j Φ i (cid:17) − S jk Φ ∂ ¯ γ jk ∂x i , (135) where the gravitational stress tensor is defined as S i Φ j = 14 π (cid:16) Φ ,i Φ ,j −
12 Φ ,k Φ ,k δ ij (cid:17) , (136)with Φ ,i = ∂ Φ /∂x i and Φ ,i = ¯ γ ij Φ ,j . We then obtain aconservative fluid momentum equation ∂F f i ∂t + 1 √ ¯ γ ∂∂x j (cid:16) √ ¯ γ h S j f i + S j Φ i i(cid:17) − (cid:16) S jk f + S jk Φ (cid:17) ∂ ¯ γ jk ∂x i = − C H i . (137)In the absence of neutrino-matter interactions, and if Carte-sian coordinates are used, Equation (137) states that the fluidmomentum is conserved.Using Equations (129), (131), and (134), we can write aconservative equation for the energy density E f + ρ Φ ∂∂t (cid:16) E f + E Φ (cid:17) + 1 √ ¯ γ ∂∂x i (cid:16) √ ¯ γ h(cid:16) E f + p (cid:17) v i + F i Φ i(cid:17) = − C E , (138)where E Φ = ρ Φ is the gravitational energy density, F i Φ = ρ Φ v i + 18 π ¯ γ ij (cid:16) Φ ∂ Φ ,t ∂x j − Φ ,t ∂ Φ ∂x j (cid:17) (139)is the gravitational energy flux density, and Φ ,t = ∂ Φ /∂t .Equation (138) states that, in the absence of neutrino-matterinteractions, the fluid energy (internal plus kinetic) plus thegravitational potential energy is conserved. Conservation Laws in Non-Relativistic, Self-GravitatingNeutrino Radiation Hydrodynamics
We discuss conservation laws for non-relativistic, self-gravitating neutrino radiation hydrodynamics in this subsec-tion. Tracking the evolution of conserved quantities is usefulwhen evaluating the physical reliability of numerical simu-lations. Equations (123) and (124) (or more precisely, theirenergy-integrated versions) combined with Equations (129)-(132) state the conservation of several physical quantities inneutrino radiation hydrodynamics. Mass conservation is triv-ially stated by Equation (129). Equations (124) and (130) re-sult in a total (fluid plus radiation) momentum equation ∂∂t (cid:16) F f i + H i (cid:17) + 1 √ ¯ γ ∂∂x j (cid:16) √ ¯ γ h S j f i + K ji i(cid:17) − (cid:16) S jk f + K jk (cid:17) ∂ ¯ γ jk ∂x i = − ρ ∂ Φ ∂x i . (140)Alternatively, combining Equations (124) and (137) results in ∂∂t (cid:16) F f i + H i (cid:17) + 1 √ ¯ γ ∂∂x j (cid:16) √ ¯ γ h S j f i + S j Φ i + K ji i(cid:17) − (cid:16) S jk f + S jk Φ + K jk (cid:17) ∂ ¯ γ jk ∂x i = 0 . (141)Equations (123) and (131) result in a conservative equa-tion for the “total” (internal plus kinetic plus radiation) energydensity of a radiating flow ∂∂t (cid:16) E f + E (cid:17) + 1 √ ¯ γ ∂∂x i (cid:16) √ ¯ γ h ( E f + p ) v i + F i i(cid:17) = − ρv i ∂ Φ ∂x i . (142)ONSERVATIVE RADIATION MOMENT EQUATIONS 17Combining Equations (123) and (138) results in a conserva-tive equation for the total (internal plus kinetic plus gravita-tional plus radiation) energy density of a self-gravitating radi-ating flow ∂∂t (cid:16) E f + E Φ + E (cid:17) + 1 √ ¯ γ ∂∂x i (cid:16) √ ¯ γ h ( E f + p ) v i + F i Φ + F i i(cid:17) = 0 . (143)In Equations (140)-(143), the neutrino energy density, mo-mentum density, and stress include contributions from all neu-trino species; e.g., n E, F i , S ji o = X s n E s , F s i , S j s i o . (144)Equations (141) (assuming Cartesian coordinates are used)and (143) are exactly (and locally) conservative. Therefore,energy and momentum conservation can be useful checksfor evaluating the physical reliability of non-relativistic neu-trino radiation hydrodynamics simulations based on the O ( v ) -minus moment equations. However, for radiation hydrody-namics based on the O ( v ) -plus moment equations, or thepseudo-Newtonian moment equations in Section 6, where thegravitational potential is obtained by solving a modified Pois-son equation (e.g., Kim et al. 2012) and the hydrodynamicsequations are promoted to include more relativistic effects, weare unable to take steps similar to those taken to obtain Equa-tions (137) and (138) from Equations (130) and (131), re-spectively. Moreover, there are additional source terms on theright-hand sides of the total momentum and energy equationsdue to changes in the radiation momentum and energy causedby the gravitational field (e.g., the ‘bending of light’ effect andgravitational redshifts). These (locally) non-vanishing grav-itational source terms (and the lack of global conservation)limit the usefulness of tracking total energy and momentum aschecks on the physical reliability of (self-gravitating) neutrinoradiation hydrodynamics simulations. (Note however that in3+1 general relativity, the AMD mass is conserved in asymp-totically flat spacetimes; e.g., Baumgarte & Shapiro 2010).Consistency with the neutrino number equation results inlepton number conservation in (self-gravitating) neutrino ra-diation hydrodynamics. Combining the energy-integratedversion of Equation (125) (for electron neutrinos and electronantineutrinos) with Equation (132) results in a conservationequation for the total electron lepton number ∂∂t (cid:16) E N,νe − E N, ¯ νe + n e (cid:17) + 1 √ ¯ γ ∂∂x i (cid:16) √ ¯ γ h ˜ F i N ,ν e − ˜ F i N , ¯ ν e + n e v i i(cid:17) = 0 . (145)Lepton number conservation in numerical simulations basedon solving Equations (123) and (124), or any of the other radi-ation moment equations presented in this paper, may dependsensitively on the chosen discretization. In particular, the dis-cretized neutrino energy and momentum equations must beconsistent (in the sense discussed in Sections 7.1 and 7.2)with a discretized version of the conservative neutrino numberequation; Equation (125) in the particular case discussed here.When discretizing the monochromatic energy and momentumequations, attention should be paid to the cancellations thatoccur when deriving the neutrino number equation from theenergy and momentum equations, which ideally also occur in the discrete limit. To achieve this, the individual terms in themoment equations for the neutrino four-momentum should bediscretized in a coordinated rather than independent fashion .Despite the lack of energy conservation due to gravitationalredshifts discussed above, the total lepton number is still con-served. Thus, lepton number conservation may serve as anextremely useful gauge on the physical consistency of numer-ical simulations of neutrino radiation hydrodynamics (e.g.,Liebend¨orfer et al. 2004; Lentz et al. 2012a). SUMMARY AND DISCUSSION
In preparation for development of numerical methods formultidimensional neutrino radiation hydrodynamics, withthe eventual goal of simulating the explosion mechanismof core-collapse supernovae, we have derived conservative,monochromatic general relativistic moment equations for theradiation four-momentum (cf. Equation (47)). The radi-ation moment equations are conservative in the sense that(modulo radiation-matter interactions and geometry sources)the time rate of change of the radiation four-momentum isgoverned by space and momentum space divergences. Wehave used the freedom to choose distinct spacetime andmomentum space coordinates (cf. Lindquist 1966; Riffert1986; Mezzacappa & Matzner 1989; Cardall & Mezzacappa2003). The evolved radiation quantities are functions of thecoordinate basis spacetime position components x µ and theradiation energy ǫ measured by an observer comoving withthe fluid. (When integrated over the comoving-frame en-ergy, the equations reduce to familiar position space conser-vation laws.) The specific choice of phase space coordinatesis motivated by our intent to develop numerical methods forcomputer simulations of multidimensional neutrino transport.Neutrino-matter interactions are on the one hand most eas-ily handled computationally in the comoving-frame. On theother hand, conservation of global quantities (e.g., energyand lepton number) is naturally expressed in the laboratoryframe. Convenient treatment of neutrino-matter interactions and global conservation are naturally handled with the cho-sen phase space coordinates.We have presented radiation moment equations valid forconformally flat spacetimes (Section 5; cf. Equations (67)and (72)). We have further specialized the radiation mo-ment equations to the pseudo-Newtonian and the Newtoniangravity, O ( v ) limits (Sections 6 and 7). Furthermore, in the O ( v ) limit, we have presented the O ( v ) -plus (Section 7.1)and O ( v ) -minus (Section 7.2) moment equations. The O ( v ) -plus radiation energy and momentum equations are givenby Equations (93) and (104), respectively. In the no grav-ity, strict O ( v ) limit, these equations are conservative for-mulations of similar, non-conservative equations presentedby other authors (e.g., Buchler 1979; Kaneko et al. 1984;Munier & Weaver 1986b). The O ( v ) -minus radiation en-ergy and momentum equations are given by Equations (123)and (124), respectively. Special relativistic moment equationsare given in Appendix A, while general relativistic momentequations for spherically symmetric spacetimes are given inAppendix B. The moment equations given in the appen-dices are also conservative versions of non-conservative equa-tions presented by other authors (e.g., Castor 1972; Mihalas1980; Munier & Weaver 1986b; M¨uller et al. 2010). (SeeShibata et al. 2011; Cardall et al. 2012, for conservative 3+1general relativistic moment equations.)We have paid special attention to the issue of neutrino num-ber (and lepton number) conservation. For numerical meth-8 Endeve et al.ods based on solving moment equations for the neutrino four-momentum, total lepton number conservation will likely serveas a very useful gauge on the physical consistency of simula-tions of neutrino radiation hydrodynamics. To this end, wehave exposed the relationship between the equations for theneutrino four-momentum and the neutrino number equationwith multiple examples (cf. Sections 4.3, 5, 6, and 7, and Ap-pendix A). The lab-frame radiation number density is relatedto the lab-frame energy and momentum densities by Equa-tion (81). The conservative neutrino number equation andthe equations for the neutrino four-momentum are similarlyrelated. The non-relativistic limits of the radiation momentequations are not uniquely defined. We obtain consistencywith the conservative neutrino number equation by adoptingdifferent orders of v for the energy and momentum equations(cf. the O ( v ) -minus and O ( v ) -plus approximations in Section7). In particular, as was detailed in Section 7.1, when car-rying out the steps to obtain the conservative neutrino num-ber equation from the conservative energy and momentumequations, we observe that terms emanating from the timeand space derivatives and the geometry sources cancel withterms emanating from the energy derivatives. The remainingterms constitute the left-hand side of the number equation.(The right-hand sides of the equations are similarly related.)The energy and momentum equations are therefore consis-tent with the conservative number equation. In the O ( v ) -plus limit, consistency with the conservative O ( v ) numberequation is obtained by adopting different orders of v and re-taining some higher-order terms in the energy and momen-tum equations (accurate to O ( v ) and O ( v ) , respectively). Inthe O ( v ) -minus limit, we obtain consistency with the con-servative O ( v ) number equation by adopting energy and mo-mentum equations accurate to O ( v ) and O (1) , respectively. Ideally, discrete representations of the radiation energy andmomentum equations can be constructed so that a conserva-tive neutrino number equation can be analogously obtained inthe discrete limit. The discretization is then consistent withneutrino number conservation (which is necessary to ensurelepton number conservation; Section 7.4). The realization ofthis consistency in a numerical method for neutrino radiationtransport based on the two-moment model derived here willbe the focus of a future study.General relativistic effects are important in the core-collapse supernova environment (e.g., Bruenn et al. 2001;M¨uller et al. 2012a), and definitive simulations elucidatingthe explosion mechanism of core-collapse supernovae musteventually be performed in full general relativity (possiblyemploying multi-energy and multi-angle neutrino transport).To this end, the pseudo-Newtonian and the Newtonian, O ( v ) -plus equations may represent useful self-consistent approxi-mations beyond the O ( v ) -minus equations. They are closelyrelated to the corresponding equations valid for conformallyflat spacetimes, which again are close to the fully general rel-ativistic equations (Shibata et al. 2011; Cardall et al. 2012).It seems plausible that numerical methods developed for theself-consistent moment equations presented in this paper—which arguably are easier to work with—can be extended insteps of increasing degree of complexity to the fully gen-eral relativistic case, as has been done in the case of con-servative methods for hydrodynamics (e.g., Font et al. 1994;Banyuls et al. 1997).This research was supported by the Office of AdvancedScientific Computing Research and the Office of NuclearPhysics, U.S. Department of Energy. APPENDIXCONSERVATIVE SPECIAL RELATIVISTIC MOMENT EQUATIONS
In this appendix we present conservative, multidimensional, monochromatic moment equations valid in special relativity. Theequations are obtained from the general relativistic moment equations presented in Sections 4 and 5 in the limit of flat spacetime.First we list sufficiently general equations to accommodate Cartesian, spherical polar, and cylindrical coordinates. Then wespecialize the special relativistic moment equations to spherical symmetry and compare with the corresponding equations givenby Mihalas & Mihalas (1999). Monochromatic, special relativistic moment equations have been presented elsewhere (e.g.,Castor 1972; Mihalas 1980; Munier & Weaver 1986b; Mihalas & Mihalas 1999). However, the equations presented in thisappendix are in fully conservative form. We also present the conservative neutrino number equation, and discuss the relationshipbetween the conservative number equation and the conservative radiation energy and momentum equations.The invariant line element is ds = g µν dx µ dx ν , (A1)where, for flat spacetime, the covariant metric tensor is diagonal and given by g µν = diag [ − , , a ( x ) , b ( x ) c ( x )] . Thecontravariant metric tensor is g µν = diag [ − , , a − , b − c − ] . The transformation from the orthonormal tetrad basis to theglobal coordinate basis may simply be written as e µ ¯ µ = diag [1 , , a − , b − c − ] , and the corresponding inverse transformationis e ¯ µµ = diag [1 , , a, b c ] ; i.e., e α ¯ µ e ¯ µµ = δ αµ . We also have √− g = √ γ = abc . The Lorentz transformation (or boost) from theorthonormal comoving basis to the orthonormal (non-comoving) tetrad basis is given by Λ ¯ µ ˆ µ = (cid:18) Λ ¯0ˆ0 Λ ¯0ˆ ı Λ ¯ ı ˆ0 Λ ¯ ı ˆ ı (cid:19) = (cid:18) W W ¯ v ˆ ı W ¯ v ¯ ı δ ¯ ı ˆ ı + W W +1 ¯ v ¯ ı ¯ v ˆ ı (cid:19) , (A2)where W = ( 1 − ¯ v ¯ ı ¯ v ¯ ı ) − / is the Lorentz factor. In Equation (A2), ¯ v ¯ ı and ¯ v ˆ ı are three-velocity components defined withrespect to the orthonormal tetrad basis, and are therefore accented with a bar ( ¯ v ¯ ı = ¯ v ¯ ı = ¯ v ˆ ı = ¯ v ˆ ı , for ¯ ı = ˆ ı )—they should not beconsidered spatial components of a four-vector. The placement of, and the accent on, the indices of the three-velocity componentsmatch those of the Lorentz transformation Λ ¯ µ ˆ µ . This results in unambiguous notation when relations between quantities in thedifferent reference frames are made explicit (see for example Equations (A16)-(A17) and (A18)-(A20) below). The Lorentztransformation from the orthonormal tetrad basis to the orthonormal comoving basis, is obtained by reversing the sign on theONSERVATIVE RADIATION MOMENT EQUATIONS 19three-velocity components; i.e., Λ ˆ µ ¯ µ = (cid:18) Λ ˆ0¯0 Λ ˆ0¯ ı Λ ˆ ı ¯0 Λ ˆ ı ¯ ı (cid:19) = (cid:18) W − W ¯ v ¯ ı − W ¯ v ˆ ı δ ˆ ı ¯ ı + W W +1 ¯ v ˆ ı ¯ v ¯ ı (cid:19) . (A3)With the specifications of the transformations between the orthonormal tetrad basis and the coordinate basis, and the Lorentztransformations between the orthonormal comoving basis and the orthonormal tetrad basis we can obtain specific expressionsfor the four-velocity of the comoving observer, Equation (58). We can also relate the Eulerian projections of the number-fluxfour-vector, the stress-energy tensor, and the third-order moments in terms of the corresponding Lagrangian projections (Section5). In the comoving-frame, the four-velocity of the comoving observer is simply u ˆ µ = ( 1 , . The four-velocity of the comovingobserver in the orthonormal tetrad basis is u ¯ µ = Λ ¯ µ ˆ µ u ˆ µ = W ( 1 , ¯ v ¯ ı ) , (A4)while the coordinate basis four-velocity of the comoving observer is u µ = e µ ¯ µ u ¯ µ = W ( 1 , v i ) , (A5)where v i = e i ¯ ı ¯ v ¯ ı . Note that ¯ v ¯ ı ¯ v ¯ ı = v i v i . Moreover, by contracting Equation (A5) with u µ , we find W = ( 1 − v i v i ) − .The Eulerian projections of the number-flux four-vector can now be expressed in terms of the corresponding Lagrangianprojections as (cf. Equations (59) and (60)) E N = W J N + v i H i N , (A6) F i N = H i N + W v i J N , (A7)where H i N is related to the comoving-frame moments by H i N = e i ¯ ı (cid:16) δ ¯ ı ˆ ı + W W + 1 ¯ v ¯ ı ¯ v ˆ ı (cid:17) H ˆ ı N . (A8)The Eulerian projections of the stress-energy tensor are related to the corresponding Lagrangian projections in a similar manner.In particular, we have (cf. Equations (61) and (62)) E = W J + 2 W v i H i + v i v j K ij , (A9) F i = W H i + W v i J + v j K ij + W v i v j H j , (A10) S ij = K ij + W (cid:0) v i H j + H i v j (cid:1) + W v i v j J . (A11)In Equations (A9)-(A11), the Lagrangian projections H i and K ij are related to the comoving-frame moments by H i = e i ¯ ı (cid:16) δ ¯ ı ˆ ı + W W + 1 ¯ v ¯ ı ¯ v ˆ ı (cid:17) H ˆ ı , (A12) K ij = e i ¯ ı e j ¯ (cid:16) δ ¯ ı ˆ ı δ ¯ ˆ + W W + 1 h ¯ v ¯ ı ¯ v ˆ ı δ ¯ ˆ + δ ¯ ı ˆ ı ¯ v ¯ ¯ v ˆ i + W ( W + 1 ) ¯ v ¯ ı ¯ v ¯ ¯ v ˆ ı ¯ v ˆ (cid:17) K ˆ ı ˆ . (A13)Note that K ˆ ı ˆ = k ˆ ı ˆ J , where the rank-two variable Eddington tensor k ˆ ı ˆ is defined in Equation (26). Finally, we express thethird-order moments in terms of the Lagrangian projections (cf. Equations (65) and (66)) ǫ Q ijk = L ijk + W (cid:16) v i K jk + v j K ik + v k K ij (cid:17) + W (cid:16) v i v j H k + v i v k H j + v j v k H i (cid:17) + W v i v j v k J , (A14)where the Lagrangian projections of the third-order tensor is related to the comoving-frame moments by L ijk = e i ¯ ı e j ¯ e k ¯ k (cid:16) δ ¯ ı ˆ ı δ ¯ ˆ δ ¯ k ˆ k + W W + 1 h ¯ v ¯ ı ¯ v ˆ ı δ ¯ ˆ δ ¯ k ˆ k + ¯ v ¯ ¯ v ˆ δ ¯ ı ˆ ı δ ¯ k ˆ k + ¯ v ¯ k ¯ v ˆ k δ ¯ ı ˆ ı δ ¯ ˆ i + W ( W + 1 ) h ¯ v ¯ ı ¯ v ˆ ı ¯ v ¯ ¯ v ˆ δ ¯ k ˆ k + ¯ v ¯ ı ¯ v ˆ ı ¯ v ¯ k ¯ v ˆ k δ ¯ ˆ + ¯ v ¯ ¯ v ˆ ¯ v ¯ k ¯ v ˆ k δ ¯ ı ˆ ı i + W ( W + 1 ) ¯ v ¯ ı ¯ v ¯ ¯ v ¯ k ¯ v ˆ ı ¯ v ˆ ¯ v ˆ k (cid:17) L ˆ ı ˆ ˆ k , (A15)and L ˆ ı ˆ ˆ k = l ˆ ı ˆ ˆ k J . The rank-three variable Eddington tensor l ˆ ı ˆ ˆ k is defined in Equation (32).Alternatively, we can write the Eulerian projections of the number-flux four-vector in terms of the number-flux four-vector inthe orthonormal tetrad basis, whose components are expressed in terms of comoving-frame angular moments by N ¯0 = W (cid:16) J N + ¯ v ˆ ı H ˆ ı N (cid:17) , (A16) N ¯ ı = δ ¯ ı ˆ ı H ˆ ı N + W ¯ v ¯ ı (cid:16) J N + WW + 1 ¯ v ˆ ı H ˆ ı N (cid:17) . (A17)0 Endeve et al.Similarly, The monochromatic radiation energy density, momentum density, and stress in the orthonormal tetrad basis can beexpressed in terms of the comoving-frame moments as T ¯0¯0 = W (cid:16) J + 2 ¯ v ˆ ı H ˆ ı + ¯ v ˆ ı ¯ v ˆ K ˆ ı ˆ (cid:17) , (A18) T ¯ ı ¯0 = W h δ ¯ ı ˆ ı H ˆ ı + W ¯ v ¯ ı J + δ ¯ ı ˆ ı ¯ v ˆ K ˆ ı ˆ + WW + 1 ¯ v ¯ ı (cid:16) [ 2 W + 1 ] ¯ v ˆ ı H ˆ ı + W ¯ v ˆ ı ¯ v ˆ K ˆ ı ˆ (cid:17)i , (A19) T ¯ ı ¯ = δ ¯ ı ˆ ı δ ¯ ˆ K ˆ ı ˆ + W (cid:16) ¯ v ¯ ı δ ¯ ˆ H ˆ + ¯ v ¯ δ ¯ ı ˆ ı H ˆ ı (cid:17) + W ¯ v ¯ ı ¯ v ¯ J + W W + 1 (cid:16)h ¯ v ¯ ı ¯ v ˆ ı δ ¯ ˆ + δ ¯ ı ˆ ı ¯ v ¯ ¯ v ˆ i K ˆ ı ˆ + 2 W ¯ v ¯ ı ¯ v ¯ ¯ v ˆ ı H ˆ ı (cid:17) + W ( W + 1 ) ¯ v ¯ ı ¯ v ¯ ¯ v ˆ ı ¯ v ˆ K ˆ ı ˆ . (A20)Equations (A18) and (A19) correspond to Equations (91.10) and (91.11) in Mihalas & Mihalas (1999) (see also Equations (182)and (183) in Munier & Weaver 1986a). Equation (A20) corresponds to Equation (91.12) in Mihalas & Mihalas (1999) (see alsoEquation (184) in Munier & Weaver 1986a). We also have ǫ U ¯ ı ¯ ¯ k = δ ¯ ı ˆ ı δ ¯ ˆ δ ¯ k ˆ k L ˆ ı ˆ ˆ k + W (cid:16) ¯ v ¯ k δ ¯ ı ˆ ı δ ¯ ˆ K ˆ ı ˆ + ¯ v ¯ δ ¯ ı ˆ ı δ ¯ k ˆ k K ˆ ı ˆ k + ¯ v ¯ ı δ ¯ ˆ δ ¯ k ˆ k K ˆ ˆ k (cid:17) + W (cid:16) ¯ v ¯ ¯ v ¯ k δ ¯ ı ˆ ı H ˆ ı + ¯ v ¯ ı ¯ v ¯ k δ ¯ ˆ H ˆ + ¯ v ¯ ı ¯ v ¯ δ ¯ k ˆ k H ˆ k (cid:17) + W W + 1 (cid:16) δ ¯ ı ˆ ı δ ¯ ˆ ¯ v ¯ k ¯ v ˆ k + δ ¯ ı ˆ ı δ ¯ k ˆ k ¯ v ¯ ¯ v ˆ + δ ¯ ˆ δ ¯ k ˆ k ¯ v ¯ ı ¯ v ˆ ı (cid:17) L ˆ ı ˆ ˆ k + W ¯ v ¯ ı ¯ v ¯ ¯ v ¯ k J + W W + 1 h(cid:16) δ ¯ ı ˆ ı ¯ v ¯ ¯ v ˆ + ¯ v ¯ ı ¯ v ˆ ı δ ¯ ˆ (cid:17) ¯ v ¯ k K ˆ ı ˆ + (cid:16) δ ¯ ı ˆ ı ¯ v ¯ k ¯ v ˆ k + ¯ v ¯ ı ¯ v ˆ ı δ ¯ k ˆ k (cid:17) ¯ v ¯ K ˆ ı ˆ k + (cid:16) δ ¯ ˆ ¯ v ¯ k ¯ v ˆ k + ¯ v ¯ ¯ v ˆ δ ¯ k ˆ k (cid:17) ¯ v ¯ ı K ˆ ˆ k i + 3 W W + 1 ¯ v ¯ ı ¯ v ¯ ¯ v ¯ k ¯ v ˆ ı H ˆ ı + W ( W + 1) (cid:16) δ ¯ ı ˆ ı ¯ v ¯ ¯ v ˆ ¯ v ¯ k ¯ v ˆ k + δ ¯ ˆ ¯ v ¯ ı ¯ v ˆ ı ¯ v ¯ k ¯ v ˆ k + δ ¯ k ˆ k ¯ v ¯ ı ¯ v ˆ ı ¯ v ¯ ¯ v ˆ (cid:17) L ˆ ı ˆ ˆ k + 3 W (1 + W ) ¯ v ¯ ı ¯ v ¯ ¯ v ¯ k ¯ v ˆ ı ¯ v ˆ K ˆ ı ˆ + W ( W + 1) ¯ v ¯ ı ¯ v ¯ ¯ v ¯ k ¯ v ˆ ı ¯ v ˆ ¯ v ˆ k L ˆ ı ˆ ˆ k . (A21)We can now list the conservative special relativistic radiation moment equations.The conservative, monochromatic special relativistic radiation moment equations are obtained from Equations (67) and (72)by setting α = ψ = 1 and β i = 0 . The radiation energy equation becomes ∂ E ∂t + 1 √ γ ∂∂x i (cid:16) √ γ F i (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ F ǫ (cid:17) = 1 ǫ Z Ω p C [ f ] d Ω , (A22)where the energy space energy flux has been defined as F ǫ = ǫ n(cid:16) F i + S ji v j + W ǫ Q jki v j v k − W E v i (cid:17) ∂v i ∂t + (cid:16) S ji + W ǫ Q jki v k − W v i F j (cid:17) ∂v i ∂x j + 12 (cid:16) S ij + W ǫ Q ijk v k (cid:17) v l ∂γ ij ∂x l o . (A23)We have written ∂ ln W/∂t = W v i ∂v i /∂t and ∂ ln W/∂x j = W v i ∂v i /∂x j . The terms proportional to derivatives of three-velocity components account for Doppler shifts of the radiation energy spectrum as measured by the comoving observer. Theradiation momentum equation becomes ∂ F i ∂t + 1 √ γ ∂∂x j (cid:16) √ γ S ji (cid:17) − S jk ∂γ jk ∂x i − ǫ ∂∂ǫ (cid:16) ǫ S ǫi (cid:17) = γ ij ǫ Z Ω p j C [ f ] d Ω , (A24)where we have defined the energy space momentum flux as S ǫi = ǫ n(cid:16) S ij + W ǫ Q kij v k − W F i v j (cid:17) ∂v j ∂t + (cid:16) W ǫ Q kij − W S ki v j (cid:17) ∂v j ∂x k + 12 W ǫ Q jki v l ∂γ jk ∂x l o . (A25)We specialize the monochromatic radiation energy and momentum equations to spherical polar coordinates; i.e., ( x , x , x ) =( r, θ, φ ) , a = b = r , and c = sin θ , and impose spherical symmetry ( ∂/∂θ = ∂/∂φ = 0 ). The only nonzero component of thefluid three-velocity is v = v ≡ v r . In spherical symmetry, the comoving-frame angular moments are { J , H , K , L } = { J , H ˆ1 , K ˆ1ˆ1 , L ˆ1ˆ1ˆ1 } = 2 π ǫ Z − f µ { , , , } dµ, (A26)ONSERVATIVE RADIATION MOMENT EQUATIONS 21with µ = cos ϑ . Then, the monochromatic radiation energy and momentum equations become ∂ E ∂t + 1 r ∂∂r (cid:16) r F (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ F ǫ (cid:17) = 2 π W (cid:16) Z − C [ f ] dµ + v r Z − C [ f ] µ dµ (cid:17) , (A27)and ∂ F ∂t + 1 r ∂∂r (cid:16) r S (cid:17) − r (cid:16) J − K (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ S ǫ (cid:17) = 2 π W (cid:16) v r Z − C [ f ] dµ + Z − C [ f ] µ dµ (cid:17) , (A28)respectively. With spherical symmetry imposed, we have E = W (cid:0) J + 2 v r H + v r K (cid:1) , (A29) F = W (cid:0) [ 1 + v r ] H + v r [ J + K ] (cid:1) , (A30) S = W (cid:0) K + 2 v r H + v r J (cid:1) . (A31)Moreover, we have S = K = K ˆ2ˆ2 = ( J − K ) and S = K = K ˆ3ˆ3 = ( J − K ) . The fluxes inside the energyderivatives can be written explicitly in terms of the comoving-frame moments as F ǫ = ǫ W n W h ( H + v r K ) + v r ( K + v r L ) i ∂v r ∂t + W h v r ( H + v r K ) + ( K + v r L ) i ∂v r ∂r + h ( J − K ) + v r ( H − L ) i v r r o , (A32) S ǫ = ǫ W n W h ( v r H + K ) + v r ( v r K + L ) i ∂v r ∂t + W h v r ( v r H + K ) + ( v r K + L ) i ∂v r ∂r + h v r ( J − K ) + (
H − L ) i v r r o , (A33)where we have used ǫ Q = W (cid:0) L + 3 v r K + 3 v r H + v r J (cid:1) , (A34) ǫ Q = 1 ǫ Q = 12 W ( [ H − L ] + v r [ J − K ] ) . (A35)Equations (A27) and (A28) are monochromatic lab-frame radiation energy and momentum equations in conservative form—validto all orders of v r . They can be expressed explicitly in terms of the comoving-frame moments, and compared with corresponding non-conservative equations in Mihalas & Mihalas (1999); their Equations (95.11) and (95.12). Equation (A27) is equivalentto Equation (95.11) plus v r times Equation (95.12) in Mihalas & Mihalas (1999). Similarly, Equation (A28) is equivalent toEquation (95.12) plus v r times Equation (95.11) in Mihalas & Mihalas (1999).The special relativistic moment equations are consistent with the conservative neutrino number equation, which we obtain byadding W ǫ − times Equation (A22) and − W ǫ − v i contracted with Equation (A24). An intermediate result is ∂∂t (cid:16) ǫ W h E − v i F i i(cid:17) + 1 √ γ ∂∂x i (cid:16) √ γ ǫ W h F i − v j S ij i(cid:17) − ǫ ∂∂ǫ (cid:16) ǫ ǫ W h F ǫ − v i S ǫi i(cid:17) + 1 ǫ W nh F i − W (cid:16) E v i − v j F j v i (cid:17)i ∂v i ∂t + h S ji − W v i (cid:16) F j − S jk v k (cid:17)i ∂v i ∂x j + S ij v k ∂γ ij ∂x k o − ǫ W (cid:16) F ǫ − v i S ǫi (cid:17) = 1 ǫ Z Ω C [ f ] d Ω . (A36)(Note that u µ p µ = − ǫ .) From Equations (A23) and (A25) we have ǫ (cid:16) F ǫ − v i S ǫi (cid:17) = h F i − W (cid:16) E v i − v j F j v i (cid:17)i ∂v i ∂t + h S ji − W v i (cid:16) F j − S jk v k (cid:17)i ∂v i ∂x j + S ij v k ∂γ ij ∂x k . (A37)Thus, by virtue of the cancellation of the terms on the second line with the last term on the left-hand side, Equation (A36) reducesto the conservative monochromatic number equation ∂ E N ∂t + 1 √ γ ∂∂x i (cid:16) √ γ F i N (cid:17) − ǫ ∂∂ǫ (cid:16) ǫ F ǫ N (cid:17) = 1 ǫ Z Ω C [ f ] d Ω , (A38)where the energy space number flux is F ǫ N = W n(cid:16) F i − W h E − F j v j i v i (cid:17) ∂v i ∂t + (cid:16) S ji − W h F j − S jk v k i v i (cid:17) ∂v i ∂x j + S ij v k ∂γ ij ∂x k o . (A39)2 Endeve et al.Equation (A38) is a conservative equation for the neutrino number density. In the absence of neutrino-matter interactions, itstates that the neutrino number is conserved. In carrying out the derivation of the number equation from the radiation energyand momentum equations, we observe that the leftover terms emanating from bringing W and W v i inside the space and timederivatives in the energy and momentum equations cancel with the leftover terms from bringing ǫ − inside the energy derivatives.The remaining terms constitute the left-hand side of the conservative number equation (Equation (A38)). These cancellations maybe key to constructing a numerical scheme for neutrino transport based on the solution of the conservative energy and momentumequations that is also consistent with the conservative number equation. Cancellations similar to those occurring in the continuumderivation of the number equation form the energy and momentum equations should also occur in the discrete limit in order toensure consistency with neutrino number conservation. Such consistency may help ensure lepton conservation in simulations ofneutrino radiation hydrodynamics. GENERAL RELATIVISTIC MOMENT EQUATIONS FOR SPHERICALLY SYMMETRIC SPACETIMES
M¨uller et al. (2010, hereafter referred to as MJD10) have recently presented numerical methods for multidimensional, generalrelativistic neutrino radiation hydrodynamics for the case where the conformal flatness condition (CFC) is imposed on the spatialmetric (e.g., Wilson et al. 1996). MJD10 employ the so-called ray-by-ray approach to multidimensional neutrino transport,where the radiation flux is assumed to be purely radial, and the radiation field is essentially obtained by solving independentspherically symmetric problems along each radial ray. (Lateral advection of neutrinos in optically thick regions and non-radialcomponents of the radiation pressure gradient are still taken into account.) The radiation moment equations for the conformallyflat spacetime become fully general relativistic when spherical symmetry is imposed. For the sake of comparing the conservativemoment equations presented here with the equations solved by MJD10, we adopt a spherically symmetric spacetime, and list—infull detail—general relativistic moment equations. The equations presented here can be obtained directly from Equations (67),(72), and (83).The invariant spacetime line element can be decomposed as ds = − α dt + γ ij (cid:0) dx i + β i dt (cid:1) (cid:0) dx j + β j dt (cid:1) . (B1)We adopt spherical polar coordinates ( x , x , x ) = ( r, θ, φ ) . With spherical symmetry imposed, non-spectral quantities areonly functions of time t and the radial coordinate r . The only nonzero component of the fluid three-velocity in the orthonormaltetrad basis is ¯ v ¯1 = ¯ v ¯1 ≡ v r . For the CFC spacetime, the spatial metric γ ij = ψ ¯ γ ij is diagonal, and the conformal metric is ¯ γ ij = diag [ 1 , r , r sin θ ] , and ψ is the conformal factor. Moreover, α is the lapse function and β i is the shift vector (only theradial component of the shift vector β r is nonzero when spherical symmetry is imposed). The lapse function, the shift vector, andthe conformal factor can be obtained by solving a system of nonlinear elliptic equations (Wilson et al. 1996). For the sphericallysymmetric spacetime, the determinant of the spatial metric is √ γ = ψ r sin θ , and we can set e µ ¯0 = α − [ 1 , − β r , , , e ı = [ 0 , , , and e i ¯ ı = ψ − diag [ 1 , r − , ( r sin θ ) − ] .The monochromatic lab-frame radiation energy and momentum equations, expressed in terms of the comoving-frame angularmoments J , H , K , and L , become α ∂∂t (cid:16) W h ˆ J + 2 v r ˆ H + v r ˆ K i(cid:17) + 1 α ∂∂r (cid:16) αW nh(cid:16) ψ − v r β r α (cid:17) + v r (cid:16) v r ψ − β r α (cid:17)i ˆ H + h v r ψ − β r α i ˆ J + v r h ψ − v r β r α i ˆ K o(cid:17) + 1 ψ ∂ ln α∂r W h(cid:16) v r (cid:17) ˆ H + v r (cid:16) ˆ J + ˆ K (cid:17)i + (cid:16) ∂ ln ψ ∂τ − α ∂β r ∂r (cid:17) W h ˆ K + 2 v r ˆ H + v r ˆ J i + (cid:16) ∂ ln ψ ∂τ − r β r α (cid:17)h ˆ J − ˆ K i − ǫ ∂∂ǫ (cid:16) ǫ W nh ψ ∂ ln α∂r + v r (cid:16) ∂ ln ψ ∂τ − α ∂β r ∂r (cid:17) + W Dv r Dτ i(cid:16) ˆ H + v r ˆ K (cid:17) + h v r ψ ∂ ln α∂r + (cid:16) ∂ ln ψ ∂τ − α ∂β r ∂r (cid:17) + W (cid:16) v r ∂v r ∂τ + 1 ψ ∂v r ∂r (cid:17)i(cid:16) ˆ K + v r ˆ L (cid:17) + h D ln ψ Dτ + 1 r (cid:16) v r ψ − β r α (cid:17)i(cid:16)h ˆ J − ˆ K i + v r h ˆ H − ˆ L i(cid:17)o(cid:17) = W (cid:16) π Z − ˆ C [ f ] dµ + v r π Z − ˆ C [ f ] µ dµ (cid:17) , (B2)ONSERVATIVE RADIATION MOMENT EQUATIONS 23and α ∂∂t (cid:16) W h(cid:16) v r (cid:17) ˆ H + v r (cid:16) ˆ J + ˆ K (cid:17)i(cid:17) + 1 α ∂∂r (cid:16) αW nh ψ − v r β r α i ˆ K + h(cid:16) v r ψ − β r α (cid:17) + v r (cid:16) ψ − v r β r α (cid:17)i ˆ H + v r h v r ψ − β r α i ˆ J o(cid:17) + 1 ψ ∂ ln α∂r W h ˆ J + 2 v r ˆ H + v r ˆ K i + (cid:16) ∂ ln ψ ∂τ − α ∂β r ∂r (cid:17) W h(cid:16) v r (cid:17) ˆ H + v r (cid:16) ˆ J + ˆ K (cid:17)i − ψ (cid:16) ∂ ln ψ ∂r + 1 r (cid:17)h ˆ J − ˆ K i − ǫ ∂∂ǫ (cid:16) ǫ W nh ψ ∂ ln α∂r + v r (cid:16) ∂ ln ψ ∂τ − α ∂β r ∂r (cid:17) + W Dv r Dτ i(cid:16) v r ˆ H + ˆ K (cid:17) + h v r ψ ∂ ln α∂r + (cid:16) ∂ ln ψ ∂τ − α ∂β r ∂r (cid:17) + W (cid:16) v r ∂v r ∂τ + 1 ψ ∂v r ∂r (cid:17)i(cid:16) v r ˆ K + ˆ L (cid:17) + h D ln ψ Dτ + 1 r (cid:16) v r ψ − β r α (cid:17)i(cid:16) v r h ˆ J − ˆ K i + h ˆ H − ˆ L i(cid:17)o(cid:17) = W (cid:16) v r π Z − ˆ C [ f ] dµ + 2 π Z − ˆ C [ f ] µ dµ (cid:17) , (B3)respectively. Equations (B2) and (B3) are conservative evolution equations for the monochromatic lab-frame radiation energydensity and momentum density, respectively. The square root of the spatial metric determinant has been absorbed in quantitiesaccented with a hat (e.g., ˆ J = √ γ J ; cf. MJD10). The expression inside the time derivative in Equation (B2) is the contributionfrom radiation to the “matter sources” in the definition of the (conserved) ADM mass (cf. MJD10). For flat spacetime (i.e., α = ψ = 1 and β r = 0 ), Equations (B2) and (B3) reduce to Equations (A27) and (A28), respectively. In Equations (B2) and(B3), we have defined the “convective derivative” and the “proper time derivative” along constant coordinate lines; DDτ = ∂∂τ + v r ψ ∂∂r and ∂∂τ = 1 α ∂∂t − β r α ∂∂r , (B4)respectively.The conservative, monochromatic lab-frame number equation, also expressed in terms of comoving-frame angular moments,can be obtained directly from Equation (83), or by adding W ǫ − times Equation (B2) and − v r W ǫ − times Equation (B3). Theresult is α ∂∂t (cid:16) W h ˆ J N + v r ˆ H N i(cid:17) + 1 α ∂∂r (cid:16) αW h(cid:16) ψ − v r β r α (cid:17) ˆ H N + (cid:16) v r ψ − β r α (cid:17) ˆ J N i(cid:17) − ǫ ∂∂ǫ (cid:16) ǫ W nh ψ ∂ ln α∂r + v r (cid:16) ∂ ln ψ ∂τ − α ∂β r ∂r (cid:17) + W Dv r Dτ i ˆ H + h v r ψ ∂ ln α∂r + (cid:16) ∂ ln ψ ∂τ − α ∂β r ∂r (cid:17) + W (cid:16) v r ∂v r ∂τ + 1 ψ ∂v r ∂r (cid:17)i ˆ K + h D ln ψ Dτ + 1 r (cid:16) v r ψ − β r α (cid:17)i(cid:16) ˆ J − ˆ K (cid:17)o(cid:17) = 2 πǫ Z − ˆ C [ f ] dµ, (B5)where the comoving-frame number density and number flux are (cid:8) J N , H N (cid:9) = 2 π Z − f µ { , } dµ. (B6)Note that { J , H } = ǫ { J N , H N } (cf. Equation (A26)). Equation (B5) corresponds to the evolution equation for the neutrinonumber density given by MJD10 (their Equation (30); corrected here for a few misprints). For flat spacetime, Equation (B5)reduces to the corresponding special relativistic equation (cf. Equation (A38)).MJD10 (cf. their Appendix B) have taken important initial steps in developing numerical methods for neutrino transport basedon moments models where simultaneous conservation of energy and lepton number is considered (see Liebend¨orfer et al. 2004,for the Boltzmann case in spherical symmetry). There may, however, be room for further improvements. Equations (B2) and(B3) differ from the energy and momentum equations solved by MJD10 (their Equations (27) and (28), respectively). We obtaintheir energy equation by adding W times Equation (B2) and − v r W times Equation (B3). Similarly, we obtain their momentumequation by adding − v r W times Equation (B2) and W times Equation (B3). The relationship between the equation for theneutrino number density and the neutrino energy density in MJD10 (their Equations (30) and (27), respectively) is relativelysimple: they differ only by a factor ǫ . Note that there is a similarly simple relationship between Equation (52) and the conservativenumber equation (Equation (41)). This simple relationship is convenient when constructing a numerical method based on theradiation energy and momentum equations, that is simultaneously consistent with the conservative neutrino number equation(which ensures number conservation). However, the energy and momentum equations solved by MJD10 (obtainable directlyfrom Equation (52)) are formulated in a non-conservative form. The “source terms” depend on time and space derivatives of the4 Endeve et al.fluid three-velocity, which do not vanish in the absence of gravity and neutrino-matter interactions. Conservation of energy isnot guaranteed when the non-conservative formulation of the energy equation is used. (Additional complications may also arisein the presence of shocks.) As a possible improvement, a numerical scheme that simultaneously conserves energy and leptonnumber may be based on solving Equations (B2) and (B3), which are conservative formulations of the radiation energy andmomentum equations, respectively. However, such a scheme is potentially much more complicated than the algorithm outlinedin Appendix B in MJD10. When carrying out the the steps for obtaining the conservative number equation from the conservativeenergy and momentum equations, we observe that terms emanating from the time derivatives, space derivatives, and the geometrysources cancel with terms emanating from the energy derivatives. The remaining terms constitute the left-hand side of the numberequation. Similar cancellations must occur in the discrete limit in order for the discretized energy and momentum to be consistentwith neutrino number conservation (which is necessary to ensure lepton number conservation).. Note that there is a similarly simple relationship between Equation (52) and the conservativenumber equation (Equation (41)). This simple relationship is convenient when constructing a numerical method based on theradiation energy and momentum equations, that is simultaneously consistent with the conservative neutrino number equation(which ensures number conservation). However, the energy and momentum equations solved by MJD10 (obtainable directlyfrom Equation (52)) are formulated in a non-conservative form. The “source terms” depend on time and space derivatives of the4 Endeve et al.fluid three-velocity, which do not vanish in the absence of gravity and neutrino-matter interactions. Conservation of energy isnot guaranteed when the non-conservative formulation of the energy equation is used. (Additional complications may also arisein the presence of shocks.) As a possible improvement, a numerical scheme that simultaneously conserves energy and leptonnumber may be based on solving Equations (B2) and (B3), which are conservative formulations of the radiation energy andmomentum equations, respectively. However, such a scheme is potentially much more complicated than the algorithm outlinedin Appendix B in MJD10. When carrying out the the steps for obtaining the conservative number equation from the conservativeenergy and momentum equations, we observe that terms emanating from the time derivatives, space derivatives, and the geometrysources cancel with terms emanating from the energy derivatives. The remaining terms constitute the left-hand side of the numberequation. Similar cancellations must occur in the discrete limit in order for the discretized energy and momentum to be consistentwith neutrino number conservation (which is necessary to ensure lepton number conservation).