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 Reference Kroupa2001; Chabrier Reference Chabrier2003, Reference Chabrier, Corbelli, Palla and Zinnecker2005) and (ii) to explain the properties of stellar multiple systems (e.g. Raghavan et al. Reference Raghavan2010; Janson et al. Reference Janson2012). 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 Reference Padoan and Nordlund2002; Hennebelle & Chabrier Reference Hennebelle and Chabrier2008, Reference Hennebelle and Chabrier2009) 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 Reference Andre, Ward-Thompson and Barsony1993, Reference Andre, Ward-Thompson and Barsony2000). This has been demonstrated in numerical simulations by Bate (Reference Bate1998, Reference Bate2000), Horton, Bate, & Bonnell (Reference Horton, Bate and Bonnell2001), Matsumoto & Hanawa (Reference Matsumoto and Hanawa2003), Goodwin & Whitworth (Reference Goodwin and Whitworth2004), Delgado-Donate, Clarke, & Bate (Reference Delgado-Donate, Clarke and Bate2004a), Delgado-Donate et al. (Reference Delgado-Donate, Clarke, Bate and Hodgkin2004b), Goodwin, Whitworth, & Ward-Thompson (Reference Goodwin, Whitworth and Ward-Thompson2004, Reference Goodwin, Whitworth and Ward-Thompson2006), Walch et al. (Reference Walch, Burkert, Whitworth, Naab and Gritschneder2009, Reference Walch, Naab, Whitworth, Burkert and Gritschneder2010), Walch, Whitworth, & Girichidis (Reference Walch, Whitworth and Girichidis2012), Lomax et al. (Reference Lomax, Whitworth, Hubber, Stamatellos and Walch2014, Reference Lomax, Whitworth, Hubber, Stamatellos and Walch2015b), and Lomax, Whitworth, & Hubber (Reference Lomax, Whitworth and Hubber2015a).
Observations of prestellar cores (e.g. Motte, Andre, & Neri Reference Motte, Andre and Neri1998; Testi & Sargent Reference Testi and Sargent1998; Johnstone et al. Reference Johnstone, Wilson, Moriarty-Schieven, Joncas, Smith, Gregersen and Fich2000; Motte et al. Reference Motte, André, Ward-Thompson and Bontemps2001; Johnstone et al. Reference Johnstone, Fich, Mitchell and Moriarty-Schieven2001; Stanke et al. Reference Stanke, Smith, Gredel and Khanzadyan2006; Enoch et al. Reference Enoch2006; Johnstone & Bally Reference Johnstone and Bally2006; Nutter & Ward-Thompson Reference Nutter and Ward-Thompson2007; Alves, Lombardi, & Lada Reference Alves, Lombardi and Lada2007; Enoch et al. Reference Enoch, Evans, Sargent, Glenn, Rosolowsky and Myers2008; Simpson, Nutter, & Ward-Thompson Reference Simpson, Nutter and Ward-Thompson2008; Rathborne et al. Reference Rathborne, Lada, Muench, Alves, Kainulainen and Lombardi2009; Könyves et al., Reference Könyves2010; Pattle et al. Reference Pattle2015) 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. (Reference Holman, Walch, Goodwin and Whitworth2013) 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 N2H+ 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 Reference Kenyon and Hartmann1995). 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,
$$\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 Reference Toomre1964). Second, the cooling time of a fragment must be shorter than the orbital period if it is to avoid being sheared part (Gammie Reference Gammie2001).
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 Reference Lodato and Rice2005).
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. Reference Tsukamoto, Takahashi, Machida and Inutsuka2015). Forgan & Rice (Reference Forgan and Rice2013) 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. Reference Lomax, Whitworth, Hubber, Stamatellos and Walch2014, Reference Lomax, Whitworth and Hubber2015a, Reference Lomax, Whitworth, Hubber, Stamatellos and Walch2015b, 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. Reference Hubber, Batty, McLeod and Whitworth2011). Gravitational forces are computed using a tree and artificial viscosity is controlled by the Morris & Monaghan (Reference Morris and Monaghan1997) 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−3 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 Reference Hubber, Walch and Whitworth2013). 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. (Reference Stamatellos, Whitworth, Bisbas and Goodwin2007b). 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,
$$\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. Reference Offner, Klein, McKee and Krumholz2009), 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 Reference Palla and Stahler1993).
We adopt the phenomenological model of episodic accretion, presented by Stamatellos, Whitworth, & Hubber (Reference Stamatellos, Whitworth and Hubber2011, Reference Stamatellos, Whitworth and Hubber2012). This is based on calculations of the magneto-rotational instability (MRI; Zhu, Hartmann, & Gammie Reference Zhu, Hartmann and Gammie2009, Reference Zhu, Hartmann and Gammie2010a; Zhu et al. Reference Zhu, Hartmann, Gammie, Book, Simon and Engelhard2010b). 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
$$\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. Reference Kenyon, Hartmann, Strom and Strom1990). Furthermore, FU Ori-type stars (e.g. Herbig Reference Herbig1977; Greene, Aspin, & Reipurth Reference Greene, Aspin and Reipurth2008; Peneva et al. Reference Peneva, Semkov, Munari and Birkle2010; Green et al. Reference Green2011; Principe et al. Reference Principe, Kastner, Grosso, Hamaguchi, Richmond, Teets and Weintraub2013) can exhibit large increases in luminosity which last ≲ 102 yr. Statistical arguments by Scholz, Froebrich, & Wood (Reference Scholz, Froebrich and Wood2013) suggest that the downtime between outbursts should be of order 104 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}))$
:
$$\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
$$\begin{equation}
\bf{ \mu }\equiv \left(\begin{array}{c}
\mu _{M} \\
\mu _{R} \\
\mu _{{\sigma _\rm {NT}}}
\end{array}\right),
\end{equation}$$
and
$$\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. Reference Motte, Andre and Neri1998; André et al. Reference André, Belloche, Motte and Peretto2007) 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.
Arithmetic means, standard deviations, and correlation coefficients of log (M), log (R), and
$\log (\sigma _\rm {NT })$
for cores in Ophiuchus.
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. (Reference Motte, Andre and Neri1998) and André et al. (Reference André, Belloche, Motte and Peretto2007). See Lomax et al. (Reference Lomax, Whitworth, Hubber, Stamatellos and Walch2014) 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:
$$\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 Reference Lomax, Whitworth and Cartwright2013). The axes are normalised to a given R, giving the dimensions of the core,
$$\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 Reference Bonnor1956; Alves, Lada, & Lada Reference Alves, Lada and Lada2001; Harvey et al. Reference Harvey, Wilner, Lada, Myers, Alves and Chen2001; Lada et al. Reference Lada, Muench, Rathborne, Alves and Lombardi2008). 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
$$\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}$$
$$\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 −4 in three dimensions. We include bulk rotation and radial excursion by modifying the amplitudes a and phases φ of the k = 1 modes:
$$\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}$$
$$\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 rx, ry, rz
, ω
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 1283 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 × 105 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. Reference André, Belloche, Motte and Peretto2007). 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}$
. Footnote
1
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 (Reference Chabrier, Corbelli, Palla and Zinnecker2005) and Kroupa (Reference Kroupa2001) fits to the observed IMF. The vertical-dashed line shows the hydrogen burning limit at M = 0.08 M⊙. See Lomax et al. (Reference Lomax, Whitworth, Hubber, Stamatellos and Walch2015b) 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 (Reference Chabrier, Corbelli, Palla and Zinnecker2005) IMF.
In star-forming regions, the observed ratio of stars to brown dwarfs is
$$\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. Reference Andersen, Meyer, Greissl and Aversa2008). This figure is best reproduced with ERF:
$\mathcal {A}_\rm {ERF}=3.9\pm 0.6$
.Footnote
2
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 Reference Stamatellos, Hubber and Whitworth2007a; Stamatellos & Whitworth Reference Stamatellos and Whitworth2009).
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 (Reference Reipurth and Zinnecker1993), we use the multiplicity frequency and the pairing factor. The multiplicity frequency is fraction of systems which is multiple
$$\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
$$\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
$$\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 primariesFootnote
3
. For comparison, we have also included the same figure for stars in the field (See Duchêne & Kraus Reference Duchêne and Kraus2013, 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 Reference Leinert, Zinnecker, Weitzel, Christou, Ridgway, Jameson, Haas and Lenzen1993; Ratzka, Köhler, & Leinert Reference Ratzka, Köhler and Leinert2005; Kraus et al. Reference Kraus, Ireland, Martinache and Hillenbrand2011). 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 Reference Kroupa1995; Parker et al. Reference Parker, Goodwin, Kroupa and Kouwenhoven2009; Parker & Goodwin Reference Parker and Goodwin2011, Reference Parker and Goodwin2012). 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.
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. Reference Kraus, Ireland, Martinache and Hillenbrand2011) and the field Eggleton & Tokovinin (Reference Eggleton and Tokovinin2008); Tokovinin (Reference Tokovinin2008).
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−2. 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 × 104 yr (b) t = 2.7 × 104 yr (c) t = 3.2 × 104 yr (d) t = 3.7 × 104 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. Reference King, Parker, Patience and Goodwin2012a, Reference King, Goodwin, Parker and Patience2012b). 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. Reference Kouwenhoven, Goodwin, Parker, Davies, Malmberg and Kroupa2010).
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. Reference Motte, Andre and Neri1998; André et al. Reference André, Belloche, Motte and Peretto2007). 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:
$$\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
$$\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}$
.
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−2. Sink particles are represented by black dots. See Lomax et al. (Reference Lomax, Whitworth and Hubber2015a) 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}$
.
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. (Reference Lomax, Whitworth and Hubber2015a) 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 (Reference Chabrier, Corbelli, Palla and Zinnecker2005) IMF.
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 (Reference Chabrier, Corbelli, Palla and Zinnecker2005) IMF. See Lomax et al. (Reference Lomax, Whitworth and Hubber2015a) 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. Reference Stamatellos, Hubber and Whitworth2007a; Stamatellos & Whitworth Reference Stamatellos and Whitworth2009) 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.
ACKNOWLEDGEMENTS
OL and APW gratefully acknowledge the support of a consolidated grant (ST/K00926/1) from the UK STFC. This work was performed using the computational facilities of the Advanced Research Computing @ Cardiff (ARCCA) Division, Cardiff University. All false-colour images have been rendered with splash (Price Reference Price2007).
