Disc truncation in embedded star clusters:
Dynamical encounters versus faceon accretion
Key Words.:
accretion, accretion discs âprotoplanetary discs â planetary systems: formation – stars: formation – galaxies: star clusters: generalObservations indicate that the dispersal of protoplanetary discs in star clusters occurs on time scales of about 5 Myr. Several processes are thought to be responsible for this disc dispersal. Here we compare two of these processes: Dynamical encounters and interaction with the interstellar medium, which includes faceon accretion and ram pressure stripping. We perform simulations of embedded star clusters with parameterisations for both processes to determine the environment in which either of these processes is dominant. We find that faceon accretion, including ram pressure stripping, is the dominant disc truncation process if the fraction of the total cluster mass in stars is regardless of the cluster mass and radius. Dynamical encounters require stellar densities pc combined with a mass fraction in stars of to become the dominant process. Our results show that during the embedded phase of the cluster, the truncation of the discs is dominated by faceon accretion and dynamical encounters become dominant when the intracluster gas has been expelled. As a result of faceon accretion, the protoplanetary discs become compact and their surface density increases. In contrast, dynamical encounters lead to discs that are less massive and remain larger.
1 Introduction
The clustered environments in which stars are generally born (Lada & Lada, 2003) can have a severe impact on the protoplanetary discs of their newborn members. The occurrence frequency of protoplanetary discs in young star clusters is observed to decrease with an efolding time scale of the order of 5 Myr (Cloutier et al., 2014). Both internal and external processes influence the dispersal of protoplanetary discs in dense star clusters. Internal processes, such as depletion by accretion onto the central star (e.g. Hartmann et al., 1998) and subsequent photoevaporation by the central star (e.g. Bally & Scoville, 1982; Shu et al., 1993; Clarke et al., 2001; Alexander et al., 2006a, b) are believed to dissipate the discs on time scales Myr (Dullemond et al., 2006; Williams & Cieza, 2011).
Furthermore, observations indicate that external photoevaporation by nearby O stars can reduce the fraction of stars that have a disc by a factor of two (e.g. Balog et al., 2007; Guarcello et al., 2007, 2009; Fang et al., 2012, 2013), although these results can, in part, be explained by sample incompleteness (Richert et al., 2015). Vicente & Alves (2005) find no obvious correlation between disc sizes and their proximity to OB stars in the Trapezium cluster. The massloss rate from the outer edge of the protoplanetary disc due to external photoevaporation can be up to (Facchini et al., 2016). However, to truncate disc radii to values smaller than 100 AU within a time scale of 10 Myr, ultraviolet fluxes that are at least an order of magnitude higher than typical values in a moderatesized cluster () are required (Adams, 2010). Mann & Williams (2010) showed that despite the hostile environment around the massive star Ori C, the discs in the Orion Nebula cluster (ONC) have similar properties to those observed in isolated lowmass starforming regions. In this work, we focus on embedded clusters, in which the presence of the pristine gas is expected to absorb radiation and reduce the effect of external photoevaporation.
Apart from (external) photoevaporation, the effects of stellar flybys on protoplanetary discs have also been investigated (e.g. Clarke & Pringle, 1993; Ostriker, 1994; Heller, 1995; Kobayashi & Ida, 2001). Scally & Clarke (2001) modelled the ONC, using fixed disc radii of 10 and 100 AU, and find that stellar encounters probably do not play a significant role in the destruction of protoplanetary discs. However, Olczak et al. (2006) also modelled the ONC and their study suggests that up to 15 % of the discs may be destroyed by stellar encounters. Moreover, Portegies Zwart (2016) finds that the observed distribution of disc sizes in the Orion Trapezium cluster can be reproduced by only considering truncation due to dynamical encounters, assuming initial disc radii of 400 AU and an initial fractal distribution for the stars. Both hydrodynamical and gravitational body studies show that truncation due to dynamical encounters is important in environments with stellar densities (Rosotti et al., 2014; Vincke et al., 2015). Observational studies also show that the observed size distribution of protoplanetary discs depends on the ambient stellar density (de Juan Ovelar et al., 2012).
In recent studies, Wijnen et al. (2016, 2017) suggest that the ambient gas in embedded starforming regions can also reduce the sizes of protoplanetary discs. On the one hand, the ram pressure exerted by the interstellar medium (ISM) can strip the outer parts of the protoplanetary disc. On the other hand, the continuous accretion of ISM with little to no azimuthal angular momentum with respect to the disc material, which we call faceon accretion, causes the disc to contract. Wijnen et al. (2017) provide a parameterisation for the evolution of the disc as a function of the ambient gas density and relative velocity. Here we use this parameterisation in body simulations with a static gas potential and we compare its effect with that of dynamical encounters. Vincke & Pfalzner (2016) find that disc truncation due to dynamical encounters is most effective in the embedded phase of the cluster. We investigate in which part of the parameterspace each of these two effects is expected to determine the discsize distribution. We also investigate whether or not, and how, the two processes affect the disc masses. To allow for a clear comparison, we do not take external photoevaporation into account. As discussed above, this may be a fair assumption to first order.
2 Methods
Our simulations are carried out within the amuse framework (Portegies Zwart et al., 2013; Pelupessy et al., 2013)^{1}^{1}1http://www.amusecode.org. The amuse framework is written in python and provides a userfriendly way to combine different astrophysical packages and software, for example body integrators and stellar evolution codes. Our setup is practically identical to the setup used by Portegies Zwart (2016) to model the disc size distribution as a result of dynamical encounters in the Trapezium cluster. We perform body simulations of star clusters in a static gas potential. Each star is given a parameterised disc and we follow the evolution of the radii and masses of the discs under the influence of faceon accretion and dynamical encounters using parameterisations for both processes. Each of these two processes is considered independently in our simulations. Below, after summarising the method used to resolve dynamical encounters, we describe our implementation of the faceon accretion model of Wijnen et al. (2017) and the addition of a static gas potential.
We assume the surface density profile of the protoplanetary discs initially follows a powerlaw, , with , corresponding to a minimummass solar nebula (Weidenschilling, 1977; Hayashi, 1981). This allows us to set the inner radius of the disc to 0 AU, instead of a value (110 AU), without realising an unrealistic mass distribution over the disc which would occur for . Following Portegies Zwart (2016), we assume the initial outer radius of all protoplanetary discs is 400 AU and their mass is equal to 10 % of the mass of their host star. We verify that the discs satisfy the Toomre (1964) criterion for stability against selfgravity.
2.1 Resolving dynamical encounters
We use the same method as described in Sects. 2.1 and 2.2 of Portegies Zwart (2016) to resolve dynamical encounters. That is, when stars come within a certain encounter radius (initially 0.02 parsec) from each other, the dynamical encounter is resolved by calculating the pericentre distance, , with Kepler’s equation (using the kepler module from the starlab package; Portegies Zwart et al. 2001). We calculate the new disc radius, , for a star with mass from the estimate found by Breslau et al. (2014) for coplanar, prograde and parabolic encounters^{2}^{2}2Recently, Bhandare et al. (2016) derived a parameterisation by averaging over all inclinations. Their parameterisation is less restrictive because coplanar, prograde encounters have the strongest effect on the disc. We choose the most restrictive prescription to model the strongest influence dynamical encounters can have on protoplanetary discs.:
(1) 
with being the mass of the other star. This parameterisation has been derived using test particles, that is, it does not consider hydrodynamical forces. Previous work shows that excluding pressure, selfgravity, and viscous effects in simulations of discs in dynamical encounters provides results that are consistent with simulations that include these effects (Pfalzner et al., 2005). Only for very close encounters, , should viscous forces be accounted for. It is questionable whether or not a disclike structure survives after such a destructive encounter, in which even stellar capture may occur (Muñoz et al., 2015). For disc radii , the parameterisation may no longer be valid and such small radii must therefore be interpreted with care.
The mass that is lost from the disc during the encounter is calculated using:
(2) 
which follows from the assumption that all material beyond is stripped. The disc of the other star accretes a fraction of the mass that is lost, which we compute as:
(3) 
Eqs. 1, 2, and 3 are all applied symmetrically in a dynamical encounter.
2.2 Parameterising faceon accretion
To parameterise faceon accretion of ambient gas onto protoplanetary discs, we use the theoretical model derived in Wijnen et al. (2017), which has been shown to adequately describe the evolution of the radius, mass, and surface density profile of the disc. This model is an extension of the thin accretion disc theory, assuming that the mass flux onto the disc is uniform over the surface area of the disc and completely accreted by each surface element. In this model, the disc is considered as consisting of concentric rings. The accretion of ISM with no azimuthal angular momentum decreases the specific angular momentum of each ring, which then contracts to a smaller radius that corresponds to its new specific angular momentum. The differential equations for the radial drift velocity and for the surface density of the disc then obtain an extra ISM accretion term with respect to the thin accretion disc theory. When viscous effects are neglected, these differential equations can be simplified to a scalefree theoretical model. This model depends on a dimensionless parameter , which is proportional to the timeintegrated mass flux and is defined as:
(4) 
where is the initial surface density of the disc at an arbitrary but fixed scaling radius . The mass flux onto the disc is determined by the density of the ambient medium and velocity relative to this medium . The factor , where is the inclination between the angular momentum vector of the disc and the flow direction of the ISM, accounts for the situation in which the disc is not aligned perpendicular to the mass flux. Assuming an initial surface density profile, , the surface density of the disc at any moment in time is then given by:
(5) 
Using the parameter , the differential equation for any radius within the disc, , can be expressed in terms of a dimensionless radius via . The differential equation for can be written as (Eq. 16 from Wijnen et al., 2017):
(6) 
We calculate the solution to Eq. 6 for an initial value and call it . The evolution of any other annulus in the disc with initial value can then be derived from:
(7) 
where we have taken advantage of the selfsimilarity of the solutions. Rather than integrating Eq. 6 for every arbitrary starting radius , we obtain its evolution from using Eq. 7. In our simulations, we follow the evolution of the outer radius of each disc, , via its dimensionless equivalent . For each disc, the evolution of is determined by its initial surface density and its timeintegrated mass flux. The solutions to Eq. 7 therefore describe the radius of any disc, at any time, via its individual parameter . In our simulations, the integration in Eq. 4 is replaced by a summation over all time steps. At each time step we use the velocity of the star and the density at its position. The time step we use is defined in Sect. 2.5.
The disc may also be truncated by the ram pressure, , of the ambient medium. We compute the truncation radius, , by solving the equation derived in Sect. 2.1 of Wijnen et al. (2017) but modified to account for the perpendicular velocity component of the flow:
(8) 
At every time step we check if the truncation radius is smaller than the disc radius using the current parameters and solving for the radius with the NewtonRaphson method. Namouni & Guzzo (2007) derived an exact solution for the truncation radius for the purely dynamical problem, in which a particle is ejected from its orbit by a constant force. This solution differs by a factor of in the denominator of Eq. 8. Neither solution takes the gasdynamical and viscous effects in the disc into account. We have verified that Eq. 8 provides a better approximation of the results in the hydrodynamical simulations of Wijnen et al. (2017) and we therefore choose to use this approximation.
Let us suppose a disc is truncated due to ram pressure stripping at a certain moment, , corresponding to via Eq. 4; after which its new radius is . Then, according to Eq. 7, we can write:
(9) 
In order to compute its further evolution with Eq. 7, we must solve Eq. 9 for the unknown new initial value . We do this in the following manner: We define an inverse function, , that for any returns the initial value for which . We have calculated these initial values for by solving the differential Eq. 6 up to the moment for a large number of initial values . This yields unique pairs of and corresponding values. These pairs form an interpolation table, which is our function . We can rewrite Eq. 7 in terms of by replacing with for which and :
(10) 
which is true for any . Using , we can also rewrite Eq. 9 to obtain:
Eqs. 10 and 2.2 can only be equal if and . Therefore, after truncation, we can obtain the new initial from the truncated disc radius and via:
(12) 
Thus when ram pressure truncates the disc, we use Eq. 12 to calculate the new initial value and then use Eq. 7 to compute the further evolution of the disc radius. If the disc is truncated by ram pressure stripping, we also assume all material beyond the truncation radius is removed from the disc as was done in Wijnen et al. (2017) in their comparison with hydrodynamical simulations.
For faceon accretion, we calculate the mass of the disc at any moment by integrating over the surface density profile (Eq. 5):
(13) 
(cf. Eq. 18 from Wijnen et al., 2017). We give the discs a random orientation vector that represents the angular momentum vector of the disc. At every time step, we calculate the absolute value of the dot product between this orientation vector and the velocity vector to obtain , necessary in Eq. 4. We take the absolute value because it does not matter whether the disc is rotating prograde or retrograde. Wijnen et al. (2017b, submitted) have shown that when the rotation axis of the disc is inclined with respect to the velocity vector, the two vectors tend to align as the disc experiences a net torque from the flow (even if the disc is initially symmetric). The mechanism that causes the tilting of the disc is known but the dissipative processes that determine the time scale of the process are not fully understood. The time scale estimate derived by Wijnen et al. (2017b) for this process suggests that a change in the inclination generally takes at least yr for the conditions encountered in our simulations, except in the simulations with the most dense initial conditions. Since this is much longer than the crossing time of the cluster, we do not take this tilt process into account.
The model described above is only valid if the accretion of ISM is sufficiently rapid that viscous effects can be neglected; that is, , where is the accretion time scale and the viscous time scale (see Sect. A.2 of Wijnen et al., 2017). This is valid as long as:
(14) 
where is the initial disc radius (to which the disc was truncated in case of ram pressure stripping). We define the accretion time scale as (cf. Eq. 8 in Wijnen et al., 2017):
(15) 
In the last equality, we have taken the timeaveraged value of , using Eq. 4, instead of the instantaneous value. For the viscous time scale, we use Eq. 7 from Wijnen et al. (2017), adopting the same viscosity parameter and sound speed km/s as in that study. At the end of each simulation we check whether every disc meets the condition in Eq. 14, which is the case in all our simulations.
Over long time scales, the radius and mass of the disc scale as simple powerlaw functions of (Wijnen et al., 2017). For , the disc radius approaches , as follows from Eq. 6. On the same time scale, the mass of the disc is dominated by the second term inside the brackets in Eq. 13 and . We use these relations in the discussion of our results.
2.3 Gas density distribution
We assume the gas in the cluster follows a Plummer profile (Plummer, 1911), which is characterised by a total mass and a characteristic radius . This distribution allows us to model the initial distribution of the stars and gas consistently: the initial positions of the stars are distributed following the same Plummer profile. A different initial stellar distribution, for example a fractal distribution, cannot be combined with a consistent analytic (equilibrium) gas distribution. On the other hand, modelling the gas with a hydrodynamics code makes the simulations computationally expensive and prevents us from exploring a large parameter space. In this work we explore the regimes in which the process of dynamical encounters and faceon accretion are relevant. In future studies, specific parts of the parameter space can be explored in more detail.
As in Portegies Zwart (2016) we simulate our simplified starforming region for 1 Myr, which is motivated by observations of the ONC. The mean stellar age of the ONC is Myr and 80 % of the stars in the Trapezium cluster have an inferred age less than 1 Myr (Prosser et al., 1994; Hillenbrand, 1997). Moreover, young clusters are observed to be embedded during the first 13 Myr of their evolution (Lada & Lada, 2003; Portegies Zwart et al., 2010). A simulation time of 1 Myr allows us to keep the density distribution timeindependent, which reduces the integration error.
2.4 Initial stellar conditions
The stars are randomly drawn from a Kroupa (2001) initial mass function between 0.01 and 100 ; the mean mass of this distribution is . The mass of the protoplanetary disc () is added to the mass of the star in the Nbody integrator, which is the fourthorder hermite Nbody code ph4 (McMillan, publicly available in amuse) and for which we assume a timestep parameter and a softening of 100 AU. The initial positions of the stars are assigned according to a Plummer distribution that follows the gas distribution. This simplified assumption may not be correct as studies have shown that the star formation efficiency increases with gas density (e.g. Burkert & Hartmann, 2013). This effect can be partly taken into account by comparing simulations with different fractions of the total cluster mass in stars. We adopt four different values of this fraction, given in Sect. 2.6. The highest value of could represent either an extreme case of the aforementioned scaling of the star formation efficiency with gas density or a starforming region at the end of its embedded phase. The cluster is set up without primordial mass segregation. The Plummer profile may underestimate both the stellar and gas density in the centre of the cluster but we show in the results section that as long as both are underestimated consistently our results still hold.
2.5 Coupling the Nbody integrator to the gas potential
The Nbody integrator and analytic potential are coupled using the AMUSE Bridge integrator (Fujii et al., 2007; Pelupessy et al., 2013), which uses Hamiltonian splitting and a second order leapfrog integration scheme. We choose a Bridge time step of , where we have assumed with the total cluster mass. This time step is small enough to ensure that the amount of sweptup ISM in one time step never exceeds . The energy conservation in the Bridge system is better than a few parts in , which is sufficient for a reliable result (Portegies Zwart & Boekholt, 2014).
2.6 Varying model parameters
Our simulations are defined by four parameters: The total cluster mass, , the virial radius of the cluster , within which 64 % of is contained and which is related to the characteristic Plummer radius via , the number of stars, , and the virial ratio , defined as , where is the kinetic energy of the stellar component if the stars are in virial equilibrium with the total cluster mass distribution. We checked that in case of virial equilibrium, the radial distribution of the stars remains roughly constant during the simulation. For a given value of we choose four fixed values for such that the fraction of mass in stars (a proxy for the star formation efficiency) is approximately 3, 10, 30 and 90 %. We calculate the total stellar mass, , and the total gas mass is then given by . The realisation of the initial stellar mass function determines and the remaining gas mass therefore varies for a given . Each of these simulations is performed with a virial ratio of and 1.0, representing subvirial, virial, and supervirial initial conditions respectively. Moreover we perform each simulation ten times with a different random seed to reduce the Poisson noise in the final results. We choose small cluster radii such that our simulations can be considered as a local substructure within a larger starforming region. Our simulated clusters are on the dense end of what is generally observed. In the results section it becomes clear that the most important implications of our simulations can be easily extended to clusters that are less dense. We address this in Discussion Sect. 4.6. The initial conditions are listed in Table 1, which also provides labels for several simulations that we discuss in detail. Two of our sets of initial conditions correspond to the early and late evolution of the Trapezium cluster, respectively, and have been labelled according to the assumed mass fraction in stars.
300  [0.1, 0.25, 0.5]  [20, 70, 205, 610]  [0.1, 0.5, 1.0] 

1000  [0.1, 0.25, 0.5]  [70, 230, 680, 2045]  [0.1, 0.5, 1.0] 
3000  [0.1, 0.25, 0.5]  [205, 680, 2045, 6135]  [0.1, 0.5, 1.0] 
Model  

Lowmass (LM)  300  [0.1, 0.25, 0.5]  610  0.5 
Intermediate mass (IM)  1000  0.25  [70, 230, 680, 2045]  0.5 
Highmass (HM)  3000  [0.1, 0.25, 0.5]  2045  0.5 
Trap30  1000  0.25  680  0.5 
Trap90  300  0.25  610  0.5 

3 Results
When we refer to the disc radius and disc mass resulting from faceon accretion, we mean the combined effect of faceon accretion and ram pressure stripping.
3.1 Disc radii
In Sects. 3.1.1 and 3.1.2 we discuss the results for the disc radius in the set of IM simulations in detail. The simulation with and an approximate mass fraction in stars of 30 % can be considered representative for the early evolution of the Trapezium cluster (Lada & Lada, 2003) and we therefore label it Trap30. All four simulations have a stellar density pc meaning that truncation due to dynamical encounters should be relevant.
3.1.1 Characterising the disc radius distribution
Because we compare the effect of dynamical encounters and faceon accretion for a large parameter space of initial conditions, it is convenient to characterise the resulting disc radius distribution by, preferably, one value. In Fig. 1 we show the distribution of disc radii for the Trap30 simulations. For dynamical encounters the distribution is bimodal, because the process is stochastic and 20 % of the stars have not experienced a flyby that was close enough for truncation. This fraction increases to 56 % in the IM simulation with . The stars that do not experience a dynamical encounter reside in the outskirts of the cluster, where the stellar density, and hence the number of encounters, is low. Since the distribution of disc radii resulting from dynamical encounters is bimodal, we choose to use the median disc radius, , instead of the mean disc radius, as a characteristic of the distribution. This provides a better indicator for the influence of dynamical encounters when few discs experience truncation due to a flyby. For faceon accretion, only a small fraction of the stars  less than 6 % in the IM simulations discussed in this section  have a disc radius larger than 380 AU. This fraction does not depend on the number of stars in these simulations. The stars with radii AU reside at large distances from the cluster centre, where not only the stellar density but also the ambient gas density and velocities are low. Fig. 1 shows the standard deviation of the distributions as a function of their median radius. It shows that for each process the standard deviation does not depend strongly on the number of stars. The standard deviation is larger for dynamical encounters due to the skewed distribution of truncated discs. The median disc radius for faceon accretion varies strongly between the ten IM simulations with . The large number of stars drawn from the initial mass function in these simulations leads to large fluctuations in the resulting gas density. We discuss this in the following section.
Based on the results illustrated by Figs. 1 and 1, we conclude that the median disc radius is a good indicator for the effect of each process on the discsize distribution. For dynamical encounters, the median disc radius does not depend on the initial disc radii, that is, the tail of the distribution. We therefore use the median disc radius in the remainder of this work to quantify the effects of the two processes.
3.1.2 Comparing stellar mass fractions
Because the total cluster mass is set, the main difference between simulations with the same initial conditions but different stellar mass fractions is the amount of mass in gas and in stars. The velocity dispersion of the stars isa determined by the total cluster mass and it is therefore practically the same for the IM simulations discussed in this section. Furthermore, the distribution of stars throughout the cluster also remains roughly the same because we assumed the initial cluster to be in virial equilibrium. The discriminating parameters between the simulations with different numbers of stars are the gas and stellar density. We therefore plot the mean disc radius versus these two quantities in Figs. 2 and 2, respectively. We determine the stellar density by counting all stars within the cluster radius of 0.25 pc at the end of the simulation and calculate the average gas density within the cluster radius analytically. By definition, the gas density decreases as the number of stars increases. There are some notable outliers in Fig. 2, particularly for the highest number of stars. These outliers are caused by the sampling of the initial stellar mass function with a fixed number of stars, that continues beyond the intended if is not yet reached. The gas that remains to form stars varies between 5 and 240 among simulations that have the same initial conditions otherwise.
The median disc radius resulting from faceon accretion decreases with increasing gas density. This is expected theoretically because, as discussed at the end of Sect. 2.2, the longterm evolution of an individual disc radius scales as . Of the quantities that determine the value of in Eq. 4, the typical velocities and surface densities are similar in each simulation. Averaged over the whole population, only the mean gas density varies between IM simulations with different , so that we expect the approximate relation . In Fig. 2 the curve of the median disc radii resulting from faceon accretion follows a slightly steeper trend with density, approximately . The difference is due to ram pressure stripping becoming more prominent at higher gas densities and causes additional truncation.
On the other hand, when only dynamical encounters are accounted for, the median disc radius increases strongly with the ambient gas density, that is, with decreasing stellar density as can be seen in Fig. 2. In simulations with the same total mass and radius, the change in stellar mass fraction affects the number of encounters only via the stellar density, because the velocity dispersion remains almost the same. At a stellar density of 1000 pc the effect of dynamical encounters is minor, AU, and at a stellar density of pc the effect of faceon accretion on the disc sizes is still larger than the effect of dynamical encounters. Only for a stellar mass fraction of 90 % and corresponding stellar density , are dynamical encounters the dominant disctruncation process. Even at the lowest gas densities, faceon accretion reduces the averaged median disc radius substantially, to 133 AU. This average is determined by three outliers and the rest of the simulations with the same initial conditions indicate that the effect of both processes is comparable. Due to their low gasmass fractions, these three outlying clusters are more dynamically evolved than the other seven clusters with the same parameters, as can be inferred from the low stellar density within the cluster radius at the end of the simulation.
These results show that the effect of dynamical encounters relative to faceon accretion becomes important once the gas is (partly) expelled or star formation is very efficient.
3.1.3 Disc radii as a function of cluster radius and mass
In Fig. 3 we show the median disc radius as a function of the stellar density for our HM simulations with three different cluster radii. For a stellar mass fraction of 30 %, these simulations have the highest stellar density at each cluster radius. Fig. 3 shows that even at stellar densities , the process of faceon accretion leads to smaller disc sizes than dynamical encounters. As the stellar density increases with decreasing cluster radius, both the gas density and velocity dispersion increase as well, which enhances the effect of faceon accretion and ram pressure stripping. The highest stellar density in Fig. 3 corresponds to a gas density of . This is at the high end of what may be expected in the cores of starforming regions if the star formation efficiency increases with gas density. These conditions illustrate that as long as the gas to stellar mass ratio is , faceon accretion is the dominant disctruncation process regardless of the stellar density. Only by expelling gas or in case of more efficient star formation can dynamical encounters become the dominant process if the stellar density is sufficiently high. We find that for stellar densities pc, dynamical encounters do not alter the disc radii, as was also shown in Rosotti et al. (2014) and Vincke et al. (2015).
In Fig. 3 we show the results of our LM simulations with three different cluster radii. These simulations correspond to a stellar mass fraction of about 90 %. The scatter in the gas density between simulations with the same initial conditions is caused by the sampling of the initial mass function and the resulting variation in the remaining gas mass (as discussed in Sect. 3.1.2). The simulations with pc can be interpreted as representing either the end of the embedded phase of the Trapezium cluster or an efficiently formed Trapezium cluster. In these Trapezium analogue simulations, labelled Trap90, the stellar density is . The contributions of faceon accretion and dynamical encounters to the decrease in disc radii are about equal in this case. In the most dense conditions, that is, and a stellar mass fraction of 90 %, dynamical encounters are mildly dominant, but the contributions of both processes are still very similar and severe: AU. If dynamical encounters result in a median disc radius of AU, they are very destructive and the parameterisation we use for their radii may no longer apply. The LM simulations with pc in Fig. 3 show that even at gas densities the process of faceon accretion is still relevant, because AU.
Our simulations show that the stellar density has to exceed a few times and simultaneously the stellar mass fraction has to be well above 30 % for dynamical encounters to be at least equally important as faceon accretion in the disc truncation process. Even if dynamical encounters are the dominant process, at stellar densities exceeding and a stellar mass fraction of 90 %, the influence of faceon accretion is still important. In these compact clusters, the corresponding velocity dispersion also enhances the effects of faceon accretion and ram pressure stripping. Generally, dynamical encounters become more important than faceon accretion in the regime where flybys are very destructive, resulting in .
3.2 Disc masses
In Fig. 4 we show the distribution of disc masses in our Trap30 simulations. The net effect of a dynamical encounter is that the disc loses more mass than it accretes from the disc of the other star. In the process of faceon accretion, on the other hand, the discs generally gain mass as more mass is swept up than is stripped by the ram pressure. The distribution of disc masses resulting from dynamical encounters is therefore expected to be shifted to lower masses compared to the initial distribution and for faceon accretion the distribution should be shifted towards higher masses. This is demonstrated in Fig. 4. In all our simulations the median disc mass resulting from faceon accretion is higher than the median disc mass produced by dynamical encounters, as we discuss in Sect. 3.2.1. As for the radius distribution, we choose the median disc mass, , to characterise the disc mass distribution. The median disc mass provides a better correspondence to the peak in the distribution than the mean disc mass, because the latter is determined by the high end of the disc mass spectrum.
3.2.1 Disc masses as a function of the parameter space
Figs. 5 and 5 show the median disc mass as a function of the gas and stellar density, respectively, for our IM simulations. These results correspond to the median disc radii shown in Figs. 2 and 2. The median disc masses resulting from faceon accretion increase slightly with increasing gas density. Applying similar reasoning as for the median disc radius in Sect. 3.1.2, the median disc mass is expected to scale as , again assuming that averaged over the population only the mean gas density varies between IM simulations with a different number of stars. When the curve for the median disc masses in Fig. 5 is fitted as a powerlaw of , this gives an exponent of 0.07. The increase in the median disc mass with density is slower than theoretically expected, because ram pressure stripping is more effective at high gas density, removing mass from the discs and decreasing the median disc mass. For dynamical encounters, the median disc mass decreases with increasing stellar density, which is also expected theoretically. Yet the increasing stellar density has only a very small effect on the median disc mass once the stellar density .
Figs. 6 and 6 show the median disc masses for the same simulations as shown in Figs. 3 and 3. These figures illustrate that for either process, the cluster parameters have a small influence on the median disc mass. In the HM simulations (Fig 6), the median disc mass decreases with decreasing cluster radius for both processes because both the ram pressure and stellar density increase as the cluster radius decreases. In the HM simulations with pc, ram pressure stripping is so strong that the median disc mass resulting from faceon accretion is 1.5 % lower than the initial median disc mass. In this particular case, the mass loss due to ram pressure stripping cannot be compensated by the mass gain from faceon accretion. In the LM simulations in Fig. 6, the ram pressure is much lower and the median disc mass resulting from faceon accretion increases for more compact clusters that have a higher density and velocity dispersion.
In general, we find that dynamical encounters lead to disc masses that are lower than those predicted by faceon accretion. Faceon accretion produces compact discs with a high surface density, while dynamical encounters generally result in larger and less massive discs.
3.3 Virial ratios
In the simulations discussed in the previous section, we assume that the stars and gas are initially in virial equilibrium. We also performed simulations with ‘cold’, that is, , and ‘warm’, , initial conditions. The simulations that start dynamically cold go through a phase of contraction. This increases the stellar density and velocity dispersion. On the other hand, the simulations that start dynamically warm expand from the beginning, which causes the stellar density to decrease. In Fig. 7 we show the mean disc radii and masses in our HM simulations with pc for the three virial ratios. For the stellar density is substantially lower at the end of the simulations and as a consequence the median disc radius resulting from dynamical encounters is larger than in the simulations in virial equilibrium. A virial ratio of leads to a higher stellar density at the end of the simulation and hence the median disc radius due to dynamical encounters is smaller. The median disc mass for dynamical encounters follows the same trend as a function of virial ratio.
For faceon accretion, the median disc radius does not depend strongly on the virial ratio but there is a slight decreasing trend with decreasing virial ratio, that is, with increasing velocity dispersion. Simultaneously, the median disc mass increases with decreasing virial ratio. If the cluster contracts, the discs generally move through regions with higher density with a higher velocity compared to a cluster that expands. The discs are therefore able to sweep up more gas when the initial conditions are cold.
In essence, if the initial conditions are supervirial (, then faceon accretion is the dominant process in truncating the discs regardless of the other cluster characteristics we use in this work. Faceon accretion is also the dominant process if the initial conditions are subvirial and the stellar mass fraction is %. Dynamical encounters become dominant if the stellar mass fraction is well above 30 % and the cluster goes through a phase of contraction.
4 Discussion
4.1 Assumptions on accretion and mass loss
We kept the mass of the stars and the gas potential constant in our simulations. Here we discuss the consistency of our simulations with this and our other assumptions on accretion and mass loss.
We assume that the total gas mass in the cluster remains constant in our simulations. In most simulations the net effect of faceon accretion does not change the total gas mass by %, except in the simulations with the highest stellarmass fraction. At most 20.9 % of the gas mass is accreted in one of these simulations (with , pc, and ). Dynamical encounters are the dominant cause of disc truncation in this regime, even if the effect of faceon accretion may have been overestimated due to the presumed constant gas density distribution. We therefore conclude that our results are not affected by the assumption that the amount of gas remains constant in our simulations.
Furthermore, the masses of the particles in the body simulation remain constant during the simulation. Except for the most compact simulations with the lowest stellar mass fraction, the fraction of stars, including their discs, that gain more than 10 % in mass is always less than 1 % in our simulations that start in virial equilibrium. This small fraction of stars is similar to what is found by Throop & Bally (2008) who simulated BondiHoyle accretion for 4 Myr with similar cluster parameters and a star formation efficiency of 33 %. In the compact virialised simulations with the lowest stellar mass fraction, at most 11.6 % of the stars gained more than 10 % in mass in the simulations with , pc, and . This fraction increases to a maximum of 38.7 % in the initially subvirial simulation with , pc, . A stellar mass fraction of roughly 3 % is probably not realistic in such dense conditions (e.g. Bonnell et al., 2003; Burkert & Hartmann, 2013) and we have therefore not discussed these simulations in detail. We merely use them to investigate how either process scales as a function of the initial conditions.
We do not account for BondiHoyle accretion in our simulations. Previous studies have modelled BondiHoyle accretion onto stars in clusters (e.g. Throop & Bally, 2008; Scicluna et al., 2014; BallesterosParedes et al., 2015). They find that the amount of accreted material scales with and for solarmass stars it is typically not more than a few per cent of their initial mass, depending on the gas density and velocity dispersion. If the BondiHoyle radius of a star is larger than its disc radius, the gravitational focusing of gas by the star enhances the effective surface area for accretion onto the disc. This increases the amount of gas that is swept up by the disc. We tested this effect on our faceon accretion model by correcting the mass flux for the gravitational focusing if the BondiHoyle radius is larger than the disc radius (BisnovatyiKogan et al., 1979; Edgar, 2004) but this did not change our distribution of disc radii and masses. Generally, the BondiHoyle radius is smaller than the disc radius in our simulations. In case the BondiHoyle radius is larger, the disc initially contracts faster and the net effect on the disc radius and mass is similar after 1 Myr.
The amount of mass lost in a dynamical encounter is probably overestimated by Eq. 2. Breslau et al. (2014) show that the surface density profile is affected by a dynamical encounter. Their Fig. 1b shows that the surface density within the cutoff radius can become higher than the initial surface density profile. Unfortunately, the effect of dynamical encounters on the surface density profile cannot be easily parameterised. For encounters between debris discs, Jílková et al. (2016) showed that the minimum disc radius beyond which particles can be unbound, , is proportional to the encounter pericenter and the mass of the accreting star divided by the total mass of both stars, in case of lowinclination and loweccentricity encounters. However, not all material beyond this radius is lost from the disc. Their parameterisation for is similar to what we use for the new disc radius in Eq. 1. The redundant mass lost from the disc in our simulations is partly compensated by accretion of material that is lost from the other disc in the encounter. As long as there is net mass loss during an encounter, the relative influence of dynamical encounters and faceon accretion on the median disc mass does not change in our simulations.
4.2 Viscous forces
The viscous spreading of a protoplanetary disc makes it more vulnerable to dynamical encounters and ram pressure stripping. When a disc moves through an ambient medium, the disc size is determined by the process of faceon accretion and ram pressure stripping at the outer edge, that is, the disc cannot spread beyond a certain radius (Wijnen et al., 2016, 2017). For faceon accretion, neglecting viscous spreading is therefore a fair assumption, as long as the accretion time scale is shorter than the viscous time scale. We verified that for the accretion time scales in our simulations, this criterion is satisfied. In case of dynamical encounters, the increasing crosssection of the disc resulting from viscous spreading would increase the probability of a truncating encounter. However, during the embedded phase, the viscous spreading would be inhibited as outlined above. At the end of the embedded phase, viscous spreading may become important and might increase the probability of dynamical encounters. Our simulations with a stellar mass fraction of 90 % and stellar densities show that, by that time, the process of dynamical encounters is dominant, without accounting for the viscous spreading. Including viscous spreading at this stage may ease the constraint on the stellar density. However, initial results by Concha Ramírez & Vaher (in prep.) find that gas expulsion leaves the cluster in a supervirial state and including viscous spreading of the disc (according to the selfsimilarity solutions of LyndenBell & Pringle, 1974) does not increase the number of truncations in an expanding cluster.
4.3 Drag forces
As the discs move through an ambient medium they experience a drag force exerted by the ram pressure, which we have not taken into account. The part of the disc that is tightly bound to the star, which is determined by in Eq. 8, can decelerate the star. We can estimate the time scale on which the drag force operates as . For an star in our Trap30 simulation, assuming an average disc radius of 200 AU, this corresponds to a velocity change of 5 % after 1 Myr. This time scale decreases with stellar mass but less massive discs contract faster. The average effect on the population is probably not larger than this 5 %, also considering that the disc of an is expected to be AU within 0.2 Myr in the centre of the Trap30 simulation.
4.4 Initial disc masses and radii
We have assumed relatively high disc masses of . The masses of protoplanetary discs are observed to be of the order of (Andrews & Williams, 2005, 2007). These estimates may be off by an order of magnitude but, in general, disc masses are probably not much higher than . Assuming lower disc masses in our simulations would not affect the relative influence of dynamical encounters on disc sizes in this study (cf. Eq. 1). A similar reasoning applies to the initial disc radii: Although the probability of experiencing a dynamical encounter depends on the initial disc radius, the resulting disc radius does not. Assuming initially larger disc radii would affect the tail of our distribution but not the median disc radius (e.g. Fig. 2 in Vincke et al., 2015).
In case of faceon accretion, distributing a lower disc mass over the same surface area, or the same disc mass over a larger surface area, decreases the surface density of the disc. As a result, the disc is more vulnerable to ram pressure stripping and the disc contracts faster. If all other conditions remain the same, then over long time scales the disc radius scales as , where and are the initial disc mass and radius (Sect. 2.2). For a given disc mass, starting with an initial radius of 1000 AU instead of 400 AU thus causes the disc to contract about twice as fast. The same reasoning applies to assuming lower disc masses, except that the scaling with disc mass in the above reasoning is less strong. The truncation radius due to ram pressure stripping also decreases with decreasing initial surface density, although the dependence is weaker. The relative importance of faceon accretion with respect to dynamical encounters would therefore increase if the initial disc radii were larger and/or disc masses were smaller.
4.5 Photoevaporation and winds
As discussed in the introduction, neglecting the effect of external photoevaporation on the discs on a time scale of 1 Myr in an embedded cluster is probably a fair assumption. Winds from massive stars also affect the ambient gas density. They can locally clear the gas and create lowdensity cavities. The density in these regions can be orders of magnitude lower than the average gas density in the cluster. The wind velocity is three orders of magnitude higher than the velocity dispersion of the clusters in our simulations. Although the integrated mass flux on the disc would be similar, the high temperature difference between these ejecta and the disc may hamper efficient accretion. Twodimensional simulations of supernova ejecta interacting with a protoplanetary disc show that only 1 % of the hot gaseous ejecta is intercepted by the disc, contrary to cold gas and dust that is accreted very efficiently onto a protoplanetary disc (Ouellette et al., 2007, 2009, 2010; Wijnen et al., 2016, 2017). Hence the integrated mass flux inside a wind cavity could be smaller by a factor of 100. The time scale on which these winds create cavities depends on the most massive star in the cluster and thus on the realisation of the initial mass function. Simulations of embedded clusters with that take feedback into account find that the impact of winds from massive stars on the gas distribution becomes relevant on a time scale of 1 Myr (Pelupessy et al., 2012), which is the duration of our simulations. The same authors find that supernovae are not expected to play a role until an age of roughly 10 Myr. The contraction of the disc caused by faceon accretion is faster in the beginning when stellar feedback may not have affected the gas distribution substantially yet. The winds could also strip the discs from their outer regions because of their high velocities.
We also neglect the influence of photoevaporation by the host star on both processes. Dynamical encounters are not affected by this process but the photoevaporative winds from the disc and jets may influence the faceon accretion process. Interaction between inflowing ambient gas and material from the star itself or evaporated disc material could slow down faceon accretion onto the disc. The winds from the star and the disc depend on the mass of the star and its spectral energy distribution. Extremeultraviolet radiation irradiates the inner few AU of the protoplanetary disc, while Xrays and farultraviolet radiation induce photoevaporation at disc radii of tens of AU and AU, respectively (e.g. Armitage, 2011; Williams & Cieza, 2011). Irradiation by the central star may become dominant when the disc has shrunk to radii AU. Photoevaporation can further erode the disc from the inside out on time scales of yr after disc lifetimes of several Myr (we refer to e.g. Alexander et al., 2014, and references therein). The increased surface density due to faceon accretion could make the disc more resistant against this erosion. Currently, these processes cannot be included in our models but they should be accounted for in future studies.
4.6 Observational constraints
The process of dynamical encounters predicts a larger spread in observed disc radii than faceon accretion. For example, in the Trap30 simulations, the radius distribution resulting from faceon accretion is concentrated around a clear peak (cf. Fig. 1), whereas the distribution produced by dynamical encounters is much broader. In Fig. 8 we compare the radius distributions of the Trap30 simulations to observations of disc radii in the ONC. Vicente & Alves (2005) assembled a sample of 149 protoplanetary discs out of 300 young stellar objects, which they claim is complete down to a disc radius of 50 AU. To compare to the observations, we only show the simulation results for stars with a mass , which roughly corresponds to the latest spectral type of the stars to which the discs could be matched by Vicente & Alves (2005). The distribution resulting from faceon accretion is in better agreement with the observed distributions. Based on the resolution limit of their data, Vicente & Alves (2005) estimate that 40 to 45 % of the discs in the Trapezium cluster have radii larger than 50 AU. This agrees with what we find for the whole simulated population in Fig. 1. Excluding the very lowmass stars () increases the median disc radius for both processes by roughly 40 AU, because the discs of these stars generally have small radii. This is also demonstrated by the absence of the bin with the smallest radii for faceon accretion in Fig. 8, when compared to Fig. 1. The relative influence of both processes is not affected by the inclusion or exclusion of these stars in our analysis.
The clusters simulated in this work are generally more dense than the observed band of young stellar clusters in the massradius plane (e.g. Pfalzner et al., 2016) and radius plane (e.g. Kuhn et al., 2015). We have shown in this work that lowering the stellar and gas density consistently, that is, decreasing the cluster mass for a given radius or increasing the cluster radius for a given mass, increases the influence of faceon accretion relative to dynamical encounters. The effect of faceon accretion can therefore be expected to be even stronger with respect to dynamical encounters in the parameter space covered by the majority of observed embedded clusters.
The initial conditions in this work are idealised but our simulations illustrate that the relative importance of both processes can be constrained by observations of disc radii. Furthermore, the combination of such observations with observations of the stellar and ambient gas density might put constraints on the duration of the embedded phase of clusters and their star formation efficiencies. A comparison based on disc masses is less straightforward. As discussed in Sect. 4.4, disc masses can only be estimated from observations to within a factor of 10. According to our simulations this is not enough to distinguish between either process based on the resulting disc masses.
4.7 Combined effect of dynamical encounters and faceon accretion
The prescriptions for dynamical encounters and faceon accretion cannot be coupled consistently, because it is not evident how the mass that is accreted in a dynamical encounter should be accounted for in the faceon accretion model. We performed test simulations that included the combined effect of both processes by neglecting the mass accreted in a dynamical encounter. These simulations show that when faceon accretion is the dominant process, the additional truncation caused by dynamical encounters on the disc radius is marginal, unless the influence of both processes is comparable. On the other hand, if dynamical encounters are the dominant truncation mechanism then faceon accretion further decreases the disc size.
The median disc mass that results from combining both processes is always intermediate between the median disc masses predicted by either process independently, whether the mass that is accreted in a dynamical encounter is added at the end of the simulation or not. The value of the median disc mass resulting from the combined effect of both processes leans towards the value predicted by the process that dominates the truncation of the disc radii.
5 Conclusions
By including the processes of faceon accretion and dynamical encounters in body simulations of embedded clusters, we find that faceon accretion, including the effect of ram pressure stripping, is dominant in truncating protoplanetary discs if the fraction of mass in stars is , regardless of other cluster parameters. The stellar mass fraction has to be well above 30 % and simultaneously the stellar density has to exceed a few times pc in order for dynamical encounters to have a comparable effect on the truncation of the discs. Even at stellar densities pc and for gas mass fractions , faceon accretion has a comparable effect to dynamical encounters, but in these circumstances both processes are destructive, resulting in median disc radii AU. We confirm the results from previous studies that dynamical encounters require stellar densities pc to be effective at all (Rosotti et al., 2014; Vincke et al., 2015). This implies that dynamical encounters only become relevant compared to faceon accretion either (1) in clusters with stellar densities pc and extremely high stellar mass fractions, or (2) at the end of the embedded phase of the cluster if the stellar density is similar and discs have not yet shrunk sufficiently due to faceon accretion. Our modelling of the faceon accretion process does not take viscous spreading of the disc or photoevaporative winds from the star and disc into account, which may both reduce the contraction of the disc due to faceon accretion. On the other hand, we have assumed massive discs with , while lowermass discs that are subject to faceon accretion contract faster and are more vulnerable to ram pressure stripping. Of the available prescriptions for the effect of dynamical encounters on the disc radius, we use the one that results in the smallest radius.
Faceon accretion leads to discs that are compact and have a relatively high surface density. On the other hand, dynamical encounters generally result in larger and less massive discs.
Acknowledgements.
We are thankful to Lucie Jílková, Francisca Concha Ramírez, Eero Vaher and Vincent HénaultBrunet for valuable discussions. We also thank the referee for his/her comments. This research is funded by the Netherlands Organisation for Scientific Research (NWO) under grant 614.001.202.References
 Adams (2010) Adams, F. C. 2010, Annual Review of Astronomy and Astrophysics, 48, 47
 Alexander et al. (2014) Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, Protostars and Planets VI, 475
 Alexander et al. (2006a) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006a, Monthly Notices of the Royal Astronomical Society, 369, 216
 Alexander et al. (2006b) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006b, Monthly Notices of the Royal Astronomical Society, 369, 229
 Andrews & Williams (2005) Andrews, S. M. & Williams, J. P. 2005, The Astrophysical Journal, 631, 1134
 Andrews & Williams (2007) Andrews, S. M. & Williams, J. P. 2007, The Astrophysical Journal, 671, 1800
 Armitage (2011) Armitage, P. J. 2011, Annual Review of Astronomy and Astrophysics, 49, 195
 BallesterosParedes et al. (2015) BallesterosParedes, J., Hartmann, L. W., PerezGoytia, N., & Kuznetsova, A. 2015, arXiv:1506.02591 [astroph], arXiv: 1506.02591
 Bally & Scoville (1982) Bally, J. & Scoville, N. Z. 1982, The Astrophysical Journal, 255, 497
 Balog et al. (2007) Balog, Z., Muzerolle, J., Rieke, G. H., et al. 2007, The Astrophysical Journal, 660, 1532
 Bhandare et al. (2016) Bhandare, A., Breslau, A., & Pfalzner, S. 2016, Astronomy and Astrophysics, 594, A53
 BisnovatyiKogan et al. (1979) BisnovatyiKogan, G. S., Kazhdan, Y. M., Klypin, A. A., Lutskii, A. E., & Shakura, N. I. 1979, Soviet Astronomy, 23, 201
 Bonnell et al. (2003) Bonnell, I. A., Bate, M. R., & Vine, S. G. 2003, Monthly Notices of the Royal Astronomical Society, 343, 413
 Breslau et al. (2014) Breslau, A., Steinhausen, M., Vincke, K., & Pfalzner, S. 2014, Astronomy and Astrophysics, 565, A130
 Burkert & Hartmann (2013) Burkert, A. & Hartmann, L. 2013, The Astrophysical Journal, 773, 48
 Clarke et al. (2001) Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, Monthly Notices of the Royal Astronomical Society, 328, 485
 Clarke & Pringle (1993) Clarke, C. J. & Pringle, J. E. 1993, Monthly Notices of the Royal Astronomical Society, 261, 190
 Cloutier et al. (2014) Cloutier, R., Currie, T., Rieke, G. H., et al. 2014, The Astrophysical Journal, 796, 127
 de Juan Ovelar et al. (2012) de Juan Ovelar, M., Kruijssen, J. M. D., Bressert, E., et al. 2012, Astronomy & Astrophysics, 546, L1
 Dullemond et al. (2006) Dullemond, C. P., Hollenbach, D., Kamp, I., & D’Alessio, P. 2006, arXiv preprint astroph/0602619
 Edgar (2004) Edgar, R. 2004, New Astronomy Reviews, 48, 843
 Facchini et al. (2016) Facchini, S., Clarke, C. J., & Bisbas, T. G. 2016, Monthly Notices of the Royal Astronomical Society, 457, 3593
 Fang et al. (2013) Fang, M., Boekel, R. v., Bouwman, J., et al. 2013, A&A, 549, A15
 Fang et al. (2012) Fang, M., van Boekel, R., King, R. R., et al. 2012, Astronomy and Astrophysics, 539, A119
 Fujii et al. (2007) Fujii, M., Iwasawa, M., Funato, Y., & Makino, J. 2007, Publications of the Astronomical Society of Japan, 59, 1095
 Guarcello et al. (2009) Guarcello, M. G., Micela, G., Damiani, F., et al. 2009, Astronomy and Astrophysics, 496, 453
 Guarcello et al. (2007) Guarcello, M. G., Prisinzano, L., Micela, G., et al. 2007, Astronomy and Astrophysics, 462, 245
 Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, The Astrophysical Journal, 495, 385
 Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
 Heller (1995) Heller, C. H. 1995, The Astrophysical Journal, 455, 252
 Hillenbrand (1997) Hillenbrand, L. A. 1997, The Astronomical Journal, 113, 1733
 Jílková et al. (2016) Jílková, L., Hamers, A. S., Hammer, M., & Portegies Zwart, S. 2016, Monthly Notices of the Royal Astronomical Society, 457, 4218
 Kobayashi & Ida (2001) Kobayashi, H. & Ida, S. 2001, Icarus, 153, 416
 Kroupa (2001) Kroupa, P. 2001, Monthly Notices of the Royal Astronomical Society, 322, 231
 Kuhn et al. (2015) Kuhn, M. A., Feigelson, E. D., Getman, K. V., et al. 2015, The Astrophysical Journal, 812, 131
 Lada & Lada (2003) Lada, C. J. & Lada, E. A. 2003, Annual Review of Astronomy and Astrophysics, 41, 57
 LyndenBell & Pringle (1974) LyndenBell, D. & Pringle, J. E. 1974, Monthly Notices of the Royal Astronomical Society, 168, 603
 Mann & Williams (2010) Mann, R. K. & Williams, J. P. 2010, The Astrophysical Journal, 725, 430
 Muñoz et al. (2015) Muñoz, D. J., Kratter, K., Vogelsberger, M., Hernquist, L., & Springel, V. 2015, Monthly Notices of the Royal Astronomical Society, 446, 2010
 Namouni & Guzzo (2007) Namouni, F. & Guzzo, M. 2007, Celestial Mech Dyn Astr, 99, 31
 Olczak et al. (2006) Olczak, C., Pfalzner, S., & Spurzem, R. 2006, The Astrophysical Journal, 642, 1140
 Ostriker (1994) Ostriker, E. C. 1994, The Astrophysical Journal, 424, 292
 Ouellette et al. (2009) Ouellette, N., Desch, S. J., Bizzarro, M., et al. 2009, Geochimica et Cosmochimica Acta, 73, 4946
 Ouellette et al. (2007) Ouellette, N., Desch, S. J., & Hester, J. J. 2007, The Astrophysical Journal, 662, 1268
 Ouellette et al. (2010) Ouellette, N., Desch, S. J., & Hester, J. J. 2010, ApJ, 711, 597
 Pelupessy et al. (2012) Pelupessy, F. I., Jänes, J., & Portegies Zwart, S. 2012, New Astronomy, 17, 711
 Pelupessy et al. (2013) Pelupessy, F. I., van Elteren, A., de Vries, N., et al. 2013, Astronomy and Astrophysics, 557, 84
 Pfalzner et al. (2016) Pfalzner, S., Kirk, H., Sills, A., et al. 2016, Astronomy and Astrophysics, 586, A68
 Pfalzner et al. (2005) Pfalzner, S., Umbreit, S., & Henning, T. 2005, The Astrophysical Journal, 629, 526
 Plummer (1911) Plummer, H. C. 1911, Monthly Notices of the Royal Astronomical Society, 71, 460
 Portegies Zwart & Boekholt (2014) Portegies Zwart, S. & Boekholt, T. 2014, The Astrophysical Journal Letters, 785, L3
 Portegies Zwart (2016) Portegies Zwart, S. F. 2016, Monthly Notices of the Royal Astronomical Society, 457, 313
 Portegies Zwart et al. (2010) Portegies Zwart, S. F., McMillan, S. L., & Gieles, M. 2010, Annual Review of Astronomy and Astrophysics, 48, 431
 Portegies Zwart et al. (2001) Portegies Zwart, S. F., McMillan, S. L. W., Hut, P., & Makino, J. 2001, Monthly Notices of the Royal Astronomical Society, 321, 199
 Portegies Zwart et al. (2013) Portegies Zwart, S. F., McMillan, S. L. W., van Elteren, A., Pelupessy, F. I., & de Vries, N. 2013, Computer Physics Communications, 184, 456
 Prosser et al. (1994) Prosser, C. F., Stauffer, J. R., Hartmann, L., et al. 1994, The Astrophysical Journal, 421, 517
 Richert et al. (2015) Richert, A. J. W., Feigelson, E. D., Getman, K. V., & Kuhn, M. A. 2015, The Astrophysical Journal, 811, 10
 Rosotti et al. (2014) Rosotti, G. P., Dale, J. E., de Juan Ovelar, M., et al. 2014, Monthly Notices of the Royal Astronomical Society, 441, 2094
 Scally & Clarke (2001) Scally, A. & Clarke, C. 2001, Monthly Notices of the Royal Astronomical Society, 325, 449
 Scicluna et al. (2014) Scicluna, P., Rosotti, G., Dale, J. E., & Testi, L. 2014, Astronomy and Astrophysics, 566, L3
 Shu et al. (1993) Shu, F. H., Johnstone, D., & Hollenbach, D. 1993, Icarus, 106, 92
 Throop & Bally (2008) Throop, H. B. & Bally, J. 2008, The Astronomical Journal, 135, 2380
 Toomre (1964) Toomre, A. 1964, The Astrophysical Journal, 139, 1217
 Vicente & Alves (2005) Vicente, S. M. & Alves, J. 2005, Astronomy and Astrophysics, 441, 195
 Vincke et al. (2015) Vincke, K., Breslau, A., & Pfalzner, S. 2015, Astronomy and Astrophysics, 577, A115
 Vincke & Pfalzner (2016) Vincke, K. & Pfalzner, S. 2016, The Astrophysical Journal, 828, 48
 Weidenschilling (1977) Weidenschilling, S. J. 1977, Monthly Notices of the Royal Astronomical Society, 180, 57
 Wijnen et al. (2016) Wijnen, T. P. G., Pols, O. R., Pelupessy, F. I., & Zwart, S. P. 2016, A&A, 594, A30
 Wijnen et al. (2017) Wijnen, T. P. G., Pols, O. R., Pelupessy, F. I., & Zwart, S. P. 2017, A&A, 602, A52
 Williams & Cieza (2011) Williams, J. P. & Cieza, L. A. 2011, Annual Review of Astronomy and Astrophysics, 49, 67