Spiraldriven accretion in protoplanetary discs
I. 2D models
^{1} Univ. Grenoble Alpes, IPAG, 38000, Grenoble, France
email: geoffroy.lesur@ujfgrenoble.fr
^{2} CNRS, IPAG, 38000 Grenoble, France
^{3} Laboratoire AIM, CEA/DSM–CNRS–Université Paris 7, Irfu/Service d’Astrophysique, CEASaclay, 91191 GifsurYvette, France
Received: 12 June 2015
Accepted: 12 September 2015
We numerically investigate the dynamics of a 2D nonmagnetised protoplanetary disc surrounded by an inflow coming from an external envelope. We find that the accretion shock between the disc and the inflow is unstable, leading to the generation of largeamplitude spiral density waves. These spiral waves propagate over long distances, down to radii at least ten times smaller than the accretion shock radius. We measure spiraldriven outward angular momentum transport with 10^{4} ≲ α < 10^{2} for an inflow accretion rate Ṁ_{inf} ≳ 10^{8} M_{⊙} yr^{1}. We conclude that the interaction of the disc with its envelope leads to longlived spiral density waves and radial angular momentum transport with rates that cannot be neglected in young nonmagnetised protostellar discs.
Key words: accretion, accretion disks / hydrodynamics / waves
© ESO, 2015
1. Introduction
Accretion is an essential phenomenon in astrophysics as it is the way through which gravitational energy is transformed into heat and therefore observable radiation. However, the physical process responsible for accretion is still debated. Although the magnetorotational instability (MRI, Balbus & Hawley 1991) provides an efficient accretion mechanism, its applicability to protoplanetary discs is questionable since these objects are very weakly ionised and might quench the MRI through various nonideal magnetohydrodynamical (MHD) processes (Turner et al. 2014).
The usual approach to protoplanetary disc dynamical modelling is to assume that these objects are isolated in space. In this context, most of the models rely on processes such as MRIdriven turbulence, selfgravity, or winds to drive accretion. However, these discs are not isolated systems, and the question is whether surrounding material, albeit much less dense than the disc, could perturb it sufficiently to drive accretion. This possibility was proposed by Padoan et al. (2005) to explain accretion rates that scale like the central protostar mass squared. This scenario was later refined by Throop & Bally (2008), Klessen & Hennebelle (2010), and Padoan et al. (2014) and interpreted as a consequence of BondiHoyle accretion onto the protostar.
In this Letter, we explore the importance of external accretion using a very simplified model. We consider a 2D hydrodynamic Keplerian disc onto which falls gas coming from a surrounding envelope. The disc is inviscid and hydrodynamically stable so that no accretion can occur without this external inflow. We first present in detail the model and the numerical method used to solve the equations of motion. We then explore some of the results obtained using this model and finally discuss the limitations and implications of our findings.
This work shares some similarities with that of Bae et al. (2015): both consider accretion discs subject to an external inflow using 2D hydrodynamical models. However, we focus here on the transport of angular momentum far from the perturbed region using extended disc models, whereas Bae et al. (2015) concentrated on the local response of the disc in the presence of viscous angular momentum transport. We discuss the differences between the outcome of these models in Sect. 4.
2. Model
2.1. Physical model
We study the dynamics of an accretion disc close to Keplerian rotation. We consider a 2D frame (R,θ) and integrate the equations of motion in the vertical direction. In this initial study, we neglect the role played by magnetic fields and selfgravity to concentrate on a purely hydrodynamical model. The equations of motion read where v is the flow velocity, Σ the surface density, P the vertically averaged pressure, G the gravitational constant, and M_{⊙} the mass of the central object, which is assumed to be a solarmass star.
Instead of the usual minimum mass solar nebula (MMSN) density profile, we use a powerlaw surface density profile with an exponential outer cutoff (3)where we choose Σ_{0} = 1700 g cm^{2} and R_{c} = 40 au, which corresponds to typical protoplanetary discs observed in Ophiuchus (Andrews et al. 2009) with a mass of 4.8 × 10^{2} M_{⊙}. The temperature profile is chosen to be locally isothermal so that c_{s}/v_{K} ≡ ε = 0.1 or 0.05, where c_{s} is the sound speed and is the Keplerian velocity. The equation of state is then simply . Assuming a vertically isothermal hydrostatic profile, this model leads to an accretion disc with a constant aspect ratio H/R = ε, where H is the typical vertical disc scale height.
Fig. 1 Density maps as a function of time for run S71. From left to right: t = 2T_{100}, t = 2.8 T_{100}, t = 6 T_{100}. 

Open with DEXTER 
2.2. Boundary conditions
The outer radial boundary condition is located at R_{out} = 400 au. To mimic matter falling onto the disc, we inject material at the outer radial boundary condition that then falls onto the disc. The material is injected over an azimuthal extent 0 <θ<θ_{inj} with a sonic radial velocity. v_{rout} = −c_{s}. We vary three parameters for the injected material: its specific angular momentum ℒ_{out} ≡ v_{θ}(R_{out})R_{out}, its accretion rate Ṁ_{inf} ≡ −∫dθ Σ_{out}v_{rout}, Σ_{out} being deduced from the desired accretion rate, and the azimuthal width of the inflow θ_{inj}.
The inner radial boundary is located at R_{in} = 1 au. The inner boundary condition is forced to the initial density and velocity profile. To avoid spurious reflexion of spiral density waves, we add a damping zone in the region 1 au <R< 3 au. In this region, we exponentially relax nonaxisymmetric velocity fluctuations on a fixed timescale set to eight local orbital periods. To avoid mass accumulation in the damping zone, we also relax the density to the initial density profile on the same timescale.
2.3. Numerical algorithm
We use Pluto 4.0 (Mignone et al. 2007) to solve the equations of motion on a cylindrical grid, using logspaced grid cells in the radial direction. Parabolic reconstruction is used in each cell to cope with the nonuniform radial grid, and a thirdorder RungeKutta algorithm is used to evolve the system in time. We use the HLL Riemann solver^{1} at cell boundaries combined with orbital advection (Mignone et al. 2012) to allow for long integration times. We have tested that starting with sonic white noise perturbation on v and without any inflow, the disc was hydrodynamically stable with α ~ 10^{7} after a few local orbits, as expected with our density profile (3). The resolution of each run was (N_{R},N_{θ}) = (512,512) for ε = 0.1 or (N_{R},N_{θ}) = (1024,1024) for ε = 0.05, which corresponds to a locally isotropic resolution of eight points per H. We checked that doubling the resolution of run S71 did not quantitatively affect our results. We are therefore confident that our simulations are numerically converged.
2.4. Units, diagnostics, and notations
Length units are given in au, surface densities in g cm^{2}, time units in orbital periods at 100 au = T_{100} (or equivalently in units of 10^{3} yr), and accretion rates in M_{⊙} yr^{1}.
In the following, several diagnostics are used to measure the behaviour of the disc coupled to the inflow. We first introduce two averages, an azimuthal average and a temporal average ⟨·⟩. We start the temporal average from t = 10 T_{100} to allow the system to relax from the initial condition (this corresponds to 10^{4} orbits at the inner boundary).
These averaging procedures allow us to define fluctuating quantities . The first diagnostic is the Shakura & Sunyaev (1973)α parameter . We also use the disc radius R_{d} defined as the location where the disc rotation velocity drops below 90% of the Keplerian velocity v_{K}. This definition might look rather arbitrary, but it nevertheless leads to a welldefined disc radius since the rotation velocity drops very rapidly with radius at the transition region between the disc and the inflow. Each simulation is integrated for 3 × 10^{4} orbits at the inner boundary, which corresponds to 30 T_{100}.
Our runs are labelled “XYZ” where X = S or A denotes symmetric (θ_{inj} = 2π) or asymmetric (θ_{inj} = π/ 2) inflows, Y is the accretion rate at R_{out} and Z is the specific angular momentum ℒ_{out}. Our simulations all have ε = 0.1, except for runs labelled “thin”, which assume ε = 0.05.
3. Results
We first concentrate on symmetric inflows. We then explore the impact of asymmetric inflows and of the disc thickness in a second part.
3.1. Symmetric accretion
3.1.1. Accretion shock
The initial evolution of the symmetric accretion run S71 is presented in Fig. 1. The falling material initially forms an accretion shock that propagates inward. At t ≃ 2 T_{100}, the shock stalls at R = R_{d}. The stationary shock then develops an instability, which appears as a shortwavelength nonaxisymmetric ondulation of the shock front (t = 2.8 T_{100}). This instability then quickly saturates and generates strong spiral waves that propagate inward on the long term (t> 3.5T_{100}).
The origin of this instability is tightly linked to the structure of the shock. Spirals are continuously produced from R ≲ R_{d}, and the instability survives for the entire duration of the simulation. The instability mechanism becomes clear when the average specific angular momentum profile ⟨ ℒ ⟩ and the average potential vorticity are plotted as a function of radius (Fig. 2). We find that at the location of the stationary shock a strong velocity gradient develops, which leads to an extremely peaked ζ. This sharp structure appears because the specific angular momentum of the falling material is lower than that of the disc material at R = 200. A potential vorticity layer is therefore unavoidable if matter keeps falling onto the disc.
As is well known, such a hydrodynamical structure naturally leads to a KelvinHelmholtz instability (KHI, also known as Rossby wave instability (RWI) in the astrophysical context^{2}). We note, however, that the shear rate is locally so strong that it also locally violates the Rayleigh criterion: in the shock region, the specific angular momentum is decreasing outward (Fig. 2). As stated above, this configuration is unavoidable given that the inflow has a lower angular momentum than the disc. The instability driving the spiral is therefore a mixture of KHI (RWI) and Rayleigh centrifugal instability. It has a growth rate close to the orbital frequency at R_{d} and is sustained, for reasons stated above, during the entire duration of the simulation.
Fig. 2 Averaged potential vorticity ⟨ ζ ⟩ (plain blue) and specific angular momentum ℒ (red dashed) of run S71. Note the sharp negative angular momentum gradient associated with a strong potential vorticity peak at R ≃ 200. 

Open with DEXTER 
Fig. 3 Profiles of α as a function of radius, measured from R_{d}. The spiral stress increases with increasing Ṁ_{inf} and with decreasing initial angular momentum ℒ_{out}. 

Open with DEXTER 
3.1.2. Angular momentum transport
Spiral waves generate a positive Reynolds stress through the entire disc, which we quantify as an α parameter (Fig. 3). We find that these spirals produce a very strong transport close to R ≃ R_{d} with α ≃ 1. They propagate inward from R_{d} and dissipate through shocks, reaching α ≳ 10^{3} at R_{d}/ 2 and α ≳ 10^{4} at R_{d}/ 10. A generation region (0.6R_{d}<R<R_{d}) can be defined where the inflow mixes with the disc, spiral modes are excited, and α values are similar for all of our simulations. Likewise, there is a propagation region (R< 0.6R_{d}) where the disc density structure is unaffected by the inflow except for propagating spiral waves.
We see that α in the generation region is not affected by Ṁ_{inf}, but α in the propagation region (R< 0.6R_{d}) clearly is. Moreover, for our strongest Ṁ_{inf} (run S61 for instance), we observe a bump of transport for R ≃ 0.6R_{d}, directly at the transition between these two regions. This is because the mass of the generation region has increased significantly due to the inflow while that of the propagation region has not, forming a jump in surface density at the interface. As a result, the generation region has a higher inertia and excites stronger waves at the interface with the propagation region. We note that smaller angular momentum in the inflow also leads to somewhat stronger spirals. These findings can be summarised by representing α at R = R_{d}/ 10 as a function of Ṁ_{inf} (Fig. 4).
Fig. 4 α measured at R = R_{d}/ 10 as a function of Ṁ_{inf}. Colours represent different ℒ_{out} or H/R, diamonds correspond to symmetric inflow, and stars correspond to asymmetric accretion. The dashed line are fits given by Eq. (4). 

Open with DEXTER 
List of simulations.
From these results, we can estimate that in the symmetric inflow case (4)For inflows with significant rotation (ℒ_{out} = 1) we have α_{0} = 1.3 × 10^{4} and γ = 0.4, whereas for inflows without initial rotation we obtain α_{0} = 3.0 × 10^{4} and γ = 0.5. These values constitute a lower bound for the average alpha found in the disc. Typical values of α at R = R_{d}/ 2 are 10 to 100 times higher than these estimates (Table 1).
As is well known, standard accretion disc theory allows deriving a theoretical mass accretion rate from the profile of α computed above (Balbus & Papaloizou 1999, Eq. (28)). We have checked that the measured accretion rates are consistent with this theoretical expectation, provided that Ṁ is computed on long enough averages (typically more than ten local orbital periods). This shows that standard accretion disc theory also applies to spiraldriven transport of angular momentum. This result does not imply, however, that “turbulent” heating follows Shakura & Sunyaev (1973)αdisc model (Balbus & Papaloizou 1999). This point cannot be tested in our model since we used a locally isothermal equation of state.
3.2. Parameter survey
Since spiral waves are primarily sound waves, a dependence on the sound speed and therefore on H/R = ε is expected. To test this effect, we reduce H/R by a factor 2 with ε = 0.05. These simulations, labelled “thin”, show that the amplitude of the spirals is significantly reduced when H/R is reduced. We find a typical attenuation of a factor 3 between ε = 0.1 and ε = 0.05. This is similar to the case of tidally induced spiral waves (Savonije et al. 1994), despite an excitation mechanism different from ours. This similitude indicates that the wave propagation is most likely mainly responsible for this dependence.
All the results presented above were obtained in a situation where gas is falling onto the disc in an axisymmetric way. To test the effect of this assumed symmetry on the results, we consider the case of accretion coming from a narrower stream (runs labelled A*). We find that the disc shape is largely unaffected for the inflow rates we used and produces spiral waves similar to the one found in the symmetric case (Fig. 1). In particular, we do not find any nonaxisymmetric feature close to the location where the stream hits the disc. In addition, the azimuthally averaged profiles of ⟨ ζ ⟩ and ⟨ ℒ ⟩ are very similar to the one found in the symmetric accretion case. The average flow is still unstable to the same instability as in the symmetric case. It is therefore expected to find that asymmetric accretion leads to α(R) profiles that are qualitatively comparable to those found in the symmetric case.
By comparing α(r) obtained in the asymmetric case to the symmetric case, we find that asymmetric accretion in general drives a slightly more efficient angular momentum transport, probably thanks to the excitation of low m modes of the instability by the stream that propagates more easily. At R = R_{d}/ 10, we obtain an α about two to four times larger than the one obtained from symmetric accretion simulations (Fig. 4). This suggests that the symmetric accretion scenario given by Eq. (4) constitutes a lower bound to spiraldriven angular momentum transport.
4. Discussion
We explored the impact of an external inflow on a protoplanetary disc using 2D hydrodynamical simulations. We found that this interaction leads to the generation of strong spiral waves that propagate on long distances (typically down to radii smaller that R_{d}/ 10). The resulting α at R_{d}/ 2 are larger than 10^{3} and reach a few 10^{4} at R_{d}/ 10 for inflow mass rates higher than 10^{8} M_{⊙} yr^{1}. This stress varies strongly with Ṁ_{inf} and H/R, but is not significantly affected by the geometry of the inflow.
As shown in the introduction, this study is similar in nature to the work of Bae et al. (2015). However, several important differences should be emphasised. First, our simulations only assumed a radial inflow and no vertical inflow. Second, the angular momentum of the falling material is always smaller than that of the outer disc so that a radial shear layer is always present, a possibility only explored in one simulation by Bae et al. (2015). Because of these differences, we do not observe the formation of vortices, most probably because the instability at R_{d} involves a mixture of the RWI and of Rayleigh centrifugal instability thanks to the strong shear layer. We also obtained higher α values than did Bae et al. (2015) in the generation region, probably due to the same shear layer.
Finally, spiral waves such as we discussed here are now detectable with infrared continuum observations. Several authors have already reported direct observations of spiral patterns at the outer edge of evolved protostellar discs (e.g. Muto et al. 2012; Benisty et al. 2015). Whether an external inflow might be at the origin of these spiral pattern is an open question that we defer to a future publication.
We also ran simulations with the HLLC Riemann solver. In this case, however, the noise level of the background spirals in the absence of inflow is α ~ 10^{4} because there is almost no numerical dissipation. We therefore only present results obtained with HLL for which the background noise level is much lower.
In our case, the RWI is not due to a density bump but to a strong and localised shear. Nevertheless, the general RWI criterion, ζ^{1} having a maximum (Lovelace et al. 1999), is satisfied.
Acknowledgments
The computations presented here were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujfgrenoble.fr). S.F. acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007−2013)/ERC Grant agreement no. 258729.
References
 Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502 [NASA ADS] [CrossRef] [Google Scholar]
 Bae, J., Hartmann, L., & Zhu, Z. 2015, ApJ, 805, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Flock, M., Stute, M., Kolb, S. M., & Muscianisi, G. 2012, A&A, 545, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJ, 748, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., Kritsuk, A., Norman, M. L., & Nordlund, Å. 2005, ApJ, 622, L61 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., Haugbølle, T., & Nordlund, Å. 2014, ApJ, 797, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Throop, H. B., & Bally, J. 2008, AJ, 135, 2380 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI, 411 [Google Scholar]
All Tables
All Figures
Fig. 1 Density maps as a function of time for run S71. From left to right: t = 2T_{100}, t = 2.8 T_{100}, t = 6 T_{100}. 

Open with DEXTER  
In the text 
Fig. 2 Averaged potential vorticity ⟨ ζ ⟩ (plain blue) and specific angular momentum ℒ (red dashed) of run S71. Note the sharp negative angular momentum gradient associated with a strong potential vorticity peak at R ≃ 200. 

Open with DEXTER  
In the text 
Fig. 3 Profiles of α as a function of radius, measured from R_{d}. The spiral stress increases with increasing Ṁ_{inf} and with decreasing initial angular momentum ℒ_{out}. 

Open with DEXTER  
In the text 
Fig. 4 α measured at R = R_{d}/ 10 as a function of Ṁ_{inf}. Colours represent different ℒ_{out} or H/R, diamonds correspond to symmetric inflow, and stars correspond to asymmetric accretion. The dashed line are fits given by Eq. (4). 

Open with DEXTER  
In the text 