Self Gravity

Physics

Equations of Self Gravity

The inclusion of self gravity involves an additional force to the momentum equation

\frac{d \rho \mathbf{v}}{dt} = \mathbf{f}_g = -\rho \nabla \phi

and the energy equation now includes the work done by gravity

\frac{d E}{dt} = \mathbf{v} \cdot \mathbf{f}_g = -\rho \mathbf{v} \cdot \nabla \phi

and requires the solution for the gas potential

\nabla^2 \phi = 4 \pi G \left (\rho - \rho_0\right)

where

• \rho_0 = \bar{\rho} for periodic boundary conditions
• \rho_0 = 0 otherwise.

Making the equations conservative

AstroBEAR uses hypre's linear solver to solve Poisson's equation for the potential. However, there are two methods for including the source terms for momentum and energy. The non-conservative approach simply includes the terms during a source update, while the conservative approach recasts the RHS of the momentum and energy equations as total divergences - so they can be differenced conservatively.

Substituting the Poisson equation into the momentum equation we have after some manipulations

\frac{d \rho \mathbf{v}}{dt} = (\frac{\nabla^2\phi}{4 \pi G}+\rho_0) \nabla \phi = \nabla \cdot \left [ \frac{1}{4 \pi G} \left ( \nabla \phi \nabla \phi - \frac{1}{2} \mathbf{I} \left ( \nabla \phi \cdot \nabla \phi \right) \right ) + \rho_0 \phi \mathbf{I} \right ]

where we have the momentum flux tensor

\mathbf{T} = \frac{1}{4 \pi G} \left ( \nabla \phi \nabla \phi - \frac{1}{2}\mathbf{I} \left ( \nabla \phi \cdot \nabla \phi \right) \right) + \rho_0 \phi \mathbf{I}

For the energy equation, we can only find a conservative approach - if we consider the evolution of the total combined energy (hydro + gravitational).

\frac{d \left(E + \frac{1}{2} \left (\rho - \rho_0 \right) \phi \right)}{dt} = -\nabla \cdot \mathbf{F_g}

Combining this with the our previous energy equation we arrive at

\nabla \cdot \mathbf{F_g} = -\frac{1}{2} \frac{d \left (\rho \phi \right)}{dt} + \frac{1}{2} \rho_0 \frac{d \phi}{dt}+ \rho \mathbf{v} \cdot \nabla \phi

Expanding the time derivative and using the continuity equation, we arrive at

\nabla \cdot \mathbf{F_g} = \frac{1}{2} \phi \nabla \cdot \left ( \rho \mathbf{v} \right) - \frac{1}{2} \left (\rho - \rho_0 \right) \dot{\phi} + \rho \mathbf{v} \cdot \nabla \phi

Now from time differencing the poisson equation, we arrive at

\nabla^2 \dot{\phi} = 4 \pi G \dot{\rho} = - 4 \pi G \nabla \cdot (\rho \mathbf{v})

and multiplying both sides by \frac{\phi}{8 \pi G} gives

\frac{1}{2} \nabla \cdot (\rho \mathbf{v}) \phi = -\frac{1}{8 \pi G} \phi \nabla^2 \dot{\phi}

Adding the LHS and subtracting the RHS and substituting

\rho - \rho_0 = \frac{1}{4 \pi G} \nabla^2 \phi

into the divergence equation gives us

\nabla \cdot \mathbf{F_g} = \nabla \cdot \left (\phi \rho \mathbf{v} \right) + \frac{1}{8 \pi G} \left (\phi \nabla^2 \dot{\phi} - \dot{\phi}\nabla^2\phi \right)

Now we can use the relationship

\nabla \cdot \left [ \phi \nabla \dot{\phi} - \dot{\phi} \nabla \phi \right ]= \phi \nabla^2 \dot{\phi} - \dot{\phi} \nabla^2 \phi

to arrive at

\mathbf{F_g} = \phi \rho \mathbf{v} + \frac{1}{8 \pi G} \left (\phi \nabla \dot{\phi} - \dot{\phi}\nabla \phi \right)

Solving for \dot{\phi}

Now to solve for \dot{\phi} we can solve directly

\nabla^2 \dot{\phi} = 4 \pi G \dot{\rho} = 4 \pi G \left ( \frac{\rho^{n+1} - \rho^n}{\Delta t} \right )

or we can simply wait until we have solved for the new potential based on the updated density

\nabla^2 \phi^{n+1} = 4 \pi G \rho^{n+1}

and then use

\dot{\phi} = \frac{\phi^{n+1} - \phi^n}{\Delta t}

to apply the energy correction.

Currently what is done in AstroBEAR

AstroBEAR's momentum flux is first calculated using

\mathbf{T^n} = \frac{1}{4 \pi G} \left ( \nabla \phi^n \nabla \phi^n - \frac{1}{2} \mathbf{I} \left ( \nabla \phi^n \cdot \nabla \phi^n \right ) + \rho_0 \phi^n \mathbf{I} \right )

Where \phi^{n} is used and everything is calculated at cell faces.

Then the energy is updated using

E^{n+1} = E^n - \rho \mathbf{v} \cdot \nabla \phi

where each component in the sum \rho \mathbf{v} \cdot \nabla \phi is averaged at the left and right face for each dimension.

Then following the second elliptic solve, the momentum is updated again by calculating the difference between the original momentum flux tensor \mathbf{T} and the time centered one \mathbf{T}^{n+1/2} = \mathbf{T}^n + \frac{1}{2} \left ( \mathbf{T}^{n+1} - \mathbf{T}^n \right)

And a corrective momentum flux tensor is calculated

\mathbf{T'} = \frac{1}{2} \left (T^{n+1} - T^n \right )

and applied.

Likewise, a energy source term correction is applied

E^{n+1} = E^{n+1} - \frac{1}{2} \rho \mathbf{v} \cdot \nabla \left (\phi^{n+1} - \phi^n \right)

Notes

The electric field generated by a point chargeQis

E=\frac{1}{4\pi\epsilon_{0}}\frac{Q}{r^{2}}

On the other hand, the gravitational acceleration by a point mass of M is

g=-G\frac{M}{r^{2}}

The Poisson equation for electrostatic potential \phi_{e} is

\nabla^{2}\phi=-\frac{\rho}{\epsilon_{0}}

Analytically we can get the Poisson equation for the gravity

\nabla^{2}\phi=4\pi G\rho

where \rho is the mass density. This equation describes how the potential \phi is determined by the mass density distribution.

For example, consider the uniform density distribution

begin{equation}\label{eq:uniDensity}\tag{1} \rho= \begin{cases} \rho_0 & r < R \\ 0 & r \ge R \end{cases} \end{equation}

Use spherical coordinates, the Poisson equation for the density distribution Equation (1) can be written as

begin{equation}\label{eq:poissonSphCoord}\tag{2} \frac{1}{r^{2}}\frac{\partial }{\partial r}\left(r^{2}\frac{\partial\phi}{\partial r} \right)=4\pi G\rho \end{equation}

Or

begin{equation}\label{eq:poissonSphCoord1}\tag{3} \frac{\partial ^{2}\phi }{\partial r^{2}}+\frac{2}{r}\frac{\partial\phi}{\partial r}=4\pi G\rho \end{equation}

Solving Equation (3) with the density Equation (2), we obtain the solution for the potential for the uniform density sphere:

begin{equation}\label{eq:solUniSphere}\tag{4} \phi(r)=\begin{cases} -\frac{4\pi}{3} G\rho \frac{R^{3}}{r} & r>R\\ -\frac{2\pi G\rho}{3}\left[R^{2}-r^{2}\right]-\frac{4\pi}{3}G\rho R^{2} &r\ge R \end{cases} \end{equation}

Free-fall Time

The free-fall time is the characteristic time that would take a body to collapse under its own gravitational attraction, if no other forces existed to oppose the collapse. Using Gauss's theorem, the gravitational acceleration g_{r} is given by

begin{equation}\label{eq:graAcce}\tag{5} g_{r}=-\frac{GM_{r}}{r^{2}} \end{equation}

where M_{r} is the mass inside r

begin{equation}\label{eq:insideMass}\tag{6} M_{r}=\int_{0}^{r}\rho4\pi r^{2}dr \end{equation}

So the equation of motion for a gas molecule under a control of the self-gravity can be written as

begin{equation}\label{eq:motion}\tag{7} \frac{d^{2}r}{dt^{2}}=-\frac{GM_{r}}{r^{2}} \end{equation}

Equation (7) can be written as the form of the conservation of mechanical energy

begin{equation}\label{eq:conserE}\tag{8} \frac{1}{2}\left(\frac{dr}{dt}\right)^{2}-\frac{GM_{r}}{r}=E \end{equation}

Suppose the gas molecule is initially located at r=R. So the total energy

begin{equation}\label{eq:totalE}\tag{9} E=-\frac{GM_{r}(R)}{R} \end{equation}

Now consider gas collapse. So M_{r}= M_{r}(R). Equation (8) gives

begin{equation}\label{eq:approx1}\tag{10} \frac{dr}{dt}=-[2GM_{r}(R)]^{\frac{1}{2}}\left(\frac{1}{r}-\frac{1}{R} \right)^{\frac{1}{2}} \end{equation}

Let

begin{equation}\label{eq:Etadef}\tag{11} r=R\cos^{2}\eta \end{equation}

From the initial state we have \eta(t=0)=0. So Equation (10) reduces to

begin{equation}\label{eq:reduced}\tag{12} \cos^{2}\eta\frac{d\eta}{dt}=\left(\frac{GM_{r}(R)}{2R^{3}}\right)^{\frac{1}{2}} \end{equation}

Solve Equation (12) we have

begin{equation}\label{eq:sol}\tag{13} t=\left(\frac{R^{3}}{2GM_{r}(R)} \right)^{\frac{1}{2}}\left(\eta+\frac{\sin2\eta}{2} \right) \end{equation}

Let r=0 or \eta=\frac{\pi}{2}. we have the free-fall time

begin{equation}\ref{eq:freeFall}\tag{14} \begin{aligned} t_{ff} & =\left(\frac{R^{3}}{2GM_{r}(R)}\right)^{\frac{1}{2}}\frac{\pi}{2}\\ & =\left(\frac{3\pi}{32G\rho_{0}}\right)^{\frac{1}{2}} \end{aligned} \end{equation}

where

begin{equation}\label{rhoBar}\tag{15} \rho_{0}=\frac{M_{r}(R)}{\frac{4\pi R^{3}}{3}} \end{equation}

is the average density of the gas.

Implementation

Uniform Density Cloud

In this section we consider the problem of the collapse of a pressureless (T=0), uniform, initially motionless cloud which has an analytic solution.

The definition of \eta in Equation (11) can also be written as

begin{equation}\tag{16} \cos\eta=\left(\frac{\rho}{\rho_{0}} \right)^{-\frac{1}{6}} \end{equation}

From Equations (13) and (14) we have

begin{equation}\label{eq:trho}\tag{17} \frac{t_{\rho}}{t_{ff}}=\frac{2}{\pi}\left[\eta + \frac{1}{2}\sin(2\eta) \right] \end{equation}

Equation (17) describes the time it takes for the cloud to collapse to a density \rho.

Boundary Conditions

Present Considerations

Some notes on periodic BC's and particles

With periodic boundary conditions, we would like the total potential to be unchanged under the conversion of gas to particle. This requires adjusting the particle potential so that it has the same constant of integration as well as being periodic. When hypre solves for the potential, the constant of integration is a free parameter and I believe hypre adjusts it so that <\phi>=0. To calculate the potential of a point charge, we need the Greens function corresponding to

\nabla_D^2 G(x,y) = \delta^D(x-y)

where D is the dimenesion of the problem. For 3D,

G(x,y)=\frac{-1}{4\pi|x-y|}

and for 2D,

G(x,y)=\frac{\ln(|x-y|)}{2\pi}

and for 1D,

G(x,y)=\frac{|x-y|}{2}

Now with periodic BC's, the solution to the gas potential is given by

\nabla^2 \phi = 4\pi G (\rho-\bar{\rho}), so the solution to the potential would actually be given by \phi(y)=-\int_V{\frac{M\delta(x-y)-M/V}{|x-y|}dx} where M is the mass of the particle. Note this is not the same as using mirror versions of the particle… An easier way to solve this however, is to use fourier transforms.

k^2\hat{\phi}(k)=4\pi G\hat{\rho}(k) with \hat{\phi}(0)=0. Now \hat{\rho}(k) is just the fourier transform of a delta function times the mass, which is M, so \hat\phi(k)=\frac{4\pi GM}{k^2} and \phi(x)=4 \pi G M\int{\frac{1}{k^2}e^{ikx} dk}

This can be discretized as

\phi(l)=\frac{4 \pi GM}{\Delta k^2} \displaystyle \sum_{j=-N/2}^{N/2}{\frac{1}{j^2}e^{2\pi ijl/N} \frac{\Delta k}{2\pi}}=\frac{4 \pi GM}{(2 \pi/L)^2}\left(\frac{2 \pi}{L}\right) \displaystyle \sum_{j=1}^{N/2}{\frac{1}{j^2}\cos(2\pi jl/N)}

This still requires summation over many wave numbers - although the higher wave numbers have less impact because of the k^{-2} dependence. However in 2D this is lessened because there are more wave vectors in each anuli, and in 3D there is equal power in each shell. This then, becomes an order N^{6 operation}

Solutions:

• Don't ever use sink potential?
• Don't use sink potential to update gas potential
• Solve for gas potential after accretion?