Momentum Conserving Self Gravity Source Terms in AstroBEAR

• q - initial state
• qLx - left edge of x - interface in predictor
• qLy - lower edge of y - interface in predictor
• qRy - upper edge of y - interface in predictor
• q2 - intermediate cell centered quantities
• q2Lx - updated left edge of x - interface used for final riemann solve
• fx - predictor x flux
• f2x - final x flux
• q3 - new values

Normal CTU:

• qLx = Reconstruct(q), etc… (for all 6 qxx states)
• fx=Riemann_Solve(qLx,qRx), etc… (for all 3 edges)
• q2 = q+fx+fy+fz
• q2Lx=qLx+fy+fz etc… (for all 6 q2xx states)
• f2x=Riemann_Solve(q2Lx,q2Rx) (for all 3 edges)
• q3=q+f2x+f2y+f2z

Strang split self-gravity:

• q = q + SG(q, phi)*hdt
• qLx = Reconstruct(q), etc… (for all 6 qxx states)
• fx=Riemann_Solve(qLx,qRx), etc… (for all 3 edges)
• q2 = q+fx+fy+fz
• q2Lx=qLx+fy+fz etc… (for all 6 q2xx states)
• f2x=Riemann_Solve(q2Lx,q2Rx) (for all 3 edges)
• q3=q+f2x+f2y+f2z
• q3=q3+SG(q3, phi)*hdt

Momentum conserving self-gravity

• qLx = Reconstruct(q), etc… (for all 6 qxx states)
• qLx(px) =qLx(px)+SG_x(q,phi)*hdt
• qLy(py) = qLy(py)+SG_y(q,phi)*hdt
• qLy(pz) = qLy(pz)+SG_z(q,phi)*hdt
• fx=Riemann_Solve(qLx,qRx), etc… (for all 3 edges)
• q2 = q+fx+fy+fz + SG(q,phi)
• q2Lx=qLx+fy+fz etc… (for all 6 q2xx states)
• q2Lx=q2Lx+SG_y(q,phi)+SG_z(q,phi) (for all 6 q2xx states)
• f2x=Riemann_Solve(q2Lx,q2Rx) (for all 3 edges)
• f2x(p)=f2x(p)+SG_PFlux_x(phi) (same for y and z)
• q3=q+f2x+f2y+f2z+SG_EFlux_(f2x(rho),f2y(rho), f2z(rho), phi)

SG_x has two non-zero terms…

[ \left [ \begin{array}{c} 0 \\ \rho \nabla_x \phi \\ 0 \\ 0 \\ \rho v_x \nabla_x \phi \end{array} \right ] \]

while SG_PFlux_x(phi) has only 1 non-zero term

[ \left [ \begin{array}{c} 0 \\ \frac{1}{8G\pi} \left ( (\nabla_x \phi)^2 - (\nabla_y\phi)^2 - (\nabla_z\phi)^2 \right ) + \bar{\rho}\phi \ \\ \frac{1}{4G\pi} \nabla_x \phi \nabla_y\phi \\ \frac{1}{4G\pi} \nabla_x \phi \nabla_z\phi \\ 0 \end{array} \right ] \]

and SG_EFlux(f2x, f2y, f2z, phi) is

[ \left [ \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ \rho v_x \nabla_x \phi + \rho v_y \nabla_y \phi + \rho v_z \nabla_z \phi \end{array} \right ] \]

Making the momentum be conserved requires recasting the gravitational force as a total differential

\dot p=-\rho \nabla \phi = -\nabla \cdot F so we need to find F such that

\nabla \cdot F = \rho \nabla \phi

Since \nabla^2 \phi = 4\pi G (\rho-\bar{\rho})

we can substitute for \rho and we have

\nabla \cdot F = (\frac{\nabla^2\phi}{4 \pi G}+\bar{\rho}) \nabla \phi which is equivalent (in 1D) to \nabla \cdot \left [\frac{1}{2} \frac{\left (\nabla \phi\right )^2}{4 \pi G} + \bar{\rho} \phi \right ]

where we can identify the equivalent momentum flux tensor as F=\frac{1}{2} \frac{\left (\nabla \phi\right )^2}{4 \pi G} + \bar{\rho} \phi

• In more than 1D we have

begin{align*} \partial_t p_i&=\rho \partial_i \phi \\ & = \left (\frac{\partial_j(\partial_j \phi)}{4 \pi G}+\bar{\rho} \right ) \partial_i \phi \\ & = \frac{\partial_j \left ( \partial_j \phi \partial_i \phi \right) - (\partial_j \phi) \partial_j (\partial_i \phi )}{4 \pi G}+\partial_i (\bar{\rho}\phi) \\ & = \frac{\partial_j \left ( \partial_j \phi \partial_i \phi \right) - (\partial_j \phi) \partial_i (\partial_j \phi )}{4 \pi G}+\partial_i (\bar{\rho}\phi) \\ & = \frac{\partial_j \left ( \partial_j \phi \partial_i \phi \right)}{4 \pi G} - \partial_i\frac{(\partial_j \phi \partial_j \phi)}{8 \pi G}+\partial_i (\bar{\rho}\phi) \\ & = \frac{\partial_j \left ( \partial_j \phi \partial_i \phi \right)}{4 \pi G} - \partial_i\left ( \frac{ \partial_j \phi \partial_j \phi}{8 \pi G} +\bar{\rho}\phi \right) \\ & = \partial_j \frac{\left ( \partial_j \phi \partial_i \phi \right)}{4 \pi G} - \partial_j \left ( \delta^i_j \left ( \frac{\partial_k \phi \partial_k \phi}{8 \pi G} + \bar{\rho}\phi \right) \right ) \\ & = \partial_j T_{ij} \mbox{ where } T_{ij} = \frac{\left ( \partial_j \phi \partial_i \phi \right)}{4 \pi G} - \delta^i_j\left(\frac{\partial_k \phi \partial_k \phi}{8 \pi G}+\bar{\rho}\phi \right) \\ \end{align*}