The effect of external environment on the evolution of protostellar disks

# The effect of external environment on the evolution of protostellar disks

###### Key Words.:
protoplanetary disks–stars:formation–hydrodynamics–stars:protostars

## Abstract

Context:

Aims:Using numerical hydrodynamics simulations we studied the gravitational collapse of pre-stellar cores of sub-solar mass embedded into a low-density external environment.

Methods:Four models with different magnitude and direction of rotation of the external environment with respect to the central core were studied and compared with an isolated model.

Results: We found that the infall of matter from the external environment can significantly alter the disk properties as compared to those seen in the isolated model. Depending on the magnitude and direction of rotation of the external environment, a variety of disks can form including compact ( AU) ones shrinking in size due to infall of external matter with low angular momentum, as well as extended disks forming due to infall of external matter with high angular momentum. The former are usually stable against gravitational fragmentation, while the latter are prone to fragmentation and formation of stellar systems with sub-stellar/very-low-mass companions. In the case of counterrotating external environment, very compact ( AU) and short-lived ( yr) disks can form when infalling material has low angular momentum. The most interesting case is found for the infall of counterrotating external material with high angular momentum, leading to the formation of counterrotating inner and outer disks separated by a deep gap at a few tens AU. The gap migrates inward due to accretion of the inner disk onto the protostar, turns into a central hole, and finally disappears giving way to the outer strongly gravitationally unstable disk. This model may lead to the emergence of a transient stellar system with sub-stellar/very-low-mass components counterrotating with respect to that of the star.

Conclusions:

## 1 Introduction

Solar-type stars form from the gravitational collapse of dense molecular cloud cores, which in turn can be the product of supersonic turbulent compression (e.g. Padoan et al., 2001) or slow expulsion of magnetic field due to ambipolar diffusion (e.g. Shu et al., 1987; Basu & Mouschovias, 1994). A large fraction of the collapsing core passes through a protostellar disk, formed due to conservation of the net angular momentum of the core, before finally landing on a protostar. The main accretion phase lasts until most of the parental core is accreted onto the disk plus protostar system or dispersed via feedback from stellar irradiation and outflows. In this ”classic” scenario, the core is usually described by an isolated, initially unstable Bonnor-Ebert-type sphere or perturbed singular isothermal sphere.

However, recent observations indicate that dense molecular cores constitute part of a complicated filamentary structure in the spatial and velocity space, with individual filaments resembling a network of twisted and intersecting threads (Hacar et al., 2013). A general agreement between the measured core masses and those predicted from the linear stability analysis suggests that the cores are likely to form from the gravitational fragmentation of these filaments (Hacar & Tafalla, 2011). In this picture, prestellar cores may not be described by a perfectly spherical shape and may be submerged into a varying density environment characterized by gas flows with varying angular momentum. As a net result, it is likely that the protostar plus disk system may continue accreting from the external environment even after the parental core has completely dissipated.

This conjecture seems to be supported by recent observations showing that large star-forming regions contain a population of relatively old stars (with age  10 Myr), which still have signatures of accretion disks (Beccari et al., 2010; De Marchi et al., 2013a, b). Scicluna et al. (2014) explain this phenomenon by invoking Bondi-Hoyle accretion to rebuild a new disk around these stars during passage through a clumpy molecular cloud. While it is difficult to determine the proper stellar age of accreting stellar objects (Baraffe et al., 2012) and it is not clear if Bondi-Hoyle accretion is the right mechanism to describe accretion onto the disk, the fact that there may be quite old disks (with age much greater than the canonical disk-dispersal -folding time of 3 Myr) supports the idea of prolonged accretion.

Theoretical and numerical studies of clustered star formation also suggest that accretion may continue from larger (than the initial extent of the parental core) spatial scales until the gas reservoir within the cluster is completely exhausted, leading to a constant accretion time that is proportional to the global freefall time of the cluster (e.g. McKee & Offner, 2010; Dunham et al., 2014). In particular, numerical simulations of, e.g., Bate et al. (2010), Fielding et al. (2014), and Padoan et al. (2014) demonstrate that molecular clouds are turbulent and inherently chaotic environments, within which individual collapsing cores may accrete mass and angular at a highly non-steady rate from a network of filaments interconnecting the cores. It is not impossible that the angular momentum vector of the accreted material can undergo significant changes both in magnitude and direction.

Finally, observations seem to indicate longer lifetimes of the main accretion phase, 0.42 Myr (Evans et al., 2009), than those derived from numerical simulations of isolated gravitationally unstable cores, 0.12 Myr (Vorobyov, 2010; Dunham & Vorobyov, 2012), again implying that mass infall onto the disk/star system may continue after the initial parental core has dissipated.

In this paper, we explore the effect of prolonged gas infall onto the properties of circumstellar disks. In particular, we consider prestellar cores submerged into a low-density, external environment with varying angular momentum. As a first approximation, we assume that the external material has constant density and angular velocity. More realistic initial conditions taken from numerical hydrodynamics simulations of clustered star formation will be considered in a follow-up study.

## 2 Model description

Our numerical hydrodynamics model for the formation and evolution of a star plus disk system in the thin-disk limit is described in detail in Vorobyov & Basu (2010) and Vorobyov et al. (2013b) and is briefly reviewed below for the reader’s convenience. We start our numerical simulations from the gravitational collapse of a starless cloud core submerged into a low-density external environment, continue into the embedded phase of star formation, during which a star, disk, and envelope are formed, and terminate our simulations after approximately 0.5 Myr of disk evolution. The protostellar disk occupies the inner part of the numerical polar grid and is exposed to intense mass loading from the infalling envelope and external environment.

To avoid too small time steps, we introduce a “sink cell” at  AU and impose a free inflow inner boundary condition and a free outflow outer boundary condition so that that the matter is allowed to flow out of the computational domain but is prevented from flowing in. The sink cell is dynamically inactive; it contributes only to the total gravitational potential and secures a smooth behaviour of the gravity force down to the stellar surface. During the early stages of the core collapse, we monitor the gas surface density in the sink cell and when its value exceeds a critical value for the transition from isothermal to adiabatic evolution, we introduce a central point-mass object. In the subsequent evolution, 90% of the gas that crosses the inner boundary is assumed to land on the central object. The other 10% of the accreted gas is assumed to be carried away with protostellar jets.

The basic equations of mass, momentum, and energy transport in the thin-disk limit are

 ∂Σ∂t=−∇p⋅(Σ% \boldmathvp), (1)
 ∂∂t(Σ\boldmathvp) + [∇⋅(Σ\boldmathvp⊗% \boldmathvp)]p=−∇pP+Σ% \boldmathgp+ (2) + (∇⋅Π)p,
 Missing or unrecognized delimiter for \left (3)

where subscripts and refer to the planar components in polar coordinates, is the mass surface density, is the internal energy per surface area, is the vertically integrated gas pressure calculated via the ideal equation of state as , is the radially and azimuthally varying vertical scale height determined in each computational cell using an assumption of local hydrostatic equilibrium, is the velocity in the disk plane, and is the gradient along the planar coordinates of the disk. The dependence of on temperature is parameterized based on the work of Masunaga & Inutsuka (2000) with for the gas temperature less than 100 K and otherwise1. The gravitational acceleration in the disk plane, , takes into account self-gravity of the disk, found by solving for the Poisson integral (see details in Vorobyov & Basu, 2010), and the gravity of the central protostar when formed. Turbulent viscosity is taken into account via the viscous stress tensor , the expression for which is provided in Vorobyov & Basu (2010). We parameterize the magnitude of kinematic viscosity using the -prescription with a spatially and temporally uniform .

The radiative cooling in equation (2) is determined using the diffusion approximation of the vertical radiation transport in a one-zone model of the vertical disc structure (Johnson & Gammie, 2003)

 Λ=FcσT4mpτ1+τ2, (4)

where is the Stefan-Boltzmann constant, is the midplane temperature of gas2, is the mean molecular weight, is the universal gas constant, and is a function that secures a correct transition between the optically thick and optically thin regimes. We use frequency-integrated opacities of Bell & Lin (1994). The heating function is expressed as

 Γ=FcσT4irrτ1+τ2, (5)

where is the irradiation temperature at the disc surface determined by the stellar and background black-body irradiation as

 T4irr=T4bg+Firr(r)σ, (6)

where is the uniform background temperature (in our model set to the initial temperature of the natal cloud core) and is the radiation flux (energy per unit time per unit surface area) absorbed by the disc surface at radial distance from the central star. The latter quantity is calculated as

 Firr(r)=L∗4πr2cosγirr, (7)

where is the incidence angle of radiation arriving at the disc surface (with respect to the normal) at radial distance .

The stellar luminosity is the sum of the accretion luminosity arising from the gravitational energy of accreted gas and the photospheric luminosity due to gravitational compression and deuterium burning in the stellar interior. The stellar mass and accretion rate onto the star are determined using the amount of gas passing through the sink cell. The properties of the forming protostar are calculated using a stellar evolution code described in Baraffe & Chabrier (2010). This code is coupled with the main hydrodynamical code in real time. However, due to heavy computational load the stellar evolution code is invoked to update the properties of the protostar only every 20 yr, while the hydrodynamical time step may be as small as a few weeks and the entire duration of numerical simulations may exceed 1.0 Myr. The input parameter for the stellar evolution code is the accretion rate , while the output are the radius and photospheric luminosity of the protostar. This coupling of the disk and protostar evolution allows for a self-consistent determination of the radiative input of the protostar into the disk thermal balance, which is important for the accurate study of disk instability and fragmentation.

## 3 Initial conditions

As the initial setup, we take a pre-stellar core submerged into a constant-density external environment. For the initial surface density profile of the core, we adopt a simplified form of a vertically integrated Bonnor-Ebert sphere (Dapp & Basu, 2009) The resulting initial distribution of the gas surface density takes the following form

 Σ=⎧⎪⎨⎪⎩r0Σ0√r2+r20for r≤Rcore,Σextotherwise, (8)

where is the gas surface density at the center of the core, is the radius of the central plateau of the core, is the radius of the core, is the initial sound speed, and is the density of the external environment, the value of which is set equal to the gas surface density at the outer edge of the core (). In all models the value of is set to 1.2 and the initial temperature is set to 10 K.

To study the effect of infall from the external environment characterized by varying angular momentum, we adopt the following form for the initial radial profile of angular velocity

 Unknown environment '% (9)

where is the angular velocity of the core and is a dimensionless scaling factor defining the angular velocity in the external environment. The angular velocity profile between the core and the external environment is smoothly joined at the outer edge of the core.

We have considered 5 models, the initial setup for which is shown in Figure 1. In particular, the left panels present the radial distribution of , while the right panels show as a function of radial distance. Model 1 represents a uniformly rotating isolated core (no external environment) with mass and ratio of rotational to gravitational energy . Models 2–5 consider the same core submerged into the external environment with different angular velocities. In particular, model 2 has a corotating external environment with the same angular velocity as that of the core, while in model 3 the angular velocity of the external environment is 10 times smaller than that of the core. Models 4 and 5 are characterized by a counterrotating external environment. In particular, in model 4 the external angular velocity is ten times smaller by absolute value than that of the core, while in model 5 the external angular velocity is equal to that of the core by absolute value. In models 2-5 the radial distribution of gas surface density is identical. The parameters of every model are provided in Table 1.

## 4 Results

In this section, we start with considering the properties of a protostellar disk formed from the gravitational collapse of an isolated pre-stellar core and continue by considering models that take into account the effect of external environment.

### 4.1 The isolated model

As a prototype model, we consider an isolated prestellar core with mass and ratio of rotational to gravitational energy (hereafter, model 1). This choice is motivated by our previous numerical hydrodynamics simulations indicating that protostellar disks formed from cores with similar characteristics are expected to be gravitationally unstable and prone to fragmentation in the early embedded stage of stellar evolution (Vorobyov, 2011). Figure 2 presents the time evolution of the gas surface density in model 1 in the inner  AU box. The whole computational region extends to 10000 AU. The time elapsed from the formation of the central star is indicated in each panel. The scale bar is in log g cm and the lower surface density threshold is set to  g cm so that the black region is not empty but filled with a low density material falling onto the disk from the parental core.

The protostellar disk in model 1 becomes gravitationally unstable at the very early stage of evolution (as manifested by a spiral structure already present at  Myr) due to a continuing mass loading from the infalling core. The first fragment forms in the disk after  Myr, when the disk mass and radius begin to exceed 0.1  and 100 AU, respectively. This is in agreement with many previous numerical and theoretical studies, indicating that the disk needs to acquire a critical mass () and grow to a critical size ( AU) before gravitational fragmentation can ensue (Stamatellos & Whitworth, 2008; Rice et al., 2010; Vorobyov, 2011; Meru & Bate, 2012). We used the fragment tracking algorithm described in Vorobyov (2013) to calculate the masses of the fragments. The fragments were identified based on two conditions: 1) the fragment must be pressure-supported, with a negative pressure gradient with respect to the center of the fragment and 2) the fragment must be kept together by gravity, with the potential well being deepest at the center of the fragment. The resulting masses of the fragments lie in the 2.5–5.0  range.

One can see in Figure 2 that the process of disk fragmentation in model 1 is not continuous but intermittent – time intervals showing fragments are alternated by intervals showing no signs of disk fragmentation. The timescale with which the disk fragments is determined by the ratio (Vorobyov et al., 2013a)

 Tfr=Md˙Minfall, (10)

where is the disk mass and is the rate of mass infall from the parental core onto the disk. Taking and  yr for the typical values of and (see Figs. 10 and 11), we obtain  Myr, which appears to agree (within a factor of unity) with a visually estimated value of  Myr (see Fig. 2). The lifetime of fragments is usually shorter than 0.1 Myr, either due to migration on the central protostar or destruction via tidal torques (Stamatellos et al., 2011; Vorobyov et al., 2013a), resulting in prolonged periods of disk evolution without fragments. Finally, it is worth noting that the disk forms only one fragment at a time, which implies that ejection of fragments into the intercluster medium due to gravitational many-body interaction in such a system is impossible (Basu & Vorobyov, 2012).

As many numerical studies indicate, a fragment can be driven into the disk inner regions and probably onto the star (Vorobyov & Basu, 2006, 2010; Machida et al., 2011a; Cha & Nayakshin, 2011), ejected into the intracluster medium (Basu & Vorobyov, 2012), dispersed by tidal torques (Boley et al., 2010; Zhu et al., 2012) or even survive and settle onto quasi-stable, wide-separation orbits (Vorobyov, 2013). In model 1, there are no accretion bursts that could be triggered by fragments migrating onto the star (see Fig. 11). There are also no ejections or survival of the fragments. This leaves us only with one outcome – the fragments in model 1 have been tidally dispersed.

### 4.2 The effect of external environment: corotating infall

In this section we consider models 2 and 3 in which a prestellar core is submerged in an external environment corotating with the core. Figure 3 presents the gas surface density maps in model 2 at the same evolution times as in Figure 2 for the isolated core. The top panels have a size of  AU centered on the star, while the other panels have a twice larger size. The early evolution in model 2 ( Myr) is similar to that in model 1— the disk is gravitationally unstable showing a flocculent spiral structure and no signs of fragmentation. This is explained by the fact that during this early period the disk accretes the matter from the parental core rather than from the external environment.

The effect of the external environment becomes evident at  Myr when its matter starts reaching the disk surface. In this phase, the disk evolution in model 2 becomes strikingly different from that in the isolated model 1. Violent gravitational instability and fragmentation leads to almost complete disintegration of the disk. The system now resembles a multi-component stellar system connected with dense filaments, within which individual components tend to form their own mini-disks. The masses of the fragments span a range from a few Jupiters to very-low-mass stars, with the most massive fragment having a mass of .

As the next step, we consider model 3 characterized by the external environment rotating at a ten times smaller angular velocity than that of the core. Figure 4 presents the gas surface density maps in model 3 in the same spatial box of  AU and at similar evolution times as for the isolated model 1 (see Fig. 2). The early evolution at  Myr is similar in both models 1 and 3. The disk is weakly gravitationally unstable and develops a flocculent spiral structure. A significant difference in the disk appearance between models 1 and 3 is seen after  Myr, when infall from the external environment triggers a vigourous gravitational instability and fragmentation in the disk. Unlike model 1 however, the process of disk fragmentation in model 3 seems to be continuous (in the time interval of  Myr) and is characterized by the presence of multiple fragments. This is explained by the fact that the time period between disk fragmentation episodes  yr, calculated using equation (10) and typical disk masses and infall rates ( and  yr), appears to be comparable to the longest migration times of fragments onto the central protostar (Vorobyov et al., 2013a). As a net result, some of the fragments are always present in the disk. The masses of the fragments range from 3.0  to 15 .

The phase of vigourous disk fragmentation ends after  Myr and no fragments are seen in the disk afterwards. The disk shrinks notably in size and its subsequent evolution is characterized by the presence of a weak flocculent spiral structure. This drastic change in the disk evolution as compared to that in model 2 (the latter showing vigorous fragmentation at these evolutionary times) can be attributed to low angular momentum of the infalling external environment, which exerts a negative torque onto the disk. This causes the disk to shrink in size. Although the disk mass remains substantial at this phase (see Section 5), a small size ( AU) and high disk temperature ( K) make gravitational fragmentation unlikely.

### 4.3 The effect of external environment: counterrotating infall

In this section we consider the effect of infall from the external environment counterrotating with respect to the prestellar core. As the first step, we consider model 4 with external environment counterrotating (with respect to the core) at an angular velocity ten times smaller than that of the core. Figure 5 presents the gas surface density maps for model 4 for the same regions and evolution times as in the isolated model 1. The early evolution of the disk in model 4 ( Myr) bears some similarity with that in the isolated model 1—the disk rotates counterclockwise, it is gravitationally unstable and shows clear signs of fragmentation. The fragments have masses around 4.0 . The drastic change with respect to all models considered so far is seen at  Myr when the disk starts quickly shrinking in size and virtually disappears after  Myr due to a strong negative torque exerted by the infalling external environment conterrotating with respect to the disk. The disk looses its rotational support and accretes onto the star, leading to a steep increase in the stellar mass as evident in Figure 10. This example demonstrates how the counterrotating external environment can significantly shorten the lifetime of circumstellar disks.

As the second step, we consider model 5 with external environment counterrotating (with respect to the core) with the same angular velocity at that of the core. Figure 6 presents the gas surface density in model 5 in the inner  AU box (top row) and in the inner  AU box (middle and bottom rows). The early evolution ( Myr) in models 4 and 5 is very similar – the disk in both models is gravitationally unstable and shows signs of fragmentation. The masses of the fragments lie in the 4.0–7.0  range.

The subsequent evolution of the disk in model 5 is however quite different from any model considered so far. First, the disk shrinks in size to just a few tens AU at  Myr due to a strong negative torque exerted by the counterotating infalling material. Then, an outer disk forms rotating in the opposite direction with respect to that of the inner, heavily reduced disk. During this transformation, a notable cavity develops in the gas surface density at a radial distance of the order of a few tens of AU. This process is illustrated in Figure 7 showing the gas velocity field (yellow arrows) superimposed on the gas surface density in the inner  AU box during a narrow time period of 8 kyr. The infall of counterrotating envelope onto the disk is already evident at  Myr, leading to a sharp reduction in the disk size by  Myr. The outer counterrotating disk starts forming at t=0.181 Myr and becomes fully developed by t=0.183 Myr. The inner edge of the cavity is initially located at  AU.

The time evolution of the cavity can be best described with a series of azimuthally averaged gas surface density () profiles shown in Figure 8. The radial profile of in the early evolution is characterized by a near-constant-density central region ( AU) and a steeply declining tail at larger radii. The cavity appears at  Myr and become fully developed by  Myr, with the maximum contrast in between the lowest value in the cavity and the maximum values in the inner and outer disks on the order of 900 and 40, respectively. As time progresses, the cavity moves closer to the central star because the matter from the outer disk flows through the cavity onto the inner disk, exerting a negative torque on to the latter. By  Myr the cavity turns into a central hole when the inner disk becomes smaller than the size of our sink cell. Finally, by  Myr the central hole fills in and the subsequent profile of becomes typical for massive, strongly unstable disks with local sharp maxima representing fragments embedded in the disk.

After t=0.2 Myr the disk experiences a second vigorous episode of disk fragmentation leading in the end to its almost complete disintegration and formation of a multiple system of fragments connected with dense filaments. The masses of the fragments at this stage reach a much higher value, ranging from 7.0 to 82 Jupiter masses.

The difference in the properties of circumstellar disks formed in the considered models can be understood by analyzing the centrifugal radius of matter initially located at a distance from the star:

 Rcf=J2(r)GM(r), (11)

where is the specific angular momentum at a radial distance , is the gravitational constant, and it the mass enclosed within distance . Figure 9 presents the centrifugal radius as a function of distance for various initial distributions of mass and angular momentum in models 1–5. More specifically, the thick black line shows for the isolated pre-stellar core in model 1, while the solid red and dashed red lines provide the centrifugal radius in models 2 and 3 characterized by a corotating external environment. The solid blue and dashed blue lines show for models 4 and 5 characterized by a counterrotating external environment.

The centrifugal radius in model 1 is gradually increasing with distance and reaches a maximum value of 132 AU for the material located at the core outer edge (shown by the arrow). We note however that the disk in model 1 has in fact grown to a size greater than suggested by the maximum centrifugal radius of  AU (see Fig. 2) due to viscous spreading.

The radial profile of in model 2 (solid red line) is similar to that of model 1, except that the matter of the external environment is characterized by a gradually increasing reaching a maximum value of 791 AU. As a net result, the disk in model 2 grows to a larger size and becomes more massive than in model 1 (see Fig. 10), breaking finally into a series of fragments due to violent gravitational instability. In model 3, however, the external environment has a much smaller angular velocity than that of the core. As a consequence, the centrifugal radius of the external environment gradually declines to just several AU (dashed red line). The infalling material from the external environment with a substentially sub-Keplerian rotation exerts an ever increasing negative torque onto the disk outer regions. This causes the disk to shrink in size with time as is evident in Figures 4 and 10.

The last two models 4 and 5 are characterized by external environment rotating in the opposite direction to that of the core. In model 4, characterized by angular velocity of the external environment 10 times smaller than that of the core, the centrifugal radius first grows with distance (as in model 1) but then shows a sharp drop near the core outer edge where the angular velocity changes its sign. of the external material always stays below several AU. As a result, the infalling external material with a counterrotating sub-Keplerian velocity exerts a strong negative torque onto the disk outer regions. The disk shrinks to a size smaller than the size of our sink cell (5 AU) and virtually disappears in Figure 5.

In model 5 the angular velocity of the counterrotating external environment is the same as that of the core. As a consequence, first drops to a small value near the core outer edge (where the rotation changes its sign) but then quickly grows to a maximum value of 685 AU. This causes the disk first to shrink in size, when the low- material from the external environment hits the disk and starts extracting disk’s angular momentum. Then, the infalling external material characterized by an ever growing hits the centrifugal barrier just outside the heavily reduced inner disk and the outer, counterrotating disks begins to form. The transition region between the inner and outer disks, where rotation changes its direction and the matter lacks centrifugal support, is manifested by a density gap clearly visible at  Myr in Figure 6.

## 5 Comparison of disk properties

In this section, we compare various properties of our model disks including the disk masses and sizes, which were computed using the disk tracking algorithm. The disk mass and radius are calculated at each time step by disentangling the disk and infalling envelope on the computational mesh. We do this in practice by adopting a surface density threshold of =0.5 g cm between the disk and envelope and also using the radial gas velocity profile. This method is described in detail by Vorobyov (2011, see fig. 1 and the pertaining text). The values of the disk radius and mass both depend somewhat on the adopted threshold. Nevertheless, a lower value of =0.1 g cm results only in a few percent increase in the net disk mass (see Fig. 2 in Dunham et al., 2014).

We found that our algorithm generally works well, but may somewhat overestimate the disk mass in the early Class 0 stage (when the inner part of the infalling envelope may have densities exceeding ) and somewhat underestimate the disk mass in the late Class I and Class II stages (when the disk spreads out and its density drops below in the outer parts). The envelope mass includes both the matter in the core and external environment (when present).

The left column in Figure 10 presents the integrated disk masses (solid lines), stellar masses (dashed lines) and envelope masses (dash-dotted lines) as a function of time passed since the formation of the central protostar, while the right column shows the corresponding disk radii . Each row corresponds to a particular model as indicated in the figure. In model 1, the disk mass steadily grows and reaches a maximum value of at  Myr, gradually declining afterwards. The growth rate of the protostar exceeds that of the disk and the final mass of the star is expected to be about . The envelope mass steadily decreases and the model enters the class II stage, defined as the time when less than 10% of the initial mass reservoir is left in the envelope, at  Myr. The disk grows in radius showing small radial pulsations and reaches a maximum value of  AU at  Myr. All in all, this model demonstrates a standard behaviour as has been found in previous numerical simulations of isolated collapsing cores (e.g. Vorobyov, 2011).

The properties of the forming star plus disk system in model 2 are strikingly different from those of model 1. First, the disk mass reaches a maximum value of , three times greater than that in model 1. Second, the disk radius exhibits large variations, which are indicative of significant radial motions within the disk caused by the gravitational interaction between the fragments and elements of spiral arms3. Third, the central star experiences several episodes of sharp increase in mass caused by accretion of infalling fragments, the effect discussed in more details in Section 6 below. The final stellar mass is expected to be around . After  Myr, there is still 30% of the total mass available for accretion onto the disk+star system, indicating that the system is still in the embedded phase. We note that the mean duration of the embedded phase obtained in numerical hydrodynamics simulations of collapsing isolated cores ( Myr, albeit with a large scatter depending on the core mass) appears to be shorter by a factor of 3–4 than that obtained by analyzing the number of embedded and Class II sources in star-forming regions (Vorobyov, 2010; Dunham & Vorobyov, 2012). A prolonged accretion from the external environment may provide a possible solution to this problem.

Model 3 with a slowly corotating external environment is peculiar in several aspects. On the one hand, the disk radius stays around 100 AU during most of the evolution, exceeding 200 AU only for a limited period. On the other hand, the disk mass becomes at least a factor of 2 greater than that in model 1, indicating the formation of a dense and compact disk. Due to an elevated mass transport through the disk caused by infall of low-angular-momentum material from the external environment, the final stellar mass is likely to exceed , forming in the long run a substantially more massive object than in models 1 and 2. The envelope mass after  Myr is 5% of the total available mass, indicating that the system has entered the Class II stage. To summarize, model 3 is capable of forming compact and dense disks around solar-mass stars.

Model 4 with a slowly counterrotating external environment is unique in its kind. The negative torque exerted onto the disk by infalling counterotating material leads to complete accretion of the disk onto the protostar by  Myr. A small circumstellar disk may still exist on sub-AU scales, as suggested by the centrifugal radius in Figure 9, but we cannot confirm its presence in our numerical simulations due to the use of a sink cell at  AU. In the subsequent evolution, the forming star is likely to accrete material either from a small sub-AU disk or directly from the infalling envelope. The fraction of stars with detectable near-infrared disk emission seems to approach 80-90% for objects with an age smaller than 1.0 Myr (Hernandez et al., 2007), suggesting that systems like those explored in model 4 are rare.

Finally, model 5 presents a very curious case showing the reversal of the disk rotation due to the infall of counterrotating external material. The initial disk growth is followed by a sharp decline when the disk mass and size drop to just and  AU at  Myr, manifesting the beginning of the disk reversal phase. In the subsequent stage, the disk evolution is similar to that in model 2 showing vigorous gravitational instability and fragmentation. The final stellar mass is expected to be slightly below that of the Sun. The formation of systems with disks counterrotating to that of the star is rare, with RW Aur A being one possible example showing both the disk and jet counterrotating with respect to the rotation of the host binary system (Woitas et al., 2005; Bisikalo et al., 2012).

We want to mention that the gas temperature in models 3–5, which feature compact disks during the evolution, becomes comparable to or higher than 200 K at =10-20 AU, a value which is at least a factor of two greater than that in the isolated model 1. We cannot calculate the gas temperature at the sub-AU scales, but by extrapolating the profile we speculate that the temperature at =0.5 AU can reach values on the order of 1200 K, which is close to what is needed for the formation of calcium-aluminium-rich inclusions (CAIs). Numerical simulations with a smaller sink cell are planned for the future to investigate this possibility.

Table 2 summarizes the main properties of the forming objects in models 1-5. The mean values are calculated starting from the formation of the disk and ending at  Myr after the formation of the protostar. In model 4, the mean values were taken over the time period covering the actual existence of the disk.

## 6 Accretion and infall rates

In this section we compare the mass accretion rates through the sink cell and mass infall rates at a distance of 2000 AU from the star in our five models. Figure 11 presents (solid lines) and (red dashed lines) as a function of time elapsed since the formation of the first hydrostatic core (FHC) in models 1–5 (from top to bottom).

The accretion rate in model 1 shows a time behavior typical for weakly gravitationally unstable disks (Vorobyov & Basu, 2010). is negligible in the pre-stellar phase and rises to a maximum value of  yr when the FHC forms at t = 0 Myr. The subsequent short period of evolution is characterized by a gradually declining when the material from the core lands directly onto the forming star. A sharp drop in the accretion rate at  Myr manifests the beginning of the disk formation phase when the infalling core material hits the centrifugal barrier near the sink cell and the accretion rate drops to a negligible value. However, the process of mass loading from the core continues, the disk grows in mass and size, and the disk quickly becomes gravitationally unstable. In the subsequent evolution shows order-of-magnitude flickering and gradually declines with time to a value of  yr at  Myr. We note that model 1 exhibits no accretion bursts with magnitude  yr, as might be expected for models with fragmenting disks whereby fragments migrate onto the star due to gravitational interaction with spiral arms (Vorobyov & Basu, 2006, 2010). This indicates that the fragments in model 1 (see Fig. 2) were tidally destroyed before reaching the inner sink cell. We also note that in model 1 exceeds only during the early evolution of the disk ( Myr) and both values are similar in magnitude in the subsequent evolution. The inequality is one of the conditions for disk gravitational fragmentation (e.g. Vorobyov & Basu, 2010; Kratter et al., 2010), explaining the presence of fragments in the disk of model 1.

The early evolution of in model 2 is similar to that of model 1. A notable difference comes about after  Myr when begins to show variations by several orders of magnitude. Accretion bursts, some of which exceeding in magnitude  yr, are caused by disk gravitational fragmentation followed by migration of the fragments onto the forming protostar – a phenomenon known as the burst mode of accretion (Vorobyov & Basu, 2006, 2010). The number of bursts in this particular model amounts to a few tens and the burst activity subsides in the early Class II stage. The mass infall rate is greater on average than , thus leading to the accumulation of mass in the disk and strengthening its gravitational instability and propensity to fragment. To check this supposition, we calculated the spatial distribution of the Gammie parameter , where is the local cooling time, and the Toomre parameter . Disk fragmentation is supposed to occur in regions where both quantities become smaller than unity (e.g. Gammie, 2001; Rice et al., 2003). We found that during the initial disk evolution, unaffected by infall from the external environment, the -parameter becomes smaller than unity only in a few spots in the disk, while the regions with cover a much larger area, indicating that the insufficient strength of gravitational instability (rather than slow disk cooling) is the main reason hindering fragmentation. On the other hand, during the external infall stage, the regions with occupy a notably larger area, while the spatial area with remains similar to what was found during the early stage, meaning that the infall from the external environment is the dominant effect driving the disk to the fragmentation boundary.

We note that sometimes shows large variations caused by the perturbing influence of the disk. More specifically, the gravitational interaction between fragments and spiral arms within the disk sometimes leads to ejections of individual fragments or elements of spiral arms to distances on the order of thousands AU.

The mass accretion rate in model 3 with a slowly corotating external environment exhibits a much smaller accretion variability than that of model 2. There are only two strong accretion bursts around  Myr and the later evolution of is characterized only by small-scale flickering. On average, is higher in model 3 than in models 1 and 2, staying during the initial 0.5 Myr of evolution in the  yr range. The matter from the external environment falling in onto the disk at a substantially sub-Keplerian speed acts to reduce the disk rotational support, thus increasing the mass flow through the disk onto the star. We note that is greater than at  Myr, causing disk fragmentation in the early disk evolution (see Fig. 4), but both values become comparable in the later evolution, indicating the presence of efficient mass transport through the compact and hot disk due to the increasing strength of viscous torques.

The effect of the external environment in model 4 is manifested by a gradual increase of the mass accretion rate after  Myr to  yr (on average). This increase is even more pronounced than in model 3 due to the fact that the infalling external material in model 4 is counterrotating with respect to the disk, thus exerting a stronger negative torque onto the disk outer regions than in model 3. As a net result, the disk loses rotational support and the mass transport rate through the disk increases leading to a complete disappearance of the disk after  Myr. In the subsequent evolution, both and are nearly identical and slowly decline with time.

Finally, the early evolution of in model 5 is very similar to that in model 4. After an initial period of near constant accretion at a value of  yr, the mass accretion rate shows a transient increase by an order of magnitude during  Myr. During this period, the matter is accreted from the inner disk, which is gradually losing its rotational support and shrinking in size due to the infall of matter from the counterrotating external environment. At  Myr, the inner disk reduces to a size smaller than the sink cell and drops substantially. In the subsequent evolution, the matter is accreted from the outer, counterrotating disk. Because the infall rate is greater on average than the mass accretion rate onto the star , the outer disk quickly becomes gravitationally unstable, leading high-amplitude variations in . It is interesting to note that the later evolution shows only one moderate accretion burst at  Myr, implying that fragments forming in the disk (see Fig. 6) are either tidally dispersed or may stay in the system with a possibility of forming in the long run a multi-component stellar/planetary system.

## 7 Discussion: the variety of disk properties

Observations of star-forming regions indicate that the fraction of (sub-)solar-mass stars with circumstellar disks approaches 80%-90% for objects with an age equal to or smaller than 1.0 Myr (Hernandez et al., 2007). The young disks seem to have diverse properties with masses ranging from a few Jupiters to a fraction of solar mass and sizes from a few to hundreds AU (Andrews et al., 2009; Jorgensen et al., 2009; Eisner, 2012; Tobin et al., 2012).

The observed variety of disk properties can in principle be reproduced by variations in the mass, angular momentum, and magnetic field strength of isolated parental cores. For instance, numerical hydrodynamics simulations of the gravitational collapse of non-magnetized isolated cores (Vorobyov, 2011) yield the following near-linear relation between the disk masses and the masses of the central sub-stellar/low-mass object

 Md,CO = (0.73+0.11−0.09)M1.05±0.07∗,CO, (12) Md,CI = (0.65+0.04−0.05)M1.0±0.04∗,CO, (13)

where the indices CO and CI correspond to the class 0 and class I phase, respectively. The minimum and maximum disk masses span a range from to . At the same time, observations of embedded sources infer disks with the lower limit extending to just one Jupiter mass (Jorgensen et al., 2009; Tobin et al., 2012) and, in some cases, fail to detect disks at all (Maury et al., 2010). Systems with very-low-mass disks () can be the result of efficient magnetic braking operating in the embedded phase of star formation (Machida et al., 2011b; Seifried et al., 2012) or can simply be formed from cores with very low angular momentum and, subsequently, having a very small centrifugal radius. The situation is similar for disk radii – numerical simulations of isolated non-magnetized cores (Vorobyov, 2011) seem to yield disk sizes that are a factor of several greater than inferred from observations (Tobin et al., 2012; Eisner, 2012).

In this paper, we have explored the evolution of pre-stellar cores submerged into an external environment and found that its effect can add another dimension to the variety of disk properties around (sub-)solar-mass stars. Depending on the value of angular velocity and the direction of rotation of the external environment with respect to the core the resulting disks can be characterized by vastly different physical parameters such as masses, sizes, lifetimes, and accretion rates onto the host star. The external environment can affect the strength of gravitational instability in the disk and hence the likelihood of giant planet and brown dwarf formation. The effect can actually work both ways: suppressing or promoting disk gravitational fragmentation depending on the amount of rotation of the external material.

Models 3 and 4 with a slowly rotating external environment predict the formation of compact disks with radii  AU and, in some cases, even smaller than 5 AU (the size of our central sink cell). Even model 2 with a fast rotating environment is characterized by a more compact disk (than the isolated model 1) in the late evolution, which is likely explained by the loss of angular momentum carried away by fragments accreted onto the star4. Observational estimates of disk radii seem to indicate that the majority of embedded and T Tauri disks are indeed rather compact,  AU (Vicente & Alves, 2005; Tobin et al., 2012; Murillo et al., 2013).

Perhaps, the most interesting finding of our numerical modelling is the ability to form disks with cavities or holes in the case of counterrotating external environment. These features are transient with the lifetime of a few tens kyr and are characterized by conterrotating inner and outer disks. In our model, the gap is located at  AU, which is similar with the position of the gap in, e.g., LkH 330 (Brown et al., 2008). The gap migrates inward with time gradually transforming into a central hole. Though being a short-lived phenomenon, this mechanism of the gap/hole formation can present an interesting alternative to other mechanisms such as disk photoevaporation by the central star (Clarke et al., 2001; Owen et al., 2010) or gravitational clearing due to a massive planet (e.g. Kley et al., 2001).

After the inner, corotating disk vanishes due to accretion of its material onto the star, a stellar system with a disk counterotating with respect to that of the star emerges. If the disk is sufficiently massive to fragment, this process may lead to the formation of a system with sub-stellar/very-low-mass components counterrotating with respect to the host star. This phenomenon may, however, be transient in nature. The matter is that angular momentum in the external environment is usually much greater than that of the star and the latter may change stellar rotation to match that of the infalling material. Although the efficiency of this process should depend on how much angular momentum is lost via jets/outflows and whether or not a binary/multiple stellar system forms instead of a single star, such counterotating systems are nevertheless expected to be rare, with RW Aur A being one possible example (Woitas et al., 2005; Bisikalo et al., 2012).

We note that we have considered a rather simplified initial setup, treating the external environment as a homogeneous medium of constant density and angular velocity. In reality, however, the external environment is likely to have a complicated velocity and density structure as often seen in numerical hydrodynamics simulations of clustered star formation. For instance Bate et al. (2010), Fielding et al. (2014), and Padoan et al. (2014) demonstrated that individual cores may accrete mass and angular momentum at a highly non-steady rate due to a filamentary and chaotic structure of the intracluster medium, leading to the formation of misaligned star-disk systems (the misalignment in some transient cases exceeds ). It is therefore possible that the angular momentum vector of the accreted material can undergo significant changes, including its orientation with respect to that of the core. Moreover, the dynamics of dust particles needs to be considered as well because the holes/gaps in the disk are often identified through the lack of dust emission there. Finally, we acknowledge that the external mass in our models () is greater than that of the core (0.63 ), which may have amplified the effect of infall from the external environment. Models with a smaller need to be also investigated in order to evaluate the dependence of our results on the available mass in the external reservoir. We plan to consider more realistic initial configurations, smaller external mass, and the dynamics of dust in a subsequent study.

## 8 Conclusions

In this paper we have performed numerical hydrodynamics simulations of the gravitational collapse of a pre-stellar core with sub-solar mass submerged into a low-density external environment characterized by different magnitude and direction of rotation with respect to the core. We followed the evolution of our models through the star and disk formation stage and terminated the simulations when the age of the forming star started to exceed 0.5 Myr. In all models, the mass of the core and external environment were fixed at and , respectively, in order to exclude systematic differences due to variations in these parameters. We found that depending on the value of angular velocity and the direction of rotation of the external environment with respect to the core the resulting disks can be characterized by vastly different physical properties.

Our most interesting finding is that infall of material from the external environment counterrotating with respect to the core may lead to the formation of the inner and outer disks rotating in the opposite directions to each other. The transitional region between the disks manifests itself as a deep gap in the gas density, which migrates inward due to accretion of the inner disk onto the protostar. After about several tens of kyr, the gap turns into a central hole, which fills in after another few tens of kyr. The formation of a transient gap/hole due to infall from counterrotating external environment presents an interesting alternative to other mechanisms of gap/hole formation such as disk photoevaporation (Clarke et al., 2001; Owen et al., 2010) or gravitational clearing due to a massive planet (e.g. Kley et al., 2001). Our main results can be summarized as follows.

• Pre-stellar cores submerged into an external environment corotating at an angular velocity similar in magnitude to that of the core tend to form extended and strongly unstable disks, which fragment into several objects connected by dense filaments and may form in the long run systems with sub-stellar/very-low-mass companions. The mass accretion rate onto the primary shows a highly variable character with episodic bursts caused by the infall of some of the fragments onto the star.

• Pre-stellar cores embedded into a slowly corotating external environment (with angular velocity ten times smaller than that of the core) are likely to form compact ( AU) disks, stable against gravitational fragmentation due to their small size. These disks can drive an order of magnitude higher mass accretion rates onto the star than those formed from isolated cores.

• Pre-stellar cores submerged into a slowly counterrotating external environment are expected to form compact ( AU) and short-lived ( yr) disks due to strong negative torques exerted by the infalling external material.

• Disks formed from pre-stellar cores embedded into an external environment counterotating with respect to the core at a similar angular velocity are likely to change their initial direction of rotation. In the process of this transformation, the inner and outer disks form counterotating with respect to each other and separated by a transient gap in the gas density. In the long run, the inner disk accretes onto the star and a transient stellar system with sub-stellar/very-low-mass components counterrotating with respect to that of the star may emerge due to gravitational fragmentation of the outer, strongly unstable disk.

We thank the referee for constructive comments that helped to improve the manuscript. This project was partly supported by the Russian Ministry of Education and Science Grant (state assignment) 3.961.2014/K. MG acknowledges support by the Austrian FWF through the NFN project grant S116 Pathways to Habitability: From Disks to Active Stars, Planets and Life, and the related subproject S116 604-N16 Radiation & Wind Evolution from T Tauri Phase to ZAMS and Beyond. The simulations were performed on the Shared Hierarchical Academic Research Computing Network (SHARCNET), on the Atlantic Computational Excellence Network (ACEnet), and on the Vienna Scientific Cluster (VSC-2). This publication is supported by the Austrian Science Fund (FWF).

### Footnotes

1. In the case that the gas temperature exceeds 2000 K, is set to 1.1, but a limited numerical resolution does not allow us to follow the second collapse of hot fragments
2. This definition of the midplane temperature is accurate within a factor of unity (Zhu et al., 2012)
3. These numbers should be taken with caution because our disk tracking algorithm may be inaccurate on systems with violently fragmenting disks.
4. We have not calculated the angular momentum flux through the sink cell. This study will be presented in a follow-up paper.

### References

1. Andrews, S. M., Wilner, D. J., Huges, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502
2. Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 2002, A&A, 382, 563
3. Baraffe, I., & Chabrier, G. 2010, A&A, 521, 44
4. Baraffe, I., Vorobyov, E. I., & Chabrier, G. 2012, ApJ, 756, 118
5. Basu, S., & Mouschovias, T. Ch. 1994, ApJ, 432, 720
6. Basu, S., & Vorobyov, E. I. 2012, ApJ, 750, 30
7. Bate, M., Lodato, G., Pringle, J. E., 2014, MNRAS, 401, 1505
8. Beccari, G., Spezzi, L., De Marchi, G., et al. 2010, ApJ, 720, 1108
9. Bell, K. R., Lin, D. N. C., 1994, ApJ, 427, 987
10. Bisikalo, D. V., Dodin, A. V., Kaigorodov, P. V., Lamzin, S. A., Malogolovets, E. V., Fateeva, A. M. 2012, Astron. Reports, 56, 686
11. Boley, A. C., Hayfield, T., Mayer, L., & Durisen, R. H. 2010, Icarus, 207, 509
12. Brown, J. M., Blake, G. A., Qi, C., Dullemond, C. P., & Wilner, D. J., ApJL, 675, 109
13. Cha S.-H., Nayakshin S., 2011, MNRAS, 415, 3319
14. Chabrier, G., Baraffe I. 2000, ARA&A, 38, 337
15. Clarke C. J., Gendrin A., & Sotomayor M., 2001, MNRAS, 328, 485
16. Dapp, W. B., & Basu, S. 2009, MNRAS, 395, 1092
17. De Marchi, G., Beccari, G., & Panagia, N. 2013a, ApJ, 775, 68
18. De Marchi, G., Panagia, N., Guarcello, M. G., & Bonito, R. 2013b, MNRAS, 435, 3058
19. Dunham, M. M., Vorobyov, E. I., Arce, H. G. 2014, MNRAS, 444, 887
20. Dunham, M. M., & Vorobyov, E. I. 2012, ApJ, 747, 52
21. Dunham, S., et al. in Protostars & Planets VI, 2014, arXiv1401.1809
22. Eisner, J. A. 2012, ApJ, 755, 23
23. Evans, N. J., II, Dunham, M. M., Jorgensen, J. K., et al. 2009, ApJS, 181, 321
24. Fielding, D. B., McKee, C. F., Socratis, A., Cunningham, A. J., Klein, R. I. 2014, arXiv:1409:5148
25. Gammie, C. F. 2001, ApJ, 553, 174
26. Hernandez, J., Hartmann, L., Megeath, T., et al. 2007, ApJ, 662, 1067
27. Hacar, A. & Tafalla, M. 2011, A&A, 533, 34
28. Hacar, A., Tafalla, M., Kauffmann, J., Kovacs, A. 2013, A&A, 554, 55
29. Johnson, B. M., & Gammie C. F. 2003, ApJ, 597, 131
30. Jorgensen, J. K., van Dishoeck, E. F., Visser, R., Bourke, T. L., Wilner, D. J., Lommen, D., Hogerheijde, M. R., & Myers, P. C. 2009, A&A, 507, 861
31. Kley, W., D’Angelo, G., Henning, T. 2001, ApJ, 547, 457
32. Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 708, 1585
33. Masunaga, H., & Inutsuka, S.-I. 2000, ApJ, 531, 350
34. Machida, M. N., Inutsuka, S., & Matsumoto, T. 2011, ApJ, 729, 42
35. Machida, M. N., Inutsuka, S., & Matsumoto, T. 2011, PASJ, 63, 555
36. Maury A. J., Andre, Ph., Hennebelle, P., Motte, F., Stamatellos, D., Bate, M., Belloche, A., Duchene, G., & Whitworth, A. 2010, A&A, 512, 40
37. McKee, C. F., & Offner S. S. R. 2010, ApJ, 716, 167
38. Meru, F., & Bate, M. R. 2012, MNRAS, 427, 2022
39. Murillo, N. M., Lai, S.-P., Bruderer, S., Harsono, D., & van Dishoeck, E. F. 2013, A&A, 560, 103
40. Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415
41. Padoan, P. Juvela, M. Goodman, A. A., Nordlund, A. 2001, ApJ, 553, 227
42. Padoan, P., Haugbolle, T., Nordlung A. 2014, arXiv:1407.1452
43. Rice, W. K. M., Armitage, P. J., Bate, M. R., & Bonnell, I. A. 2003, MNRAS, 339, 1025
44. Rice, W. K. M., Mayo, J. N., & Armitage, P. J. 2010, MNRAS, 402, 1740
45. Scicluna, P., Rosotti, G. Dale, J. E. & Testi, L. 2014, A&A, 566, 3
46. Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2012, MNRAS, 423, L40
47. Shu, F. H., Adams, F. C., Lizano, S. 1987, ARA&A, 25, 23
48. Stamatellos, D., & Whitworth, A. P. 2008, A&A 480, 879
49. Stamatellos D., Maury A., Whitworth A., André P., 2011, MNRAS, 413, 1787
50. Tobin, J. J., et al. 2012, Nature, 492, 83
51. Vicente, S. M., & Alves, J. 2005, A&A, 441, 195
52. Vorobyov, E. I., Basu, S., 2006, ApJ, 650, 956
53. Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896
54. Vorobyov, E. I. 2010, ApJ, 713, 1059
55. Vorobyov, E. I. 2011, ApJ, 729, 146
56. Vorobyov, E. I. 2013, A&A, 552, 129
57. Vorobyov, E. I., Zakhozhay, O., Dunham, M. M. 2013, MNRAS, 433, 3256
58. Vorobyov, E. I., Baraffe, I., Harries, T., Chabrier, G. 2013, A&A, 557, 35
59. Woitas, J., Bacciotti, F., Ray, T. P. et al. 2005, A&A, 432, 149
60. Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F. 2012, ApJ, 746, 110
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters