diff --git a/Docs/source/index.rst b/Docs/source/index.rst
index c0e9fc762..b4ca4e7f1 100644
--- a/Docs/source/index.rst
+++ b/Docs/source/index.rst
@@ -60,6 +60,7 @@ of state routines.
    networks
    templated_networks
    screening
+   neutrinos
 
 .. toctree::
    :maxdepth: 1
@@ -71,6 +72,13 @@ of state routines.
    nse
    sdc
 
+.. toctree::
+   :maxdepth: 1
+   :caption: Util / external libraries
+   :hidden:
+
+   util
+
 .. toctree::
    :maxdepth: 1
    :caption: Unit tests
diff --git a/Docs/source/integrators.rst b/Docs/source/integrators.rst
index 80fa229d0..5d110b9bb 100644
--- a/Docs/source/integrators.rst
+++ b/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::
 
diff --git a/Docs/source/networks.rst b/Docs/source/networks.rst
index 660054d3a..908a6ee20 100644
--- a/Docs/source/networks.rst
+++ b/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
diff --git a/Docs/source/neutrinos.rst b/Docs/source/neutrinos.rst
new file mode 100644
index 000000000..28deab6ab
--- /dev/null
+++ b/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}$
diff --git a/Docs/source/refs.bib b/Docs/source/refs.bib
index 9e8cd982b..67eb16774 100644
--- a/Docs/source/refs.bib
+++ b/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}
+}
diff --git a/Docs/source/sdc.rst b/Docs/source/sdc.rst
index 38821a2a1..1d9b52e29 100644
--- a/Docs/source/sdc.rst
+++ b/Docs/source/sdc.rst
@@ -147,7 +147,7 @@ The reactions don’t modify the total density, :math:`\rho`, or momentum,
       \Adv{\rho X_k}^{n+1/2} \\ \Adv{\rho e}^{n+1/2} \\
    \end{array} \right ) +
    \left (
-      \begin{array}{c} \rho \omegadot_k \\ \rho \Sdot \end{array}
+      \begin{array}{c} \rho \omegadot_k \\ \rho \epsilon \end{array}
    \right )
 
 Here the advective terms are piecewise-constant (in time)
diff --git a/Docs/source/util.rst b/Docs/source/util.rst
new file mode 100644
index 000000000..1ac479819
--- /dev/null
+++ b/Docs/source/util.rst
@@ -0,0 +1,46 @@
+******************************
+Helper Functions and Libraries
+******************************
+
+The ``util/`` directory contains a number of external libraries and simple
+utilities that are used by the different components of Microphysics.
+
+* ``approx_math/`` : these are a set of headers that implement
+  approximations to ``atan()``, ``exp()``, ``log()``, and ``pow()``.
+  These can be much faster than the C++ library versions, especially
+  on GPUs.
+
+* ``autodiff/`` : this is a clone of the C++ autodiff library from
+  https://github.com/autodiff/autodiff
+
+  The header ``microphysics_autodiff.H`` provides a set of interfaces
+  for working with the AMReX datatypes and interfacing with the
+  autodiff library.
+
+* ``build_scripts/`` : a set of python scripts used during the build
+  process to parse the runtime parameters.
+
+* ``cj_detonation/`` : a simple routine to compute the Chapman-Jouguet
+  detonation speed for one of our networks.
+
+* ``esum.H`` : an implementation of the exact sum algorithm based on the
+  msum algorithm by Raymond Hettinger.  It is generated automatically
+  by the ``esum_cxx.py`` script and creates implementations for exact
+  numbers of terms (``esum3()``, ``esum4()``, ...)
+
+* ``gcem/`` : a templated math library that provides implementations of
+  the standard library math functions that can be used in ``constexpr``
+  expressions.  This is from https://github.com/kthohr/gcem
+
+  Some of the constants are redefined in 64-bit floating point in
+  ``microphysics_math.H`` to avoid ``long double`` issues on some
+  architectures.
+
+* ``hybrj/`` : a C++ port of the MINPACK hybrid Powell minimization function
+  to zero a set of functions.
+
+* ``linpack.H`` : a C++ port of the LINPACK ``dgesl`` and ``dgefa`` LU
+  decomposition Gaussian elimination routines.
+
+* ``microphysics_sort.H`` : a set of sorting routines for
+  ``amrex::Array1D`` data.
diff --git a/EOS/breakout/actual_eos.H b/EOS/breakout/actual_eos.H
index ac0ba24bc..a4b511388 100644
--- a/EOS/breakout/actual_eos.H
+++ b/EOS/breakout/actual_eos.H
@@ -21,6 +21,10 @@ void actual_eos_init ()
         gamma_const = 5.0_rt / 3.0_rt;
     }
 
+    // this EOS assumes that has a 1/mu entry -- make sure that is valid
+    if (aux_names_cxx[AuxZero::iinvmu] != "invmu") {
+        amrex::Error("invalid aux state for breakout EOS");
+    }
 }
 
 template <typename I>
@@ -52,7 +56,7 @@ void actual_eos (I input, T& state)
 
     // Calculate mu. This is the only difference between
     // this EOS and gamma_law.
-    state.mu = 1.0_rt / state.aux[1];
+    state.mu = 1.0_rt / state.aux[AuxZero::iinvmu];
 
     switch (input) {