Protostellar Feedback in Turbulent Fragmentation: Consequences for Stellar Clustering and Multiplicity
Abstract
Stars are strongly clustered on both large () and small (binary) scales, but there are few analytic or even semianalytic theories for the correlation function and multiplicity of stars. In this paper we present such a theory, based on our recentlydeveloped semianalytic framework called MISFIT, which models gravitoturbulent fragmentation, including the suppression of fragmentation by protostellar radiation feedback. We compare the results including feedback to a control model in which it is omitted. We show that both classes of models robustly reproduce the stellar correlation function at scales, which is well approximated by a powerlaw that follows generally from scalefree physics (turbulence plus gravity) on large scales. On smaller scales protostellar disk fragmentation becomes dominant over common core fragmentation, leading to a steepening of the correlation function. Multiplicity is more sensitive to feedback: we found that a model with the protostellar heating reproduces the observed multiplicity fractions and mass ratio distributions for both Solar and subSolar mass stars (in particular the brown dwarf desert), while a model without feedback fails to do so. The model with feedback also produces an atformation period distribution consistent with the one inferred from observations. However, it is unable to produce shortrange binaries below the length scale of protostellar disks. We suggest that such close binaries are produced primarily by disk fragmentation and further decrease their separation through orbital decay.
keywords:
stars: formation – stars: binaries: general – galaxies: star clusters: general – turbulence – galaxies: star formation – cosmology: theory1 Introduction
Star formation (SF) is complex problem that involves nonlinear physics (turbulence, chemistry, gravity, radiation, etc.) on a vast dynamic range. To achieve a deeper understanding of this process a number of simplified theoretical models have been proposed that try to pinpoint the physical processes responsible for individual qualitative features. The most common test of these models is a comparison to the initial mass function (IMF), but that is just one aspect of star formation. It has been long proposed that the stellar clustering and multiplicity properties carry the imprints of the physical processes of star formation (Kuiper 1935), making them an ideal secondary test for different star formation models.
It is well known that starforming regions are highly structured, with stellar positions correlated on a wide range of scales Lada & Lada (2003); Portegies Zwart et al. (2010); Bressert et al. (2010); Gouliermis et al. (2015). The stellar correlation function has been measured in a wide range (about 5 orders of magnitude in radius) and is found to be rising monotonically on smaller scales in all star clusters (Simon 1997; Nakajima et al. 1998; Hartmann 2002; Hennekemper et al. 2008; Kraus & Hillenbrand 2008). Despite the overwhelming observational data and statistical analysis (Bate et al., 1998; Cartwright & Whitworth, 2004) there has been little effort to formulate a theoretical understanding of why star formation is clustered. A number of authors have measured the clustering of the stars produced in numerical simulations (e.g., Klessen & Burkert 2000; Hansen et al. 2012; see the review by Krumholz 2014 for further references) and found reasonable agreement with observations, but the physical origin of the result was not completely clear. Hopkins (2013a) (henceforth referred to as H13) was the first to provide a quantitative explanation in terms of the statistics of turbulence. Using the excursion set formalism H13 showed that the correlation function of “last crossing objects”^{1}^{1}1Smallest selfgravitating structures in a fully developed turbulent medium at a fixed time. They are considered to be the analogues of protostellar cores. is remarkably similar to that of observed cores, which itself is similar to the correlation function of stars (Stanke et al. 2006). However, this model has a significant limitation in that it is calculated at a fixed time, so the collapse and further fragmentation of cores is not taken into account; it cannot therefore predict the correlation function of stars, nor their multiplicity statistics.
There is similarly an abundance of observational data about the multiplicity properties of stars (e.g. Raghavan et al. 2010 for Solartype stars, Burgasser et al. 2007 for brown dwarfs; see Duchêne & Kraus 2013 for a more detailed review). It is generally understood that most multiple star systems either form during the star formation phase through common core fragmentation and protostellar disk fragmentation (Tohline 2002) or during the cluster dissolution phase (Kouwenhoven et al., 2010; Parker & Meyer, 2014). Most theoretical work is focused on modeling these processes in detailed numerical studies. Hydrodynamical simulations (e.g. Bate 2009a, 2012; Offner et al. 2010; Krumholz et al. 2012) have shown good agreement with observed multiplicity statistics and found that radiation feedback is essential. However, these simulations necessarily have limited dynamical range and statistics, of key importance for highmass stars and long range binaries, and pinpointing the key physics in them is quite challenging.
There has also been significant effort to infer atformation multiplicity properties from observations. Both observations (Duchêne 1999; Kraus et al. 2008, 2011) and simulations have shown that stars are born in complex, multiple systems that are broken up by dynamical effects (e.g., ejection of stars) causing multiplicity to drop (Goodwin et al. 2007; Kaczmarek et al. 2011) and the period distribution (commonly referred to as the binary distribution function) to shift to shorter periods (Kroupa 1995; Marks et al. 2011). This can be understood as the result of long range binaries being preferentially broken up by ejection events, which also increase the binding between leftover stars (“hardening”). This means that to reproduce the present day multiplicity and binary distribution functions the atformation multiplicity should be of order unity for massive stars, and their period distribution should be flat. These findings, however, have recently been called into question. Parker (2014) showed that the densities of star forming regions are constant or increasing with time, while Parker & Meyer (2014) found that an initial distribution of stars with unit multiplicity and an excess of wide binaries will not evolve through Nbody processes into a distribution consistent with that observed in field stars.
The aim of this paper is to expand upon the work of H13 by investigating the features imprinted by isothermal fragmentation and protostellar heating through common core fragmentation in the stellar correlation and multiplicities. This is accomplished by utilizing the MISFIT (Minimalistic Star Formation Including Turbulence) semianalytical framework described by Guszejnov et al. (2016) (hereafter referred to as GKH16), which combines the fragmentation formalism of Guszejnov & Hopkins (2015) and the protostellar heating model of Krumholz (2011) (henceforth referred to as GH15 and K11) to follow the evolution and collapse of a statistical ensemble of giant molecular clouds (GMCs) down to the protostellar size scale.
The remainder of this paper is organized as follows. First, in Sec. 2 we briefly outline the MISFIT framework that we use. In Sec. 3.1 we show that the stellar correlation function is insensitive to both initial conditions and underlying physics and that the predicted 2D correlation function agrees well with observations. In Sec. 3.2 we show that for low mass stars, turbulent fragmentation mediated by radiation feedback can roughly reproduce the observed multiplicities and mass ratio distribution, and provides qualitative agreement with the expected binary distribution function. However, we also show that protostellar disk fragmentation is necessary to explain the short period tail of the distribution. Finally, in Sec. 4 we summarize our findings and conclude.
2 Model and Methodology
2.1 SemiAnalytic Framework
In this study we use an improved version of the MISFIT semianalytical framework introduced in GH15 and GKH16 (see Appendix A for a detailed description of all changes from the previously published version) which allows us to simulate the evolution and fragmentation of GMC sized clouds at a modest computational cost (compared to full radiationhydro simulations). This not only allows the rapid exploration of different initial conditions and underlying physics but also enables a statistical analysis as we are able to simulate an ensemble of clouds.
This, of course, comes at the cost of some simple approximations. The main assumption of MISFIT is that density fluctuations in collapsing GMCs are created by turbulence and thus obey “random walk” statistics (see e.g. Hopkins 2013c). As the cloud collapses it pumps energy into turbulence (so that virial equilibrium is maintained) as motivated by Robertson & Goldreich (2012) and Murray et al. (2015). Unlike most analytical models MISFIT preserves spatial and temporal information and can be easily expanded with additional physics (e.g. equation of state, angular momentum etc.). We show in Appendix B that, despite these strong assumptions, our results are roughly in agreement with the detailed radiation hydrodynamical simulation of Bate (2012).
The simulation starts from a GMC with fully developed turbulence and follows its collapse. The density field of the cloud is resolved on a grid with points and is evolved in Fourier space following a FokkerPlanck equation (see Hopkins 2013b and GH15). For the bulk of this paper we use , and in Appendix C we show that this is sufficient to achieve convergence. Every time a new selfgravitating substructure appears (i.e., the cloud fragments) the code is run recursively for each substructure. When the cloud size reaches the predefined relative size scale (the relative termination scale) the simulation stops. This termination scale represents the length scale where the initial assumptions break down and the selfsimilar fragmentation cascade is terminated.
The primary effect that breaks selfsimilarity and imposes a scale in our calculations is angular momentum, which leads to the formation of a disk once the object has contracted a certain amount. Disk formation is the natural termination scale. In our model the source of angular momentum is random turbulent motion, which in the supersonic limit means that the distribution of the rotational kinetic energy fraction is strongly peaked around a few percent (Burkert & Bodenheimer, 2000), consistent with the observed distribution of protostellar core rotation rates (Goodman et al., 1993). If one translates this into an angular momentum and assumes that the specific angular momentum of fluid elements is conserved during collapse, the characteristic radius of disk formation is . In this paper we adopt as our fiducial value for most calculations, and we explore the sensitivity of the results to our choice in Appendix C.
The initial conditions of the parent clouds are defined by their mass , the sonic length (the scale at which the turbulent velocity dispersion is equal to the sound speed), the sonic mass (the minimum selfgravitating mass at the sonic scale), and the termination scale . All other parameters (e.g. temperature, Mach number) can be derived from these. Moreover, the total mass only affects our results by changing the outer scale of the turbulent cascade, a result we demonstrate in Appendix C, so we shall not discuss it further here. For details about initial conditions and the detailed algorithm see GH15 in which a detailed stepbystep guide to the model is provided in Appendix A.
The final output of the simulation is a list of protostars and their initial properties (e.g. mass, velocity, position). As we are not accounting for later dynamical processes, our results only apply at the time of formation. The leftover unbound material is assumed to escape.
2.2 Implementation of Stellar Feedback
In this paper we investigate the clustering properties of two classes of models: the case of pure isothermal fragmentation and a model with feedback from protostellar heating based on K11. Isothermal turbulence is scalefree (McKee et al. 2010; Krumholz 2014), so we expect no absolute scales in any results (although scales from initial conditions may appear), while the heated model imprints a mass scale that is insensitive to initial conditions (hence there is a peak in the IMF, as shown in GKH16. For comparison we also include some runs where in addition to protostellar heating the gas has a “stiffening” equation of state (EOS). This means that the gas reacts to compression as a subisothermal medium at very large scales, isothermally at intermediate scales and transitions to an adiabatic behavior after reaching a threshold volume density where it becomes opaque to its own cooling radiation. We model this effect using a physically motivated EOS based on Masunaga & Inutsuka (2000) and Glover & Mac Low (2007). In this case the effective polytropic index is depends on the local volume density as
(1) 
where we set and corresponding to and . See GKH16 for more details.
Our treatment of protostellar radiative feedback is a fairly crude approach motivated by K11, and supported numerically by Krumholz et al. (2016). We assume that the center of selfgravitating clouds collapses first, forming a protostellar seed, then the rest of the cloud accretes onto it. The energy of the matter accreted by this seed is radiated within the optically thick core. The temperature of the material depends on the accretion rate onto the protostar (and thus the mass and dynamical time of the gas around it), and on the energy yield per unit mass from accretion, which we denote . The value of is set by the protostellar massradius relation, and K11 shows that it is determined primarily by the effects of deuterium burning, which regulates the central temperatures of protostars. Because deuterium burning begins when protostars are only a few , and, for low mass protostars continues for Myr, it is the dominant factor in setting during the bulk of a molecular cloud’s starforming history. Comparing with detailed protostellar evolution calculations, K11 finds that J kg to better than half a dex accuracy for all protostellar masses in the range , and to better than a dex accuracy from . We therefore adopt this value of throughout the remainder of this paper. If we assume a spherically symmetric system then ,following K11, the gas at distance from an accreting protostar is heated up to a temperature of
(2) 
where is the mass enclosed in radius , while , are the gravitational and StefanBoltzmann constant respectively. Crudely, this scaling reflects energy conservation as for the opaque cloud (see K11 for more details). Combined, internal heating and the physical processes captured by the EOS of the model set the temperature as
(3) 
Note that in the feedback only case we use an isothermal EOS which means that where is the initial temperature of the cloud.
It is important to note that protostellar feedback is not scalefree. By using Eq. 2 and assuming virial equilibrium we can find the length scale around a protostar below which heating becomes important ():
(4) 
where is the mean molecular weight measured in units of and is the Boltzmann constant. For our numerical calculations we adopt , appropriate for fully molecular H with 1 He per 10 H nuclei. We can similarly find the characteristic mass scale
(5) 
that sets the peak of the IMF (see K11 for a more detailed calculation that leads to ).
To easily identify the results for different models and parameters we use the labels shown in Table 1. The STD label refers to initial conditions similar to Milky Way GMCs, while ULIRG runs have the very high temperature and strong turbulence characteristic to the clouds of Ultra Luminous Infrared Galaxies (ULIRGs). There are also a number of runs where the physical parameters are not varied but the numerical ones are, so that we can identify numerical artifacts in our results.
Label  Input Parameters  Derived Parameters  Thermodynamics  
[]  N  []  []  [K]  []  
Isohermal  MW  32  1.6  0.1  10  11.1  10.5  Isothermal  
Isothermal_SmallR  32  1.6  0.1  10  11.1  10.5  Isothermal  
Isohermal  ULIRG  32  0.31  0.0026  75  0.66  13.1  Isothermal  
Heating  MW  32  1.6  0.1  10  11.1  10.5  Protostellar Heating  
Heating  ULIRG  32  0.31  0.0026  75  0.66  13.1  Protostellar Heating  
Heating_N16  16  1.6  0.1  10  11.1  10.5  Protostellar Heating  
Heating_N64  64  1.6  0.1  10  11.1  10.5  Protostellar Heating  
Heating_M1E3  32  1.6  0.1  10  3.5  5.2  Protostellar Heating  
Heating_M1E5  32  1.6  0.1  10  35.4  16  Protostellar Heating  
Heating_SmallRmin  32  1.6  0.1  10  11.1  10.5  Protostellar Heating  
Heating_LargeRmin  32  1.6  0.1  10  11.1  10.5  Protostellar Heating  
Heating+EOS  MW  32  1.6  0.1  10  11.1  10.5  Heating+EOS  
Heating+EOS  ULIRG  32  0.31  0.0026  75  0.66  13.1  Heating+EOS 
2.3 Clustering and Multiplicity Statistics
For each simulation we have as output a list of stellar masses and positions. From these, we compute several statistical quantities describing the stellar distribution. Our first quantity of interest is the correlation function. In this paper we adopt the usual definition of the 3D correlation function ,
(6) 
where is the number of objects whose distance is , is the density of objects, and denotes ensemble averaging.
We can similarly define the 2D correlation function , which is identical to except that one computes the distance only in 2 of the 3 orthogonal directions. Unlike , the stellar is measurable and it is easy to show that , where is the mean surface density of stars measured in an annulus at distance from other stars. For the purpose of generating quantities that can be readily compared to observations, we must also account for sensitivity limits, which make it difficult to detect low mass objects. Since studies of stellar correlation have been performed with a wide range of sensitivities, we simply choose a roughly representative limiting mass , and compute our correlation function using only stars more massive than this limit.
Our second characteristic of interest is the multiplicity properties of the stars – both the multiplicity fractions and the distribution of periods and mass ratios. Since our calculations involve no dynamical evolution after the protostars are formed, deriving these statistics is not trivial, as the newly formed stars form a fractallike structure where each star is bound to a number of other stars. Such a configuration is expected for young star clusters based on simulations, and is completely consistent with the observed distribution of newlyformed stars (e.g., Kruijssen, 2009; Bate, 2009b; Krumholz et al., 2012). However, it makes identification of distinct, bound systems difficult, and leads to structures which are very unlikely to survive for even a single cluster crossing time (e.g. nonhierarchical quadruple systems orbiting each other). Thus it is important that we try to correct for this behavior. In this paper we use the hierarchical algorithm introduced by Bate (2009b), which has the following steps:

Calculate the binding energy between all pairs of stars.

Find the most bound pair and replace it with a single point mass with the same total mass and momentum, located at the center of mass of the removed pair.

Recursively repeat steps 1 and 2 until no more bound stars are left, with the exception that we do not combine pairs of objects if the resulting bound aggregate would consist of more than 4 individual stars. If such an aggregate is the most bound pair at any point, proceed to the next most bound pair, terminating if no other bound pair exists. Also, in order to get stable, hierarchical multiples we require that the period of a newly assigned star is at least ten times higher than that of the original aggregate.
This algorithm provides a list of single, binary, triple and quadruple star systems with which we can calculate the multiplicity fraction , defined as
(7) 
where are the number of single, binary, triple, quadruple systems within which the most massive star (primary star) has mass . This definition has the advantage that it can be observed fairly robustly (Hubber & Whitworth 2005; Bate 2009b), as this does not differentiate between the classes of multiple star systems, so does not change if a new companion star is discovered in a binary system. Note that to account for the decreased sensitivity of observations to very low mass stars we neglect companions below .
3 Results
3.1 The Stellar Correlation Function
Figure 1 shows the predicted stellar correlation function for a selection of our models, computed using an ensemble average of the GMCs we have run for each case. The shape is close to a powerlaw, with , with a cutoff at the size scale of the parent GMC. These properties are remarkably robust to changes in initial conditions and even to changes in the underlying small scale physics.
Qualitatively, the isothermal purepowerlaw behavior can be understood with a simple toy model: small objects form after significant contraction and a number of fragmentation events for which the physics is selfsimilar in the isothermal case. So imagine that a cloud contracts by a factor of , then fragments into two equalmass fragments. Then each of these two fragments contracts and produces two more fragments, and so on. This prescription is similar to a wellstudied fractal, the Cantor Set (Cantor Dust more specifically). For this, the correlation function is a powerlaw with slope for (see Appendix D).
Figure 2 compares our predictions to the observed surface density of stars (proportional to the 2D correlation function). In examining this plot, note that the absolute values of the stellar surface densities are not meaningful, since these are just dictated by a our choice of sonic length, and thus can be tuned freely by considering slightly different physical scalings, exactly as one might expect when considering a range of starforming regions of widely varying density, mass, and velocity dispersion. Instead, the meaningful comparison is the shapes of the functions. In this regard, we see that the simulated correlation functions have a slope quite similar to the observations at scales larger than . Below this scale our models cannot reproduce the significant steepening of the correlation function. We show below that this directly manifests in the distribution of short period binaries where the simulation fall short of observations at the same scale. This is the length scale of the largest protostellar disks, which suggests that disk physics (which are neglected in these models) is responsible for the steepening. However, we must stress that dynamical relaxation does affect the observed, finite age systems and is probably responsible for outlier systems like Trapezium and Upper Sco. Both of these systems are dynamically older, in the sense that they have existed for more crossing times, than the other systems shown, which supports this conjecture. One should be careful not to draw the false conclusion that this model fully explains the observed stellar spatial distribution simply because it reproduces the correlation function. It has been shown that very different geometries (e.g., fractal vs spherical) can lead to similar correlation function slopes (Gouliermis et al., 2014). Nevertheless we can say that this model is at least consistent with the observed stellar correlation function in the large scale, fractallike regime.
3.2 Multiplicity
After grouping stars into bound systems following the procedure described in Section 2.3, we assign each star one of the following labels:

Single: The star is not bound to any other stars.

Multiple: The star is the most massive (primary) star of a multiple star system.

Nonprimary: The star is part of a multiple star system, but it is not the primary star.
Not all of the companion stars that emerge from our analysis would be detectable by current techniques. In particular, brown dwarf companions to main sequence stars are quite hard to detect. Therefore we must correct for completeness before comparing to observations. As a guide to the current observational capabilities, we follow the summary given in Table 8 of De Rosa et al. (2014). Based on this summary, we apply the following cuts to our data:

For primaries with mass , we discard any companions with masses below

For primaries with mass , we discard companions for which the secondary to primary mass ratio is .
While these cuts are only an approximate representation to the diversity of observational survey selection functions in the literature, they provide a reasonable approximation to the capabilities of the current state of the art.
Fig. 3 shows the fraction of stars in each of the three classes as a function of stellar mass before and after applying the observational bias. Since the isothermal model has no inherent physical scale, in an ideal case we would not expect any mass dependence. However, the finite initial mass and the cutoff imposed by observational selection affect the results. The former leads to finite size effects at larger masses. Specifically, since there is a finite total mass, there must be a single most massive star, and for obvious reasons it cannot be nonprimary. Similarly, other stars that are near the most massive are also biased against being nonprimary. This effect is responsible for the decline in the nonprimary curve at high masses. At the other end of the mass spectrum, the fact that brown dwarf companions to hydrogenburning primaries are difficult to detect explains the sharp decline in the nonprimary fraction and sharp rise in the multiple fraction for the lowest mass bin. The sharp change in behaviour above and below has a simple explanation: for hydrogenburning stars, multiplicity surveys are primarily conducted in the field, while for brown dwarfs they are mainly conducted in young clusters. Since brown dwarfs are easier to detect in young clusters than in the field, surveys of brown dwarf primaries are much more complete in finding brown dwarf companions than surveys of hydrogenburning primaries.
The case including stellar radiation feedback looks qualitatively similar to the isothermal one. In both cases we recover the simple rule that more massive stars tend to be the primary stars of systems while smaller stars tend to be their companions. However, most of the stars in the isothermal model were born in systems of multiple stars, while there is a significant number of single stars in the radiative heating case. The transition from where it becomes more common for stars to be the primaries of multiple systems than to be single is which is a result of the peak in the stellar mass function imposed by heating, which suppresses the formation of stars below the IMF peak at . Smaller stars are unlikely to be primaries mainly because there are increasingly few lower mass stars available to be their companions.
We can also compare the results of our models to observations. In Fig. 3 we plot in the right panel the results of Moe & Di Stefano (2016), based on analysis of the observations of Raghavan et al. (2010). Compared to these observations, our model slightly overpredicts the multiplicity of solar type stars and underpredicts the fraction that are single.
We compare the stellar multiplicity as a function of primary mass with observations in Figure 4. We find that the observed results for the heated and isothermal cases are similar, and both qualitatively reproduce the observational result that the multiplicity fraction is near unity for primaries substantially above , dropping to tens of percent for or smaller primaries. However, the apparent similarity between the observed distributions for the isothermal and heated cases is primarily an illusion due to observational completeness effects. In the isothermal case, essentially every primary has an undetected brown dwarf companion, and thus the true multiplicity fraction for primaries of this mass is close to unity. It is only our inability to detect these brown dwarfs that makes the predicted distribution in the isothermal case at all compatible with the observations.
It is also worth investigating how these results would be affected by protostellar disk fragmentation. To do so, we construct a toy model for protostellar disk fragmentation that can be used to postprocess our simulation results, based on the works of Kratter et al. (2010) and Offner et al. (2010). They define the thermal parameter of the disk as , where is the infall accretion rate onto the disk and is the sound speed of the disk. Both find that is the main parameter in determining whether a protostellar disk fragments or not. Physically, is just the ratio of the accretion rate into the disk to the maximum rate at which a gravitationally stable disk with dimensionless viscosity can deliver mass to the central star, which is .
Using our protostellar heating prescription we can express where we have used that in an opaque medium. As noted above, the outer edge of the disk should be found at a radius , where is the rotational kinetic energy divided by the binding energy. Using Eq. 2 to evaluate at this radius yields
(8) 
where we used with as the molecular weight of the gas. has a weak dependence on the radius so we can safely use the protostellar disk size scale. This leads to . Based on Fig. 2 of Kratter et al. (2010) fragmentation is very likely if , which corresponds to collapsing final fragment masses in our model. Thus, in this crude approximation, the only effect protostellar disk fragmentation would have on our multiplicity fraction in Fig. 4 is that it would reach unity at a somewhat lower stellar mass. Our conclusion that low mass disks are for the most part too warm to fragment, but that disk fragmentation should be common for somewhat superSolar and larger stars, is consistent with the numerical results of Offner et al. (2010).
3.3 Demographics of the Binary Population
3.3.1 Mass Ratios and the Brown Dwarf Desert
One of the key observed properties of binaries is the apparent flat mass distribution of companion masses with a cutoff at very low masses (the socalled “brown dwarf desert”). In Fig. 5 we test to what extent our models can reproduce this observation by comparing the mass distribution of the most massive companions in our simulations with observations of this quantity for Solar type and very low mass (VLM) stars (). Although in principle we could compute other mass ratios (e.g., the mass ratios of all pairs of stars in multiple systems, c.f. Raghavan et al. 2010), we focus on the most massive companions because these are the most robustly determined from observations. It is extremely challenging observationally to identify secondary and tertiary companions of a star in a triple or quadruple system. As a result, observations are most likely to discover the most massive companion rather than all companions, making the primary to secondary mass ratio the most welldetermined. This also has the advantage that the most massive companion is the least likely to be ejected by dynamical processes. For our heated models, however, in practice it makes relatively little difference whether we include all companions or just the most massive one.
Examining Figure 5, it is clear that the isothermal and heated models are both roughly consistent with observations after the observational bias is applied. The primary exception is that both models somewhat underpredict the frequency of nearequal mass companions; such companions can plausibly be attributed to disk fragmentation, which tends to produce mass ratios close to unity (Bate, 2000). It is important to note that without the observational bias the isothermal model predicts an overwhelming number of very low mass ratio companions. Meanwhile the results for Solar type stars in the heated case is only slightly affected by observational bias, which means that the brown dwarf desert is not an observational bias.
To gain further insight into why our heated models are able to reproduce the brown dwarf desert, while our isothermal models fail, let us compare these companion mass distributions with the null hypothesis that that companion masses are randomly drawn from the IMF^{2}^{2}2Binaries forming from randomly sampling the IMF has been ruled out (Reggiani & Meyer, 2011), making it an important test for theoretical models.. Fig. 6 compares our measured companion mass ratio distribution with that we would expect under the null hypothesis, again considering only the most massive pair in a given star system^{3}^{3}3Computing the null hypothesis distribution requires some care, because for systems with stars, even if all companions are drawn randomly from the IMF, the mass distribution for the most massive companion does not follow the IMF. Specifically, suppose we have an IMF , so that the cumulative distribution function (CDF) of masses (i.e., the probability that a randomly chosen star has mass ) is . Now consider a system where the primary has mass . Since we require companion masses to be smaller than , they follow the conditional CDF , which for has the same shape as the CDF for single stars. However, now consider a system consisting of stars. The most massive companion has a mass only if all companions have mass , and if the companion masses are independent the probability of this is . This does not have the same shape as the single star CDF. For the purposes of Figure 6, we account for this effect by generating our null hypothesis lines as a weighted sum , where both the single star CDFs and the relative frequencies of multiplicity are measured directly from the simulations for each primary mass .. The figure shows that in no case are the results consistent with random sampling of companions from the IMF. In the isothermal case the companion distribution for both Solar and VLM primaries follows the IMF for very low mass companions, but that there is a significant excess of companions at mass ratios . In the heated case the situation is qualitatively similar, in that mass ratios near unity are overrepresented compared to the null hypothesis.
Now let us consider the implications of this finding for the brown dwarf desert. The companion mass ratio distribution is a product of two factors: the underlying IMF of all stars, and any biases imposed by the fact that the stars whose mass distribution we are computing are nonprimaries. With or without heating, we find that mass ratios near unity are favoured compared to a null hypothesis of random IMF sampling. That is, if we collect two samples of stars with the same upper mass limit, and for one sample we randomly select only nonprimary stars and for the other we randomly select stars without regard to multiplicity characteristics, the nonprimary sample will typically be more massive. For Solartype primaries, the combination of a bias towards higher mass companions and the overall negative slope of the IMF near 1 (so that lower mass stars are more probable overall) yields a relatively flat mass ratio distribution – the IMF shape and the bias nearly cancel.
Now let us consider VLM stars. For VLM primaries, the bias towards equal mass companions is qualitatively similar to that for Solartype stars. For our isothermal case, and unlike in reality, the IMF slope near is also about the same as that near , due to the overall scalefree nature of isothermal fragmentation. Because both the IMF slope and the bias are about the same for Solar and VLM stars, the distribution of companion mass ratios is also qualitatively similar. For our heated case, as in reality, we have a very different situation. The slope of the IMF is negative near , but positive (or at least close to flat) near . As a result, for VLM primaries both the bias towards massive companions and the IMF itself favour more massive objects as companions. The result is a companion mass ratio distribution that is sharply biased towards stellar companions and away from brown dwarfs, producing the observed brown dwarf desert. We therefore find that the brown dwarf desert is a result of the change in the IMF slope between and , which in turn is imposed by thermal feedback causing a deviation from scalefree behaviour during gas collapse and fragmentation.
3.3.2 Binary Separations
In addition to the mass ratio distribution, our spatiallyresolved model allows us to examine the predicted semimajor axis distribution of binaries. We do so in Fig. 7 for Solartype stars. The distribution appears peaked which comes from the peak of the companion mass ratio distribution (Fig. 6) with the corresponding length scale of .
Comparing with the observations from Marks et al. (2011) we can see that on large scales our model of common core fragmentation seems to very roughly reproduce the present day observations. Although our results only give the “atformation” period distribution, the comparison is still meaningful because, as explained in Sec. 2, we have attempted to limit the systems we count in our model to hierarchical systems that should be dynamically stable. In any case, it is clear that, similar to the case of the 2D correlation function (Fig. 2), turbulent fragmentation is unable to reproduce the observations on small scales. Another source for such binaries is required at , for which protostellar disk fragmentation is a good candidate^{4}^{4}4 It should be noted that disc fragmentation simulations also fail to produce extremely close binaries (). These are likely to have either formed from a wider binary whose separation decreased due to orbital decay (e.g. Stahler 2010; Korntreff et al. 2012), or from exchange interactions in star clusters.. Note that decreasing would technically improve the fit (see Appendix C), but this would require unphysically low values.
4 Conclusions
The aim of this paper is to investigate the origin of the stellar correlation function and multiplicity statistics, and in particular to understand which features of these distributions result from pure scalefree isothermal fragmentation, and which bear the imprints of scaledependent stellar feedback. Using the MISFIT semianalytical turbulent fragmentation framework of GH15 and GKH16 we find that the shape of the correlation function is almost entirely set by isothermal turbulence. Stellar feedback, which operates primarily on small scales, has little effect. On smaller scales () both a purely isothermal model and one including stellar radiation feedback underpredict the stellar correlation, suggesting that our turbulent fragmentation models lack certain small scale physics (likely protostellar disk fragmentation). As with the correlation function, we find that our models provide a reasonable match to the observed the binary period distribution at large separations regardless of whether we include protostellar heating or not, but that we fail to produce enough very close binaries. We again conjecture that these close binaries are a result of disk fragmentation and Nbody interaction, which our model does not include.
The situation for the mass ratios and multiplicity fraction of binaries is quite different. Isothermal fragmentation produces far too many multiple stars compared to what is observed, with even stars predicted to have multiplicities near unity. Adding protostellar heating substantially improves the situation, though the multiplicity fraction is still somewhat too high, likely because our models do not include dynamical evolution that will disrupt unstable systems. These differences, however, are almost completely washed out by the observational biases.
Most interestingly, if we neglect the observational bias we find that while turbulent fragmentation with or without protostellar heating can adequately reproduce the observed companion mass distribution for Solar type stars (except for very low mas companions), but only when protostellar heating is included can we reproduce the mass distribution of companions for lowmass primaries. In particular, only our models including radiative feedback reproduce the “brown dwarf desert”, whereby the companions to low mass stars ( ) are overwhelmingly stellar objects (i.e., close to a mass ratio of unity) rather than brown dwarfs. Models that include only scalefree physics predict a companion mass ratio distribution for low mass stars that is qualitatively similar to that for Solartype stars, a direct consequence of the scalefree nature of these models. In contrast, protostellar heating suppresses the number of brown dwarfs relative to stars, so that the companion mass ratio distribution is very different for Solartype stars that lie above the IMF peak and lowmass stars that lie at ore below it. We therefore conclude that the brown dwarf desert is a consequence of the physical mass scale imprinted by protostellar heating into the otherwise scale free star formation process.
Acknowledgments
Support for PFH and DG was provided by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1411920 and CAREER grant #1455342. MRK is supported by NASA ATP grant NNX15AT06G and Australian Research Council grant DP160100695. Numerical calculations were run on the Caltech compute cluster “Zwicky” (NSF MRI award #PHY0960291) and allocation TGAST130039 granted by the Extreme Science and Engineering Discovery Environment (XSEDE) supported by the NSF.
References
 Bate (2000) Bate M. R., 2000, MNRAS, 314, 33
 Bate (2009a) Bate M. R., 2009a, MNRAS, 392, 590
 Bate (2009b) Bate M. R., 2009b, MNRAS, 392, 590
 Bate (2012) Bate M. R., 2012, MNRAS, 419, 3115
 Bate et al. (1998) Bate M. R., Clarke C. J., McCaughrean M. J., 1998, MNRAS, 297, 1163
 Bressert et al. (2010) Bressert E., et al., 2010, MNRAS, 409, L54
 Burgasser et al. (2007) Burgasser A. J., Reid I. N., Siegler N., Close L., Allen P., Lowrance P., Gizis J., 2007, Protostars and Planets V, pp 427–441
 Burkert & Bodenheimer (2000) Burkert A., Bodenheimer P., 2000, ApJ, 543, 822
 Cartwright & Whitworth (2004) Cartwright A., Whitworth A. P., 2004, MNRAS, 348, 589
 De Rosa et al. (2014) De Rosa R. J., et al., 2014, MNRAS, 437, 1216
 DelgadoDonate et al. (2004) DelgadoDonate E. J., Clarke C. J., Bate M. R., Hodgkin S. T., 2004, MNRAS, 351, 617
 Duchêne (1999) Duchêne G., 1999, A&A, 341, 547
 Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, ARA&A, 51, 269
 Glover & Mac Low (2007) Glover S. C. O., Mac Low M.M., 2007, ApJS, 169, 239
 Goodman et al. (1993) Goodman A. A., Benson P. J., Fuller G. A., Myers P. C., 1993, ApJ, 406, 528
 Goodwin et al. (2004) Goodwin S. P., Whitworth A. P., WardThompson D., 2004, A&A, 414, 633
 Goodwin et al. (2007) Goodwin S. P., Kroupa P., Goodman A., Burkert A., 2007, Protostars and Planets V, pp 133–147
 Gouliermis et al. (2014) Gouliermis D. A., Hony S., Klessen R. S., 2014, MNRAS, 439, 3775
 Gouliermis et al. (2015) Gouliermis D. A., et al., 2015, MNRAS, 452, 3508
 Guszejnov & Hopkins (2015) Guszejnov D., Hopkins P. F., 2015, preprint, (arXiv:1507.06678)
 Guszejnov et al. (2016) Guszejnov D., Krumholz M. R., Hopkins P. F., 2016, MNRAS, 458, 673
 Hansen et al. (2012) Hansen C. E., Klein R. I., McKee C. F., Fisher R. T., 2012, ApJ, 747, 22
 Hartmann (2002) Hartmann L., 2002, ApJ, 578, 914
 Hennekemper et al. (2008) Hennekemper E., Gouliermis D. A., Henning T., Brandner W., Dolphin A. E., 2008, ApJ, 672, 914
 Hopkins (2013a) Hopkins P. F., 2013a, MNRAS, 428, 1950
 Hopkins (2013b) Hopkins P. F., 2013b, MNRAS, 430, 1653
 Hopkins (2013c) Hopkins P. F., 2013c, MNRAS, 430, 1880
 Hubber & Whitworth (2005) Hubber D. A., Whitworth A. P., 2005, A&A, 437, 113
 Kaczmarek et al. (2011) Kaczmarek T., Olczak C., Pfalzner S., 2011, A&A, 528, A144
 Klessen & Burkert (2000) Klessen R. S., Burkert A., 2000, ApJS, 128, 287
 Korntreff et al. (2012) Korntreff C., Kaczmarek T., Pfalzner S., 2012, A&A, 543, A126
 Kouwenhoven et al. (2010) Kouwenhoven M. B. N., Goodwin S. P., Parker R. J., Davies M. B., Malmberg D., Kroupa P., 2010, MNRAS, 404, 1835
 Kratter et al. (2010) Kratter K. M., Matzner C. D., Krumholz M. R., Klein R. I., 2010, ApJ, 708, 1585
 Kraus & Hillenbrand (2008) Kraus A. L., Hillenbrand L. A., 2008, ApJ, 686, L111
 Kraus et al. (2008) Kraus A. L., Ireland M. J., Martinache F., Lloyd J. P., 2008, ApJ, 679, 762
 Kraus et al. (2011) Kraus A. L., Ireland M. J., Martinache F., Hillenbrand L. A., 2011, ApJ, 731, 8
 Kroupa (1995) Kroupa P., 1995, MNRAS, 277
 Kruijssen (2009) Kruijssen J. M. D., 2009, A&A, 507, 1409
 Krumholz (2011) Krumholz M. R., 2011, ApJ, 743, 110
 Krumholz (2014) Krumholz M. R., 2014, Phys. Rep., 539, 49
 Krumholz et al. (2012) Krumholz M. R., Klein R. I., McKee C. F., 2012, ApJ, 754, 71
 Krumholz et al. (2016) Krumholz M. R., Myers A. T., Klein R. I., McKee C. F., 2016, MNRAS, 460, 3272
 Kuiper (1935) Kuiper G. P., 1935, PASP, 47, 15
 Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
 Low & LyndenBell (1976) Low C., LyndenBell D., 1976, MNRAS, 176, 367
 Marks et al. (2011) Marks M., Kroupa P., Oh S., 2011, MNRAS, 417, 1684
 Masunaga & Inutsuka (2000) Masunaga H., Inutsuka S.i., 2000, ApJ, 531, 350
 McKee et al. (2010) McKee C. F., Li P. S., Klein R. I., 2010, ApJ, 720, 1612
 Moe & Di Stefano (2016) Moe M., Di Stefano R., 2016, preprint, (arXiv:1606.05347)
 Murray et al. (2015) Murray D. W., Chang P., Murray N. W., Pittman J., 2015, preprint, (arXiv:1509.05910)
 Nakajima et al. (1998) Nakajima Y., Tachihara K., Hanawa T., Nakano M., 1998, ApJ, 497, 721
 Offner et al. (2010) Offner S. S. R., Kratter K. M., Matzner C. D., Krumholz M. R., Klein R. I., 2010, ApJ, 725, 1485
 Parker (2014) Parker R. J., 2014, MNRAS, 445, 4037
 Parker & Meyer (2014) Parker R. J., Meyer M. R., 2014, MNRAS, 442, 3722
 Portegies Zwart et al. (2010) Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431
 Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
 Reggiani & Meyer (2011) Reggiani M. M., Meyer M. R., 2011, ApJ, 738, 60
 Robertson & Goldreich (2012) Robertson B., Goldreich P., 2012, ApJ, 750, L31
 Simon (1997) Simon M., 1997, ApJ, 482, L81
 Stahler (2010) Stahler S. W., 2010, MNRAS, 402, 1758
 Stanke et al. (2006) Stanke T., Smith M. D., Gredel R., Khanzadyan T., 2006, A&A, 447, 609
 Tohline (2002) Tohline J. E., 2002, ARA&A, 40, 349
Appendix A Improvements to Previous Model
Two papers (Guszejnov & Hopkins 2015; Guszejnov et al. 2016) have been published so far using the MISFIT (Minimalist Star Formation Including Turbulence) semianalytical star formation framework as this paper. Since the publication of those results, several improvements have been made to the algorithm, all of which are implemented for this paper. These do not change any of the published qualitative IMF results (e.g. general shape, sensitivity to initial conditions). They include:

Correction of a bug that suppressed fragmentation at the end of the cloud evolution, violating selfsimilarity at a weak level. The effects on our previous work are small, but are substantial on the statistics on low mass companions. This is now fixed.

Fragments are properly tracked and taken into account for the evolution of their parent (e.g. their contribution to the gravitational potential is taken into account as long as the parent has not yet contracted beyond their position). This causes no qualitative difference.

Instead of using an absolute termination scale (taken to be in the previous papers, roughly the size of protostellar disks), the collapse of clouds now terminates once clouds have contracted to a fixed fraction of their initial radius, chosen to be roughly when angular momentum support becomes dominant. This assumes that the source of angular momentum for clouds is from random turbulent motion. The resulting distribution for is strongly peaked around a few percent (Burkert & Bodenheimer 2000). If collapse happens at constant virial parameter than the size scale where angular momentum starts dominating is . This is the scale where the cloud flattens and forms a disk, which we choose as our termination point.^{5}^{5}5The choice of a relative termination scale instead of an absolute value has the added benefit of imprinting no absolute length scale into the problem, preserving selfsimilarity.

We set a lower limit of on fragment masses based on the opacity arguments of Low & LyndenBell (1976). This is in fact equivalent to a simplified EOS model, where we terminate the fragmentation once the cloud reaches the adiabatic limit. This provides a natural termination for the fragmentation cascade in our “isothermal” models, otherwise our results would not converge (see GKH16)
Appendix B Comparison with detailed hydrodynamic simulations
There have been a number of hydrodynamical simulations attempting to find the multiplicity statistics and separation distribution of newly formed stars (e.g. Goodwin et al. 2004; DelgadoDonate et al. 2004; Bate 2009a, 2012; Offner et al. 2010; Krumholz et al. 2012). The semianalytical approach we present in this paper has several advantages over these (e.g. faster, no absolute resolution limit, starts from GMC) but at the cost of several strong assumptions, so it is crucial that we compare our results with theirs. We choose to compare with the simulations of Bate (2012), as these have the largest sample of multiple systems. Since we know the full binary distribution from the simulations, we make this comparison without the observational completeness correction that we apply when comparing to observations in the main text.
Figures 810 show that our results are qualitatively, and in many cases quantitatively, consistent with the simulations of Bate (2012). An important difference between our current model and the traditional simulations is that MISFIT does not have a finite resolution limit, but it does neglect disk physics. This leads to the discrepancy at small separations shown in Fig. 10.
Appendix C Numerical Tests and Convergence
In this section we show how the GMC mass, the termination scale , and the resolution parameter affect our results. To explore these questions, we have repeated our fiducial Heating  MW run with different resolutions (, GMC masses ( ), and different termination scales (). The full list of runs performed in given in Table 1. Note that unlike the results in the main text these are not modified to account for observational biases.
We first examine how our results affect the same of the IMF produced by our models. Fig. 11 shows that the shape of the IMF is robust to changes in , reaching convergence around . The location of the IMF peak and the high mass slope are also essentially insensitive to the GMC mass; the location of the peak does shift by an extremely small amount as we vary the GMC mass, as a result of its dependence on the Mach number of the turbulence; the two are related thanks to our assumption that clouds have virial ratios . However this shift is only over a plausible range of GMC masses. The IMF shows its greatest sensitivity to the relative termination scale , particularly the abundance of brown dwarfs beyond the peak. In reality our choice of a single value is an oversimplification, since real turbulence fields produce a distribution of rotational kinetic energies , and thus a distribution of parameters; the real IMF should therefore resemble a weighted average of the curves shown in Fig. 11.
Fig. 12 shows how variation of our three parameters affects the stellar correlation function. As with the IMF, we find that the correlation function is insensitive to changes in both the resolution parameter and the initial GMC mass – the former produces no noticeable differences past , while the latter mostly rescales the outer cutoff/size scale. Also, as with the IMF, larger initial masses lead to verly slightly shallower slopes.This is consistent with the discussion in Section 3.1 and Appendix D: larger masses mean stronger initial turbulence which in turn means easier fragmentation. However, as with the IMF, the effect is extremely modest. Finally, the bottom panel of Fig.12 shows that the relative termination scale has no effect on the correlation function apart from introducing a smallscale cutoff.
Fig. 13 shows the peak of the separation distribution converges above , and the parent GMC mass has little effect on it. The relative termination scale sets the width of the peak; its position is set by , as discussed in Sec. 3.2. Decreasing increases the abundance of binaries at small separations, since it pushes the transition between common core fragmentation (which we are modeling) and disc fragmentation (which are are not) to smaller scales. However, as noted in the main text, only a value of that is unphysically small would produce enough shortperiod binaries to be consistent with the observations.
Appendix D Cantorlike Model of Fragmentation
One of the most important properties of isothermal fragmentation is that it is scalefree, so we expect selfsimilar, fractallike structures to emerge. We can formulate a simple toy model to describe this process where selfgravitating clouds contract to about relative scale before breaking into two (see Fig. 14) along a random axis. The distance of the two fragments is uniformly chosen between 0 and . The fragments then rearrange themselves into spheres at the same density as their parent (meaning their radius is ). This model is very similar to the generalized Cantor dust (3D analogue of the generalized 1D Cantor set) that have the fractal dimensions of and , leading to a 3D correlation function of .
We can analytically calculate the fractal dimension of our Cantorlike model if we take the separation between fragments to be the mean value of . If we take the initial radius of the first sphere to be unity then, after iterations, the number of the objects is while their size is . If we choose a random fragment then the number of fragments within an radius is , thus
(9) 
Fig. 15 shows that this result is actually exact. Since isothermal fragmentation is quite similar to this toy model, we expect the predicted stellar correlation function to have a slope between 1 and 3 (for reasonable values). This also implies that if some additional physics makes fragmentation harder (increasing the density threshold and thus decreasing ) then the correlation function becomes steeper.
This model has the free parameter which we can restrict by assuming that the fragmentation criteria is set by the Jeansinstability. It is known that the mass of fragments would be of the Jeansmass . The number of fragments is which leads to . For this value Eq. 9 gives leading to a power law, in perfect agreement with our results from Fig. 1.