Computing the Jeans Length

The first step in setting up this problem was to consider the Jeans Length for a given ambient density and temperature. Recall,

\lambda_J = C_s(\frac{\pi} { G \rho_0})^{1/2}

I used the ambient density and temperature of the BP case for the Bonnor Ebert runs. This gave

\rho0 = 3.34 \times 10^{-23} g/cm^3

\ T = 100K

Using the isothermal sound speed of,

C_s = (\frac{K_B*T}{m_H})^{1/2}

where Kb is the boltzmann constant and mH is the mass of hydrogen, I calculated

\lambda_J = 1.08 \times 10^{20} cm \approx 35 pc

Determining the form of the perturbation function

Recall, in our linear perturbative analysis on the gravito-hydrodynamic equations, we had:

\rho = \rho_0 + \rho_1

v = v_1

P = P_0 + P_1

\phi = \phi_0 + \phi_1

This led to a 2nd order ODE in the density perturbation which had the following solution:

\rho_1 = C \cos(kx) e^{[4 \pi G \rho_0 (1-\frac{\lambda_J^2}{\lambda^2})]^{1/2}t}

In addition to specifying the sound speed and ambient density for this perturbation (implicit to the Jeans length), we need to specify the wave length of the perturbation. With the desire to keep the ratio \lambda_J / \lambda small, I chose

\lambda = 10 \lambda_J = 350 ~ pc

Using this, we can set the initial (t=0) density inside of the mesh as:

\rho = \rho_0 + \rho_1 = \rho_0 + cos(kx)

To get the form of the corresponding velocity perturbation, we use the linearized continuity equation:

\frac{\partial \rho_1}{\partial t} + \rho_0 \nabla \cdot \vec{v_1} = 0

Rearranging,

\frac{d}{dx} v_1 = -\frac{1}{\rho_0} \frac{\partial \rho_1}{\partial t} = - \frac{1}{\rho_0} e^{ikx}\Gamma e^{\Gamma t}

where \Gamma is the growth rate given by

\Gamma = [4 \pi G \rho_0 (1- (\frac{\lambda_J}{\lambda})^2)]^{1/2}

(cf. ​growth rate)

\frac{d}{dx} v_1 = - \frac{\Gamma e^{\Gamma t}}{\rho_0} cos(kx)

and so,

v_1 = (\frac{-\Gamma e^{\Gamma t}}{\rho_0})(\frac{1}{k})sin(kx)

Thus, we have the form of the initial (t=0) velocity perturbation:

v = v_1 = \frac{- \Gamma}{\rho_0 k} sin(kx)

Finally, to get the form of the pressure perturbation for an isothermal simulation, we use the Ideal Gas Equation:

P = \frac{\rho K_B T}{m_H}

which gives

P = P_0 + P_1 = \frac{K_B T}{m_H}(\rho_0 + \rho_1)

Thus, we now have the form of all perturbations of the fluid variables.

Setting up the problem domain

I decided on a 1D grid (a string of cells in x) that was much greater than the Jeans length:

L = 1750 ~ pc

The Jeans Length needed to be resolved to prevent artificial fragmentation, so

\triangle x < \frac{1}{4} \lambda_J \Rightarrow

\triangle x < 9 ~pc

Choosing 1750 cells in the x direction (small computational cost given 1D sim) satisfies this by having

\triangle x = 1~pc

so there are about 35 cells/ Jeans length.

Given the simulation is 1D, I set the ylower/yupper and zlower/zupper bounds to 0 in global.data, and set the number of cells in y = 1 and z = 1.

Problem Scales

lscale	\lambda = 350 ~pc = 1.08*10^{21}~cm
rscale	\rho_0 = 3.34*10^{23}~g/cm^3
tempscale	1

A meaningful simulation time seemed to be the e-folding time for the perturbation function. That is,

t_{final} = \tau \approx (4 \pi G \rho_0)^{-1/2}

In order to calculate this in computational units, I first ran astrobear to get the output file scales.data, which lists timescale.

After computing

t_{final} = \tau \approx (4 \pi G \rho_0)^{-1/2} \approx 6 ~Myr

It was easy to solve for tfinal,

t_{final} = 6~Myr / tscale

The code

The following is for a 1D isothermal sim of the Jeans Instability.

In ProblemModuleInit, I read in ambient density and temperature. These are used to uniformly set a corresponding Ambient Object.

READ ~ rho

Info \% q(i, :, :, 1) = rho

Ambient\%density = info \% q(i, :, :, 1)

In ProblemGridInit, I loop over the cells in x (indexed by i) adding the following density perturbation:

info \% q(i, :, :, 1) = info\%q(i, :, :, 1) + 0.1*info\%q(i, :, :, 1)*Cos(kx)*e^0

I am not 100% sure, but I think if I wanted to place the analytical function in astrobear, I would then add to ProblemBeforeStep a loop over cells in i that continues to perturb the density function in time t:

info \% q(i, :, :, 1) = info\%q(i, :, :, 1, t=0)*e^{\omega t}

Thus we should see the characteristic growth rate of \omega.

For further info, you can download the Mathematica calculations ​here

and the code ​here

Results

Analytic function -

Here is the plotted analytical function for the initial density perturbation and the perturbation after an e-folding time in Mathematica. The e-folding time in physical units is ~6 Myr, and in computational units corresponds to t = 0.0016.

t=0

t=1/omega

These next plots show AstroBEAR's output using density, velocity, and pressure perturbations described earlier:

t=0

Vx is in red, rho is in green. The density peaks at the zero's of Cos(kx). E.g., a local max is at 1. The vx perturbation is a Sin(x) function, and it has minima that corresponds with the density maxima. This produces converging flows at the max in density, as you would expect you need to develop the Jeans instability.

t=5/omega

By the end of the sim, the density perturbation has grown a little, but velocity has grown a lot.

(These are using periodic physical and elliptical BCs, and an ideal EOS with gamma = 1.0001).

Now, a query-over-time of rho at r = (1,0,0) produces the following output when compared to the analytical result:

Pretty tight fit!