The Role of Discs in the Collapse and Fragmentation of Prestellar Cores

Abstract

Disc fragmentation provides an important mechanism for producing low-mass stars in prestellar cores. Here, we describe smoothed particle hydrodynamics simulations which show how populations of prestellar cores evolve into stars. We find the observed masses and multiplicities of stars can be recovered under certain conditions.

First, protostellar feedback from a star must be episodic. The continuous accretion of disc material on to a central protostar results in local temperatures which are too high for disc fragmentation. If, however, the accretion occurs in intense outbursts, separated by a downtime of ~ 10⁴ yr, gravitational instabilities can develop and the disc can fragment.

Second, a significant amount of the cores’ internal kinetic energy should be in solenoidal turbulent modes. Cores with less than a third of their kinetic energy in solenoidal modes have insufficient angular momentum to form fragmenting discs. In the absence of discs, cores can fragment but results in a top-heavy distribution of masses with very few low-mass objects.

1 INTRODUCTION

Two of the main goals of star formation theory are (i) to understand the origin of the stellar initial mass function (IMF; e.g. Kroupa 2001; Chabrier 2003, 2005) and (ii) to explain the properties of stellar multiple systems (e.g. Raghavan et al. 2010; Janson et al. 2012). One possible solution to this problem is the turbulent fragmentation of giant molecular clouds. Here, turbulent flows within molecular clouds produce dense cores of gas (e.g. Padoan & Nordlund 2002; Hennebelle & Chabrier 2008, 2009) of order 0.01 to 0.1 pc across. These may be Jeans unstable, in which case they collapse to form stars (e.g. Andre, Ward-Thompson, & Barsony 1993, 2000). This has been demonstrated in numerical simulations by Bate (1998, 2000), Horton, Bate, & Bonnell (2001), Matsumoto & Hanawa (2003), Goodwin & Whitworth (2004), Delgado-Donate, Clarke, & Bate (2004a), Delgado-Donate et al. (2004b), Goodwin, Whitworth, & Ward-Thompson (2004, 2006), Walch et al. (2009, 2010), Walch, Whitworth, & Girichidis (2012), Lomax et al. (2014, 2015b), and Lomax, Whitworth, & Hubber (2015a).

Observations of prestellar cores (e.g. Motte, Andre, & Neri 1998; Testi & Sargent 1998; Johnstone et al. 2000; Motte et al. 2001; Johnstone et al. 2001; Stanke et al. 2006; Enoch et al. 2006; Johnstone & Bally 2006; Nutter & Ward-Thompson 2007; Alves, Lombardi, & Lada 2007; Enoch et al. 2008; Simpson, Nutter, & Ward-Thompson 2008; Rathborne et al. 2009; Könyves et al., 2010; Pattle et al. 2015) show that the core mass function (CMF) is very similar in shape to the IMF (i.e. a lognormal distribution with a power-law tail at high mass), albeit shifted upwards in mass by a factor of 3 to 5. This had led to the suggestion that there is a self-similar mapping of the CMF onto the IMF. Statistical analysis by Holman et al. (2013) suggests that a core should spawn on average four to five stars in order to explain the observed abundance of multiple systems.

A core may fragment into multiple objects via either turbulent fragmentation—similar to the molecular cloud—or disc fragmentation. Observed N₂H⁺ line widths in cores indicate that the internal velocity dispersion in most cores is sub to trans-sonic. This suggests that a typical core is unlikely to collapse and fragment into more than one or two objects through turbulence alone. However, the first protostars to form are usually attended by accretion discs (Kenyon & Hartmann 1995). These discs may fragment if two criteria are fulfilled. First, the disc must have a sufficiently large surface density Σ(R) so that fragments can overcome thermal and centrifugal support,

(1) $$\begin{eqnarray} \Sigma (R)\gtrsim \frac{c(R)\,\kappa (R)}{\pi \,G}, \end{eqnarray}$$

where κ(R) is the epicyclic frequency and c(R) is the sound speed (Toomre 1964). Second, the cooling time of a fragment must be shorter than the orbital period if it is to avoid being sheared part (Gammie 2001).

The above criteria apply well to low-mass discs. However, disc dynamics are more complex when the disc and the central protostar have similar masses. In cases where the disc is marginally unstable (i.e. $t_\rm {COOL}\sim t_\rm {ORBIT}$ and π G Σ ~ c κ), instabilities develop which bolster the accretion rate onto the central protostar. This lowers the mass (and hence surface density) of the disc, restabilising it. This suggests that an otherwise unstable disc may be able to remain stable by undergoing episodic accretion events (Lodato & Rice 2005).

When non-local effects such as radiative transfer dominate the disc temperature, local cooling timescales become very difficult—if not impossible—to calculate. Here, again, a seemingly unstable disc may be able to resist fragmentation (e.g. Tsukamoto et al. 2015). Forgan & Rice (2013) also show that irradiation (from a central protostar and/or interstellar radiation field) increases the jeans mass within disc spiral structures. They deduce that disc fragmentation is far more likely to result in low-mass stars and brown dwarfs than gas giant planets.

In this paper, we analyse simulations which follow the evolution of an ensemble of synthetic cores, based on the properties of those in the Ophiuchus star-forming region (See Lomax et al. 2014, 2015a, 2015b, for more details). We examine how protostellar feedback affects the fragmentation of discs, and how the results compare with observed stars. We also examine how the ratio of solenoidal to compressive modes in the turbulent velocity field affects the formation of discs and filaments within a prestellar core. In Section 2, we describe the numerical method used to evolve the cores. In Section 3, we describe (i) how to generate realistic core initial conditions and (ii) how they evolve under different prescriptions of radiative feedback. In Section 4, we show how changing the ratio of solenoidal to compressive modes in the velocity field affects core fragmentation. We summarise and conclude in Section 5.

2 NUMERICAL METHOD

2.1. Smoothed particle hydrodynamics

Core evolution is simulated using the seren ∇h-SPH code (Hubber et al. 2011). Gravitational forces are computed using a tree and artificial viscosity is controlled by the Morris & Monaghan (1997) prescription. In all simulations, the smoothed particle hydrodynamics (SPH) particles have mass $m_\rm {SPH}=10^{-5}\,{\rm M}_{\odot }$ , so that the opacity limit (~ 3 × 10⁻³ M_⊙) is resolved with ~ 300 particles. Gravitationally bound regions with density higher than $\rho _\rm {SINK}=10^{-9}\,{\rm g}\,{\rm cm}^{-3}$ are replaced with sink particles (Hubber, Walch, & Whitworth 2013). Sink particles have radius $r_\rm {SINK}\simeq 0.2\,{\rm au}$ , corresponding to the smoothing length of an SPH particle with density equal to $\rho _\rm {SINK}$ . The equation of state and the energy equation are treated with the algorithm described in Stamatellos et al. (2007b). Magenetic fields and mechanical feedback (e.g. stellar winds) is not included in these simulations.

2.2. Accretion feedback

Radiative feedback from the protostars (i.e. sink particles formed in the simulations) is included. The dominant contribution to the luminosity of a protostar is usually from accretion,

(2) $$\begin{eqnarray} L_\star \simeq \frac{f\,G\,M_\star \,\dot{M}_\star }{R_\star }, \end{eqnarray}$$

where, f = 0.75 is the fraction of the accreted material’s gravitational energy that is radiated from the surface of the protostar (the rest is presumed to be removed by bipolar jets and outflows; Offner et al. 2009), M _⋆ is the mass of the protostar, $\dot{M}_\star$ is the rate of accretion onto the protostar, and R _⋆ = 3 R_⊙ is the approximate radius of a protostar (Palla & Stahler 1993).

We adopt the phenomenological model of episodic accretion, presented by Stamatellos, Whitworth, & Hubber (2011, 2012). This is based on calculations of the magneto-rotational instability (MRI; Zhu, Hartmann, & Gammie 2009, 2010a; Zhu et al. 2010b). In the outer disc of a protostar (outside the sink radius, and therefore resolved by the simulation), angular momentum is redistributed by gravitational torques and material spirals inwards towards the sink. At distances within the sink radius (unresolved by the simulation), the inner disc is so hot that it is gravitationally stable, and unable to fragment. Material continues to accrete onto this inner disc (i.e. the sink particle) until the gas is hot enough to thermally ionise and couple with the local magnetic field. The MRI cuts in and magnetic torques allow material in the inner disc to rapidly accrete onto the protostar. This results in extended periods of very low accretion luminosity, punctuated by intense, episodic outbursts. The length of the downtime between outbursts—during which disc fragmentation may occur—is given by

(3) $$\begin{eqnarray} && t_\rm {DOWN}\sim 1.3\times 10^4\text{\,yr}\nonumber\\ &&\quad\times \left(\frac{M_\star }{0.2\,{\rm M}_{\odot }}\right)^{2/3}\left(\frac{\dot{M}_\rm {SINK}}{10^{-5}\,{\rm M}_{\odot }{\rm yr}^{-1}}\right)^{-8/9}, \end{eqnarray}$$

where $\dot{M}_\rm {SINK}$ is the rate at which material flows into the sink.

There is also observational motivation for adopting an episodic model. The luminosities of young stars are about an order of magnitude lower than expected from continuous accretion (this is the luminosity problem, first noted by Kenyon et al. 1990). Furthermore, FU Ori-type stars (e.g. Herbig 1977; Greene, Aspin, & Reipurth 2008; Peneva et al. 2010; Green et al. 2011; Principe et al. 2013) can exhibit large increases in luminosity which last ≲ 10² yr. Statistical arguments by Scholz, Froebrich, & Wood (2013) suggest that the downtime between outbursts should be of order 10⁴ yr, similar to the timescale given in Equation (3).

3 PRESTELLAR CORES IN OPHIUCHUS

3.1. Initial conditions

Using the observed properties of cores as a basis for numerical simulations presents a difficult inverse problem. The mass, temperature, projected area, and projected aspect ratio of a core can be reasonably inferred from bolometric measurements of the dust in prestellar cores. In addition, the line-of-sight velocity dispersion within the core can be inferred from the width of molecular lines. However, the initial boundary conditions of a core simulation must represent the full spatial and velocity structure of the system. This occupies six dimensions, whereas observational data can only provide information on three (i.e. two spatial and one velocity).

Rather than trying to emulate a specific core—which is arguably an impossibility—we can relatively easily define a distribution of cores which have the same, or at least very similar, statistical properties to those in a given region. We based the synthetic core initial conditions on Ophiuchus. This is a well-studied region, for which many of the aforementioned core properties have been measured.

3.1.1. Mass, size, and velocity dispersion

Only some of the measured core masses in Ophiuchus have both an associated size and velocity dispersion. In order to make the most of the data, we define the following lognormal probability distribution of ${\bf{x}}\equiv (\log (M),\log (R),\log (\sigma _\rm {NT}))$ :

(4) $$\begin{equation} P(\bf{ x})=\frac{1}{(2\pi )^{3/2}|\bf{ \Sigma }|}\exp \left(-\frac{1}{2}(\bf{ x}-\bf{ \mu })^{\rm T}\bf{ \Sigma }^{-1}(\bf{ x}-\bf{ \mu }) \right), \end{equation}$$

where

(5) $$\begin{equation} \bf{ \mu }\equiv \left(\begin{array}{c} \mu _{M} \\ \mu _{R} \\ \mu _{{\sigma _\rm {NT}}} \end{array}\right), \end{equation}$$

and

(6) $$\begin{equation} \bf{ \Sigma }\equiv \left( \begin{array}{c@{\hspace*{8pt}}c@{\hspace*{8pt}}c} \sigma _{M}^2 & \rho _{{M,R}}\,\sigma _{M}\sigma _{R} & \rho _{{M,\sigma _\rm {NT}}}\,\sigma _{M}\sigma _{{\sigma _\rm {NT}}} \\ \rho _{{M,R}}\,\sigma _{M}\sigma _{R} & \sigma _{R}^2 & \rho _{{R,\sigma _\rm {NT}}}\,\sigma _{R}\sigma _{{\sigma _\rm {NT}}} \\ \rho _{{M,\sigma _\rm {NT}}}\,\sigma _{M}\sigma _{{\sigma _\rm {NT}}} & \rho _{{R,\sigma _\rm {NT}}}\,\sigma _{R}\sigma _{{\sigma _\rm {NT}}} & \sigma _{{\sigma _\rm {NT}}}^2 \end{array}\right). \end{equation}$$

The coefficients of μ and Σ are calculated from the observed Ophiuchus data (Motte et al. 1998; André et al. 2007) and are given in Table 1. From P( x ), we are able to draw any number of masses, sizes, and velocity dispersions, all of which were statistically similar to those in Ophiuchus. The distribution of P(x) is shown in Figure 1.

Table 1. Arithmetic means, standard deviations, and correlation coefficients of log (M), log (R), and $\log (\sigma _\rm {NT })$ for cores in Ophiuchus.

Figure 1. The multivariate lognormal distribution, P( x ) where $\bf{ x}=(\log (M),\log (R),\log (\sigma _\rm {NT}))$ . The figure shows the projections through $\log (\sigma _\rm {NT})$ , through log (M), and through log (R). The concentric ellipses show the 1σ, 2σ, and 3σ regions of the distribution. The green circles are randomly drawn points from P( x ). The red squares are the observational data from Motte et al. (1998) and André et al. (2007). See Lomax et al. (2014) for the original version of this figure.

3.1.2. Shape

Molecular cloud cores often have elongated, irregular shapes. We include this in the simulations by assuming that each intrinsic core shape can be drawn from a family of triaxial ellipsoids. Each ellipsoid had axes:

(7) $$\begin{eqnarray} A&=&1,\nonumber\\ B&=&\exp (\tau \mathcal {G}_\rm {B}),\\ C&=&\exp (\tau \mathcal {G}_\rm {C}), \nonumber \end{eqnarray}$$

where $\mathcal {G}_\rm {B}$ and $\mathcal {G}_\rm {C}$ are random numbers drawn from a Gaussian distribution with zero mean and unit standard deviation. The scale-parameter τ ≈ 0.6 is a fit to the distribution of projected aspect ratios in Ophiuchus (Lomax, Whitworth, & Cartwright 2013). The axes are normalised to a given R, giving the dimensions of the core,

(8) $$\begin{eqnarray} A_\rm {CORE}&=&\frac{R}{(BC)^{{1}/{3}}},\nonumber\\ B_\rm {CORE}&=&BA_\rm {CORE},\\ C_\rm {CORE}&=&CA_\rm {CORE}. \nonumber \end{eqnarray}$$

3.1.3. Density profile

Density profiles of cores are often well fitted by those of critical Bonnor–Ebert spheres (e.g. Bonnor 1956; Alves, Lada, & Lada 2001; Harvey et al. 2001; Lada et al. 2008). We use such a profile for the ellipsoidal cores. Here, ρ = ρ_Ce^−ψ(ξ), where ρ_C is the central density, ψ is the Isothermal Function, ξ is the dimensionless radius, ξ_B = 6.451 is the boundary of the sphere. The density at any given point (x, y, z), where (0, 0, 0) is the centre of the core, is given by

(9) $$\begin{eqnarray} \xi &=&\xi _{{\rm FWHM}}\left(\frac{x^2}{A_\rm {CORE}^2}+\frac{y^2}{B_\rm {CORE}^2}+\frac{z^2}{C_\rm {CORE}^2}\right)^{1/2}, \end{eqnarray}$$

(10) $$\begin{eqnarray} \rho (\xi )&=&\frac{M_\rm {CORE}\xi _{{\rm B}}{\rm e}^{-\psi (\xi )}}{4\pi A_\rm {CORE}B_\rm {CORE}C_\rm {CORE}\psi ^{\prime }(\xi _{{\rm B}})},\hspace{14.22636pt}\xi <\xi _{{\rm B}}, \end{eqnarray}$$

where ψ′ is the first derivative of ψ and ξ_FWHM = 2.424 is the full width at half maximum of the column density through a critical Bonnor–Ebert sphere.

3.1.4. Velocity field

Each core is given a turbulent velocity field with power spectrum P∝k ⁻⁴ in three dimensions. We include bulk rotation and radial excursion by modifying the amplitudes a and phases φ of the k = 1 modes:

(11) $$\begin{eqnarray} \left[\begin{array}{c}\bf{ a}(1,0,0) \\ \bf{ a}(0,1,0) \\ \bf{ a}(0,0,1) \end{array}\right]=\left[\begin{array}{c@{\hspace*{8pt}}c@{\hspace*{8pt}}c}r_x & \omega _z & -\omega _y \\ -\omega _z & r_y & \omega _x \\ \omega _y & -\omega _x & r_z \end{array}\right], \end{eqnarray}$$

(12) $$\begin{eqnarray} \left[ \begin{array}{c} \bf{ \varphi }(1,0,0) \\ \bf{ \varphi }(0,1,0) \\ \bf{ \varphi }(0,0,1) \end{array}\right]= \left[\begin{array}{c@{\hspace*{8pt}}c@{\hspace*{8pt}}c}\pi /2 & \pi /2 & \pi /2 \\ \pi /2 & \pi /2 & \pi /2 \\ \pi /2 & \pi /2 & \pi /2 \end{array}\right]. \end{eqnarray}$$

The amplitude components r_x, r_y, r_z , ω _x , ω _y , and ω _z are drawn independently from a Gaussian distribution with zero mean and unit variance. The r terms define the amount of excursion along a given axis. The ω terms define the amount rotation about a given axis. The fields are generated on a 128³ grid and interpolated onto the SPH particles. The velocity dispersion of the particles is normalised to a given value of $\sigma _\rm {NT}$ .

3.2. Results

One hundred synthetic cores have been evolved for 2 × 10⁵ yr. This is roughly an order of magnitude greater than the average core free-fall time and of the same order as the estimated core–core collisional timescale (André et al. 2007). Of the one hundred cores, sixty are prestellar. Each simulation has been performed three times: once with no radiative feedback from accretion (NRF); again with episodic feedback from accretion (ERF); and finally with continuous feedback from accretion (CRF). With CRF, the protostellar luminosity was calculated using Equation (2), with $\dot{M}_\star$ set to $\dot{M}_\rm {SINK}$ .

3.2.1. Stellar masses

We find two very noticeable trends when progressing along the series NRF → ERF → CRF. First, the median number stars formed per core decreases: $\mathcal {N}_\rm {S/C}^\rm {NRF}=6.0^{+4.0}_{-2.0}$ , $\mathcal {N}_\rm {S/C}^\rm {ERF}=3.5^{+3.5}_{-2.5}$ , $\mathcal {N}_\rm {S/C}^\rm {CRF}=1.0^{+0.0}_{-0.0}$ . Second, the median protostellar mass (see Figure 2) shifts upwards: $\log _{10} M_\star ^\rm {NRF}=-1.1^{+0.4}_{-0.6}$ , $\log _{10} M_\star ^\rm {ERF}=-0.8^{+0.2}_{-0.4}$ , $\log _{10} M_\star ^\rm {CRF}=-0.5^{+0.2}_{-0.2}$ . ¹

Figure 2. The black histograms show stellar mass functions for (a) NRF, (b) ERF, and (c) CRF. The blue-dotted straight lines, and the red-dashed lognormal curve, show, respectively, the Chabrier (2005) and Kroupa (2001) fits to the observed IMF. The vertical-dashed line shows the hydrogen burning limit at M = 0.08 M_⊙. See Lomax et al. (2015b) for the original version of this figure.

These trends occur because disc fragmentation is strongly affected by protostellar feedback. In the case with NRF, discs are relatively cold and easily fragment. Recall from Equation (1) that fragmentation can occur if the disc has a sufficiently large column density or the sound speed is sufficiently low. With ERF, fragmentation is occasionally interrupted by the episodic outbursts. With CRF, the discs are constantly heated, and fragmentation becomes difficult. As a consequence, the central protostar or binary usually accretes the entire mass of its disc. Protostars with NRF, and to a slightly lesser extent ERF, are often attended by multiple low-mass companions which partially starve the primary protostar of the remaining gas. Of the three sets of simulations, those with ERF best reproduce the Chabrier (2005) IMF.

In star-forming regions, the observed ratio of stars to brown dwarfs is

(13) $$\begin{equation} \mathcal {A}=\frac{N(0.08\,{\rm M}_{\odot }<M\le 1.0\,{\rm M}_{\odot })}{N(0.03\,{\rm M}_{\odot }<M\le 0.08\,{\rm M}_{\odot })}=4.3\pm 1.6, \end{equation}$$

(See Andersen et al. 2008). This figure is best reproduced with ERF: $\mathcal {A}_\rm {ERF}=3.9\pm 0.6$ . ² Simulations with NRF and CRF yield $\mathcal {A}_\rm {NRF}=2.2\pm 0.3$ and $\mathcal {A}_\rm {CRF}=17\pm 8$ . While the figure for NRF is not completely incompatible with observations, the value with CRF is far too high. This is because brown dwarfs are unable to form via disc fragmentation (e.g. Stamatellos, Hubber, & Whitworth 2007a; Stamatellos & Whitworth 2009).

3.2.2. Stellar multiplicities

The core simulations also produce a wide variety of multiple systems. These are either simple binary systems or high order (N ⩾ 3) hierarchical multiples. A high-order system can be thought of as a binary where one or both components is another binary system. Multiple systems of protostars are extracted from the end state simulations if they have been tidally stable for at least one orbital period.

There are many ways to statistically describe the multiplicity of a population of stellar systems. Here, following Reipurth & Zinnecker (1993), we use the multiplicity frequency and the pairing factor. The multiplicity frequency is fraction of systems which is multiple

(14) $$\begin{equation} {\rm mf}=\frac{B+T+Q\dots }{S+B+T+Q\dots }, \end{equation}$$

where S is the number of single systems, B is the number of binaries, T is the number of triples, etc. The pairing factor is the average number of orbits per system

(15) $$\begin{equation} {\rm pf}=\frac{B+2T+3Q\dots }{S+B+T+Q\dots }. \end{equation}$$

These two quantities are particularly useful in conjunction, as their ratio gives the average number of objects per multiple system

(16) $$\begin{equation} \mathcal {O}_\rm {SYS}=1+\frac{{\rm pf}}{{\rm mf}}. \end{equation}$$

By definition, $\mathcal {O}_\rm {SYS}\ge 2$ .

Figure 3 shows mf, pf, and $\mathcal {O}_\rm {SYS}$ for systems with very low mass or brown dwarf primaries, M-dwarf primaries, and solar-type primaries ³ . For comparison, we have also included the same figure for stars in the field (See Duchêne & Kraus 2013, and references therein). Observations of pre-main-sequence stars in Ophiuchus and Taurus suggest that multiplicity is high when stars are young (e.g. Leinert, Zinnecker, Weitzel, Christou, Ridgway, Jameson, Haas 1993; Ratzka, Köhler, & Leinert 2005; Kraus et al. 2011). As stars age, particularly if they are in a clustered environment, dynamical interactions erode these systems and the multiplicity frequency drops to that observed in the field (e.g. Kroupa 1995; Parker et al. 2009; Parker & Goodwin 2011, 2012). The multiplicity of protostars (i.e. the objects formed in these simulations) should therefore be higher than the multiplicity of field stars. This requirement is met for the simulations with NRF and ERF. However, with CRF, the multiplicity frequency is too low for solar-type primaries.

Figure 3. Multiplicity frequency (a), pairing factor (b), and mean system order (c) for systems with very low mass, M-dwarf and solar-type primaries. The red boxes give the values for the NRF simulations, blue for the ERF, and green for the CRF. The black points give the values observed in field main-sequence stars. In all cases, the width of a box shows the extent of the mass bin, and the height shows the uncertainty.

Very high order systems, e.g. sextuples, are also form in these simulations. Figure 4 shows a sextuple system formed from a solar mass core with ERF. Initially the core fragments into two objects, which form a binary system with a circumbinary disc. This disc fragments into a further five objects. It finally settles into a sextuple system (a binary system orbiting a quadruple) with separations ranging from ~ 600 AU for the outer orbit to ~ 0.1 AU for the innermost orbits. Sextuple systems similar to this are found in star-forming regions (Kraus et al. 2011) and the field Eggleton & Tokovinin (2008); Tokovinin (2008).

Figure 4. A sequence of column density maps of a core during disc fragmentation. The initial core has M = 1.3 M_⊙, R = 3000 AU, $\sigma _\rm {NT}=0.3\,{\rm km\,s^{-1}}$ and is evolved with ERF. The colour-bar gives shows the column density in units of g cm⁻². The black dots show the positions of sink particles, i.e. protostars. Figure 4(a) shows gravitational instabilities developing in a circumbinary disc. Figure 4(b) shows seven protostars in an unstable configuration. In Figure 4(c), the protostars are configured in a quadruple system (right) and a binary (left); a single protostar is being ejected (centre). Figure 4(d) shows a stable sextuple system which lasts until the end of the simulation. (a) t = 2.2 × 10⁴ yr (b) t = 2.7 × 10⁴ yr (c) t = 3.2 × 10⁴ yr (d) t = 3.7 × 10⁴ yr

All three sets of simulations produce systems with semi-major axes ranging between roughly 0.1 AU (the resolution limit of the simulation) and 1000 AU. This upper limit corresponds to that observed in star-forming regions (e.g. King et al. 2012a, 2012b). Systems with wider orbits ( $\gtrsim 10^3\,{\rm \text{AU}}$ ) are probably assembled later through dynamical interactions in clustered environments (e.g. Kouwenhoven et al. 2010).

3.3. Summary

We find that the mass distribution and multiplicity statistics of young stars can be recovered from simulations if radiative feedback from protostellar accretion is episodic. The periods of low luminosity provide a window of opportunity, during which protostellar discs can fragment. When feedback is continuous, the disc is too warm to permit disc fragmentation. In this instance, the simulations result in (i) a top-heavy IMF, and (ii) too few stars per core to satisfactorily reproduce observed multiplicity statistics. When there is no radiative feedback, good multiplicity statistics are recovered, but the resultant IMF has too many brown dwarfs.

4 TURBULENCE: DISCS AND FILAMENTS

To examine how the structure of a core’s velocity field affects the star formation process, we take a single core set-up and vary the structure of the velocity field. The core is spherical, with M = 3 M_⊙, R = 3000 AU, and $\sigma _\rm {NT}=0.44\,{\rm km\,s^{-1}}$ . These dimensions are similar to those of SM1 in Oph-A region (Motte et al. 1998; André et al. 2007). We alter the partition of kinetic energy in solenoidal modes (i.e. shear and rotation) and compressive modes (i.e. compression and rarefaction).

4.1. Initial conditions

The amplitude of a turbulent mode $\bf {a}(\bf{ k})$ can be split into its longitudinal (compressive) component $\bf {a_\rm {L}}(\bf{ k})$ and transverse (solenoidal) component $\bf {a_\rm {T}}(\bf{ k})$ using Helmholtz decomposition:

(17) $$\begin{eqnarray} \bf{ a_\rm {L}}(\bf{ k})&=&\hat{\bf{k}}(\bf{ a}(\bf{ k})\cdot \hat{\bf{k}}),\nonumber\\[-6pt] &&\\[-6pt] \bf{ a_\rm {T}}(\bf{ k})&=&\bf{ a}(\bf{ k})-\hat{\bf{k}}(\bf{ a}(\bf{ k})\cdot \hat{\bf{k}}).\nonumber \end{eqnarray}$$

When the magnitude and direction of the amplitude is random, there is on average twice as much energy in transverse modes as there is in solenoidal modes.

We define the parameter $\delta _\rm {SOL}$ as the average fraction of solenoidal kinetic energy in a velocity field. We modify the field to have given $\delta _\rm {SOL}$ by performing the transformation

(18) $$\begin{equation} \bf{ a}(\bf{ k})\rightarrow \sqrt{3\,(1-\delta _\rm {SOL})}\,\bf{ a_\rm {L}}(\bf{ k})+\sqrt{\frac{3}{2}\,\delta _\rm {SOL}}\,\bf{ a_\rm {T}}(\bf{ k}). \end{equation}$$

We generate ten initial cores, each with a unique random velocity seed. For each core, we apply the transformation in Equation (18) with values $\delta _\rm {SOL}=0,\,\frac{1}{9},\,\frac{1}{3},\,\frac{2}{3},\,1$ , yielding a total fifty core set-ups.

4.2. Results

Figure 5 shows a montage of simulation snapshots where the random seed is fixed and $\delta _\rm {SOL}$ is varied from 0 to 1. When the field is purely compressive ( $\delta _\rm {SOL}=0$ ), protostars form within a network of filaments. Due to the low angular momentum of the system, the resultant protostars are only attended by small discs. In contrast, when the field is purely solenoidal ( $\delta _\rm {SOL}=1)$ , a single protostar forms and is attended by a coherent disc structure. The disc proceeds to fragment into multiple objects. A smooth transition between filament fragmentation and disc fragmentation is seen in the snapshots with intermediate values of $\delta _\rm {SOL}$ .

Figure 5. Column density maps of the central 820 AU × 820 AU of the (x, y)-plane, from the simulations with fixed random seed and different values of $\delta _\rm {SOL}$ , at times $t=1.00,1.25,1.50,\text{ and }1.75\times 10^4\,{\rm yr}$ . The colour scale gives the logarithmic column density in units of g cm⁻². Sink particles are represented by black dots. See Lomax et al. (2015a) for the original version of this figure. (a) $\delta _\rm {SOL}=0$ . (b) $\delta _\rm {SOL}=1/9$ . (c) $\delta _\rm {SOL}=1/3$ . (d) $\delta _\rm {SOL}=2/3$ . (e) $\delta _\rm {SOL}=1$ .

Figure 6 shows the fraction of protostars formed by filament fragmentation (i.e. in relative isolation) and disc fragmentation (i.e. in discs around more massive protostars) as a function of $\delta _\rm {SOL}$ . These values are averaged over all random seeds. Here, we see that the occurrence of disc fragmentation increases monotonically with $\delta _\rm {SOL}$ . Filament fragmentation, therefore, decreases monotonically with $\delta _\rm {SOL}$ .

Figure 6. The fraction of stars formed by filament fragmentation (red crosses) and disc fragmentation (green boxes) for different values of $\delta _\rm {SOL}$ , averaged over all random seeds. The error bars show the Poison counting uncertainties. See Lomax et al. (2015a) for the original version of this figure.

On average, the number of protostars spawned per core ranges from 5.4 ± 7 when $\delta _\rm {SOL}=0$ to 8.1 ± 9 when $\delta _\rm {SOL}=1$ . The median and interquartile range of mass is shown in Figure 7. Here, we see that increasing $\delta _\rm {SOL}$ pushes the median sink mass down from roughly 0.6 M_⊙ when $\delta _\rm {SOL}=0$ to 0.3 M_⊙ when $\delta _\rm {SOL}=1$ . We also find that purely compressive fields form very few brown dwarfs. We note that—even for this limited set of initial conditions—the interquartile range of protostellar masses when $\delta _\rm {SOL}\gtrsim 2/3$ is very similar to that of the Chabrier (2005) IMF.

Figure 7. The black points show the median stellar mass, and the vertical black bars show interquartile range of mass, for different values of $\delta _\rm {SOL}$ , averaged over all random seeds. The solid and dashed horizontal red lines show the median and interquartile range for the Chabrier (2005) IMF. See Lomax et al. (2015a) for the original version of this figure.

4.3. Summary

We show that the collapse and fragmentation of prestellar cores is strongly influenced by the structure of the velocity field. Disc formation and fragmentation dominates the star formation process when $\delta _\rm {SOL}\gtrsim 1/3$ . At values below this, stars form mostly through the fragmentation of filamentary structures. The value of $\delta _\rm {SOL}$ also affects distribution of protostellar masses. The distribution most resembles the observed IMF when $\delta _\rm {SOL}\gtrsim 2/3$ . Reducing $\delta _\rm {SOL}$ reduces the level of disc fragmentation, resulting in smaller number of objects with greater average mass. In extreme cases (i.e. $\delta _\rm {SOL}<1/9$ ), the formation of low-mass stars and brown dwarfs is heavily suppressed. This suggests that disc fragmentation may be a requirement for forming these objects.

5 SUMMARY AND CONCLUSIONS

Previous numerical work (e.g. Stamatellos et al. 2007a; Stamatellos & Whitworth 2009) shows that disc fragmentation is important mechanism for reproducing the properties of low-mass stars. We demonstrate that disc fragmentation also plays an important role in the conversion of prestellar cores into stars. Importantly, the observed masses and multiplicities of stars can be recovered from simulations if the following criteria are satisfied:

• • Radiative feedback from accretion onto protostars is episodic. Simulations with episodic radiative feedback produce both an IMF and multiplicity statistics in good agreement with those observed. Simulations with continuous radiative feedback fail to produce both the observed number of brown dwarfs and the multiplicity statistics associated with young objects. Furthermore, continuous radiative feedback produces protostellar luminosities greater than those observed in young stars. Simulations with no radiative feedback can produce good multiplicity statistics and a reasonable—albeit bottom-heavy—IMF, but are unrealistic.

• • A significant proportion of the core’s internal kinetic energy is in solenoidal turbulent modes. Cores with more than a third of their kinetic energy in solenoidal modes are able to easily produce stars via disc fragmentation. Decreasing this fraction results in more stars forming in fragmenting filaments. These objects tend to be of greater mass, resulting in a top-heavy mass distribution relative to the observed IMF. Furthermore, filament fragmentation struggles to produce brown dwarfs and very low mass stars.