A HYBRID SCENARIO FOR THE FORMATION OF BROWN DWARFS AND VERY LOW MASS STARS

The formation of brown dwarfs (BDs) is usually attributed to either of two distinct mechanisms. One model is the direct collapse of a very low mass cloud fragment whose mass straddles the substellar mass limit. Formation mechanisms for such low-mass (and consequently high-density) cores include colliding flows in a turbulent magnetic medium (Padoan & Nordlund 2004) and photoerosion of more massive cores by the ionizing UV radiation from nearby OB stars (Whitworth & Zinnecker 2004). Since the collapse of a prestellar core has been shown by numerous studies to lead to a very focused collapse with a power-law density profile, the probability of fragmentation during the runaway collapse phase (i.e., before the formation of a central hydrostatic object) is low. This leads to the scenario that the stellar and substellar mass function is simply a scaled-down version of the core mass function.

However, during the protostellar phase (i.e., after the formation of a central hydrostatic object), and unlike in the prestellar phase, the relative scaling of flow variables in the presence of angular momentum can lead to the formation of a centrifugally balanced disk (Basu & Mouschovias 1995a). Such a disk has a high density and evolves on a timescale much longer than its dynamical time. Therefore, the protostellar disk has time for fragmentation if conditions are suitable, e.g., the criteria for gravitational instability and fragmentation are satisfied (Toomre 1981; Gammie 2001). Furthermore, fragmentation may be induced by close encounters with members of a stellar cluster in otherwise gravitationally stable disks (e.g., Thies et al. 2010). This leads to the second model: that disk fragmentation is responsible for the formation of low-mass and substellar objects in the vicinity of a forming star, even if the mass of the parent core is far greater than the substellar limit. In this scenario, free-floating BDs would be the result of ejections from the disk, which is a plausible outcome based on extrapolation from three-body calculations of interacting point masses and can account for peculiar properties of BDs as compared to those of low-mass stars (Kroupa & Bouvier 2003). Bate (2009) obtained BD ejection from smoothed particle hydrodynamics (SPH) modeling of molecular cloud fragmentation resolved from parsec scales down to the opacity limit that occurs on AU scales. In that model, multiple collapsing fragments each form a disk. Some disks undergo fragmentation, and low-mass objects, treated as point masses in the SPH approximation, can be ejected from multiple systems. In the models of Bate (2009), a sink cell is introduced on AU scales and gas is cleared away from a region around them of size up to 10 AU. Furthermore, Stamatellos & Whitworth (2009) have modeled the formation and ejection of BDs after disk fragmentation, also using SPH and sink cells, and starting from an initial condition of a pre-existing massive disk that has a mass equal to that of the central star.

In this paper, we pursue models that solve the self-consistent formation of a centrifugal disk from the collapse of a larger cloud core, employing the thin-disk approximation. Computational efficiency allows us to run a large number of models, with various initial conditions of cloud mass and angular momentum, as well as study the long-term evolution of disks. We employ a central sink cell of size 6 AU, representing a central star and some circumstellar material, but the surrounding disk that forms is followed without recourse to sinks, even when fragmentation occurs. Our main result is that ejection of protostellar-/proto-BD-mass clumps occurs through the interaction of multiple gravitationally bound gaseous clumps, which are essentially the first Larson cores formed via gravitational fragmentation of protostellar disks.

Our results support a new hybrid paradigm of BD and low-mass star formation by disk fragmentation followed by clump ejection (rather than by ejection of finished low-mass stars/BDs). This scenario naturally accounts for the presence of disks around BDs and the relatively low ejection speeds in this scenario can account for the observed velocity dispersion and physical location of BDs relative to young stellar objects (YSOs). It is also consistent with the so-called brown dwarf desert, i.e., the lack of BD companions to primary stars at distances less than several tens of AU (McCarthy & Zuckerman 2004; Marcy et al. 2005). A disk fragmentation origin for BDs can explain the brown dwarf desert (see, e.g., Stamatellos & Whitworth 2011) since cooling constraints in disks (Rafikov 2007) favor fragmentation at radii ≳ 50 AU from the primary object. Furthermore, fragments that do form at closer range or are torqued inward from wider orbits are often subject to a runaway inward migration toward merger with the central object (Vorobyov & Basu 2005, 2006), so that low-mass gaseous companions are not typically expected to settle into stable inner orbits. However, it has been suggested that some of the eventually tidally disrupted clumps may have enough time to build up solid cores. These may then be released into inner orbits after the disruption, providing thermally processed solids or even rocky planetary cores (Boley et al. 2010; Nayakshin 2010; Cha & Nayakshin 2011; Vorobyov 2011). This scenario may account for planet formation in the inner disk but presumably not something as massive as a BD. Detailed simulations of this possible mechanism remain to be worked out.

In contrast, the direct core collapse paradigm for BD formation also suffers from the need to explain very low mass gravitationally bound objects that are one to two orders of magnitude less than the ambient Jeans mass in molecular clouds. This can be accomplished with very high bounding pressure, whether due to very strong magnetic fields (Basu & Mouschovias 1995b) or high turbulent pressure (Padoan & Nordlund 2002; Chabrier & Hennebelle 2011). Both the ideas of extremely high magnetic or turbulent bounding pressure on a fragment suffer from a lack of direct observational evidence at the moment (see André et al. 2009). Turbulent pressure may also not be present in such a focused manner as to act like an isotropic bounding pressure.

Our model and method of solution are the ones presented in detail by Vorobyov & Basu (2010), which is an update to the model presented by Vorobyov & Basu (2006). We solve the mass, momentum, and energy transport equations in the thin-disk approximation, using vertically integrated quantities. Disk self-gravity is calculated using a two-dimensional Fourier convolution theorem for polar coordinates. The energy equation includes compressional heating, heating due to stellar and background irradiation, radiative cooling from the disk surfaces, heating due to turbulent viscosity (parameterized via the usual Shakura & Sunyaev prescription with α = 0.005), and a diffusion approximation to link the effective surface temperature of the disk with the midplane temperature. Frequency-integrated opacities are adopted from the calculations of Bell & Lin (1994) and a smooth transition is introduced between the optically thin and optically thick regimes. In the opaque regime, gravitationally bound clumps that form in the disk are essentially "first hydrostatic cores" or "first Larson cores" in the parlance of star formation (Larson 1969; Masunaga & Inutsuka 2000). A "second collapse" down to stellar dimensions can take place when the central temperature reaches ∼2000 K and H₂ is dissociated. This possibility is not allowed in our energy equation so that we may continue to track the evolution without recourse to sink cells;⁴ however, clumps in the disk typically do not reach this temperature before disruption or ejection (see Section 4). The stellar irradiation is based on a luminosity that is a combination of accretion luminosity and photospheric luminosity calculated from the pre-main-sequence tracks of D'Antona & Mazitelli (1997). The pressure is related to the internal energy through an ideal gas equation of state.

Models presented in this paper are run on a polar coordinate (r, ϕ) grid with 512 × 512 zones. A central sink cell of radius 6 AU is employed and the potential of a central point mass is added to the disk self-gravity once a central hydrostatic core is formed. The radial points are logarithmically spaced, with the innermost cell outside the central sink having size 0.07–0.1 AU depending on the cloud core size (i.e., the radius of the computational region). The latter varies in the 0.025–0.12 pc (5000–24,000 AU) limits. The radial and azimuthal resolution are about 1 AU at a radius of 100 AU. The inner and outer boundary conditions are set to allow for free outflow from the computational domain. The initial conditions for the prestellar core correspond to profiles of column density Σ and angular velocity Ω of the form

$$\begin{eqnarray} \Sigma & = & {r_0 \Sigma _0 \over \sqrt{r^2+r_0^2}}\:, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \Omega & = &2\Omega _0 \left({r_0\over r}\right)^2 \left[\sqrt{1+\left({r\over r_0}\right)^2 } -1\right]. \end{eqnarray} \tag{ 2 }$$

These profiles have a small near-uniform central region of size r₀ and then transition to an r⁻¹ profile; they are representative of a wide class of observations and theoretical models (Basu 1997; André et al. 2009; Dapp & Basu 2009). Our previous papers have also shown that the qualitative features of the ultimate centrifugal disk evolution are independent of the specific profiles adopted for the initial prestellar core. Our reference model presented in this paper has rotational support in the initial state characterized by the dimensionless parameter η ≡ Ω²₀r₀²/c²_s = 1.4 × 10⁻², where c_s = 0.188 km s⁻¹ (corresponding to a temperature T = 10 K) in the initial state but varies spatially and temporally during the evolution. Alternately, the ratio of rotational energy to the magnitude of gravitational potential energy is β = 1.3 × 10⁻². Table 1 shows a list of parameters and some results from models we have run for this study. The reference model of our study (model 5) has r₀ = 2006 AU, a total core mass M_c = 0.92 M_☉, and β = 1.3 × 10⁻².

Our recent numerical hydrodynamics simulations with realistic cooling have shown that fragment formation is enhanced at radii ≳ 50 AU (in agreement with many previous studies on disk fragmentation), and that episodic accretion due to clump formation and inspiralling of the clumps is a robust result (Vorobyov & Basu 2010). In this paper, we report on a qualitatively new kind of result—an occasional ejection of fragments from the disk into the intracluster medium. Long integration times (∼1 Myr) need to be run in order to capture the occasional ejection events. Below, we discuss in detail the results for our prototype model.

Figure 1 shows a sequence of column density images for our reference model (model 5), at various times after the formation of a central object, and in a region of size 2400 AU on each side. The circumstellar disk forms at t ≈ 0.01 Myr and already at t = 0.05 Myr it undergoes fragmentation. Within the first ≈0.25 Myr, multiple fragments are formed in the relatively massive disk, at distances of ≳  50 AU and ≲ few hundred AU. They are generally torqued inward through gravitational interaction with trailing spiral arms, as first found by Vorobyov & Basu (2005, 2006). Others located at large radii may eventually disperse. However, under sometime-favorable conditions, a clump within a multi-clump environment can be ejected through many-body interaction. The ejection is also aided by the non-axisymmetric potential of the relatively massive disk. The arrow in the image in the second row of the middle column points to a clump that is subsequently ejected from the system. The velocity of the ejected clump (at the moment when it leaves the computational boundary at r = 12, 000 AU) is about 0.9 km s⁻¹, which is a factor of three greater than the escape velocity v_esc = (2GM_sd/r)^1/2, where the total mass of the central star and its disk is M_sd = 0.6 M_☉. This means that the clump will truly be lost to the system. The ejection event is transient and the clump is not seen in the subsequent image at t = 0.27 Myr. The total mass of the clump, calculated as the mass passing the computational boundary during the ejection event, is 0.15 M_☉. It should be noted, however, that this value includes not only the compact core but also a diffuse envelope and even fragments of spiral arms (see Figure 2). We speculate that upon contraction this clump may form a substellar object, given that a significant fraction of mass remains in the disk until it is ejected due to outflows and/or dispersed due to photoevaporation. No further major fragmentation episodes in the protostellar disk are seen after t = 0.27 Myr, likely due to a sharp drop in the total disk mass caused by the ejection. The disk then gradually evolves toward a nearly axisymmetric state.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To illustrate the ejection mechanism in more detail, Figure 2 shows a sequence of column density images on a much larger spatial scale of the simulation box, 20,000 AU on each side, and on a much shorter timescale, t = 0.26–0.31 Myr. The arrows in the images identify the same clump that was identified in Figure 1, and show its ejection from the disk and indeed the entire simulation box. As the clump moves nearly radially outward, it grows in size and diffuses somewhat but this is an artifact of the decreasing resolution of our numerical grid at larger radii. Future higher resolution simulations will be needed to verify the evolution of the clumps after they undergo an ejection event.

Figure 3 shows the mass accretion rate onto the central sink as a function of time since the beginning of the simulation. The time t = 0.12 Myr represents the moment that a central object is formed. After a brief period of smooth accretion, a disk forms outside the sink cell, and the subsequent accretion is highly variable, with the highest rates exceeding 10⁻⁴ M_☉ yr⁻¹ and the lowest rates plunging below 10⁻⁹ M_☉ yr⁻¹. Interestingly, there is a late burst at 0.38 Myr (0.26 Myr after the formation of the central star) that is greater in magnitude than all previous bursts. This burst corresponds to the inward flow of material that accompanies the ejection of the clump. The effects are twinned due to the need for overall angular momentum conservation—one clump is ejected while the other one is driven into the center. After this last major burst, the disk settles down to a more quiescent phase of evolution with a much lower accretion variability, characterized by flickering in the accretion rate but no bursts, and a gradual decline in accretion rate. The residual accretion rate of a few × 10⁻⁸ M_☉ yr⁻¹ at 1.0 Myr is consistent with observed accretion rates for T Tauri stars (e.g., Calvet et al. 2004). This behavior agrees well with a gradual gravitational stabilization of the disk that is caused the ejection event and the associated substantial mass loss. From this time onward, viscous torques start to dominate gravitational torques in the disk mass transport (Vorobyov & Basu 2009).

Zoom In Zoom Out Reset image size

Figure 4 illustrates the evolution of mass in various components of the system, as well as in the overall simulation region. The total mass in the system suffers a drop as the ejected clump exits through the outflow outer boundary applied to the simulation region. The mass loss equals the clump mass: ≈0.15 M_☉. The clump maintains its form all the way to the edge of our simulation box notwithstanding a constantly decreasing resolution on the logarithmic grid, and is gravitationally bound, with the magnitude of its gravitational potential energy being greater than the sum of its rotational and thermal energies. Hence, we expect it to form a low-mass star or substellar object depending on the efficiency with which it converts into a collapsed object. Figure 4 also shows a notable sharp drop in the disk mass and an equivalent increase in the envelope mass at t = 0.26 Myr. These effects are the manifestation of the ejection event when the clump leaves the parent disk and joins the envelope. At about the same time, the stellar mass increases sharply due to consumption of the twin clump that has lost its angular momentum due to the many-body interaction.

Zoom In Zoom Out Reset image size

We have considered 20 models that exhibit disk fragmentation. The corresponding initial prestellar core masses lie in the range M_c = 0.3–1.8 M_☉, and ratios of rotational to gravitational energy are in the range β = (0.25–2.3) × 10⁻². Only 11 models have shown ejection events and the parameters of these models are described in Table 1. A star symbol next to the model number means that the ejected clump was initially a binary clump, and a dagger symbol means that the clump is smeared out before reaching the computational boundary. Both these phenomena are discussed later in this section. The last five columns present the mass of the central star at the moment of ejection M_s, the mass M_eject of the ejected clump, the time t_eject after the formation of the star when the ejection event takes place, the ejection velocity v_eject, and the escape velocity v_esc, respectively. Both v_eject and v_esc are calculated at the time when the clump passes through the outer computational boundary. A diamond symbol indicates that the corresponding value is calculated at a smaller distance, since the ejected clump has been smeared out before reaching the outer computational boundary. When more than one ejection event is present in the same model, the different ejection speeds are separated by a slash. In all cases, v_eject > v_esc, indicating that the ejected clumps will be lost to the system and become freely floating objects. The mean and median ejection velocities in our models are 0.7 km s⁻¹ and 0.6 km s⁻¹, and the minimum and maximum values are 0.38 km s⁻¹ and 1.2 km s⁻¹, respectively. These are consistent with the velocities of low-mass stars and BDs in the Chamaeleon I and Taurus star-forming regions (e.g., Joergens 2006). There is practically no dependence of the ejection velocities on the mass of the ejected clumps.

Figure 5 presents the total mass in the system and illustrates the ejection events in 11 models of Table 1. A sharp decrease in the total mass indicates the time when the clump escapes the computational region. However, approximately half of the models (marked with a dagger symbol in Table 1) demonstrate a rather smooth decline in the total mass, extended over (1–2) × 10⁵ yr. In these models, clumps have dispersed before reaching the outer computational boundary. They may become gravitationally unbound due to strong tidal torques during the close encounters that lead to ejection, and/or due to decreasing numerical resolution on a logarithmically spaced grid as they propagate radially outward.

Zoom In Zoom Out Reset image size

Table 1 demonstrates that some of the ejected clump masses M_eject are below the substellar mass limit, giving us more confidence to propose this as a potential BD formation mechanism. Note that because this process occurs due to the interactions of gaseous clumps and a gaseous disk, it is somewhat different from the scenario of ejections by multi-body interactions of point masses, i.e., objects that have already fully collapsed down to compact equilibrium configurations.

Analysis of the ejection times in Table 1 and Figure 5 shows that most ejection events take place between 0.2 Myr and 0.8 Myr after the formation of the central star and no ejections are seen after t = 1 Myr. This is due to the fact that gravitational instability and fragmentation are limited exclusively to the early stages of disk evolution.

It is interesting to note that models 6, 7, and 11 experience multiple ejection events. The same may be true for model 10 but the ejected masses at t ≈ 0.25 Myr and t ≈ 0.4 Myr are too small for the events to be resolved. In addition, in models 2, 10, and 11 (marked with the star symbol), the ejected clumps were initially close binaries with separations of order a few AU, but later merged as they propagated toward the outer boundary, most likely due to decreasing resolution at large radii. The possibility of ejection of close binary clumps is a particularly appealing feature of the clump ejection scenario for BD formation since the separation of binary BDs seems to peak at around 1–4 AU (Whitworth et al. 2007).

Figure 6 shows the β − M_c phase space covered in our modeling. The dark-shaded area highlights the region where both disk fragmentation and clump ejection are observed, while the light-shaded area defines the region where only disk fragmentation takes place. There is a increasing tendency for ejections to take place as β and M_c increase. Models that fall into the light-shaded domain usually form just one fragment at a time, in which case the ejection mechanism cannot work even in principle. Furthermore, not all models that are in the ejection domain may undergo ejection events due to the highly stochastic nature of this mechanism. This means that the dark-shaded area outlines the region where ejection is possible in principle, but in reality it might not be observed if a chance arrangement of clumps would not favor their close encounters. The robustness of Figure 6 is also affected by a rather coarse grid of models explored in the paper. For instance, we have considered only four values of β and about ten models for each β. The uncertainty associated with a small number of models is illustrated in Figure 6 by arrows showing where the fragmentation (and ejection) regions could potentially extend, had we considered a finer grid of models. We restricted our modeling with M_c ⩽ 2.0 M_☉ and ratios β ⩽ 0.025 due to the fact that models with higher initial core masses and angular momenta tend to form wide-separation (>10 AU) binary/multiple stellar systems for which our numerical code is not well suited. Disks around such systems may be strongly truncated due to tidal effects, which may impede disk fragmentation and ejection. On the other hand, close encounters with other members of a stellar cluster may promote disk fragmentation and ejection (Thies et al. 2010). A dedicated study using suitable numerical codes is needed to address these issues. We also note that a steep cutoff of the fragmentation region at small β ≲ 0.001 is a physical effect related to the fact that such models form disks of too low a mass and size for gravitational fragmentation to take place (Vorobyov 2011).

Zoom In Zoom Out Reset image size

We find that high-resolution models with a realistic treatment of heating and cooling and in which a disk forms consistently from prestellar collapse can sometimes lead to the ejection of gaseous clumps from the disk and eventually from the system (defined by the extent of the original prestellar core). These are the first confirmed examples of ejections due to interactions of compact gaseous clumps and their parent gas disks and host star, as opposed to being caused by the interaction of point particles. Previous numerical simulations of ejections from disks by Bate (2009) and Stamatellos & Whitworth (2009) have used sink particles with a smoothing length as a free parameter. While the smoothing length can in principle be used to mimic the size of a clump, our numerical simulations are free from such uncertainty and the size of the clump is instead determined self-consistently by the disk thermodynamics. We show explicitly that the clumps formed via disk fragmentation need not contract to stellar densities in order to be ejected into the intracluster medium. Moreover, some of the clumps may be tidally disrupted during the ejection event and disperse as they leave the system, a phenomenon that is impossible to capture using sink particles. In fact, most of the clumps in the disk never achieve temperatures required for a second collapse before they are disrupted, as pointed out previously (Nayakshin 2010; Boley et al. 2010; Vorobyov 2011). For instance, Figure 2 (bottom-left panel) in Vorobyov (2011) shows that most of the clumps have maximum temperatures below 1000 K and only one clump achieves a temperature higher than 2000 K before being disrupted or torqued onto the star. This makes the phase of interaction of multiple first cores, as studied in this paper, all the more relevant to disk evolution.

An important argument usually in favor of the direct collapse scenario is the observation of some very low mass isolated proto-BD clumps (Luhman et al. 2007). However, by allowing for the possibility of ejected clumps rather than finished BDs, our hybrid scenario also allows for this possibility. Clumps that survive can cool and contract to form a (sub-)stellar core and a disk, thus also explaining the presence of accretion disks around BDs. Furthermore, the ejection speeds of the clumps are <1 km s⁻¹ ≈1 pc Myr⁻¹, reflecting the fact that ejections take place in a much more diffuse situation than interacting point masses in close proximity, i.e., with impact parameters of tens of AU rather than several radii of finished BDs or low-mass stars. With these ejection speeds, the velocity dispersion of BDs is not expected to differ much from that of YSOs and they should remain spatially co-located within cluster forming clumps of parsec scale for at least the typical 1–2 Myr age of YSO clusters. Our models can also explain qualitatively why the BDs are less numerous than stars, since the ejections only occur for some region of parameter space. While a quantitative prediction of BD fraction remains difficult without more insight into core initial conditions, we can use a stellar initial mass function (IMF) to obtain some estimates. Assuming that about half the mass of each clump accretes onto a star (i.e., 50% formation efficiency), that only clumps with mass M_c ⩾ 0.8 M_☉ undergo ejection events, and using the Kroupa stellar IMF (Kroupa 2001) with lower and upper cutoff masses 0.08 M_☉ and 100 M_☉, respectively, we estimate that the number of ejected clumps is about one for every 10 stars. Here we take into account that clumps survive ejection in only about half the realizations, but this decrease is mostly offset by models showing multiple ejection events. The net ejection rate is then about one clump per model. This estimate is consistent with the estimated empirical ratio of 5–10 between stars and BDs (Luhman et al. 2007). We must keep in mind that our estimate is rather tentative and we expect it to increase somewhat if more models are run at higher resolution. There are also hints in our simulations that ejected clumps with substantial rotation may form close binary objects, in agreement with the prevalence of close binaries in BD–BD systems. All in all, our hybrid scenario agrees, at least qualitatively, with nearly all of the key observational features (see Luhman et al. 2007) of BDs. The only observation not easily accounted for are the wide BD binaries, but these are very rare and they may indeed form from direct collapse since the total mass required is likely not substellar.

Our results reveal a fundamental process that is intrinsic to gaseous disks with moderate to high levels of angular momentum. This clump ejection mechanism should be robust to environmental variations of metallicity and temperature, as it depends more directly on angular momentum and multi-body dynamics. Therefore, we believe that low-mass stars and substellar objects should form in a variety of situations in different cosmic locales and epochs.

The authors thank the anonymous referee, as well as Pavel Kroupa and Alexander DeSouza, for helpful comments on the manuscript. E.I.V. gratefully acknowledges support from the RFBR grants 10-02-00278 and 11-02-92601 and also from a Lise Meitner Fellowship (FWF, Austria). S.B. was supported by a grant from NSERC. We thank the Atlantic Computational Excellence Network (ACEnet) and the SHARCNET consortium for access to computational facilities.