We carry around three different 1D radial quantities: \(\etarhoec\) (edge-centered), \(\etarhocc\) (cell-centered), and \(\divetarho\) (cell-centered). These notes discuss when each of these is used, and how they are computed, in both plane-parallel and spherical.

# The Mixing Term, \(\etarho\)

The base state evolves in response to heating and mixing in the star. The density evolution is governed by

with

designed to keep the average value of the full density, \(\rho\), over a layer of constant radius in the star equal to \(\rho_0\). To complete the update of the base state, we need evolution equations for the pressure, \(p_0\), and velocity, \(w_0\). For spherical geometry, the derivation of \(w_0\) constraint equation is shown in the multilevel paper, resulting in the following system

In paper III, we introduced a mixing term, \(\etarho\), to the density evolution equation (Eq.116), with the objective of keeping the base state density equal to the average of the density over a layer. For a spherical base state, it is best to define the average in terms of spherical coordinates,

where \(\int_{\Omega_H} d\Omega = 4\pi\) represents the integral over the spherical \(\theta\) and \(\phi\) angles at constant radius.

Recall from Paper III that if initially \(\overline{\rho'} = 0\), there is no guarantee that \(\overline{\rho'} = 0\) will hold at later time. To see this, start with the equation for the perturbational density

written in a slightly different form:

We integrate this over a spherical shell of thickness \(2h\) at radius \(r_0\), i.e., \(\Omega_H \times (r_0-h, r_0+h)\), and normalize by the integration volume, which is \(\sim 4\pi r_0^2 2h\) for small \(h\), to obtain:

where we have used the divergence theorem in spherical coordinates to transform the volume integral on the right hand side into an area integral over \(\Omega_H\). We see that the right hand side disappears since \(\int_{\Omega_H}\rho_0(\Ubt\cdot\eb_r)d\Omega=0\), which follows from the definition of \(\Ubt\). Now, expanding the remaining terms and taking the limit as \(h\rightarrow 0,\) we can write

again using the divergence theorem, extracting the time derivative from the spatial integral, and switching the order of operations as appropriate.

In short,

and thus, \(\etarho = \overline{\left(\rho'\Ub\cdot\eb_r\right)}\).

We need both \(\etarho\) alone and its divergence for the various terms in the construction of \(w_0\) and the correction to \(\rho_0\). The quantity \(\etarho\) is edge-centered on our grid, and for Cartesian geometries, we constructed it by averaging the appropriate fluxes through the grid boundaries. For a spherical base state, this does not work, since the spherical shells do not line up with the Cartesian grid boundaries.

Therefore, we take a different approach. We compute \(\etarho\) by constructing the quantity \(\rhop \Ubt \cdot \er\) in each cell, and then use our average routine to construct a 1-d, cell-centered \(\eta_{\rho,r}\) (this is essentially numerically solving the integral in Eq.117. The edge-centered values of \(\etarho\), \(\eta_{\rho,r+1/2}\) are then constructed by simple averaging:

Instead of differencing \(\eta_{\rho,r+1/2}\) to construct the divergence, we instead use Eq.124 directly, by writing:

where we have made use of the fact that \(\overline{\rhop^n} = 0\) by construction.

# \(\eta\) Flow Chart

Enter

`advance_timestep`

with \([\etarhoec, \etarhocc]^{n-\myhalf}\).Call

`make_w0`

. The spherical version uses uses \(\etarho^{{\rm ec},n-\myhalf}\) and \(\etarho^{{\rm cc},n-\myhalf}\).Call

`density_advance`

. The plane-parallel version computes \(\etarho^{{\rm flux},n+\myhalf,*}\).Call

`make_etarho`

to compute \([\etarhoec, \etarhocc]^{n+\myhalf,*}\). The plane-parallel version uses \(\etarho^{{\rm flux},n+\myhalf,*}\).Call

`make_psi`

. The plane-parallel version uses \(\etarho^{{\rm cc},n+\myhalf,*}\).Call

`make_w0`

. The spherical version uses uses \(\etarho^{{\rm ec},n+\myhalf,*}\) and \(\etarho^{{\rm cc},n+\myhalf,*}\).Call

`density_advance`

. The plane-parallel version computes \(\etarho^{{\rm flux},n+\myhalf}\).Call

`make_etarho`

to compute \([\etarhoec, \etarhocc]^{n+\myhalf}\). The plane-parallel version uses \(\etarho^{{\rm flux},n+\myhalf}\).Call

`make_psi`

. The plane-parallel version uses \(\etarho^{{\rm cc},n+\myhalf}\).

# Computing \(\etarhoec\) and \(\etarhocc\)

This is done in `make_eta.f90`

.

## Plane-Parallel

We first compute a radial edge-centered multifab, \(\etarho^{\rm flux}\), using

\(\etarhoec\) is the edge-centered “average” value of \(\eta_{\rho}^{\rm flux}\),

\(\etarhocc\) is a cell-centered average of \(\etarhoec\),

## Spherical

First, construct \(\eta_{\rho}^{\rm cart} = [\rho'(\Ubt\cdot\eb_r)]^{n+\myhalf}\) using:

Then, \(\etarhocc\) is the cell-centered average of \(\eta_{\rho}^{\rm cart}\),

On interior faces, \(\etarhoec\) is the average of \(\etarhocc\),

At the upper and lower boundaries, we use

# Using \(\etarhoec\)

## Plane-Parallel

NOT USED.

## Spherical

In `make_w0`

, \(\etarhoec\) is used in the construction of the RHS
for the \(\delta w_0\) equation.

# Using \(\etarhocc\)

## Plane-Parallel

In `make_psi`

, \(\psi = \etarhocc g\).

## Spherical

In `make_w0`

, \(\etarhocc\) is used in the construction of the RHS
for the \(\delta w_0\) equation.