EDP Sciences
Free Access
Issue
A&A
Volume 594, October 2016
Article Number A30
Number of page(s) 14
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201527886
Published online 12 October 2016

© ESO, 2016

1. Introduction

Stars generally form in clusters (Lada & Lada 2003). These dense environments affect the formation and evolution of the stars they host. For example, globular clusters (GCs) were once considered the archetype of coeval, simple stellar populations, but during the last two decades they have been shown to harbour multiple stellar populations. Observations imply that a considerable fraction (up to 70%) of the stars currently in GCs have a very different chemical composition from the initial population (see e.g. Gratton et al. 2012). They indicate that a second (and in some cases even higher order, e.g. Piotto et al. 2007) population1 of stars has formed from material enriched by ejecta from first generation stars.

To explain the formation of these enriched stellar populations, several scenarios have been proposed (see e.g. Decressin et al. 2007; D’Ercole et al. 2008; de Mink et al. 2009; Bastian et al. 2013b). One of the recently proposed scenarios applies in particular to star formation in GCs, but could, in theory, also leave its signature on stellar systems that are less dense. Bastian et al. (2013b; hereafter B13), have suggested a scenario in which the enriched population is not formed by a second star formation event, but rather by the accretion of enriched material, which was expelled by the high-mass stars of the initial population, onto the low-mass stars of the same (initial) population. Because Bondi-Hoyle accretion, i.e. gravitational focusing of material onto the star, is unlikely to be efficient in a GC environment with a high velocity dispersion, they suggest that the protoplanetary discs of low-mass stars sweep up the enriched matter. To account for the observed abundances in the enriched population, the low-mass stars have to accrete of the order of their own mass, i.e. a 0.25 M star has to accrete about 0.25 M of enriched material in the most extreme case (as inferred from, for example, the main sequence of NGC 2808 (Piotto et al. 2007)). The timescale of this scenario is limited by the lifetimes and sizes of protoplanetary discs. B13 assume that the protoplanetary discs can accrete material for up to 20 Myr. Current disc lifetimes are observed as being 5–15 Myr, but B13 argue that their lifetimes may have been considerably longer in GCs. The accretion rate averaged over 20 Myr therefore has to be about 10-8M/ yr. In their scenario, they assume that the accretion rate is proportional to the size of the disc, , density of the interstellar medium (ISM), ρISM, and the velocity, vISM, of the disc with respect to the ISM, i.e. . Furthermore, they assume an average and constant disc radius of 100 AU. In this work, we test several of these assumptions of the early disc accretion scenario.

A similar scenario has been studied before by Moeckel & Throop (2009; M09 hereafter). They performed smoothed particle hydrodynamics simulations of a protoplanetary disc that is embedded in a flow of ISM with a velocity of 3 km s-1. They found that the mean accretion rate onto the star equals the rate expected from Bondi-Hoyle theory, whether a disc is present or not. We note that for the parameters they assumed, the theoretical Bondi-Hoyle radius is almost twice the radius of the disc. Here we follow up on the work by M09 by simulating the accretion process onto a protoplanetary disc for the typical conditions expected in a GC environment. We directly compare the outcome of two different smoothed particle hydrodynamics codes for the same set of initial conditions and different particle resolutions. We first discuss both the physical and numerical effects in our reference model and, subsequently, we compare the results of the different codes and particle numbers.

2. Expected physical effects

In dense stellar systems, where both the velocity dispersion and the possibility of close stellar encounters are high, we expect the following physical effects to play an important role in the process of accretion onto protoplanetary discs: ram pressure, angular momentum transport, dynamical encounters, and external photo-evaporation.

2.1. Ram pressure stripping

As the protoplanetary disc moves through the ISM, it experiences ram pressure, . This drag force can truncate the disc, depending on the gravitational force of the disc that keeps the latter bound to the star. By equating the gravitational force per unit area, or “pressure”, Pgrav = GMΣ(r)r-2, of the disc to Pram, we determine beyond which radius the ram pressure dominates and the disc is expected to be truncated. This truncation radius is given by (1)with M the mass of the star, and we have assumed that the surface density profile of the disc can be written as Σ(r) = Σ0(r/r0)n, where r0 is an arbitrary but constant radius to which Σ0 is scaled. After the material at radii larger than Rtrunc has been stripped from the disc, we expect the further evolution of the disc to be determined by the redistribution of angular momentum owing to accretion and viscous evolution. The pressure in the mid-plane of the disc is at least one order of magnitude smaller than the gravitational “pressure” and therefore does not play a significant role in resisting the ram pressure.

2.2. Redistribution of angular momentum

The redistribution of angular momentum in the disc is governed by two processes: (1) the accretion of material with no angular momentum material with respect to the disc and (2) the viscous evolution of the disc and consequent redistribution of its mass and angular momentum.

2.2.1. Accretion of ISM

The accretion of material with no azimuthal angular momentum lowers the specific angular momentum of the disc. Since the total angular momentum of the disc has to be conserved, mass and angular momentum will be redistributed within the disc. We can estimate this redistribution to first order, if we consider the disc consists of concentric rings with a thickness dr and mass mring(r) = 2πrΣ(r)dr. In a time interval dt, the ring will accrete an amount of mass from the ISM equal to maccr(r) = 2πrρISMvISMdtdr. The specific angular momentum, h, in the ring will decrease by a factor Σ(r) / (Σ(r) + ρISMvISMdt) and the ring will migrate to a smaller radius which corresponds to its new specific angular momentum. This process causes inward migration of material that belongs to the disc and leads to a contraction of the disc. This derivation does not take into account any interaction between adjacent rings, which is determined by viscous processes.

2.2.2. Viscous redistribution

Although the nature of the viscous processes that occur in accretion discs are still debated (see, for example, Armitage 2011), we do know that they are responsible for transporting material inwards through the disc. When this happens, some material has to move outward to conserve the total angular momentum of the disc, Jdisc. Both the viscous evolution and accretion of ISM cause material to drift inwards where it is eventually accreted onto the star and lost from the disc together with the angular momentum it carried. The outward spreading of material at the outer edge of the disc could make it more vulnerable to ram pressure stripping, which in turn also robs the disc of its angular momentum.

2.3. Dynamical encounters

In dense stellar systems, disc radii have been shown to be truncated because of close stellar encounters (Breslau et al. 2014; Rosotti et al. 2014; Vincke et al. 2015). The last study shows that, in dense clusters (), almost 40% of the discs is smaller than 100 AU and the median disc radius could be as small as 20 AU in the core. Close stellar encounters thus affect the surface area of the disc and, thereby, also the rate at which the disc can sweep up ISM. In this work, we only focus on the hydrodynamic aspects of accretion onto protoplanetary discs, in particular the accretion rate. The disc radius we find in our simulations is of similar size to the 20 AU found by Vincke et al. (2015). The question whether the process of ram pressure or dynamical encounters dominates the truncation of the disc in embedded dense stellar systems is beyond the scope of this paper.

2.4. External photo-evaporation

Globular clusters host a large number of massive stars in the early phases of their evolution and the large UV flux they produce may also strip material from the disc. Studies have shown that the fraction of stars that have discs can decrease by a factor of two close to O stars (see, for example, Balog et al. 2007; Guarcello et al. 2007, 2009; Fang et al. 2012), but Richert et al. (2015) argue that these results could be partly affected by sample incompleteness. Recently, Facchini et al. (2016) estimated that the mass-loss rate from the outer edge of a protoplanetary disc that is due to photo-evaporative winds could be of the order of 10-8−10-7M/ yr, see their Fig. 12. We do not take radiative processes into account in this work, but we note that any mass loss from the disc, additional to that found in this work, will shorten its lifetime with respect to our findings.

3. Numerical set-up

Our simulations are performed using the AMUSE environment (Portegies Zwart et al. 2013; Pelupessy et al. 2013)2. AMUSE is a Python framework in which a variety of astrophysical codes that have been published and well tested by the community, i.e. “community codes” can be used and combined. The simulation is set up in Python and AMUSE takes care of the communication between Python and the code(s) desired for use. This way, it is very easy to use the same set-up with different community codes and compare their outcomes.

Currently, two different smoothed particle hydrodynamics (SPH) codes are included in AMUSE, e.g. Fi (Pelupessy et al. 2004) and Gadget2 (Springel 2005; Springel et al. 2001). These two codes solve the dynamics of a self-gravitating hydrodynamic fluid using a Tree gravity-solver (Hernquist & Katz 1989). In this work, we use Gadget2 for our reference model and we verify the consistency and robustness of our results by comparing with Fi. Both SPH codes include self-gravity.

Because we want to be able to perform future simulations for a wide range of parameters, we need to find a balance between computational effort and convergence, i.e. the most efficient set-up. To test whether our particle resolution is sufficient, we perform simulations with different numbers of particles in the disc. This way, we can test whether the results converge and the numerical noise diminishes. The verification and validation are discussed in Sect. 4.3.1.

For our reference model, we use the vanilla version of Gadget2, which is freely available online3 and included in AMUSE, with the only exception that we have implemented the Morris & Monaghan (1997) viscosity formalism as is done in M09, see Sect. 3.2.

The set-up of our simulations is as follows (see Fig. 1): a stationary protoplanetary disc is positioned coaxially in a cylindrical, laminar flow of gas, representing the ISM. The protoplanetary disc is positioned in the centre of the cylinder. The inflow and outflow boundaries are each located at a distance of 500 AU from the centre of the protoplanetary disc. The radius of the cylinder is also set to 500 AU. Particles that flow outside the computational domain are removed from the simulation. The inflow consists of new particles. The parameters we have adopted are listed in Table 1.

thumbnail Fig. 1

Schematic overview of the numerical set-up.

Open with DEXTER

Table 1

Initial values of the parameters in our simulations.

3.1. Initial conditions and parameters

We adopt equal mass particles for the flow and the disc to prevent spurious forces on the interface between the ISM and the disc. We set the number of neighbours in the SPH codes to 64 ± 2 and use an isothermal equation of state. The sound speed, cs, equals 0.3 km s-1, which corresponds to a temperature of approximately 25 K for all particles. The isothermal equation of state can be justified by considering the cooling timescale (2)where n is the number density, kB the Boltzmann constant, T the temperature, and Λ(n,T) the cooling rate. For T = 25 K and n ≥ 103 cm-3, the cooling timescale τcool ≲ 16 yr, assuming Λ ≈ 10-26 erg cm3 s-1 (Neufeld et al. 1995). Any departure from the equilibrium temperature of 25 K will therefore be restored quickly with respect to our simulation timescale. We discuss the parameters and assumptions for the ISM and the disc separately below in Sects. 3.1.1 and 3.1.2, respectively.

3.1.1. Interstellar medium

By adopting a temperature of 25 K, we assume that cooling of the ISM is much more efficient than the radiative transfer of heat from nearby stars. It has been suggested that the Lyman-Werner flux plays an important role in the formation of multiple populations (Conroy & Spergel 2011), but a survey on 130 young massive clusters has shown little to no ionized gas (Bastian et al. 2013a). Young massive clusters are considered as being the modern-day counterpart of proto-GCs and this study therefore implies that heating does not play an important role in the environment where the accretion process is believed to take place.

We can also estimate the cooling timescale for stellar ejecta from Eq. (2), to justify that the expelled gas cools fast to low temperatures. As a lower limit for the density of stellar ejecta we assume n = 102 cm-3 and T = 104 K (Smith et al. 2007,see also B13). Combined with an appropriate cooling rate of Λ ≈ 10-25 erg cm3/ s (Richings et al. 2014) this results in a cooling timescale of roughly 6500 yr. We consider this an upper limit, because assuming a higher density for the stellar ejecta would result in an even shorter cooling timescale. Consequently, the cooling timescale is at least three to four orders of magnitude smaller than the 107 yr timescale on which the early disc accretion scenario is expected to take place. Our assumption that the ISM has cooled sufficiently and can be approximated with a temperature of 25 K is therefore justified.

To estimate the density of the ISM in the early disc accretion scenario we use the values given in B13 for the available processed material and the core radius of a typical GC, respectively 1.3 × 105M and 1 parsec. The average density would then be around 2 × 10-18 g cm-3. However, the density varies and can be higher locally. To provide a better comparison with M09, we adopt their assumed number density of n = 5 × 106 cm-3 and mean molecular weight μ = 2.3, which translates to a mass density, ρISM, of 1.92 × 10-17 g cm-3. In the case of GCs with multiple stellar populations, μ may be somewhat larger because the enriched populations exhibit a helium enhancement compared to primordial molecular clouds.

We assume an inflow velocity, vISM, of 20 km s-1 to approximate the typical velocity dispersion in GCs. This means that our set-up is in the supersonic regime and the treatment of the artificial viscosity is important. We discuss the artificial viscosity in Sect. 3.2. The high velocity gives a very small Bondi-Hoyle accretion radius, as we discuss in Sect. 3.1.3. The inflow is modelled by adding a slice of ISM with thickness vISMdt at the inflow boundary and a random uniform distribution of the SPH-particles within this slice. The ISM flow reaches the disc after 100 yr and the computational domain contains 6Ndisc ISM-particles when it is completely filled.

3.1.2. Disc

In the case of a disc in a steady state, the mid-plane temperature profile follows a simple power law, Tcrp, and the surface density profile, Σ, is proportional to , assuming a constant viscosity parameter αν in the Shakura & Sunyaev (1973) formalism (see e.g. Armitage 2011). Since we assume a constant temperature in the whole disc, p = 0 and , which corresponds to a minimum mass solar nebula (Weidenschilling 1977; Hayashi 1981) and is also assumed in, for example, M09 and Rosotti et al. (2014). The typical temperature of a protoplanetary disc at radii >10 AU is roughly consistent with a temperature of 25 K (20 K, see, for example, Armitage 2011). Although heating of the outer layers of the disc may cause photoevaporation, the mid-plane of the disc is shielded. At approximately 10 AU, the disc has to be heated to temperatures >103 K to effectively lose mass as a result of photoevaporation. An 0.4 M star has an effective temperature of 3 × 103 K, which is not high enough to unbind gas from the surface layer of the disc at radii >10 AU. At radii <10 AU extreme ultraviolet radiation from the star may cause mass loss from the disc in the order of 10-11−10-10M/yr (Armitage 2011; Font et al. 2004). This is two to three orders of magnitude less than the mass-loss rates found in this work and we conclude that taking heating by the star into account would not affect our results.

The stability of differentially rotating gaseous discs against self-gravity can be expressed in terms of the Toomre parameter Q (Toomre 1964), which is defined as (3)where Ω is the angular frequency. A disc becomes unstable if Q is less than unity at the outer edge of the disc. In our set-up, Q has a value of 44. Self-gravity is therefore not expected to lead to instabilities. As in M09, we assume that the initial mass of the disc is 1% of the mass of the star and we set the outer radius of the disc to 100 AU and the inner radius to 10 AU. Realistically, the disc should probably extend inwards towards the radius of the star, but it is computationally very expensive to simulate the disc on the orbital timescales at these small radii. During the simulation, particles in the disc will migrate inwards owing to viscous evolution of the disc and are accreted onto the star, resulting in a disc that extends further inwards. This set-up allows us to postpone the expensive calculations towards later times in our simulations, thereby significantly decreasing the duration of our simulations.

We position the disc perpendicular to the flow (see Fig. 1). This perfect alignment may not occur frequently in nature, but enables us to discern the relevant physical processes more clearly. In future work, we will investigate the influence of giving the disc an inclination with respect to the flow direction.

3.1.3. Star

We assume that the mass of the star is concentrated in a single point, which has a mass of 0.4 M. This corresponds to the typical mass expected in the early disc accretion scenario. The point mass is added as a collisionless particle to the SPH code and we treat it as a sink particle, i.e. it can accrete gas particles that fall within a certain (fixed) radius. The sink particle absorbs the mass and momentum carried by the gas particles it accretes.

We assume the radius of the sink particle is 5% of the Bondi-Hoyle accretion radius, RA, defined as (see also Bondi 1952; Shima et al. 1985) (4)in which M is the mass of the star, cs is the sound speed, and vISM the velocity of the ISM (with respect to the star). Using the values mentioned above gives us an accretion radius of 1.8 AU, which is significantly smaller than in M09, where RA ≈ 2Rdisc.

3.2. Artificial viscosity

Since the typical velocity dispersion in a GC is highly supersonic, it is important to resolve shocks. To do so, SPH codes introduce a numerical viscosity which is characterised by the parameters αSPH and βSPH, where βSPH has been introduced to prevent particle interpenetration in shocks with high Mach numbers (Springel 2010). Typical values for the parameters are αSPH ≃ 0.5−1 and βSPH = 2αSPH.

The numerical (shear) viscosity can also be used to model the physical viscous transport of matter in an accretion disc (Artymowicz & Lubow 1994), but this implies lower values of αSPH. This value can be derived using Artymowicz & Lubow (1994) to relate the artificial viscosity parameter αSPH to the standard viscosity parameter αν proposed by Shakura & Sunyaev (1973). Assuming that αν ≈ 0.01 (Armitage 2011) would correspond to an αSPH of roughly 0.02 in our reference model (), which is too low for numerical reliability. We therefore set αSPH initially to 0.1, as was also done in, for example, M09 and Rosotti et al. (2014). We implemented the viscosity switch proposed by Morris & Monaghan (1997) in Gadget24. In this treatment of the viscosity, every particle has its individual viscous parameter αSPH (and βSPH = 2αSPH), which is important because we simultaneously need a low value of αSPH in the disc and a high value of αSPH, up to 1, to resolve the shock in front of the disc. In the SPH code Fi, the viscosity remains constant throughout the simulation at a value of αSPH = 0.1. We have adopted βSPH = 1 for Fi. With these values, we minimize the viscous stresses in the disc which, owing to the low relative velocities, are dominated by the first order αSPH term, while preventing particle interpenetration in the strong accretion shock. We tested lowering the βSPH value below the adopted value and found that it can cause numerical artefacts.

Following Artymowicz & Lubow (1996), we calculate the timescale on which the disc undergoes significant viscous evolution, τν = Re Ω-1, assuming the Reynolds number, Re, and angular frequency, Ω, are given by (5)with r the radial distance to the star and h the (resolution-dependent) average smoothing length in the disc at that radius. For a simulation of 16 000 disc particles, and assuming the corresponding αν, this gives us a viscous timescale, τν, of roughly 10 000 yr at 10 AU. We can also calculate the physical viscous timescale at 10 AU by replacing h with the scale height of the disc, , at that radius. This leads to a timescale of 3 × 105 yr. To remain safely in the regime where the simulation is not dominated by the viscous evolution of the disc, we evolve the whole set-up for 2500 yr and start the inflow at t = 0. The numerical viscosity in the disc changes during the simulations with Gadget2 and we will discuss the numerical effects in Sects. 4.3.1 and 4.3.2.

thumbnail Fig. 2

a) Mass of the disc determined in 2 different ways: (1) by drawing a cylinder around the disc and adding the mass of all the particles in that volume (dotted black line) and (2) from the Hop algorithm we use (solid black line, see Sect. 3.3). The green dashed line shows the cumulative swept-up ISM mass, i.e. the sum of ISM in the disc and ISM that has been accreted onto the star, as derived with the Hop algorithm. b) Radius of the disc (solid black line) and the radius beyond which the ram pressure dominates, both derived from the surface density profile (see Sect. 3.3). c) Solid green horizontal line gives the mass of the volume around the disc, assuming it is completely filled with ISM. The dashed green line gives the mass of ISM that is in the volume but not in the disc, as derived with the Hop algorithm. The dotted black line gives the mass of ISM in the volume by calculating the difference between the total mass in the volume and the hop estimate for the disc, i.e. the difference between the solid black line and the dotted black line from the top figure. After ~1500 yr, the methods agree very well.

Open with DEXTER

3.3. Distinguishing the disc from the ISM

thumbnail Fig. 3

Edge-on (top row) and face-on (bottom row) snapshots from the simulations with our reference model. The column density is shown in logarithmic scale and integrated along the full computational domain. The spatial scale is indicated in the bottom left and the flow direction on the top left. The red circle (bottom row) and dotted bars (top row) indicate the size of the disc as determined from its radius (see Sect. 3.3).

Open with DEXTER

To determine how much ISM has been swept up by the disc at any moment in time, we need to differentiate between the disc and the ISM flow. This is not straightforward since there is no clear separation between the outer edge of the disc and the flow. To distinguish between the disc and the ISM flow, we use a clump-finding algorithm to identify different groups in the parameter space of the angular speed (vθ), the density (ρ), and the speed along the axis of the cylinder (vx) of each particle. In this parameter space, we expect the ISM particles to cluster around values of, respectively, 0 km s-1, ρISM, and vISM, while the disc particles will not. For the clump-finding algorithm we use Hop (Eisenstein & Hut 1998).

Figure 2a shows a comparison of how well this algorithm performs for our reference model, compared to drawing a coaxial cylinder around the disc, with a length of 100 AU and a radius of 120 AU, and adding the mass of all the SPH-particles in that volume. The difference between the two lines can be attributed almost completely to the ISM that is in the volume, but not part of the disc. To demonstrate this, Fig. 2c shows the difference between these two methods, the mass of ISM in the volume as determined with the Hop algorithm, and assuming the volume is completely filled with ISM. At the beginning of the simulation, when the disc is stripped of its outer parts by the ram pressure exerted by the flow, the algorithm has some difficulties in differentiating the disc from the ISM flow, but as soon as a stable situation is reached and the stripped outer edges of the disc have left the computational domain, it performs very well. The few spikes visible in Fig. 2 are artefacts: some particles are marked as part of the disc, while they are either being stripped from the outer edge (and already carry some angular momentum) or flow along the outer edge of the disc (and pick up some angular momentum from the disc). These artefacts do not influence our results since they are smoothed out when we linearly fit the data from the moment a steady state is reached.

Once we have determined which particles belong to the disc, we can also determine the surface density profile, Σ(r), and radius, Rdisc, of the disc. To do so, we calculate the column density of the disc, viewing the disc face-on. We then bin the obtained 2D surface density in concentric annuli, which gives us the surface density profile as a function of the radius. At t = 0, we determine Σ(100 AU), i.e. at the initial outer radius of the disc. At consecutive times, we recalculate the surface density profile and consider the radius at which Σ(r)|t = Σ(100 AU)|t = 0 is the radius of the disc at that moment in time. The snapshots of our reference model, shown in Fig. 3, illustrate that this method performs very well if the simulation has reached a steady state. Furthermore, the radius determined in this way also agrees with the estimate of the truncation radius owing to ram pressure stripping, Eq. (1), as can be seen in Fig. 2b, where we plotted both radii as a function of time. Much in the same way as we did for the disc radius, the truncation radius in Fig. 2b is derived from the surface density profile at each moment in time.

3.4. Angular momentum conservation

For TreeSPH codes, like Fi and Gadget2, angular momentum is not conserved exactly. The reasons for this are that (1) the gravity forces for a Barnes-Hut tree code are not exactly symmetric for particles in different parts of the domain and (2) interactions between particles in different levels of the time-step hierarchy are not exactly symmetric. To make sure that the loss of angular moment that we measure is not due to (the lack of) angular momentum conservation, we determined how well angular momentum is conserved in our reference model. To do so, we look at the change of total angular momentum, i.e. of all particles that have been in the computational domain over the course of the simulation, between 1500 and 2500 yr. The total angular momentum fluctuates around a mean value, without any increasing or decreasing trend over time. To establish a measure of angular momentum conservation, we determined the mean total angular momentum during the last 1000 yr of the simulation, which is 1.5 × 1044 kg m2 s-1, while the standard deviation is 1 × 1040 kg m2 s-1. Thus angular momentum is conserved up to about 1 part in 104. We compare this to the angular momentum loss of the disc in Sect. 4.2.

4. Results

We performed a number of simulations with both Gadget2 and Fi at different resolutions, as listed in Table 2. The main difference between both codes is the treatment of the artificial viscosity (see Sect. 3.2). Our reference model is a simulation with Gadget2 and 16 000 disc particles, labelled G16. We discuss our reference model in Sect. 4.2 and we use the other models for the convergence and consistency study, which we discuss in Sect. 4.3. Whenever we halve the number of disc particles, we double the number of simulations. We compare these results with runs from the SPH code Fi, which we have performed once for each resolution5. We finish the comparison between both SPH codes by discussing a run of Gadget2 with the same artificial viscosity implementation as used in Fi.

We choose simulation G16 as our reference model, so that we can compare our results with those of M09, who used the same code. We have not chosen one of the Gadget2 simulations with a higher resolution as our reference model, because we cannot compare this model to a simulation with the same number of disc particles with Fi. As a benchmark for the accretion rate onto the star and the mass loss from the outer edge of the disc, we performed a simulation of our reference model without inflow of ISM, labelled G16NF. We discuss this simulation first.

4.1. Reference model without inflow of ISM

We performed a simulation, labelled G16NF, of a protoplanetary disc without inflow of ISM to gauge the mass-loss rates at the inner and outer edge of the disc. We compare the mass and angular momentum loss of subsequent simulations to the quantities derived for this simulation. The results of this simulation are summarized in Table 3, which contains the summary of the simulations with Gadget2. The mass and angular momentum change rates are determined by plotting the quantity in question as a function of time, e.g. the mass of the disc, and fitting a linear function to it from the moment the simulation reaches a steady state, i.e. from 1500 yr onward. Every rate has been defined in this way, i.e. by a linear fit between 1500 and 2500 yr. We find an accretion rate of 4.5 × 10-8M/yr onto the star, star. This is similar to what has been observed for stars of 0.4 M (Muzerolle et al. 2005) and is theoretically expected for a disc in a steady state and constant αν, corresponding to αSPH = 0.1 (Shakura & Sunyaev 1973; Armitage 2011). We note that M09 found an accretion rate onto the star of 1.5 × 10-7M/yr in their isolated disc simulation.

Table 2

Different models for the convergence and consistency study.

Roughly half of the mass loss of the disc is due to accretion onto the star, because the inner edge of the disc migrates inwards. This can be seen in Fig. 4c, which shows the azimuthally averaged surface density at three different moments in time, i.e. at t = 0 (solid black), 1500 (dashed purple) and 2500 yr (dotted green). The other half is “lost” owing to the outward spreading of the outer edge as it leaves our computational domain, i.e. when the radial distance of the SPH particle to the star is more than 500 AU. The outward spreading can be inferred from the decreasing surface density profile at R > 60 AU in Fig. 4c. The rates and timescales that are quoted for model G16NF in Table 3 under stripping and angular momentum loss are therefore actually associated with viscous spreading of the disc. We have determined the timescales for viscous spreading by taking the total disc mass and angular momentum at the beginning of the interval and dividing it by the disc and , respectively. According to this extrapolation, the disc would viscously dissipate on a timescale of 7 × 104 yr. This is almost half the numerical viscous timescale at 100 AU, see Sect. 3.2, but should be considered as being more representative for the disc as a whole.

Table 3

Quantities of different runs with Gadget2, of which some are plotted in Fig. 5.

4.2. Reference model

Figure 3 shows snapshots from the simulations with our reference model. At t = 250 yr, we can see the ram pressure at play as the flow drags along disc material from radii larger than Rtrunc, which is 55 AU in this case. The third column shows that the simulation has reached a steady state at t = 1500 yr in which the disc is continuously accreting. The last snapshot illustrates that the disc is shrinking in size during the steady state. In Fig. 2a we plotted the mass of the disc in simulation G16 and the total amount of swept-up ISM as a function of time. The swept-up ISM is defined as the sum of the ISM that is in the disc at each moment in time and the ISM that has been accreted by the star up to that moment. The total amount of swept-up ISM increases and we can determine the ISM loading rate from the slope of this line, as described in Sect. 4.1. This gives a value of ISM = 1.03 × 10-7M/yr. We can compare this value to the rate that we expect based on a simple geometric estimate: (6)where ρISM and vISM are the density and velocity of the ISM, respectively, and Rdisc the radius of the disc. The ISM loading rate in the simulation is a factor of about two lower than the geometric rate. We show in Sect. 4.3.1 that the difference is independent of the resolution of our simulation or the code we used.

The result that the ISM loading rate is consistently a factor of two lower than given by Eq. (6)can be understood because, at the outer edge of the disc, ISM is not entrained by the disc. When the ISM flow first hits the disc, at t ≈ 100 yr, the steep increase of swept-up ISM does agree with the theoretical rate. This is because, initially, almost all ISM colliding with the surface of the disc is considered part of the disc, and this can be seen as an increase in the mass of the disc in Fig. 2a at t ≈ 100 yr. As the outer edges of the disc are dragged along with the flow, these regions are no longer associated with the disc and there is a large jump in the disc mass and a correspondingly smaller one in the swept-up ISM. From that moment on, the effective cross-section of the disc, i.e. the area with which it sweeps up ISM, is smaller than the actual surface area of the disc. This is illustrated in Fig. 4a, which shows the evolution of the surface density profile for the G16 simulation. The contribution of ISM to the surface density profile at each moment is indicated with filled regions in the corresponding colour. This illustrates that ISM is only entrained by the inner regions and not by the outer regions of the disc.

thumbnail Fig. 4

a) Azimuthally averaged surface density profile for the G16 simulation, Σ(r), at the start (solid black line), after 1500 yr (dashed purple line) and after 2500 yr (dotted green line). The contribution of swept-up ISM at each moment is plotted in the corresponding transparent colour. The vertical arrows in the same colour as the surface density profile indicate the radius of the disc at that moment. b) Same as the left figure, but for the F16 simulation. The initial inner radius at 10 AU is better resolved by Fi than by Gadget2. c) Same as a), but for the reference model without inflow of ISM.

Open with DEXTER

Although the disc continuously sweeps up ISM, it actually loses mass, see Fig. 2a, at a rate of 1.68 × 10-7M/yr, which is faster than the rate at which the disc sweeps up ISM. The disc loses mass at the inner edge owing to accretion onto the star, and at the outer edge owing to continuous stripping. The accretion rate onto the star is 1.27 × 10-7M/yr, roughly equal to the ISM loading rate and almost three times the accretion rate onto the star for an isolated disc. Sweeping up ISM thus enhances the accretion rate onto the star. The remaining mass is lost from the outer edges at a rate of 1.44 × 10-7M/yr. Not only is the ISM not entrained by the outer edge, it actually removes mass from those regions. This result can be explained in the following way. Owing to viscous torques within the disc, as discussed in Sect. 2.2, disc material migrates inwards. To conserve the total angular momentum of the disc, some disc material moves outwards. This outward diffusion lowers the surface density profile at the outer edge of the disc and material is therefore continuously stripped from the disc by the ram pressure of the ISM flow. The rate at which material is lost from the outer edge of the disc is four times higher than in the isolated case. The rate of angular momentum transport in the disc depends on the (artificial) viscosity and will be discussed in more detail in Sect. 4.3.2, where we compare our two different methods for modelling the viscosity.

The ablation at the outer edge of the disc causes both mass and angular momentum to be lost from the disc. During the last 1000 yr of the simulation, the angular momentum lost by the disc equals 5.43 × 1042 kg m2 s-1, which is more than 460σ, as determined in Sect. 3.4, and thus highly significant in comparison with errors in the angular momentum conservation in the code. Therefore, the loss of angular momentum of the disc cannot be ascribed to errors in time integration, but must be due to (the modelling of) physical processes in our simulation. The angular momentum lost owing to accretion onto the star in the same time interval, 1.4 × 1040 kg m2 s-1, is insignificant compared to the loss of angular momentum from the outer edge of the disc. In Sect. 4.3.1 we discuss that at least half of the angular momentum lost from the disc is caused by the interaction with the ISM at the outer edge of the disc in which angular momentum is transferred to the ISM and carried away in the flow. The remaining angular momentum loss is due to continuous stripping of original disc material from the outer edge of the disc. The angular momentum loss suggests a disc lifetime of 6 × 103 yr, which is an order of magnitude shorter than the timescale derived in the simulation without flow.

4.2.1. Evolution of the surface density profile

The slope of the surface density profile at a given radius in Fig. 4a increases with time as more ISM is entrained and disc material is transported inwards as a result of viscous evolution. From the start of the simulation, the artificial viscosity in the mid-plane of the disc increases and is higher than the initial value of 0.1. During the simulation, the typical value of αSPH in the mid-plane of the disc is 0.6. This means that the transport of mass and angular momentum through the disc in the simulation is faster than estimated in Sect. 3.2. In Sect. 4.3.2, we discuss how this compares to a simulation in which αSPH remains equal to 0.1 in the disc.

The purple and green filled regions show the contribution of ISM in the disc to the surface density profile at t = 1500 and 2500 yr. At each snapshot, the relative contribution of the ISM to the surface density profile is highest at small radii and decreases towards larger radii. The ISM in the disc migrates inwards owing to viscous torques and thus follows the same trend as the total surface density profile of the disc. Furthermore, new ISM with no angular momentum is continuously entrained by the disc, which also contributes to the inward migration of disc material. Both processes, continuous accretion, and viscous evolution steepen the exterior slope of both the total and ISM surface density profile. The total fraction of ISM in the disc increases with time as more material is swept-up.

4.3. Convergence and consistency

thumbnail Fig. 5

a) Average mass loss and gain rates, , as a function of the number of disc particles, Ndisc for both SPH-codes Fi (dotted) and Gadget2 (dashed). The change of the disc mass (black) consists of three components: mass is lost owing to accretion onto the sink particle (purple) and stripping at the outer edge (brown), and gained from sweeping up ISM by the disc and star (green). The rates are determined by a linear fit (see Sect. 4.3.1). b) Loss of angular momentum of the disc during the final 1000 yr interval as a fraction of its total angular momentum plotted against resolution for Fi (dotted) and Gadget2 (dashed). We only consider the angular momentum component along the (initial) rotation axis of the disc. c) Radius (black) and mass (purple) of the disc at the end of the simulation as a function of the number of disc particles, Ndisc for both SPH-codes Fi (dotted) and Gadget2 (dashed). The determination of these quantities is described in Sect. 3.3.

Open with DEXTER

In this section we compare the outcome of simulations with different numbers of disc particles and SPH codes. The runs for the convergence study are listed in Table 2. Figure 5a shows the average mass gain and loss rates that are relevant for our study: the accretion rate onto the star, the net rate at which ISM material is swept-up, the rate at which initial disc material is stripped, and the total rate of change in the disc mass. The rates are determined by considering the various mass quantities, e.g. the mass of the disc and the star, as a function of time during the last 1000 yr of the simulation and applying a linear fit to their slope. We summarize all the quantities plotted in Fig. 5 in Tables 3 and 4. We first discuss the convergence of the simulations, i.e. the effect of changing the number of disc particles, and subsequently the consistency between Gadget2 and Fi.

4.3.1. Convergence

Figure 5a shows that, except for the ISM loading rate, all rates have a decreasing trend with increasing resolution, with the exception of the accretion rate onto the star in the low resolution F4 simulation. The decreasing trends of the accretion rate onto the star, the stripping rate and the mass loss of the disc can be understood as follows. When the number of disc particles increases, the average smoothing length of the SPH particles decreases. Therefore, at lower resolution, the viscous timescale is shorter (according to Eq. (5)) and the transport of mass and angular momentum through the disc is faster than at higher resolutions. Moreover, as we discuss in Sect. 4.3.2, in our simulations with Gadget2 the artificial viscosity in the disc is higher for lower resolutions. This further increases the dependence of mass and angular momentum transport on the number of disc particles. As discussed in Sect. 4.2, there is a physical relation between at the inner edge of the disc, i.e. accretion onto the star, and at the outer edge of the disc, i.e. the ablation rate, which is the trend we observe in Fig. 5a. This implies that the transport of angular momentum in the simulations is dominated by numerical effects, since these rates are sensitive to the number of disc particles and have not yet converged.

The rate of mass change of the disc follows the same trend as the accretion and ablation rate, since the disc loses less mass at both the inner and outer edge with increasing resolution. The ISM loading rate, however, appears to be independent of the number of disc particles and is almost constant at all resolutions. This implies that the effective cross-section with which the disc entrains ISM is roughly equal for all resolutions. Figure 5c shows that the radius, i.e. total surface area, of the disc is larger at lower resolution, but the outskirts of the disc are also more diffuse in that case. As discussed above, if the outer edge of the disc is too diffuse, the ISM particles are not entrained by the disc. At higher resolution, the disc is more compact, as can be seen from both the decreasing radius and increasing mass in Fig. 5c. Again, this can be understood from the concept of mass and angular momentum transport: at high resolution, less mass is ablated and more mass remains in the disc, which in turn migrates inwards making the disc more compact. The surface density at the outer edge of the disc is therefore higher for a larger number of disc particles. The net effect across different resolutions is that the effective cross-section of the disc remains approximately the same, i.e. the effective cross-section depends on the height of the surface density profile in the outskirts of the disc. However, the ratio of the effective cross-section over the actual surface area of the disc, σdisc/Adisc, increases with resolution.

thumbnail Fig. 6

Angular momentum lost from the disc during the final 1000 yr in terms of the initial angular momentum of the disc for both Fi (dotted) and Gadget2 (dashed). The angular momentum loss associated with processes in the disc, e.g. accretion of ISM, angular momentum exchange between accreted ISM and original disc material, is plotted in green and angular momentum loss associated with ablation is shown in brown. We have omitted the angular momentum loss due to accretion onto the star, because it is insignificant (as discussed in Sect. 4.2).

Open with DEXTER

Figure 5b shows the amount of angular momentum lost from the disc during the last 1000 yr of the simulation as a fraction of its angular momentum at the start of that interval. It shows that more angular momentum is lost at lower resolution. The timescales that can be derived from this angular momentum loss are shown in Tables 3 and 4. These timescales are of the same order as, although mostly consistently smaller than, the timescales derived from the mass loss of the disc.

Table 4

Same as Table 3, but for the simulations with Fi.

To better understand how the angular momentum is lost, we split the angular momentum loss during the final 1000 yr into two components in Fig. 6: (1) the angular momentum loss associated with ablation and (2) the angular momentum lost owing to processes in the disc, i.e. the accretion of ISM and the exchange of angular momentum between the swept-up ISM and original disc material. The angular momentum loss is normalized to the initial angular momentum of the disc so that we can compare the components on an absolute scale. We do not show the angular momentum lost as a result of accretion onto the star, because it is two orders of magnitude smaller than the angular momentum loss caused by ablation (see Sect. 4.2). Ideally, the second component attributed to processes in the disc should be very small. The swept-up ISM should have no net azimuthal angular momentum and the amount of angular momentum gained by swept-up ISM in the disc should equal the amount that is lost by material that remains in the disc. This is not expected to be exactly zero, since the angular momentum that is lost from the outer edge of the disc was gained at the expense of material that remains in the disc.

The component related to stripped material can be separated into two constituents: angular momentum carried away by the original disc material and angular momentum carried away by ISM that briefly interacts with the disc, gains angular momentum and then moves along with the flow6. In our simulations with Gadget2, both constituents contribute equally to the angular momentum loss associated with stripping. However, in our simulations with Fi, 93% of the angular momentum loss associated with stripping is actually due to angular momentum being taken away by the ISM. This explains why the angular momentum loss associated with stripping differs by roughly a factor of two between both codes. Since stripping dominates the total angular momentum loss of the disc, the same factor of two difference between both codes is also seen in in Tables 3 and 4, when comparing simulations with the same resolution. We therefore argue that the angular momentum loss of the disc is dominated by frictional interactions between the ISM flow and the outer edge of the disc. The lifetime of the disc in our simulations is thus predominantly limited by angular momentum loss associated with this process. Both codes show a decreasing trend of angular momentum loss with increasing resolution and in Sect. 5.1 we discuss to what extent this angular momentum loss is physical.

4.3.2. Consistency

Although the trend of all quantities in Fig. 5 is the same for both codes, those that are determined from the simulations with Gadget2 are consistently lower than, or at most equal to, the same quantity determined from the simulations with Fi, with the exception of the disc radius at the lowest resolution. One of the fundamental differences between Gadget2 and Fi is the treatment of the artificial viscosity: in Fi it is constant, i.e. αSPH = 0.1, while in Gadget2 it can vary by an order of magnitude, depending on the local velocity gradient. As mentioned in Sect. 4.2.1, during the simulations with Gadget2, αSPH is not equal to 0.1 in the disc, as derived from accretion disc theory and observations (see Sect. 3.2), but approximately 0.6. In fact, αSPH in the mid-plane of the disc increases steeply inwards as a function of the radius between about and the star. At the outer edge of the disc αSPH is also higher, owing to the high velocity gradient caused by the shock. The value of αSPH as a function of radius is reflected by the surface density profiles of the G16 simulation in Fig. 4a. The radius at which αSPH has a minimum in the G16 simulation corresponds to the peak in the surface density profile. Thus matter accumulates at a radius in the disc where the transport of mass and angular momentum is slowest. In contrast, αSPH is constant in the F16 simulation and the surface density profile is therefore broader and less steep than in the G16 simulation (see Fig. 4b). Furthermore, αSPH depends on the number of disc particles: the higher the resolution, the lower αSPH (αSPH ≈ 0.4 at the highest resolution). Therefore, in a given simulation with Gadget2, more mass is transported inwards, and angular momentum is transported outwards accordingly than in the same simulation with Fi. Figure 6 shows, as discussed in the previous section, that in the simulations with Fi the angular momentum loss associated with ablation is consistently a factor of two lower than in the simulations with Gadget2. This implies that the angular momentum loss not only depends on resolution, but also on the modelling of physical processes. In principle, the smoothing length should be the same in simulations with the same number of particles for both codes. However, since the artificial viscosity is higher at smaller radii in Gadget2, numerical effects start dominating at an earlier stage in the simulation and at these small radii the smoothing length increases as a result of the faster inward movement of particles. At radii smaller than 2 AU, the numerical viscous timescale in the disc becomes smaller than the orbital period.

Despite these differences, the ISM loading rate in both codes and at all resolutions agree within a factor of 1.5, because we are looking at the whole star and disc system. In the case of Fi, accreted ISM resides in the disc for a longer time before it is accreted onto the star. So, to determine the ISM loading rate, it is less relevant how rapidly material is transported through the disc, as long as it has been swept up by the disc. Both SPH codes agree that the disc will always lose angular momentum, even if the disc gains mass in some simulations (i.e. the F16 and G128 simulations), and that the angular momentum is predominantly lost from the outer edge of the disc. However, the timescales derived from the angular momentum loss of the disc in the simulations with Fi differ by a factor of three from those derived with Gadget2 (see Tables 3 and 4).

4.3.3. Gadget2 and Fi with the same viscosity implementation

We performed a simulation with 16 000 disc particles using Gadget2 with the same (constant) viscosity implementation as in Fi, i.e. αSPH = 0.1 and βSPH = 1, to address the differences between the results of Gadget2 and Fi. We label this simulation G16CV and summarize the results in Table 3. The amount of swept-up ISM and the ISM loading rate are almost the same as in the F16 simulation. The main difference between this simulation and the F16 simulation is that both disc and ISM material are accreted faster onto the star in the G16CV simulation than they are in the F16 simulation. As a result, the disc is losing mass in the G16CV simulation, even though no original disc material is stripped from the outer edge of the disc, i.e. strip = 0 M/yr. The ISM loading rate in the G16CV simulation is dominated by the accretion rate of ISM onto the star; the amount of ISM in the disc remains roughly constant as a function of time when the system has reached a steady state. This also partly explains the difference with in the G16 simulation, where the accretion rate onto the star is dominated by the accretion of disc material and, therefore, a factor of two lower. As in the F16 simulation, the disc loses its angular momentum predominantly through interaction with the ISM flow at the outer edge of the disc. However, in the G16CV more angular momentum is transferred from the outer edge of the disc to the ISM flow than in the F16 simulation. Even though the viscosity prescription is the same in both simulations, the inward transport of mass occurs faster in the G16CV simulation. We have not been able to pinpoint the exact cause of the difference between the two codes. It could be that the discrepancies are attributable to other differences in, for example, the implementation of the time-stepping or the limiters on the acceleration, which are more difficult to discern.

5. Discussion

We discuss the uncertainties in the modelling of the physical processes and other caveats of our simulations below. After discussing the limitations in our simulations, we make a comparison to other work and finally discuss the implications of our results.

5.1. Angular momentum loss

As discussed in Sect. 4.3.2, the angular momentum loss due to ablation of the disc is strongly dependent on the artificial viscosity parameter αSPH in the disc and the ablation rate of the disc also shows a decreasing trend as a function of the resolution (see Sect. 4.3.1). Furthermore, in our simulations with Fi, angular momentum is predominantly lost to the ISM through the exchange of angular momentum at the outer edge of the disc, instead of being carried away by stripped disc material. The angular momentum loss associated with these processes shows a decreasing trend with increasing resolution. This suggests that the ablation in our simulation, and the associated angular momentum loss, is dominated by numerical effects and that, physically, it may not be the dominant process for the loss of angular momentum from the disc. In the simulation of M09, with the same version of Gadget2 but with a higher resolution and smaller flow velocity (i.e. lower ram pressure), the stripping of the outer edge does not play a significant role (see Sect. 5.5).

If the ablation of the disc is predominantly numerical, then the timescale on which the disc loses its mass cannot be considered a reliable indicator for the lifetime of the disc. Although the timescale of angular momentum loss is generally shorter than the timescale derived from the mass loss of the disc, we consider the former to be a more reliable estimator for the disc life time. In particular, we consider the timescale of angular momentum loss determined with Fi to be more indicative, because in those simulations the viscosity in the disc agrees better with accretion disc theory and the disc actually gains mass. The angular momentum loss timescale in the F16 simulation is similar to that in the G128 simulation, i.e. highest resolution with Gadget2, but it is still an order of magnitude smaller than the physical estimate in Sect. 3.2.

We cannot directly interpret the timescale of angular momentum loss as the lifetime of the disc, since it is dominated by numerical effects in our simulations. Physically, the decrease in specific angular momentum of the disc owing to accretion of ISM with zero azimuthal angular momentum may contribute equally to, or even dominate, the decreasing disc size. Both processes are also non-linear and we consider our estimate of the angular momentum loss timescale as a lower limit to the actual lifetime of the disc.

5.2. Viscous evolution

We have used two SPH codes that model the viscosity in a different way. It is clear that the transport of mass and angular momentum in the disc behaves differently in both codes even if the same viscosity prescription is implemented in both codes. To model the (viscous) evolution of the disc correctly, additional relevant physical processes, e.g. magnetohydrodynamics, radiative transfer and radial mixing (which are not all physically well understood and beyond the scope of this paper), should be incorporated. In this work, we are interested in how much mass a protoplanetary disc can sweep up and how this process would affect its lifetime. The result for the ISM-loading appears to be independent of the viscous modelling of the disc and we therefore conclude that this result is robust. Furthermore, the simulations indicate that the process of face-on accretion onto a protoplanetary disc will likely shorten, rather than prolong its lifetime.

5.3. Resolution

As mentioned in Sect. 3, we are trying to find a balance between computing time and convergence. Our reference model, G16, took 3.5 days to run on 32 cores, while the comparison F16 simulation took a little more than two months on 22 CPUs. When looking at the ablation rate as a function of increasing resolution, the number of particles that are stripped increases, but the mass that they carry away decreases more rapidly. We therefore argue that the ablation rate is not caused by numerical noise but rather, as discussed above, by the treatment of viscosity in the disc. This is supported by the result that the Poisson noise for the low-resolution models that we have performed multiple times is less than 5%, see Table 3. Furthermore, the ISM loading rate is robust for all our simulations and does not seem to be affected by numerical noise.

5.4. Modelling of the ISM

To be able to discern the relevant physical processes from the simulations, we have modelled the ISM in a very idealized way: as a homogeneous gas with no clumps and no turbulence. To a certain extent we have modelled the reaction of the disc to a density gradient when the flow first hits the disc. In that case, the size of the disc is determined by the ram pressure. If the disc were to encounter a region of gas with lower density, the surface area of the disc can extend to larger radii without being stripped and that would probably increase its effective cross-section. Despite the lower density, the ISM loading rate could still be of the same order, because a lower ISM density allows a larger surface area of the disc to collect ISM. In that sense, accretion onto protoplanetary discs may be a self-regulating process. In a follow-up work, we will perform simulations with different ISM densities and velocities to investigate if and how the ISM loading rate depends on these quantities.

In principle, it would be possible to add angular momentum to the disc, such that it maintains or prolongs its lifetime, if the ISM has turbulence or substructures on scales at or slightly smaller than the disc scale. However, even in that case, the net effect is unlikely to increase the overall lifetime, since the disc would quickly encounter a different part of the ISM where an adverse, disruptive configuration of substructure might exist.

5.5. Comparison to other work

A similar simulation has been performed by M09, but for a higher disc mass and lower velocity of the ISM with respect to the disc. In their work, both the Bondi-Hoyle radius and the radius beyond which ram pressure would dominate are much larger than the radius of the disc. They used a higher resolution for the disc, initially ~2.5 × 105 particles, which had 8 times the mass of their ISM particles. However, in their work the ISM also does not reside in the outer regions of the disc (see their Fig. 3), meaning that even at much higher resolution the effective surface area of the disc is smaller than its actual surface area. We interpret this as a physical effect: at the outer edge of the disc, the ISM is forced to flow around it and is not entrained.

Furthermore, the disc in the simulation of M09 also shows a decreasing radius and steepening surface density profile as a function of time. They attribute this to the redistribution, i.e. inward movement of disc material owing to the accretion of ISM with no angular momentum. The loss of disc material to the ISM in their simulation is negligible, the disc mass only decreases as a result of accretion onto the star. This agrees with the trend suggested in our convergence study, i.e. that the continuous stripping of the disc at the outer edge is a numerical artefact. However, it could also be due to the much lower velocity in their work.

A future parameter study at other, e.g. lower, densities and velocities of the ISM could provide a more decisive answer.

5.6. Implications for the early disc accretion scenario

Although we find that the ISM-loading rate corresponds relatively well to the geometric rate used in B13 for their early disc accretion scenario, they assume a disc radius that is significantly larger, i.e. 100 AU, than the radius found in our simulation. Moreover, the size of the disc decreases continuously during the process of accretion, and we have assumed the idealized case in which the disc is positioned exactly perpendicular to the flow. The ISM-loading rate derived in this work is therefore most probably an upper limit for the average rate that one would expect for a population of discs that all have different inclinations between their rotational axis and velocity vector.

Furthermore, B13 assume a disc lifetime of 107 yr from observations, while we find that angular momentum transport plays a significant role in shortening the lifetime of the disc in dense ISM environments. It therefore seems unlikely that a low-mass star can accrete of the order of its own mass via its protoplanetary disc, as required, in the early disc accretion scenario.

5.7. Implications for planet formation

The process of accretion onto protoplanetary discs may also affect planet formation. It could play a role in the formation of “hot Jupiters”, i.e. massive planets that have formed further out in the protoplanetary disc and migrated inwards. The mass loading and consequent redistribution may also increase the probability of forming a planet in the habitable zone. Ronco & de Ela (2014) find that a steeper surface density profile, i.e. more mass at smaller radii, increases the probability of forming a planet with a significant water content in the habitable zone. The constraint for an initial steep surface density profile could perhaps be eased, when the entraining of ISM onto the protoplanetary disc causes the disc material to migrate inward. These suggestions could be tested by incorporating our findings in detailed planet formation models.

6. Conclusions

We have performed simulations of accretion of interstellar material (ISM) onto a protoplanetary disc with two different smoothed particle hydrodynamics codes. We find that, as theoretically expected, when the flow of ISM first hits the disc, all disc material beyond the radius where ram pressure dominates is stripped. As ISM is being accreted and disc material migrates inwards, the disc becomes more compact and the surface density profile increases at smaller radii. We find that the ISM loading rate, i.e. the rate at which ISM is entrained by the disc and star, is approximately constant across all our simulations with both codes and is a factor of two lower than the rate expected from geometric arguments (see Eq. (6)). This difference arises because the outskirts of the disc do not entrain ISM and therefore the effective cross-section of the disc is smaller than its physical surface area. We find that, despite the accretion of ISM, the net effect is that the disc loses mass, except in the highest resolution simulations with both codes, where the disc gains mass. This decreasing trend with resolution implies that the net mass loss from the disc in our low-resolution simulations is numerical. Considering the timescale on which the disc loses all of its angular momentum rather than its mass, provides an estimate of about 104 yr. The angular momentum loss from the disc in our simulations is dominated by continuous stripping of disc material and by transfer of angular momentum to the ISM as it flows past the outer edge of the disc. Our convergence and consistency study, as well as previous work, indicate that these effects are predominantly numerical. The timescales estimated from the simulations with the highest resolution therefore provide a lower limit to the lifetime of the disc. The loss of angular momentum owing to accretion of disc material onto the star, which is governed by the (modelling of) viscous processes in the disc, is two orders of magnitude smaller than the loss associated with stripping.

Even if the disc grows in mass, the (specific) angular momentum of the disc will always decrease in this scenario, if not for the aforementioned angular momentum-loss processes then by accretion of ISM with no azimuthal angular momentum. Either way, the disc will shrink in size, thereby decreasing its effective cross-section. Although our ISM-loading rate agrees within a factor of two with the geometrically estimated rate, the lifetime and size of the disc are probably not sufficient to accrete the amount of mass required in the early disc accretion scenario.

In future work, we will extend our simulations to explore the parameter space and conditions that correspond to a broader range of stellar environments to find a parametrization for the mass-loading rate onto a protoplanetary disc system in terms of the density and velocity of the ambient medium and the size of the disc.


1

Most scenarios proposed to date imply subsequent epochs of star formation and hence refer to multiple generations of stars in a GC. Since it is still not clear whether GCs can facilitate an extended star formation history or if the enriched stars actually belong to the initial population, we will refer to multiple populations of stars.

4

This implementation is done in the source code of Gadget2, which is not in accordance with the philosophy of AMUSE, but is always possible if required.

5

We only did one simulation for each resolution with Fi, because it is computationally more expensive. 16 000 disc particles is the highest resolution we could simulate with Fi in a reasonable amount of time.

6

The Hop algorithm considers this ISM as being associated with the disc during a few time steps. This is an artefact of the Hop algorithm, but these interactions do carry away angular momentum from the outer edge of the disc nonetheless. In general, it is difficult to truly disentangle all components, also considering that angular momentum stripped from the outer edge of the disc has been gained from material that still resides in the disc. We therefore only provide the total angular momentum loss rates and timescales.

Acknowledgments

We thank Selma de Mink, Nate Bastian, and Nickolas Moeckel for valuable discussions and input. We are also grateful to the referee whose valuable comments have improved this work. This research is funded by the Netherlands Organisation for Scientific Research (NWO) under grant 614.001.202. NWO also granted computational resources on Cartesius under grant SH-295-14.

References

Online material

Animation of the simulation (Access here)

All Tables

Table 1

Initial values of the parameters in our simulations.

Table 2

Different models for the convergence and consistency study.

Table 3

Quantities of different runs with Gadget2, of which some are plotted in Fig. 5.

Table 4

Same as Table 3, but for the simulations with Fi.

All Figures

thumbnail Fig. 1

Schematic overview of the numerical set-up.

Open with DEXTER
In the text
thumbnail Fig. 2

a) Mass of the disc determined in 2 different ways: (1) by drawing a cylinder around the disc and adding the mass of all the particles in that volume (dotted black line) and (2) from the Hop algorithm we use (solid black line, see Sect. 3.3). The green dashed line shows the cumulative swept-up ISM mass, i.e. the sum of ISM in the disc and ISM that has been accreted onto the star, as derived with the Hop algorithm. b) Radius of the disc (solid black line) and the radius beyond which the ram pressure dominates, both derived from the surface density profile (see Sect. 3.3). c) Solid green horizontal line gives the mass of the volume around the disc, assuming it is completely filled with ISM. The dashed green line gives the mass of ISM that is in the volume but not in the disc, as derived with the Hop algorithm. The dotted black line gives the mass of ISM in the volume by calculating the difference between the total mass in the volume and the hop estimate for the disc, i.e. the difference between the solid black line and the dotted black line from the top figure. After ~1500 yr, the methods agree very well.

Open with DEXTER
In the text
thumbnail Fig. 3

Edge-on (top row) and face-on (bottom row) snapshots from the simulations with our reference model. The column density is shown in logarithmic scale and integrated along the full computational domain. The spatial scale is indicated in the bottom left and the flow direction on the top left. The red circle (bottom row) and dotted bars (top row) indicate the size of the disc as determined from its radius (see Sect. 3.3).

Open with DEXTER
In the text
thumbnail Fig. 4

a) Azimuthally averaged surface density profile for the G16 simulation, Σ(r), at the start (solid black line), after 1500 yr (dashed purple line) and after 2500 yr (dotted green line). The contribution of swept-up ISM at each moment is plotted in the corresponding transparent colour. The vertical arrows in the same colour as the surface density profile indicate the radius of the disc at that moment. b) Same as the left figure, but for the F16 simulation. The initial inner radius at 10 AU is better resolved by Fi than by Gadget2. c) Same as a), but for the reference model without inflow of ISM.

Open with DEXTER
In the text
thumbnail Fig. 5

a) Average mass loss and gain rates, , as a function of the number of disc particles, Ndisc for both SPH-codes Fi (dotted) and Gadget2 (dashed). The change of the disc mass (black) consists of three components: mass is lost owing to accretion onto the sink particle (purple) and stripping at the outer edge (brown), and gained from sweeping up ISM by the disc and star (green). The rates are determined by a linear fit (see Sect. 4.3.1). b) Loss of angular momentum of the disc during the final 1000 yr interval as a fraction of its total angular momentum plotted against resolution for Fi (dotted) and Gadget2 (dashed). We only consider the angular momentum component along the (initial) rotation axis of the disc. c) Radius (black) and mass (purple) of the disc at the end of the simulation as a function of the number of disc particles, Ndisc for both SPH-codes Fi (dotted) and Gadget2 (dashed). The determination of these quantities is described in Sect. 3.3.

Open with DEXTER
In the text
thumbnail Fig. 6

Angular momentum lost from the disc during the final 1000 yr in terms of the initial angular momentum of the disc for both Fi (dotted) and Gadget2 (dashed). The angular momentum loss associated with processes in the disc, e.g. accretion of ISM, angular momentum exchange between accreted ISM and original disc material, is plotted in green and angular momentum loss associated with ablation is shown in brown. We have omitted the angular momentum loss due to accretion onto the star, because it is insignificant (as discussed in Sect. 4.2).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.