# Equation of State Unit

## Introduction

The `Eos`

unit implements the equation of state needed by the
hydrodynamics and nuclear burning solvers. The function
`physics/Eos/Eos`

provides the interface for operating on a
one-dimensional vector. The same interface can be used for a single cell
by reducing the vector size to 1. Additionally, this function can be
used to find the thermodynamic quantities either from the density,
temperature, and composition or from the density, internal energy, and
composition. For user’s convenience, a wrapper function
(`physics/Eos/Eos_wrapped`

) is provided, which takes a section of a
block and translates it into the data format required by the
`physics/Eos/Eos`

function, then calls the function. Upon return from
the `physics/Eos/Eos`

function, the wrapper translates the returned
data back to the same section of the block.

Four implementations of the (`Eos`

) unit are available in the current
release of Flash-X: `Gamma`

which implements a perfect-gas equation of
state; `Gamma/RHD`

which implements a perfect-gas equation taking
relativistic effects into account; `Multigamma`

which implements a
perfect-gas equation of state with multiple fluids, each of which can
have its own adiabatic index (\(\gamma\)); and `Helmholtz`

which
uses a fast Helmholtz free-energy table interpolation to handle
degenerate/relativistic electrons/positrons and includes radiation
pressure and ions (via the perfect gas approximation).

As described in previous sections, Flash-X evolves the Euler equations for compressible, inviscid flow. This system of equations must be closed by an additional equation that provides a relation between the thermodynamic quantities of the gas. This relationship is known as the equation of state for the material, and its structure and properties depend on the composition of the gas.

It is common to call an equation of state (henceforth EOS) routine more than \(10^9\) times during a two-dimensional simulation and more than \(10^{11}\) times during the course of a three-dimensional simulation of stellar phenomena. Thus, it is very desirable to have an EOS that is as efficient as possible, yet accurately represents the relevant physics. While Flash-X is capable of using any general equation of state, we discuss here the three equation of state routines that are supplied: an ideal-gas or gamma-law EOS, an EOS for a fluid composed of multiple gamma-law gases, and a tabular Helmholtz free energy EOS appropriate for stellar interiors. The two gamma-law EOSs consist of simple analytic expressions that make for a very fast EOS routine both in the case of a single gas or for a mixture of gases. The Helmholtz EOS includes much more physics and relies on a table look-up scheme for performance.

## Gamma Law and Multigamma

Flash-X uses the method of Colella & Glaz (1985) to handle general equations of state. General equations of state contain 4 adiabatic indices (Chandrasekhar 1939), but the method of Colella & Glaz parameterizes the EOS and requires only two of the adiabatic indices The first is necessary to calculate the adiabatic sound speed and is given by

The second relates the pressure to the energy and is given by

These two adiabatic indices are stored as the mesh-based variables
`GAMC_VAR`

and `GAME_VAR`

. All EOS routines must return
\(\gamma_1\), and \(\gamma_4\) is calculated from
[Eqn:game].

The gamma-law EOS models a simple ideal gas with a constant adiabatic index \(\gamma\). Here we have dropped the subscript on \(\gamma\), because for an ideal gas, all adiabatic indices are equal. The relationship between pressure \(P\), density \(\rho\), and specific internal energy \(\epsilon\) is

We also have an expression relating pressure to the temperature \(T\)

where \(N_a\) is the Avogadro number, \(k\) is the Boltzmann constant, and \(\bar{A}\) is the average atomic mass, defined as

where \(X_i\) is the mass fraction of the \(i\)th element. Equating these expressions for pressure yields an expression for the specific internal energy as a function of temperature

The relativistic variant of the ideal gas equation is explained in more detail in .

Simulations are not restricted to a single ideal gas; the multigamma EOS simulations with several species of ideal gases each with its own value of \(\gamma\). In this case the above expressions hold, but \(\gamma\) represents the weighted average adiabatic index calculated from

We note that the analytic expressions apply to both the forward (internal energy as a function of density, temperature, and composition) and backward (temperature as a function of density, internal energy and composition) relations. Because the backward relation requires no iteration in order to obtain the temperature, this EOS is quite inexpensive to evaluate. Despite its fast performance, use of the gamma-law EOS is limited, due to its restricted range of applicability for astrophysical problems.

## Helmholtz

The Helmholtz EOS provided with the Flash-X distribution contains more physics and is appropriate for addressing astrophysical phenomena in which electrons and positrons may be relativistic and/or degenerate and in which radiation may significantly contribute to the thermodynamic state. Full details of the Helmholtz equation of state are provided in Timmes & Swesty (1999). This EOS includes contributions from radiation, completely ionized nuclei, and degenerate/relativistic electrons and positrons. The pressure and internal energy are calculated as the sum over the components

Here the subscripts “rad,” “ion,” “ele,” “pos,” and “coul” represent the contributions from radiation, nuclei, electrons, positrons, and corrections for Coulomb effects, respectively. The radiation portion assumes a blackbody in local thermodynamic equilibrium, the ion portion (nuclei) is treated as an ideal gas with \(\gamma \, = \, 5/3\), and the electrons and positrons are treated as a non-interacting Fermi gas.

The blackbody pressure and energy are calculated as

where \(a\) is related to the Stephan-Boltzmann constant \(\sigma_B \, = \, a c/4\), and \(c\) is the speed of light. The ion portion of each routine is the ideal gas of () with \(\gamma \, = \, 5/3\). The number densities of free electrons \(N_{\rm ele}\) and positrons \(N_{\rm pos}\) in the noninteracting Fermi gas formalism are given by

where \(h\) is Planck’s constant, \(m_{\rm e}\) is the electron rest mass, \(\beta \: = \: k T / (m_{\rm e} c^2)\) is the relativity parameter, \(\eta \: = \: \mu / k T\) is the normalized chemical potential energy \(\mu\) for electrons, and \(F_{k}(\eta,\beta)\) is the Fermi-Dirac integral

Because the electron rest mass is not included in the chemical potential, the positron chemical potential must have the form \(\eta_{{\rm pos}} \, = \, -\eta - 2/\beta\). For complete ionization, the number density of free electrons in the matter is

and charge neutrality requires

Solving this equation with a standard one-dimensional root-finding
algorithm determines \(\eta\). Once \(\eta\) is known, the
Fermi-Dirac integrals can be evaluated, giving the pressure, specific
thermal energy, and entropy due to the free electrons and positrons.
From these, other thermodynamic quantities such as \(\gamma_1\) and
\(\gamma_4\) are found. Full details of this formalism may be found
in Fryxell *et al.* (2000) and references therein.

The above formalism requires many complex calculations to evaluate the thermodynamic quantities, and routines for these calculations typically are designed for accuracy and thermodynamic consistency at the expense of speed. The Helmholtz EOS in Flash-X provides a table of the Helmholtz free energy (hence the name) and makes use of a thermodynamically consistent interpolation scheme obviating the need to perform the complex calculations required of the above formalism during the course of a simulation. The interpolation scheme uses a bi-quintic Hermite interpolant resulting in an accurate EOS that performs reasonably well.

The Helmholtz free energy,

is the appropriate thermodynamic potential for use when the temperature and density are the natural thermodynamic variables. The free energy table distributed with Flash-X was produced from the Timmes EOS (Timmes & Arnett 1999). The Timmes EOS evaluates the Fermi-Dirac integrals [Eqn:eos7] and their partial derivatives with respect to \(\eta\) and \(\beta\) to machine precision with the efficient quadrature schemes of Aparicio (1998) and uses a Newton-Raphson iteration to obtain the chemical potential of [Eqn:eos9]. All partial derivatives of the pressure, entropy, and internal energy are formed analytically. Searches through the free energy table are avoided by computing hash indices from the values of any given \((T,\rho \bar{Z}/\bar{A})\) pair. No computationally expensive divisions are required in interpolating from the table; all of them can be computed and stored the first time the EOS routine is called.

We note that the Helmholtz free energy table is constructed for only the
electron-positron plasma, and it is a 2-dimensional function of density
and temperature, *i.e.* \(F(\rho,{\rm T})\). It is made with
\({\bar {\rm A}} \, = \, {\bar {\rm Z}} = 1\) (pure hydrogen), with
an electron fraction \(Y_{\rm e} \, = \, 1\). One reason for not
including contributions from photons and ions in the table is that these
components of the Helmholtz EOS are very simple (), and one doesn’t need
fancy table look-up schemes to evaluate simple analytical functions. A
more important reason for only constructing an electron-positron EOS
table with \(Y_{\rm e} \, = \, 1\) is that the 2-dimensional table
is valid for *any* composition. Separate planes for each
\(Y_{\rm e}\) are not necessary (or desirable), since simple
multiplication by \(Y_{\rm
e}\) in the appropriate places gives the desired composition scaling. If
photons and ions were included in the table, then this valuable
composition independence would be lost, and a 3-dimensional table would
be necessary.

The Helmholtz EOS has been subjected to considerable analysis and testing (Timmes & Swesty 2000), and particular care was taken to reduce the numerical error introduced by the thermodynamical models below the formal accuracy of the hydrodynamics algorithm (Fryxell, et al. 2000; Timmes & Swesty 2000). The physical limits of the Helmholtz EOS are \(10^{-10}\,<\,\rho\,<\,10^{11}~({\rm g~cm}^{-3})\) and \(10^{4}\,<\,T\,<\,10^{11}\) (K). As with the gamma-law EOS, the Helmholtz EOS provides both forward and backward relations. In the case of the forward relation (\(\rho, T\), given along with the composition) the table lookup scheme and analytic formulae directly provide relevant thermodynamic quantities. In the case of the backward relation (\(\rho, \epsilon\), and composition given), the routine performs a Newton-Rhaphson iteration to determine temperature. It is possible for the input variables to be changed in the iterative modes since the solution is not exact. The returned quantities are thermodynamically consistent.

## Usage

### Initialization

The initialization function of the Eos unit `physics/Eos/Eos_init`

is
fairly simple for the two ideal gas gamma law implementations included.
It gathers the runtime parameters and the physical constants needed by
the equation of state and stores them in the data module. The Helmholtz
EOS `physics/Eos/Eos_init`

routine is a little more complex. The
`Helmholtz`

EOS requires an input file `helm_table.dat`

that
contains the lookup table for the electron contributions. This table is
currently stored in ASCII for portability purposes. When the table is
first read in, a binary version called `helm_table.bdat`

is created.
This binary format can be used for faster subsequent restarts on the
same machine but may not be portable across platforms. The `Eos_init`

routine reads in the table data on processor 0 and broadcasts it to all
other processors.

### Runtime Parameters

Runtime parameters for the `Gamma`

unit require the user to set the
thermodynamic properties for the single gas. `Eos/gamma`

,
`Eos/eos_singleSpeciesA`

, `Eos/eos_singleSpeciesZ`

set the ratio of
specific heats and the nucleon and proton numbers for the gas. In
contrast, the `Multigamma`

implementation does not set runtime
parameters to define properties of the multiple species. Instead, the
simulation `Config`

file indicates the requested species, for example
helium and oxygen can be defined as

```
SPECIES HE4
SPECIES O16
```

The properties of the gases are initialized in the file
`Simulation/Simulation_initSpecies`

`.F90`

, for example

```
subroutine Simulation_initSpecies()
use Multispecies_interface, ONLY : Multispecies_setProperty
implicit none
#include "Simulation.h"
#include "Multispecies.h"
call Multispecies_setProperty(HE4_SPEC, A, 4.)
call Multispecies_setProperty(HE4_SPEC, Z, 2.)
call Multispecies_setProperty(HE4_SPEC, GAMMA, 1.66666666667e0)
call Multispecies_setProperty(O16_SPEC, A, 16.0)
call Multispecies_setProperty(O16_SPEC, Z, 8.0)
call Multispecies_setProperty(O16_SPEC, GAMMA, 1.4)
end subroutine Simulation_initSpecies
```

For the Helmholtz equation of state, the table-lookup algorithm requires
a given temperature and density. When temperature or internal energy are
supplied as the input parameter, an iterative solution is found.
Therefore, no matter what mode is selected for `Helmholtz`

input, the
best initial value of temperature should be provided to speed
convergence of the iterations. The iterative solver is controlled by two
runtime parameters `Eos/eos_maxNewton`

and `Eos/eos_tolerance`

which
define the maximum number of iterations and convergence tolerance. An
additional runtime parameter for `Helmholtz`

, `Eos/eos_coulumbMult`

,
indicates whether or not to apply Coulomb corrections. In some regions
of the \(\rho\)-\(T\) plane, the approximations made in the
Coulomb corrections may be invalid and result in negative pressures.
When the parameter `eos_coulombMult`

is set to zero, the Coulomb
corrections are not applied.

### Direct and Wrapped Calls

The primary function in the `Eos`

unit operates on a vector, taking
density, composition, and either temperature, internal energy, or
pressure as input, and returning \(\gamma_1\), and either the
pressure, temperature or internal energy (whichever was not used as
input). This equation of state interface is useful for initializing a
problem. The user is given direct control over the input and output,
since everything is passed through the argument list. Also, the vector
data format is more efficient than calling the equation of state routine
directly on a point by point basis, since it permits pipelining and
provides better cache performance. Certain optional quantities such
electron pressure, degeneracy parameter, and thermodynamic derivatives
can be calculated by the `physics/Eos/Eos`

function if needed. These
quantities are selected for computation based upon a logical mask array
provided as an input argument. A .true. value in the mask array results
in the corresponding quantity being computed and reported back to the
calling function. Examples of calling the basic implementation `Eos`

are provided in the API description, see `physics/Eos/Eos`

.

The hydrodynamic and burning computations repeatedly call the Eos
function to update pressure and temperature during the course of their
calculation. Typically, values in all the cells of the block need of be
updated in these calls. Since the primary Eos interface requires the
data to be organized as a vector, using it directly could make the code
in the calling unit very cumbersome and error prone. The wrapper
interface, `physics/Eos/Eos_wrapped`

provides a means by which the
details of translating the data from block to vector and back are hidden
from the calling unit. The wrapper interface permits the caller to
define a section of block by giving the limiting indices along each
dimension. The `Eos_wrapped`

routine translates the block section thus
described into the vector format of the `physics/Eos/Eos`

interface,
and upon return translates the vector format back to the block section.
This wrapper routine cannot calculate the optional derivative
quantities. If they are needed, call the `Eos`

routine directly with
the optional mask set to true and space allocated for the returned
quantities.

## Unit Test

The unit test of the Eos function can exercise all three
implementations. Because the Gamma law allows only one species, the
setup required for the three implementations is specific. To invoke any
three-dimensional `Eos`

unit test, the command is:

`./setup unitTest/Eos/`

implementation`-auto -3d`

where *implementation* is one of `Gamma`

, `Multigamma`

,
`Helmholtz`

. The `Eos`

unit test works on the assumption that if the
four physical variables in question (density, pressure, energy and
temperature) are in thermal equilibrium with one another, then applying
the equation of state to any two of them should leave the other two
completely unchanged. Hence, if we initialize density and temperature
with some arbitrary values, and apply the equation of state to them in
`MODE_DENS_TEMP`

, then we should get pressure and energy values that
are thermodynamically consistent with density and temperature. Now after
saving the original temperature value, we apply the equation of state to
density and newly calculated pressure. The new value of the temperature
should be identical to the saved original value. This verifies that the
`Eos`

unit is computing correctly in `MODE_DENS_PRES`

mode. By
repeating this process for the remaining two modes, we can say with
great confidence that the `Eos`

unit is functioning normally.

In our implementation of the Eos unit test, the initial conditions applied to the domain create a gradient for density along the \(x\) axis and gradients for temperature and pressure along the \(y\) axis. If the test is being run for the Multigamma or Helmholtz implementations, then the species are initialized to have gradients along the \(z\) axis.