Radiation-hydrodynamic Simulations of Spherical Protostellar Collapse for Very Low-mass Objects

Torsten Stamer; Shu-ichiro Inutsuka

doi:10.3847/1538-4357/aaee81

1. Introduction

Stars form in molecular clouds—regions of the interstellar medium which are very dense (>10²/cm³) and cold (10 K), allowing for the hydrogen gas of which they consist to be present in molecular form (H₂). Within these clouds, there exist regions of especially high density, usually called prestellar cores (or protostellar core if a protostar is present within). It is believed that star formation occurs where such cloud cores become massive enough to collapse under their own gravity in a process called protostellar collapse, the numerical simulation of which is the subject of this paper.

How these cores form in the first place is still a matter of active research. It has been firmly established by the Herschel survey (e.g., André et al. 2014) that the internal structure of molecular clouds is dominated by filaments, and that the majority of prestellar and protostellar cores in star-forming regions are found within them. Theoretical consideration of an isothermal, self-gravitating cylinder (Inutsuka & Miyama 1992, 1997) shows that filaments are gravitationally stable only if their line mass (i.e., mass per length) does not exceed a certain critical value of ${M}_{\mathrm{line},\mathrm{crit}}\approx 16\,{M}_{\odot }\,{\mathrm{pc}}^{-1}\times ({T}_{\mathrm{gas}}/10\,{\rm{K}})$ . Intersections between different filaments may create the regions most conducive to star formation, and especially star cluster formation, due to localized higher density (Myers 2009, 2011; Schneider et al. 2012). We may therefore state that star formation occurs in regions of overdense, radially collapsing filaments. Inutsuka & Miyama (1992, 1997) studied the stability of self-gravitating cylindrical structures and found that filaments are unstable to longitudinal perturbations, which grow and cause the filament to fragment into a number of prestellar cores.

Once a prestellar core has formed, the next step on the way to creating a star is for this core to collapse. Under isothermal conditions, the only requirement for such a collapse is that the core's mass should be larger than its thermal Jeans mass. In reality, however, after sufficient compression the system ceases to be isothermal and instead behaves closer to adiabatic. We must therefore consider the importance of thermodynamics: is the increase in pressure strong enough to halt the gravitational collapse?

The following description of protostellar collapse is considerably simplified since it assumes spherical symmetry and ignores magnetic fields and rotation. We will continue to make these assumptions throughout this work, since we focus on 1D radiation-hydrodynamical simulations. There exist much more sophisticated simulations of protostellar collapse, including for instance magnetic field effects (Tomisaka 2002; Machida et al. 2006; Commerçon et al. 2010; Federrath 2015), protoplanetary disk formation (Stamatellos & Whitworth 2008; Inutsuka et al. 2010; Machida et al. 2011; Nordlund et al. 2014; González et al. 2015; Tomida et al. 2015; Seifried et al. 2016), chemical evolution (Visser et al. 2015; Dzyurkevich et al. 2016; Hincelin et al. 2016), and non-ideal MHD (Tsukamoto et al. 2015a, 2015b; Masson et al. 2016; Wurster et al. 2016). In addition, there are synthetic observations which can be compared to actual data (Commerçon et al. 2012; Frimann et al. 2016; Seifried et al. 2016). However, in this study we aim to analyze the effects of initial conditions (mass, radius, etc.) on the collapse process in a simplified environment. Our work is rooted in Masunaga et al. (1998) and Masunaga & Inutsuka (2000) and it is in many ways similar to Vaytet & Haugbølle (2017), which we will refer to often in this paper. There are two main differences between our work and theirs: first, we focus on a different mass range, namely brown dwarfs and sub-brown dwarfs. This is a mass regime which has so far been under explored, and we wish to investigate if the results for higher masses (the lowest mass is 0.2 M_⊙ in Vaytet & Haugbølle 2017 ) also hold for brown dwarfs, or whether there are any new effects. The second difference lies in the radiative transfer scheme. While Vaytet & Haugbølle (2017) employ flux-limited radiative diffusion as well as the gray approximation, we use a frequency-dependent method which we developed in our previous paper (Stamer & Inutsuka 2018) specifically for spherically symmetric systems.

In the following, we will give a brief overview of the process of protostellar collapse. We model the cloud core as a homogeneous sphere with total mass M, radius R, and density ρ. We express the pressure using a polytropic equation of state, i.e., $P\propto {\rho }^{{\gamma }_{\mathrm{eff}}}$ , where ${\gamma }_{\mathrm{eff}}$ is the effective ratio of specific heats. The gravitational and pressure forces are then given by

$\begin{eqnarray}&&{F}_{{\rm{G}}}=\displaystyle \frac{{GM}}{{R}^{2}},\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&{F}_{{\rm{P}}}=\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial P}{\partial R}\propto \displaystyle \frac{{\rho }^{{\gamma }_{\mathrm{eff}}-1}}{R},\end{eqnarray} \tag{ 2 }$

for a homologous collapse, ρ ∝ R⁻³, so the relationship between the pressure and radius is ${F}_{{\rm{P}}}\propto {R}^{2-3{\gamma }_{\mathrm{eff}}}$ . The ratio between the two forces then scales as

$\begin{eqnarray}&&\displaystyle \frac{{F}_{{\rm{P}}}}{{F}_{{\rm{G}}}}\propto {R}^{4-3{\gamma }_{\mathrm{eff}}}.\end{eqnarray} \tag{ 3 }$

It follows that there is a critical value of ${\gamma }_{\mathrm{eff}}$ , namely 4/3. If ${\gamma }_{\mathrm{eff}}$ is larger than this value, pressure becomes stronger relative to gravity as the radius decreases, which inevitably leads to a stabilization and stop of the collapse. On the other hand, if it is smaller than this value, gravity becomes ever more dominant over pressure as the radius decreases, leading to an accelerating collapse.

We describe the evolution (summarized in Figure 1) as follows. Initially, the cloud core is essentially isothermal ( ${\gamma }_{\mathrm{eff}}=1$ ), since it is optically thin so that any heat gained through compression is efficiently radiated away and the gas is kept at the temperature of the surrounding molecular cloud, which is typically around 10 K. The temperature actually drops slightly since the core becomes opaque to the high-frequency heating radiation of background stars while the low-frequency cooling radiation is still able to escape, an effect which we will discuss further in Sections 3.1 and 3.4.

As the collapse progresses, compressional heating becomes ever stronger and the efficiency of the radiative cooling diminishes. Once the heating overtakes the cooling, the system's behavior transitions from isothermal to adiabatic. ${\gamma }_{\mathrm{eff}}$ then becomes 5/3 (for temperatures below a few hundred Kelvin) and 7/5 above that. The difference is due to the fact that for low temperatures, the rotational degree of freedom of the hydrogen molecule is frozen out, preventing it from contributing to the heat capacity. In any case, both of these values are larger than the critical value of 4/3, which means that the pressure force ultimately stops the gravitational contraction and a hydrostatic object with a typical size of a few astronomical units, the so-called first hydrostatic core, is formed. Its surface is marked by a shock in which the infalling material is decelerated. At this point, roughly one free fall time has passed, which for the densities typical of a solar mass cloud core corresponds to about 10⁵ yr. The first core then continues to increase in mass and temperature as it accretes material from the surrounding gas envelope. Once the temperature reaches ≈10³ K, hydrogen molecules begin to dissociate. This is a strong cooling effect, which acts to reduce the effective ratio of specific heats to about 1.1 and thereby triggers a second collapse. Because of this, the first core is generally thought to exist for no more than a few thousand years, although theoretical considerations (Tomida et al. 2010) and the surprisingly large number of candidate objects (e.g., Chen et al. 2010; Chen & Arce 2010; Dunham et al. 2011 and Pineda et al. 2011) suggest that it may live much longer in some cases. We will discuss this question in detail in Section 3.3. The second collapse continues until all of the hydrogen has been dissociated, so that ${\gamma }_{\mathrm{eff}}$ increases again to 5/3 and a second hydrostatic core forms, which is also known as a protostar. The collapse phase ends here and the system's pre-main-sequence evolution begins.

This paper aims to explore the process of protostellar collapse for very low-mass objects, i.e., the brown dwarf and sub-brown dwarf regime. We wish to investigate which of the results that have previously been found for higher mass objects also hold in this mass range. In addition, we will attempt to interpret the differences that arise as well as their possible impact on observations.

This paper is organized as follows. In Section 2, we describe the numerical method used, including the opacities and equation of state. The results of our simulations are then shown and interpreted in Section 3 and further discussed in Section 4, before we summarize our work in Section 5.

Finally, the Appendix deals with a numerical problem related to energy conservation in a self-gravitating system. This is an issue that may be important to other researchers doing similar work, and we are unaware of it being discussed in the literature previously, so we analyze it in some detail.

2. Method

2.1. Basic Equations

Our hydrodynamics module is based on Colella & Woodward 1984), although we use a spatially second-order (piecewise linear) Godunov scheme rather than their piecewise parabolic one. The governing equations of the system in Lagrangian conservation form are

$\begin{eqnarray}&&\displaystyle \frac{{dV}}{{dt}}=\displaystyle \frac{d({r}^{2}v)}{{dm}},\end{eqnarray} \tag{ 4 }$

$\begin{eqnarray}&&\displaystyle \frac{{dv}}{{dt}}=-{r}^{2}\displaystyle \frac{{dp}}{{dm}}-{\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 5 }$

$\begin{eqnarray}&&\displaystyle \frac{{dU}}{{dt}}=-\displaystyle \frac{d({r}^{2}{vp})}{{dm}}-v{\rm{\nabla }}{\rm{\Phi }}+{\left(\displaystyle \frac{{dU}}{{dt}}\right)}_{\mathrm{rad}},\end{eqnarray} \tag{ 6 }$

$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Phi }}=\displaystyle \frac{4\pi G}{V}.\end{eqnarray} \tag{ 7 }$

The equations are formulated in terms of the mass coordinate m(r), which is equal to the enclosed mass divided by 4π. V, v, U, and Φ are the specific volume, velocity, total specific energy, and gravitational potential. The effect of radiation is included as the term ${\left(\tfrac{{dU}}{{dt}}\right)}_{\mathrm{rad}}$ in the energy equation. This term is calculated according to the radiative transfer scheme for spherically symmetric systems that we developed in Stamer & Inutsuka (2018). This scheme is frequency-dependent (we use 40 frequency bins) and it is very accurate in all regimes of optical depth since it avoids any kind of radiative diffusion approximation. It calculates the radiative quantities in the laboratory frame, but under the assumption that the source function is approximately the same as in the comoving frame. We ignore radiation pressure throughout this work, since it only becomes relevant for very high-mass stars, which we are not dealing with here.

2.2. Equation of State

We use an equation of state (EOS) table developed by Kengo Tomida (Princeton University) and Yasunori Hori (National Astronomical Observatory of Japan), which is described in detail in Tomida et al. (2013, 2015) and whose effect is visualized in Figure 2. This table was calculated under the assumptions of local thermodynamical and chemical equilibrium. The latter is only valid if the timescale for chemical reactions is much shorter than the dynamical timescale, but since we are dealing with very dense gas this condition holds. Furthermore, the EOS table assumes solar abundance (X = 0.7, Y = 0.3; heavy metals are not considered), and a constant ratio of ortho- to parahydrogen of 3:1. This ratio affects the thermodynamic properties of the gas, but its actual value in molecular cloud cores is poorly constrained. We do not investigate its effect at this point. Finally, as the authors themselves note in the Appendix of Tomida et al. (2013), the EOS becomes unphysical in very high density regions (ρ > 0.1 g cm⁻³) since the neglected non-ideal effects are important in that region. This is not an issue for us since our calculations only reach densities about an order of magnitude lower.

**Figure 2.** Visual representation of the realistic EOS by Tomida et al. (2013). The black line shows a typical evolutionary track (Run C0) for the central gas element in a protostellar collapse process in the ρ-T-plane (see also Figure 1), while the colored background shows the value of the adiabatic index γ at that point. Note the initial decrease in temperature caused by the attenuation of background stellar radiation, which is only captured in frequency-dependent calculations (Section 3.1). The region below about 100 K corresponds to γ = 5/3, the value for a monatomic gas, due to the hydrogen molecule's rotational degree of freedom being inactive at these temperatures. At higher temperatures, it decreases to γ = 7/5, the value expected for a diatomic gas. The four bands of lowered γ correspond, from bottom to top, to four endothermic reactions: the dissociation of ${H}_{2}$ , the ionization of hydrogen, and the first and second ionization of helium.
Download figure:
Standard image High-resolution image

**Figure 2.** Visual representation of the realistic EOS by Tomida et al. (2013). The black line shows a typical evolutionary track (Run C0) for the central gas element in a protostellar collapse process in the ρ-T-plane (see also Figure 1), while the colored background shows the value of the adiabatic index γ at that point. Note the initial decrease in temperature caused by the attenuation of background stellar radiation, which is only captured in frequency-dependent calculations (Section 3.1). The region below about 100 K corresponds to γ = 5/3, the value for a monatomic gas, due to the hydrogen molecule's rotational degree of freedom being inactive at these temperatures. At higher temperatures, it decreases to γ = 7/5, the value expected for a diatomic gas. The four bands of lowered γ correspond, from bottom to top, to four endothermic reactions: the dissociation of ${H}_{2}$ , the ionization of hydrogen, and the first and second ionization of helium.
Download figure:
Standard image High-resolution image

2.3. Opacity

We use the monochromatic opacities by Vaytet & Haugbølle (2017).¹ These are a combination of the Semenov et al. (2003) dust opacities for temperatures below 1500 K, molecular opacities calculated based on Ferguson et al. (2005) for the temperature range between 1500 and 3200 K, and atomic opacities from the Opacity Project (Badnell et al. 2005) above 3200 K. The number of frequency points is different in each of these three sources of opacity data. For our program, we sample each of them into 40 log-equally spaced bins between 10¹⁰ and 10¹⁸ Hz. The actual opacity data used for a given shell during program execution depends on the shell's density and temperature. For instance, if the temperature is low, the Semenov et al. (2003) data is used, but if it is very high the Badnell et al. (2005) data is used instead.

2.4. Initial and Boundary Conditions

The simplest setup for our initial conditions is to envision the cloud core as a self-gravitating, homogeneous, isothermal sphere. In this case, we simply determine an initial density and radius, which together determine the mass. Due to the constant density, this is an extremely unstable setup. We will also consider the case of a sphere that is initially in hydrostatic equilibrium, i.e., a Bonnor–Ebert sphere (Bonnor 1956; Ebert 1955). The condition of hydrostatic equilibrium can be expressed as

$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dr}}=-G\displaystyle \frac{M(r)}{{r}^{2}},\end{eqnarray} \tag{ 8 }$

where M(r) is the enclosed mass at radius r. We assume a constant temperature of 10 K and set up the system by fixing the central density and initial radius and numerically integrating this equation outward from the center until that radius. The total mass is then determined, as is ${P}_{{\rm{O}}}$ , the pressure at the outer boundary. An important derived parameter is the Bonnor–Ebert mass, ${M}_{\mathrm{BE}}=1.18{c}_{{\rm{s}}}^{4}{P}_{{\rm{O}}}^{-1/2}{G}^{-3/2}$ , which is the critical mass above which the core is unstable to gravitational collapse. Thus for all of our calculations, $M\gt {M}_{\mathrm{BE}}$ .

The boundary condition at the center is v = 0, while the outer boundary is infalling. This infall velocity is calculated from the velocity of the outermost shell, under the assumption that v in that shell varies linearly in space at a given time step. For the radiative boundary condition, we employ the setup of Masunaga & Inutsuka (2000), where the heating from the outside is a combination of three separate contributions are as follows.

1.
The first contribution is heating due to cosmic rays, for which we adopt the value of $6.4\times {10}^{-28}n({{\rm{H}}}_{2})$ erg m⁻³ s⁻¹ (Goldsmith & Langer 1978). This is implemented simply as an additional, constant term in the energy equation. For this purpose, we assume that the gas consists entirely of H₂.
2.
The second contribution is a 2.7 K blackbody radiation representing the cosmic microwave background.
3.
The final contribution is a 6000 K blackbody radiation representing background stellar photons, multiplied by a dilution factor (10/6000)^x. We set x = 4.9, which together with the other two heating sources keeps the initial cloud close to 10 K.

Following Vaytet & Haugbølle (2017), we stop each run once the density in the central shell decreases from one time step to the next for the first time after the beginning of the second collapse. There are several reasons for this: first, it allows a comparison of different runs at the same point during their evolution, namely the point where the second collapse ends and a first hydrostatic bounce occurs. Second, our EOS becomes invalid at densities above ≈0.1 g cm⁻³, so we should stop the calculation before such densities are reached. And finally, the main accretion phase that follows the second collapse is difficult to deal with in a Lagrangian scheme because the time steps become extremely short, on the order of hours and less, and the shells inside the second core become extremely thin, even though we do not really need a high spatial resolution there. For this reason, an effective treatment of the further evolution requires a change in method, such as switching to a Eulerian grid or using an implicit scheme. Finally, it should be noted that the exact density at which the hydrostatic bounce occurs depends on the resolution.

We perform simulation runs with initial central densities between 5 × 10⁻¹⁷ and 5 × 10⁻¹⁵ g cm⁻³ and initial radii between 65 and 750 au. These are very dense and compact objects, which is necessary for them to be gravitationally unstable. The masses are between 0.8% and 10% of the solar mass, covering the brown dwarf and some of the sub-brown dwarf mass regime (the deuterium burning limit, conventionally used as the lower bound for brown dwarfs, lies at ≈1.3% of solar mass).

We consider two different initial density structures: one is a simple homogeneous density, which makes for a very unstable setup. The labels of these runs begin with the letter A. The second setup is a hydrostatic, Bonnor–Ebert sphere profile. The labels of these runs begin with the letter B. In addition, we include a one solar mass run with initial conditions identical to Masunaga & Inutsuka (2000) for comparison (run C0), as well as a number of runs to test for first coreless collapses (runs D0 through D3; see Section 3.4). A complete overview of all the runs is given in Table 1.

Table 1. Overview of Our Simulation Runs

Run	Profile	${\rho }_{{\rm{C}}}$ (g cm⁻³)	${R}_{{\rm{O}}}(\mathrm{au})$	$M({M}_{\odot })$	$\tfrac{M}{{M}_{\mathrm{BE}}}$	${t}_{\mathrm{ff}}(\mathrm{yr})$	${t}_{\mathrm{FC}}(\mathrm{yr})$	${\tau }_{\mathrm{FC}}(\mathrm{yr})$	${\rho }_{\mathrm{fin}}$ (g cm⁻³)
A0	Hom	5 × 10⁻¹⁵	58	7 × 10⁻³	2.4	9.42 × 10²	1.22 × 10³	1.49 × 10⁴	1.44 × 10⁻²
A1	Hom	1 × 10⁻¹⁵	104	8 × 10⁻³	1.24	2.11 × 10³	3.06 × 10³	6.81 × 10³	1.29 × 10⁻²
A2	Hom	1 × 10⁻¹⁵	112	1.0 × 10⁻²	1.54	2.11 × 10³	2.67 × 10³	3.87 × 10³	8.09 × 10⁻³
A3	Hom	1 × 10⁻¹⁵	119	1.2 × 10⁻²	1.85	2.11 × 10³	2.52 × 10³	3.96 × 10³	7.24 × 10⁻³
A4	Hom	1 × 10⁻¹⁵	129	1.5 × 10⁻²	2.36	2.11 × 10³	2.36 × 10³	2.27 × 10³	8.17 × 10⁻³
A5	Hom	5 × 10⁻¹⁶	181	2.1 × 10⁻²	2.30	2.98 × 10³	3.25 × 10³	9.04 × 10²	4.59 × 10⁻³
A6	Hom	5 × 10⁻¹⁷	463	3.5 × 10⁻²	1.22	9.42 × 10³	1.04 × 10⁴	9.91 × 10²	7.33 × 10⁻³
A7	Hom	5 × 10⁻¹⁷	548	5.8 × 10⁻²	2.02	9.42 × 10³	9.56 × 10³	6.38 × 10²	5.83 × 10⁻³
A8	Hom	5 × 10⁻¹⁷	586	7.1 × 10⁻²	2.47	9.42 × 10³	9.44 × 10³	5.45 × 10²	5.83 × 10⁻³
A9	Hom	5 × 10⁻¹⁷	622	8.5 × 10⁻²	2.95	9.42 × 10³	9.36 × 10³	4.59 × 10²	6.56 × 10⁻³
A10	Hom	5 × 10⁻¹⁷	659	1.0 × 10⁻¹	3.51	9.42 × 10³	9.33 × 10³	3.64 × 10²	6.58 × 10⁻³
B0	BE	5 × 10⁻¹⁵	65	7 × 10⁻³	1.93	1.10 × 10³	1.54 × 10³	1.09 × 10⁴	1.28 × 10⁻²
B1	BE	1 × 10⁻¹⁵	110	8 × 10⁻³	1.04	2.31 × 10³	3.83 × 10³	1.45 × 10⁴	1.02 × 10⁻²
B2	BE	1 × 10⁻¹⁵	120	1.0 × 10⁻²	1.28	2.35 × 10³	3.31 × 10³	4.78 × 10³	9.14 × 10⁻³
B3	BE	1 × 10⁻¹⁵	130	1.2 × 10⁻²	1.53	2.39 × 10³	3.07 × 10³	2.58 × 10³	9.16 × 10⁻³
B4	BE	1 × 10⁻¹⁵	140	1.5 × 10⁻²	1.79	2.43 × 10³	2.93 × 10³	2.07 × 10³	1.15 × 10⁻²
B5	BE	5 × 10⁻¹⁶	200	2.1 × 10⁻²	1.83	3.45 × 10³	3.92 × 10³	8.21 × 10²	6.30 × 10⁻³
B6	BE	5 × 10⁻¹⁷	500	3.5 × 10⁻²	1.08	1.04 × 10⁴	1.20 × 10⁴	1.06 × 10³	8.97 × 10⁻³
B7	BE	5 × 10⁻¹⁷	600	5.8 × 10⁻²	1.64	1.08 × 10⁴	1.11 × 10⁴	6.95 × 10²	8.97 × 10⁻³
B8	BE	5 × 10⁻¹⁷	650	7.1 × 10⁻²	1.93	1.10 × 10⁴	1.10 × 10⁴	6.39 × 10²	6.59 × 10⁻³
B9	BE	5 × 10⁻¹⁷	700	8.5 × 10⁻²	2.24	1.12 × 10⁴	1.09 × 10⁴	5.50 × 10²	8.23 × 10⁻³
B10	BE	5 × 10⁻¹⁷	750	1.0 × 10⁻¹	2.55	1.15 × 10⁴	1.08 × 10⁴	5.67 × 10²	9.31 × 10⁻³
C0	Hom	1.42 × 10⁻¹⁹	10000	1	1.85	1.77 × 10⁵	1.85 × 10⁵	9.01 × 10²	1.33 × 10⁻²
D0	Hom	3.28 × 10⁻¹⁹	6000	0.5	1.41	1.16 × 10⁵	1.28 × 10⁵	8.78 × 10²	1.57 × 10⁻²
D1	Hom	3.28 × 10⁻¹⁹	8000	1.18	3.34	1.16 × 10⁵	1.28 × 10⁵	6.20 × 10²	1.13 × 10⁻²
D2	Hom	3.28 × 10⁻¹⁹	11000	3.08	8.68	1.16 × 10⁵	⋯	⋯	3.77 × 10⁻³
D3	Hom	3.28 × 10⁻¹⁹	12000	4	11.27	1.16 × 10⁵	⋯	⋯	3.80 × 10⁻³

Note. The values are, from left to right, the run label, density profile, central density at the beginning of the simulation, initial radius, mass, ratio of mass to Bonnor–Ebert mass, free fall time, time until first core formation, first core lifetime, and the central density at the end of the run. The time of first core formation is defined as the first time step where the density in the central shell decreases from one time step to the next. The first core lifetime is the time elapsed between that same time step and the end of the simulation.

Download table as: ASCII Typeset image

2.5. Numerical Resolution

Figure 3 shows the final density profiles for our run C0 for different resolutions between 100 and 1000 shells. We find that the resolution has a slight effect on the profile, especially near the first core region. For more than 500 shells, further changes are very small, but require a significant increase in computational time. One of our most significant results concerns the lifetime of the first core (Section 3.3). For very low resolutions, its value increases significantly (1550 yr for 50 shells, 1100 yr for 100 shells), but converges to about 900 yr at 200 shells and more. We use 500 shells for all of the calculations discussed in this paper.

3. Results

3.1. Core Evolution

The evolution in the temperature-density plane of a fluid element near the center of the system was already shown in Figure 2. It generally behaves as expected (see also Figure 1), with a near-isothermal phase until ≈10⁻¹² g cm⁻³, followed by an adiabatic phase until ≈10⁻⁸ g cm⁻³ and the second collapse until ≈10⁻³ g cm⁻³, when hydrogen molecule dissociation is complete and the EOS again becomes stiffer.

It is noteworthy, however, that the initial contraction is not completely isothermal. The temperature oscillations at densities below 10⁻¹⁸ g cm⁻³ are a numerical artifact, but the following temperature decrease is significant. The temperature drops from the initial value of ≈10 K to a minimum of ≈6 K before increasing again. This effect is only captured in frequency-dependent radiative transfer calculations such as this work or Masunaga & Inutsuka (2000), but not in gray cases such as Vaytet & Haugbølle (2017) and Tomida et al. (2010). The cause is due to the increased opacity for radiation from the outside: as the cloud core starts to contract, its optical thickness begins to increase. Since the opacity is generally larger at higher frequencies, the system first becomes optically thick to high-frequency photons while still remaining optically thin to low frequencies. The inner regions are thus shielded from the background stellar radiation, so that heating is supplied only by cosmic rays and the cosmic microwave background, but the cooling radiation emitted from the dust can still escape. In conclusion, a numerical simulation must satisfy two conditions to show this effect: it must be frequency-dependent, and it must include high-frequency radiation from the outside (such as stellar background radiation) as a heating source.

3.2. Radial Profiles

The evolution of run C0's radial profiles of density, temperature, and velocity is plotted in Figure 4. The formation of the first and second cores at ≈10¹ and ≤10⁻² au, respectively, is evident as a sudden increase in density and temperature and as a maximum in the infall velocity (accretion shock). Snapshot 2 also shows the initial temperature decrease mentioned in the previous section, and it can be seen that this decrease affects the whole inner region of the cloud core, but not the outermost regions.

The final radial profiles (Figure 5) show that the second core's properties (density, radius) are essentially independent of the initial conditions. They are all very similar irrespective of the cloud core mass, initial density, and density structure. This agrees with the results of Vaytet & Haugbølle (2017) for higher masses, but here we show that it holds true in the brown dwarf regime as well.

However, we do find significant differences in the first core. Specifically, the very low-mass runs A0–A4 and B0–B4 look qualitatively different, with a very small first core (<1 au) that essentially makes up the whole system. The unusual behavior of these runs is closely related to the question of first core lifetime and is discussed in the following section.

3.3. First Core Lifetimes

The first hydrostatic core is generally believed to be a very short-lived object, which combined with the low luminosity should make it very hard to observe. Nevertheless, as was mentioned in the Introduction, a relatively large number of first core candidates have been observed. While none of these are confirmed to be first cores, and it is entirely possible that some or all of them may turn out to be more evolved objects, this large number of candidate observations may also indicate that at least some first cores are more long-lived and therefore more common than typically believed. Tomida et al. (2010) investigated first core lifetimes using 3D RHD simulations and found that for low-mass cores ( $\lt 0.1\,{M}_{\odot }$ ), the lifetime can be in excess of 10⁴ yr. The reason for this is that for the second collapse to start, the first core must reach a critical temperature of about 2000 K. How quickly that temperature is reached depends on the mass accretion rate. However, in a very low-mass case the mass reservoir may essentially be used up so that the accretion onto the first core becomes very weak before the temperature has risen sufficiently. This is what we see in the radial profiles of the very low-mass runs A0–A4 and B0–B4 (Figure 5): the first core basically makes up the whole system. If—somewhat arbitrarily—we define the surrounding envelope as the region where the temperature is below 100 K, then for example in the final profile of run B10, 44% of the total mass remain in the envelope. For run B0, it is a mere 2%. In this situation, radiative cooling then plays a crucial role in removing entropy from the system and thus allowing further contraction and heating in a Kelvin–Helmholtz-like process on a much slower timescale.

As shown in Table 1, our calculations confirm that the first core lifetime rises drastically for these very low-mass runs, when it becomes comparable to or even much larger than the free fall time. However, this only becomes evident below about $0.02\,{M}_{\odot }$ , not $0.1\,{M}_{\odot }$ as suggested by Tomida et al. (2010), and a lifetime of 10⁴ yr is only reached in the sub-brown dwarf regime. This is likely to be due to differences in the numerical method. As a 3D calculation, Tomida et al. (2010) include the effect of rotation, but on the other hand they use an idealized EOS with γ = 5/3 and a much more simplified radiative transfer scheme (gray flux-limited diffusion approximation).

Vaytet & Haugbølle (2017) find an anti-correlation of the first core lifetime with an instability parameter given by the ratio of mass to Bonnor–Ebert mass, where the lifetime drops at first slowly and then significantly for very large ratios of $M/{M}_{\mathrm{BE}}\gt \approx 5$ . However, the lowest mass they consider is 0.2 ${M}_{\odot }$ , so they do not probe the very low-mass regime, where the behavior is different. Figure 6 shows that for the very low-mass end, the first core lifetime depends strongly on the mass of the cloud core rather than on the instability parameter.

**Figure 6.** Upper panel: correlation of the first core lifetime ${\tau }_{\mathrm{FC}}$ with the mass. While there is little effect at masses larger than about ${10}^{-2}{M}_{\odot }$ , a very significant increase in the lifetime can be observed at extremely low masses. Lower panel: correlation of ${\tau }_{\mathrm{FC}}$ with the degree of instability. A decreasing tendency with increasing instability is still discernible, but the correlation is much less pronounced than for the mass.
Download figure:
Standard image High-resolution image

**Figure 6.** Upper panel: correlation of the first core lifetime ${\tau }_{\mathrm{FC}}$ with the mass. While there is little effect at masses larger than about ${10}^{-2}{M}_{\odot }$ , a very significant increase in the lifetime can be observed at extremely low masses. Lower panel: correlation of ${\tau }_{\mathrm{FC}}$ with the degree of instability. A decreasing tendency with increasing instability is still discernible, but the correlation is much less pronounced than for the mass.
Download figure:
Standard image High-resolution image

In conclusion, despite the quantitative differences with Tomida et al. (2010), our results confirm their prediction that very low-mass first cores can be very long-lived due to exhaustion of the accretion reservoir, and that generally the behavior of the first core lifetime is qualitatively different from higher mass ranges. This may account for the surprisingly large number of observational candidates. Furthermore, we predict that future discoveries of first cores will show a bias toward low-mass objects due to their longer lifetimes.

3.4. Collapses without First Cores

Vaytet & Haugbølle (2017) found that for certain very unstable setups ( $M/{M}_{\mathrm{BE}}\gt \approx 10$ ), the collapse can proceed without the formation of a first core. This occurs because the vigorous collapse creates such a high ram pressure that the thermal pressure does not become strong enough to balance it before the temperature reaches 2000 K and thus directly proceeds to the second collapse.

We confirm this behavior in our runs D0 through D3 (Figure 7). The initial conditions of D3 are identical to the ones Vaytet & Haugbølle (2017) used to investigate this phenomenon. The other three runs differ from D3 only in that their initial radii are smaller, which also implies a smaller mass and a smaller degree of instability. Figure 7 shows the transition between cases with and without the first core. The first core's signature in the density profile, clearly visible at around 10 au in runs D0 and D1, is absent in runs D2 and D3. The lower panel of Figure 7 confirms the importance of the ratio between thermal and ram pressure. In the first coreless cases, this ratio is smaller than unity everywhere outside the second core. The most favorable conditions for skipping the first core phase are those in which the ram pressure is very large (high degree of instability leading to large infall velocity), while the thermal pressure is relatively low. With respect to the initial temperature of the cloud, there is one difference between our result and Vaytet & Haugbølle (2017): they report that a collapse without the first core occurs only for an initial temperature of 5 K, but not for higher values such as 10 K or above, although our calculations show that it occurs even with 10 K. The boundary condition for all of our runs is calibrated such as to keep the initial temperature at 10 K, but the frequency-dependent radiative transfer then causes a decrease to ≈6 K (see Section 3.1), so that the first core phase can be skipped despite the relatively high initial temperature. This example illustrates that, even though in most cases the overall behavior of gray and frequency-dependent calculations is similar, in certain cases such as this they may produce significantly different results. In summary, the unique condition for the formation of a protostar without the first core phase is that the initial core mass must be much larger than the Bonnor–Ebert mass (M ≥ 5 ${M}_{\mathrm{BE}}$ ).

4. Discussion

4.1. Spectral Energy Distribution

We discuss the observational signature that would be expected from a brown dwarf in the process of formation. The total radiation escaping from the system is equal to what we call the boundary term in our radiative transfer scheme (Stamer & Inutsuka 2018). By considering only a single direction for the escaping radiation, namely the direction toward the observer, we may use this to compute the spectral energy distribution (SED) that would be observed for a given beam size at a given distance. For the distance, we assume 150 pc, and for the beam size we vary between 1 and 1000 au (Figure 8). For the calculation of the SEDs only, we increase the number of frequency bins to 400 for better frequency resolution. We note that for a relatively high-mass case such as B10, the emission is dominated by the cold envelope at all stages, i.e., the hot core region, is obscured. This is also the main reason why observation of the early phases of star formation is difficult. The situation is exacerbated in the case of brown dwarfs because of their even lower luminosity compared to solar-type stars. A smaller beam size helps to better distinguish different stages of the collapse, since the effect of the contribution of the envelope is reduced. These results essentially confirm those of Masunaga et al. (1998) and Masunaga & Inutsuka (2000), although the latter also found that the SED becomes markedly different in the following main accretion phase, which we have not considered. We do not show the SEDs of the A runs or any other high-mass B runs since there are no significant differences. However, the situation is somewhat different for the very low-mass cases, as seen in the lowest panel of Figure 8. Here, the surrounding envelope has disappeared in the later stages, allowing the observer to look deeper into the warm interior. Therefore the SED is clearly warmer than at the beginning, and it is possible to look at a naked first core. Another difference is that for the very low-mass case, the initial temperature decrease is less pronounced (run B0 never goes as low as 6 K, unlike run B10). Recall that the temperature decrease is caused by the blocking of the optical background radiation while low-frequency cooling radiation is still able to escape. Comparing runs B0 and B10, the initial density of B0 is larger by 2 orders of magnitude while the initial radius is smaller by about 1 order of magnitude. In consequence, B0's optical depth is about 1 order of magnitude larger, and in fact the initial state can be described as optically intermediate. As a result, the cloud becomes fully optically thick and the cooling radiation is blocked as soon as the collapse begins, and therefore the central temperature hardly drops at all.

**Figure 8.** Spectral energy distribution at various times for run B10 (upper two panels) and run B0 (lowest panel), at an assumed distance of 150 pc and with different beamsizes. T_C is the central temperature at the time of the snapshot. The 6-K plot is missing in the lower panel since run B0 never reached such low core temperatures. The evolution becomes clearer with smaller beamsizes, and is especially visible for the very low-mass cases where the obscuring envelope has disappeared.
Download figure:
Standard image High-resolution image

Finally, we discuss the wave-like pattern that appears in most of the SEDs near log ν = 12.4. It stems from the Semenov et al. (2003) monochromatic opacity data, which shows this kind of pattern at the same frequency range. This opacity fine structure is only visible when viewing through an optically thin region, such as the initial cloud or, in later stages, the surrounding envelope. In the very low-mass case, it is mostly hidden in the later stages due to the lack of an envelope. It is even hidden in the initial state since in this run, the initial optical depth at the relevant frequencies is around unity.

4.2. Observational Implications

There are so far no confirmed observations of a first hydrostatic core, even though a number of candidate objects have been proposed. As a result of our simulations, we have found two effects that may affect future observations. First, the lifetime of the first core can be increased by about an order of magnitude for extremely low-mass objects due to depletion of the envelope. Second, the lack of an obscuring envelope for these same objects exposes the first core and allows more high-frequency radiation to escape. Both of these effects should make the observation of very low-mass first cores easier than it would be otherwise. The longer lifetime increases the number of objects that exist, and the warmer SED makes them easier to identify and to separate them from more evolved objects. While their extremely low intrinsic luminosity makes them very hard to detect, nonetheless, we propose that detection is somewhat less difficult than would be expected naively. Since our simulation is spherically symmetric and does not include the effects of rotation and magnetic fields, our quantitative results (i.e., the first core lifetime for a given mass) may not be very reliable. Qualitatively, however, we predict a bias toward very low-mass objects in future confirmed observations of first cores.

5. Conclusion

In this paper, we performed numerical simulations of spherically symmetric protostellar collapse, with a special focus on very low-mass objects, i.e., brown dwarfs. Our numerical model combines a second-order Godunov hydrodynamic module with the frequency-dependent radiative transfer scheme we developed in a previous paper (Stamer & Inutsuka 2018). We also include a realistic EOS to account for the effects of hydrogen dissociation, as well as frequency-dependent opacities for dust and gas. Overall, this is still a very simplified model due to the spherical symmetry, which ignores multi-dimensional effects such as magnetohydrodynamics and rotation. Nevertheless, it is useful in investigating certain properties of the collapse process and of the objects formed as its result.

When we were exploring numerical results with various initial conditions, we noticed a numerical problem with the hydrodynamic code that may be relevant to other investigations of self-similar gravitational collapse processes. To our knowledge, the solution for the problem has not been mentioned in the literature. Specifically, we found that in this self-gravitating system, a violation of energy conservation occurs if the effect of gravity in the energy equation is calculated according to the method of Colella & Woodward (1984). Over time, this causes a spurious increase in the total energy, which in certain cases can blow the whole system apart due to the increased pressure. We provide a method to overcome this problem by explicitly calculating the gravitational energy in each shell before and after each time step, as described in the Appendix.

Our collapse calculations showed that the properties of the protostar formed are mostly independent of the initial conditions of the cloud core, such as its mass and radius. This had previously been established for solar mass and larger stars, but we confirmed that it also holds true at the low-mass end. The first core phase, however, showed some novel behavior. For the lowest mass runs near the deuterium burning limit, the first core lifetime increases dramatically when lowering the mass and can exceed 10,000 yr. This is due to the depletion of the object's gas envelope. At such low masses, essentially the entire cloud core's mass precipitates onto the first core, causing the accretion rate to become very small. This slows the rate of temperature increase inside the first core and therefore delays the onset of the second collapse, which requires a temperature of around 2000 K.

This envelope depletion also has observational implications. We calculated the SEDs that would be observed from our simulations, and we found that the first core, which is usually obscured by the cold envelope, can be significantly more visible in the very low-mass cases where the envelope has disappeared. Combined with the much longer lifetime of these first cores, we conclude that very low-mass first cores should be somewhat easier to observe than might be expected otherwise.

On the opposite end of the spectrum of first core lifetimes, we also confirmed that for sufficiently unstable setups, the first core stage may be skipped entirely, as has been previously shown by Vaytet & Haugbølle (2017). A first core only forms if at some point the thermal pressure is larger than the ram pressure of the infalling gas. For very unstable and/or cold setups, this may never occur so the collapse directly proceeds to the formation of the second core. In contrast to Vaytet & Haugbølle (2017), however, our simulations do not require a low initial temperature, since our frequency-dependent radiative transfer always causes the central temperature to drop to ≈6 K due to blocking of the high-frequency stellar background radiation.

For our future work, we hope to expand these calculations, specifically our radiative transfer scheme, to multiple dimensions in order to obtain more realistic results. We are especially interested in the behavior of the first core lifetimes and their observational implications once effects such as rotation, magnetic fields, and outflows are included. In addition, with some modifications our code should also be applicable to the formation of planets through disk instability. Apart from the different boundary conditions that result from embedded in a protoplanetary disk, this is a situation which is very similar to protostellar collapse.

Appendix: Spurious Re-expansion After Second Core Formation

When we were testing our program, we encountered an unexpected behavior. After the formation of the second core, the system expanded again to first core densities and below. Suspecting a numerical error, we analyzed the development of the total energy, shown in the upper panel of Figure 9. In this plot, the formation of the first core takes place around time step 5000, and that of the second core around time step 27,000. In the early phase, the internal energy oscillates rapidly as the system's energetics is dominated by the radiative exchange with the background. The magnitude of these oscillations depends on the size of the time step. During the first core phase, we see a different kind of oscillation in the internal and gravitational energy. This is because the first core periodically expands and contracts slightly. During this time, the total energy slowly decreases due to radiation escaping from the system. The problem becomes apparent at the formation of the second core, when the total energy starts to increase steeply. This is obviously unphysical, since at this stage, radiative cooling should far exceed radiative and cosmic-ray heating, and there are no other diabatic processes. In fact, this spurious increase occurred even when completely deactivating radiation and cosmic rays, implying that the cause is in the hydrodynamic module.

**Figure 9.** Evolution of the system's total energy and its individual constituents in the original calculation (upper panel) and in the corrected version (lower panel). See the text for details.
Download figure:
Standard image High-resolution image

We found that the term −v∇Φ in Equation (6) (the energy equation) was responsible for this error. This term represents the energy gained by the fluid due to conversion of gravitational potential energy. Colella & Woodward (1984) include a non-specified external body force g in their scheme, in which case the numerical implementation of the resulting term vg is an average of the value before the current time step and the value after the time step, i.e., (vⁿgⁿ + vⁿ⁺¹gⁿ⁺¹)/2. However, in our self-gravitating system, where g = −∇Φ is not constant, this approach is incorrect. The increase in specific energy when moving the fluid element in a gravitational field is dU = −dr∇Φ. Taking the derivative with respect to time yields dU/dt = −v∇Φ, but only if ∇Φ is independent of time, which in our case it is not.

This inaccuracy leads to a continuous increase in the total energy of the system. In other words, the increase in the specific energy in each time step is larger than the corresponding decrease in gravitational potential energy. In the early stages of the collapse and during the first core phase, the effects of this error are too small to counteract the radiative losses, but during second core formation much more gravitational energy is converted to internal energy, and the abnormally large increase in pressure ultimately blows apart the second core.

We have solved this problem by explicitly calculating the gravitational energy in each shell before and after the time step, and inserting the difference in the energy equation in place of the term −v∇Φ. The gravitational potential energy of any shell is given by

$\begin{eqnarray}&&{E}_{\mathrm{grav}}=-G{\int }_{{m}_{{\rm{I}}}}^{{m}_{{\rm{O}}}}\displaystyle \frac{m}{r(m)}{dm},\end{eqnarray} \tag{ 9 }$

where m is equal to the enclosed mass and ${m}_{{\rm{I}}}$ and ${m}_{{\rm{O}}}$ are the values of this mass coordinate at the inner and outer boundary of the shell. Expressing m as a function of the radial coordinate r, we obtain

$\begin{eqnarray}&&{E}_{\mathrm{grav}}=-G{\int }_{{r}_{{\rm{I}}}}^{{r}_{{\rm{O}}}}\displaystyle \frac{{m}_{{\rm{I}}}+\tfrac{4\pi }{3}\rho \left({r}^{3}-{r}_{{\rm{I}}}^{3}\right)}{r}4\pi \rho {r}^{2}{dr}.\end{eqnarray} \tag{ 10 }$

The integral in the above equation can be solved analytically as

$\begin{eqnarray}&&\begin{array}{l}4\pi \rho \displaystyle \int \left({m}_{{\rm{I}}}+\displaystyle \frac{4\pi }{3}\rho \left({r}^{3}-{r}_{{\rm{I}}}^{3}\right)\right){rdr}\\ \quad =\ 4\pi \rho \left(\left(\displaystyle \frac{{m}_{{\rm{I}}}}{2}-\displaystyle \frac{2\pi }{3}\rho {r}_{{\rm{I}}}^{3}\right){r}^{2}+\displaystyle \frac{4\pi }{15}\rho {r}^{5}\right)+C.\end{array}\end{eqnarray} \tag{ 11 }$

Since all necessary values are known both before and after the time step, the change in gravitational energy can be calculated and inserted into the energy equation.

The lower panel of Figure 9 shows the behavior after applying this correction. The main difference is that there is no longer any increase in total energy at the point of second core formation, which here occurs around time step 37,000. In addition, the oscillations during the first core phase are significantly smaller.

Figure 10 shows the evolution of the central density before and after the correction. In the improved version, the second core is stable and no unphysical re-expansion occurs.

Radiation-hydrodynamic Simulations of Spherical Protostellar Collapse for Very Low-mass Objects

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction