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,

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

so in addition to specifying the sound speed and ambient density, 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

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.