Add neutrino docs

Docs/source/index.rst

@@ -60,6 +60,7 @@ of state routines.
    networks
    templated_networks
    screening
+   neutrinos
 
 .. toctree::
    :maxdepth: 1

Docs/source/integrators.rst

@@ -14,13 +14,13 @@ The equations we integrate to do a nuclear burn are:
    :label: eq:spec_integrate
 
 .. math::
-   \frac{de}{dt} = f(\rho,X_k,T)
+   \frac{de}{dt} = \epsilon(\rho,X_k,T)
    :label: eq:enuc_integrate
 
 Here, :math:`X_k` is the mass fraction of species :math:`k`, :math:`e`
 is the specific nuclear energy created through reactions. Also needed
 are density :math:`\rho`, temperature :math:`T`, and the specific
-heat. The function :math:`f` provides the energy release from
+heat. The function :math:`\epsilon` provides the energy release from
 reactions and can often be expressed in terms of the instantaneous
 reaction terms, :math:`\dot{X}_k`. As noted in the previous section,
 this is implemented in a network-specific manner.
@@ -46,7 +46,7 @@ energy. This allows us to easily call the EOS during the burn to obtain the temp
       :label: eq:spec_n_integrate
 
    .. math::
-      \frac{de}{dt} = f(\rho,n_k,T)
+      \frac{de}{dt} = \epsilon(\rho,n_k,T)
       :label: eq:enuc_n_integrate
 
    The effect of this flag in the integrators is that we don't worry
@@ -355,7 +355,7 @@ with the initial temperature, density, and composition:
 As the system is integrated, :math:`e` is updated to account for the
 nuclear energy release (and thermal neutrino losses),
 
-.. math:: e(t) = e_0 + \int_{t_0}^t f(\dot{Y}_k) dt
+.. math:: e(t) = e_0 + \int_{t_0}^t \epsilon(\dot{Y}_k) dt
 
 .. note::
 

Docs/source/networks.rst

@@ -26,10 +26,10 @@ is stored as ``mion(:)`` in the network.
 
 The energy release per gram is converted from the rates as:
 
-.. math:: \edot = -N_A c^2 \sum_k \frac{dY_k}{dt} M_k - \edotnu
+.. math:: \epsilon = -N_A c^2 \sum_k \frac{dY_k}{dt} M_k - \epsilon_\nu
 
 where :math:`N_A` is Avogadro’s number (to convert this to “per gram”)
-and :math:`\edotnu` is the neutrino loss term.
+and :math:`\edotnu` is the neutrino loss term (see :ref:`neutrino_loss`).
 
 
 ``general_null``
@@ -137,7 +137,7 @@ network is interpolated based on the density. It is stored as the
 binding energy (ergs/g) *per nucleon*, with a sign convention that
 binding energies are negative. The energy generation rate is then:
 
-.. math:: \edot = q \frac{dX(\isotm{C}{12})}{dt} = q A_{\isotm{C}{12}} \frac{dY(\isotm{C}{12})}{dt}
+.. math:: \epsilon = q \frac{dX(\isotm{C}{12})}{dt} = q A_{\isotm{C}{12}} \frac{dY(\isotm{C}{12})}{dt}
 
 (this is positive since both :math:`q` and :math:`dY/dt` are negative)
 
@@ -244,7 +244,7 @@ Finally, for the
 energy generation, we take our reaction to release a specific energy,
 :math:`[\mathrm{erg~g^{-1}}]`, of :math:`\qburn`, and our energy source is
 
-.. math:: \edot = -\qburn \frac{dX_f}{dt}
+.. math:: \epsilon = -\qburn \frac{dX_f}{dt}
 
 There are a number of parameters we use to control the constants in
 this network. This is one of the few networks that was designed

Docs/source/neutrinos.rst

@@ -0,0 +1,47 @@
+.. _neutrino_loss:
+
+***************
+Neutrino Losses
+***************
+
+We model thermal neutrino losses (plasma, photo-, pair-,
+recombination, and Bremsstrahlung) using the method described in
+:cite:`itoh:1996`.  This neutrino loss term is included in all of the
+reaction networks by default, and modifies the energy equation to have
+the form (for Strang splitting):
+
+.. math::
+
+   \frac{de}{dt} = \epsilon - \epsilon_\nu
+
+where $\epsilon_\nu$ are the thermal neutrino losses.
+
+.. note::
+
+   Thermal neutrino losses can be disabled at compile time by setting the make
+   variable ``USE_NEUTRINOS = FALSE``.
+
+The interface for the neutrino loss function is:
+
+.. code:: c++
+
+   template <int do_derivatives>
+   AMREX_GPU_HOST_DEVICE AMREX_INLINE
+   void sneut5(const amrex::Real temp, const amrex::Real den,
+               const amrex::Real abar, const amrex::Real zbar,
+               amrex::Real& snu, amrex::Real& dsnudt, amrex::Real& dsnudd,
+               amrex::Real& dsnuda, amrex::Real& dsnudz)
+
+Here, the template parameter, ``do_derivatives``, can be used to disable the code
+the computes the derivatives of the neutrino loss, for example, if a numerical Jacobian
+is used.  The output is
+
+* ``snu`` : $\epsilon_\nu$, the neutrino loss in erg/g/s
+
+* ``dsnudt`` : $d\epsilon_\nu/dT$
+
+* ``dsnudd`` : $d\epsilon_\nu/d\rho$
+
+* ``dsnuda`` : $d\epsilon_\nu/d\bar{A}$
+
+* ``dsnudz`` : $d\epsilon_\nu/d\bar{Z}$

Docs/source/refs.bib

@@ -688,3 +688,16 @@ @article{langanke:2001
 abstract = {The weak interaction rates in stellar environments are computed for pf-shell nuclei in the mass range A=45–65 using large-scale shell-model calculations. The calculated capture and decay rates take into consideration the latest experimental energy levels and log ft-values. The rates are tabulated at the same grid points of density and temperature as those used by Fuller, Fowler, and Newman for densities ρY e =10–1011 g/cm3 and temperatures T=107–1011 K, and hence are relevant for both types of supernovae (Type Ia and Type II). Effective 〈ft〉 values for capture rates and average neutrino (antineutrino) energies are also given to facilitate the use of interpolated rates in stellar evolution codes.}
 }
 
+@ARTICLE{itoh:1996,
+       author = {{Itoh}, Naoki and {Hayashi}, Hiroshi and {Nishikawa}, Akinori and {Kohyama}, Yasuharu},
+        title = "{Neutrino Energy Loss in Stellar Interiors. VII. Pair, Photo-, Plasma, Bremsstrahlung, and Recombination Neutrino Processes}",
+      journal = {\apjs},
+     keywords = {DENSE MATTER, ELEMENTARY PARTICLES, RADIATION MECHANISMS: NONTHERMAL, STARS: INTERIORS, METHODS: NUMERICAL},
+         year = 1996,
+        month = feb,
+       volume = {102},
+        pages = {411},
+          doi = {10.1086/192264},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/1996ApJS..102..411I},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}