Questions about the opening time of self-gravity

[YaoguangPei]

Hi all,

I'm currently running a simulation that includes both turbulence and self-gravity. However, I'm a bit confused about the timing for turning on self-gravity. In several simulation papers I’ve read, researchers usually run the simulation without gravity for about 2 crossing times first, and then turn on self-gravity for another 2 crossing times. The reason for this approach is to avoid extremely high density values caused by initial density perturbations, which could trigger gravitational collapse too early—before turbulence has built up sufficient pressure support.

In my own simulations, I ran without self-gravity until t < 25 to allow turbulence to reach a steady state. Then, I turned on self-gravity and continued until t = 45. However, in both cases, the final maximum density only reached about 6000 × 1.4 m_H/cm³, which is much lower than what I expected. I was hoping the density would reach values on the order of 10⁶.

By contrast, in a previous simulation where I turned on self-gravity from the very beginning, the maximum density reached extremely high values (~10⁸), which even broke the Truelove condition. Since my goal is to study star formation, I do require high density regions, but I'm unsure how to set up the simulation properly.

My question is: Should I enable self-gravity only after running the simulation without it for about 2 crossing times? Also, what would be an appropriate dedt if I want to achieve high densities and velocities around 1 km/s?

Thank you in advance for your help!
Here are my input file, when the gravity turns off, the parameter of 4_pi_G will vanish.
'''

problem = Turbulence with power-law PS
reference =
configure = --prob=turb -fft

problem_id = Turb # problem ID: basename of output filenames file_type = hst # History data dump dt = 0.01 # time increment between outputs file_type = vtk # Binary data dump variable = prim,p,E # variables to be output dt = 1 # time increment between outputs cfl_number = 0.1 # The Courant, Friedrichs, & Lewy (CFL) Number nlim =-1 # cycle limit tlim = 45 # time limit integrator = vl2 # time integration algorithm xorder = 2 # order of spatial reconstruction ncycle_out = 1 # interval for stdout summary info dt_diagnostics=0 nx1 = 512 # Number of zones in X1-direction x1min = -2 # minimum value of X1 x1max = 2 # maximum value of X1 ix1_bc = periodic # inner-X1 boundary flag ox1_bc = periodic # outer-X1 boundary flag

nx2 = 512 # Number of zones in X2-direction
x2min = -2 # minimum value of X2
x2max = 2 # maximum value of X2
ix2_bc = periodic # inner-X2 boundary flag
ox2_bc = periodic # outer-X2 boundary flag

nx3 = 512 # Number of zones in X3-direction
x3min = -2 # minimum value of X3
x3max = 2 # maximum value of X3
ix3_bc = periodic # inner-X3 boundary flag
ox3_bc = periodic # outer-X3 boundary flag

refinement = none

nx1 =512 nx2 =32 nx3 =32 iso_sound_speed = 0.19 # equivalent to sqrt(gamma*p/d) for p=0.1, d=1 turb_flag = 3 # 1 for decaying, 2 (impulsive) or 3 (continuous) for driven turbulence dedt = 3e2 # Energy injection rate (for driven) or Total energy (for decaying) nlow = 1 # cut-off wavenumber at low-k nhigh =3 # cut-off wavenumber at high-k expo = 5/3 # power-law exponent tcorr = 0.1 # correlation time for OU process (both impulsive and continuous) dtdrive = 0.1 # time interval between perturbation (impulsive) f_shear = 0.5 # the ratio of the shear component rseed = -1 # if non-negative, seed will be set by hand (slow PS generation) unit_system=ism four_pi_G=1.871e-3 '''