We carry around three different 1D radial quantities: ${{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}}$ (edge-centered), ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}$ (cell-centered), and $\mathrm{\nabla }\cdot ({\eta }_{\rho }{\mathbf{e}}_r)$ (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, ${\eta }_{\rho }$

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

(116)$$\frac{\partial {\rho }_0}{\partial t}=-\mathrm{\nabla }\cdot \left({\rho }_0{w}_0{\mathbf{e}}_r\right)-\mathrm{\nabla }\cdot \left({\eta }_{\rho }{\mathbf{e}}_r\right),$$

with

(117)$${\eta }_{\rho }(r)=\overline{({\rho }^{\prime }\stackrel{~}{\mathbf{U}}\cdot {\mathbf{e}}_r)}=\frac{1}{A({\mathrm{\Omega }}_H)}{\int }_{{\mathrm{\Omega }}_H}({\rho }^{\prime }\tilde{\mathbf{U}}\cdot {\mathbf{e}}_r)dA,$$

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

(118)$$\begin{array}{r}\begin{array}{rl}{w}_0& ={\bar{w}}_0+\delta w_0\\ \frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{\bar{w}}_0\right)& =\bar{S}\\ \frac{\partial }{\partial r}[\frac{{\bar{\mathrm{\Gamma }}}_1{p}_0}{r^2}\frac{\partial }{\partial r}(r^2\delta w_0)]& =-\frac{g}{r^2}\frac{\partial (r^2{\eta }_{\rho })}{\partial r}-\frac{4({\bar{w}}_0+\delta w_0){\rho }_{0}g}{r}-4\pi G{\rho }_0{\eta }_{\rho }\end{array}\end{array}$$

In paper III, we introduced a mixing term, ${\eta }_{\rho }$, 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,

(119)$$\bar{q}=\frac{1}{4\pi }{\int }_{{\mathrm{\Omega }}_H}q(r,\theta ,\varphi )d\mathrm{\Omega }$$

where ${\int }_{{\mathrm{\Omega }}_H}d\mathrm{\Omega }=4\pi$ represents the integral over the spherical $\theta$ and $\varphi$ angles at constant radius.

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

(120)$$\frac{\partial {\rho }^{\prime }}{\partial t}=-\mathbf{U}\cdot \mathrm{\nabla }{\rho }^{\prime }-{\rho }^{\prime }\mathrm{\nabla }\cdot \mathbf{U}-\mathrm{\nabla }\cdot \left({\rho }_0\tilde{\mathbf{U}}\right),$$

written in a slightly different form:

(121)$$\frac{\partial {\rho }^{\prime }}{\partial t}+\mathrm{\nabla }\cdot \left({\rho }^{\prime }\mathbf{U}\right)=-\mathrm{\nabla }\cdot \left({\rho }_0\tilde{\mathbf{U}}\right).$$

We integrate this over a spherical shell of thickness $2h$ at radius $r_0$, i.e., ${\mathrm{\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:

(122)$$\begin{array}{r}\begin{array}{rl}\frac{1}{4\pi r_0^{2}2h}{\int }_{r_0-h}^{r_0+h}{r}^{2}dr{\int }_{{\mathrm{\Omega }}_H}[\frac{\partial {\rho }^{\prime }}{\partial t}+\mathrm{\nabla }\cdot \left({\rho }^{\prime }\mathbf{U}\right)]d\mathrm{\Omega }=& -\frac{1}{4\pi r_0^{2}2h}{\int }_{r_0-h}^{r_0+h}{r}^{2}dr{\int }_{{\mathrm{\Omega }}_H}\mathrm{\nabla }\cdot \left({\rho }_0\tilde{\mathbf{U}}\right)d\mathrm{\Omega }\\ =& {-\frac{1}{4\pi r_0^{2}2h}{\int }_{{\mathrm{\Omega }}_H}\left[{\rho }_0(\tilde{\mathbf{U}}\cdot {\mathbf{e}}_r)\right]r^{2}d\mathrm{\Omega }|}_{r_0-h}^{r_0+h}\\ =& 0,\end{array}\end{array}$$

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 ${\mathrm{\Omega }}_H$. We see that the right hand side disappears since ${\int }_{{\mathrm{\Omega }}_H}{\rho }_0(\tilde{\mathbf{U}}\cdot {\mathbf{e}}_r)d\mathrm{\Omega }=0$, which follows from the definition of $\tilde{\mathbf{U}}$. Now, expanding the remaining terms and taking the limit as $h\to 0,$ we can write

(123)$$\begin{array}{r}\begin{array}{rl}0=& \underset{h\to 0}{lim}\frac{1}{4\pi r_0^{2}2h}{\int }_{r_0-h}^{r_0+h}{r}^{2}dr{\int }_{{\mathrm{\Omega }}_H}[\frac{\partial {\rho }^{\prime }}{\partial t}+\mathrm{\nabla }\cdot ({\rho }^{\prime }\mathbf{U})]d\mathrm{\Omega }\\ =& \frac{\partial }{\partial t}\left(\underset{h\to 0}{lim}\frac{1}{4\pi r_0^2}\frac{1}{2h}{\int }_{r_0-h}^{r_0+h}{r}^{2}dr{\int }_{{\mathrm{\Omega }}_H}{\rho }^{\prime }d\mathrm{\Omega }\right)+\underset{h\to 0}{lim}[\frac{1}{4\pi r_0^{2}2h}{\int }_{r_0-h}^{r_0+h}{r}^{2}dr{\int }_{{\mathrm{\Omega }}_H}\mathrm{\nabla }\cdot ({\rho }^{\prime }\mathbf{U})d\mathrm{\Omega }]\\ =& \frac{\partial }{\partial t}\left(\frac{1}{4\pi }{\int }_{{\mathrm{\Omega }}_H}{\rho }^{\prime }d\mathrm{\Omega }\right)+\underset{h\to 0}{lim}\left\{{\frac{1}{r_0^{2}2h}[\frac{1}{4\pi }{\int }_{{\mathrm{\Omega }}_H}{\rho }^{\prime }(\mathbf{U}\cdot {\mathbf{e}}_r)d\mathrm{\Omega }]r^2|}_{r_0-h}^{r_0+h}\right\}\\ =& \frac{\partial }{\partial t}\overline{{\rho }^{\prime }}+\underset{h\to 0}{lim}\left\{{\frac{1}{r_0^{2}2h}\left[\overline{{\rho }^{\prime }(\mathbf{U}\cdot {\mathbf{e}}_r)}\right]r^2|}_{r_0-h}^{r_0+h}\right\}\\ =& \frac{\partial }{\partial t}\overline{{\rho }^{\prime }}+\underset{h\to 0}{lim}\frac{1}{r_0^{2}2h}{\int }_{r_0-h}^{r_0+h}\mathrm{\nabla }\cdot \left[\overline{{\rho }^{\prime }(\mathbf{U}\cdot {\mathbf{e}}_r)}{\mathbf{e}}_r\right]r^{2}dr\\ =& \frac{\partial }{\partial t}\overline{{\rho }^{\prime }}+\mathrm{\nabla }\cdot \left[\overline{{\rho }^{\prime }(\mathbf{U}\cdot {\mathbf{e}}_r)}{\mathbf{e}}_r\right]\end{array}\end{array}$$

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

In short,

(124)$$\frac{\partial }{\partial t}\overline{{\rho }^{\prime }}=-\mathrm{\nabla }\cdot \left[\overline{{\rho }^{\prime }(\mathbf{U}\cdot {\mathbf{e}}_r)}{\mathbf{e}}_r\right]=-\mathrm{\nabla }\cdot \left({\eta }_{\rho }{\mathbf{e}}_r\right),$$

and thus, ${\eta }_{\rho }=\overline{({\rho }^{\prime }\mathbf{U}\cdot {\mathbf{e}}_r)}$.

We need both ${\eta }_{\rho }$ alone and its divergence for the various terms in the construction of $w_0$ and the correction to ${\rho }_0$. The quantity ${\eta }_{\rho }$ 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 ${\eta }_{\rho }$ by constructing the quantity ${\rho }^{\prime }\tilde{\mathbf{U}}\cdot {\mathbf{e}}_r$ 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 ${\eta }_{\rho }$, ${\eta }_{\rho ,r+1/2}$ are then constructed by simple averaging:

(125)$${\eta }_{\rho ,r+1/2}=\frac{{\eta }_{\rho ,r}+{\eta }_{\rho ,r+1}}{2}.$$

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

(126)$${[\mathrm{\nabla }\cdot ({\eta }_{\rho }{\mathbf{e}}_r)]}^{n+1/2}=-\frac{\overline{{{\rho }^{\prime }}^{n+1}}-\overline{{{\rho }^{\prime }}^n}}{\mathrm{\Delta }t}=-\frac{\overline{{{\rho }^{\prime }}^{n+1}}}{\mathrm{\Delta }t},$$

where we have made use of the fact that $\overline{{{\rho }^{\prime }}^n}=0$ by construction.

$\eta$ Flow Chart

1. Enter advance_timestep with $[{{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}},{{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}{]}^{n-\frac{1}{2}}$.

2. Call make_w0. The spherical version uses uses ${{\eta }_{\rho }}^{\mathrm{e}\mathrm{c},n-\frac{1}{2}}$ and ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c},n-\frac{1}{2}}$.

3. Call density_advance. The plane-parallel version computes ${{\eta }_{\rho }}^{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x},n+\frac{1}{2},\ast }$.

4. Call make_etarho to compute $[{{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}},{{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}{]}^{n+\frac{1}{2},\ast }$. The plane-parallel version uses ${{\eta }_{\rho }}^{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x},n+\frac{1}{2},\ast }$.

5. Call make_psi. The plane-parallel version uses ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c},n+\frac{1}{2},\ast }$.

6. Call make_w0. The spherical version uses uses ${{\eta }_{\rho }}^{\mathrm{e}\mathrm{c},n+\frac{1}{2},\ast }$ and ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c},n+\frac{1}{2},\ast }$.

7. Call density_advance. The plane-parallel version computes ${{\eta }_{\rho }}^{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x},n+\frac{1}{2}}$.

8. Call make_etarho to compute $[{{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}},{{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}{]}^{n+\frac{1}{2}}$. The plane-parallel version uses ${{\eta }_{\rho }}^{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x},n+\frac{1}{2}}$.

9. Call make_psi. The plane-parallel version uses ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c},n+\frac{1}{2}}$.

Computing ${{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}}$ and ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}$

This is done in make_eta.f90.

Plane-Parallel

We first compute a radial edge-centered multifab, ${{\eta }_{\rho }}^{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x}}$, using

(127)$${\eta }_{\rho ,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x}}=[({\tilde{\mathbf{U}}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{n+\frac{1}{2}}\cdot {\mathbf{e}}_r)+w_{0,r+\frac{1}{2}}^{n+\frac{1}{2}}]{\rho }_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{n+\frac{1}{2}}-w_{0,r+\frac{1}{2}}^{n+\frac{1}{2}}{\rho }_{0,r+\frac{1}{2}}^{n+\frac{1}{2},\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}$$

${{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}}$ is the edge-centered “average” value of ${\eta }_{\rho }^{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x}}$,

(128)$${\eta }_{\rho ,r+\frac{1}{2}}^{\mathrm{e}\mathrm{c}}=\stackrel{―}{{\eta }_{\rho ,\mathbf{i}+\frac{1}{2}{\mathbf{e}}_r}^{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x}}}$$

${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}$ is a cell-centered average of ${{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}}$,

(129)$${\eta }_{\rho ,r}^{\mathrm{c}\mathrm{c}}=\frac{{\eta }_{\rho ,r+\frac{1}{2}}^{\mathrm{e}\mathrm{c}}+{\eta }_{\rho ,r-\frac{1}{2}}^{\mathrm{e}\mathrm{c}}}{2}.$$

Spherical

First, construct ${\eta }_{\rho }^{\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{t}}=[{\rho }^{\prime }(\tilde{\mathbf{U}}\cdot {\mathbf{e}}_r){]}^{n+\frac{1}{2}}$ using:

(130)$$[\frac{{\rho }^n+{\rho }^{n+1}}{2}-{\left(\frac{{\rho }_0^n+{\rho }_0^{n+1}}{2}\right)}^{\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{t}}]\sum _d\left(\frac{{\tilde{\mathbf{U}}}_{\mathbf{i}+\frac{1}{2}{\mathbf{e}}_d}^{n+\frac{1}{2}}\cdot {\mathbf{e}}_d+{\tilde{\mathbf{U}}}_{\mathbf{i}-\frac{1}{2}{\mathbf{e}}_d}^{n+\frac{1}{2}}\cdot {\mathbf{e}}_d}{2}\right)n_d.$$

Then, ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}$ is the cell-centered average of ${\eta }_{\rho }^{\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{t}}$,

(131)$${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}=\overline{{\eta }_{\rho }^{\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{t}}}.$$

On interior faces, ${{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}}$ is the average of ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}$,

(132)$${\eta }_{\rho ,r-\frac{1}{2}}^{\mathrm{e}\mathrm{c}}=\frac{{\eta }_{\rho ,r-1}^{\mathrm{c}\mathrm{c}}+{\eta }_{\rho ,r}^{\mathrm{c}\mathrm{c}}}{2}.$$

At the upper and lower boundaries, we use

(133)$$\begin{array}{r}\begin{array}{rlr}{\eta }_{\rho ,-\frac{1}{2}}^{\mathrm{e}\mathrm{c}}& =& 0,\\ {\eta }_{\rho ,\mathrm{n}\mathrm{r}-\frac{1}{2}}^{\mathrm{e}\mathrm{c}}& =& {\eta }_{\rho ,\mathrm{n}\mathrm{r}-1}^{\mathrm{c}\mathrm{c}}.\end{array}\end{array}$$

Using ${{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}}$

Plane-Parallel

NOT USED.

Spherical

In make_w0, ${{\eta }_{\rho }}^{\mathrm{e}\mathrm{c}}$ is used in the construction of the RHS for the $\delta w_0$ equation.

Using ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}$

Plane-Parallel

In make_psi, $\psi ={{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}g$.

Spherical

In make_w0, ${{\eta }_{\rho }}^{\mathrm{c}\mathrm{c}}$ is used in the construction of the RHS for the $\delta w_0$ equation.