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.

Results

Analytic function -

Since I am having difficulty saving an initial density function in astrobear for use with ProblemBeforeStep, I just plotted the 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 are what I would expect AstroBEAR should give us after 1) strictly initializing a sinusoidal density perturbation, and 2) allowing it to evolve in time with self-gravity turned on. These results are as follows:

t=0

t=1/omega

These are using periodic physical and elliptical BCs.

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

At t=0.0016

And at t=0.005:

It looks like exponential growth…

Makes me wonder:

1) Scaling error?

2) Error in analytic function at t=6 Myr?

Here's the last curiosity, when I look at the frame where the density = e*rhoInit, the peaks match the 'analytical function' I plotted in Mathematica, but the valleys don't. .