© 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 non-ideal 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 MRI-driven turbulence, self-gravity, 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 Bondi-Hoyle 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 self-gravity to concentrate on a purely hydrodynamical model. The equations of motion read $\begin{array}{ccc}\frac{\mathit{\partial }\mathrm{\Sigma }}{\mathit{\partial t}}\mathrm{+}{\nabla }{·}\mathrm{\left(}\mathrm{\Sigma }{v}\mathrm{\right)}& \mathrm{=}& \\ \frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}{v}{·}{\nabla }{v}& \mathrm{=}& \end{array}$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 solar-mass star.

Instead of the usual minimum mass solar nebula (MMSN) density profile, we use a power-law surface density profile with an exponential outer cut-off $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}(\frac{\mathit{R}}{\mathrm{1}\hspace{0.17em}\mathrm{au}}{)}^{-1}\exp (\mathrm{-}\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{c}}})\mathit{,}& & \end{array}$(3)where we choose Σ₀ = 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 ${\mathit{v}}_{\mathrm{K}}\mathrm{\equiv }\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}\mathit{R}}$ is the Keplerian velocity. The equation of state is then simply $\mathit{P}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\Sigma }$. 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 S7-1. From left to right: t = 2T100, t = 2.8 T100, t = 6 T100.

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θ Σ_outv_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 non-axisymmetric 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 log-spaced grid cells in the radial direction. Parabolic reconstruction is used in each cell to cope with the non-uniform radial grid, and a third-order Runge-Kutta algorithm is used to evolve the system in time. We use the HLL Riemann solver¹ 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 S7-1 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₁₀₀ (or equivalently in units of 10³ 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 $\overline{\hspace{0.17em}\hspace{0.17em}\mathrm{·}\hspace{0.17em}\hspace{0.17em}}$ and a temporal average ⟨·⟩. We start the temporal average from t = 10 T₁₀₀ to allow the system to relax from the initial condition (this corresponds to 10⁴ orbits at the inner boundary).

These averaging procedures allow us to define fluctuating quantities $\mathit{\delta }{v}\mathrm{=}{v}\mathrm{-}\overline{\mathrm{\Sigma }{v}}\mathit{/}\overline{\mathrm{\Sigma }}$. The first diagnostic is the Shakura & Sunyaev (1973)α parameter $\mathit{\alpha }\mathrm{\equiv }\overline{\mathrm{\Sigma }\mathit{\delta }{\mathit{v}}_{\mathit{R}}\hspace{0.17em}\mathit{\delta }{\mathit{v}}_{\mathit{\theta }}}\mathit{/}\overline{\mathit{P}}$. 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 well-defined 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⁴ orbits at the inner boundary, which corresponds to 30 T₁₀₀.

Our runs are labelled “XY-Z” 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 S7-1 is presented in Fig. 1. The falling material initially forms an accretion shock that propagates inward. At t ≃ 2 T₁₀₀, the shock stalls at R = R_d. The stationary shock then develops an instability, which appears as a short-wavelength non-axisymmetric ondulation of the shock front (t = 2.8 T₁₀₀). This instability then quickly saturates and generates strong spiral waves that propagate inward on the long term (t> 3.5T₁₀₀).

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 $\mathrm{⟨}\mathit{\zeta }\mathrm{⟩}\mathrm{\equiv }\mathrm{⟨}\mathrm{(}\mathrm{2}\mathrm{\Omega }\mathrm{+}\mathit{R}{\mathrm{\Omega }}^{\mathrm{\prime }}\mathrm{)}\mathit{/}\overline{\mathrm{\Sigma }}\mathrm{⟩}$ 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 Kelvin-Helmholtz instability (KHI, also known as Rossby wave instability (RWI) in the astrophysical context²). 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 S7-1. Note the sharp negative angular momentum gradient associated with a strong potential vorticity peak at R ≃ 200.

Fig. 3 Profiles of α as a function of radius, measured from Rd. The spiral stress increases with increasing Ṁinf and with decreasing initial angular momentum ℒout.

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 S6-1 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 = Rd/ 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).

From these results, we can estimate that in the symmetric inflow case $\begin{array}{ccc}\mathit{\alpha }\mathrm{(}{\mathit{R}}_{\mathrm{d}}\mathit{/}\mathrm{10}\mathrm{)}\mathrm{\simeq }{\mathit{\alpha }}_{\mathrm{0}}\hspace{0.17em}(\frac{\mathit{Ṁ}\inf }{{\mathrm{10}}^{-7}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{yr}}^{-1}}{)}^{\mathit{\gamma }}\mathrm{·}& & \end{array}$(4)For inflows with significant rotation (ℒ_out = 1) we have α₀ = 1.3 × 10^-4 and γ = 0.4, whereas for inflows without initial rotation we obtain α₀ = 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 spiral-driven 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 non-axisymmetric 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 spiral-driven 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.

Acknowledgments

The computations presented here were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.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

All Tables

All Figures

	Fig. 1 Density maps as a function of time for run S7-1. From left to right: t = 2T100, t = 2.8 T100, t = 6 T100.
In the text

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

	Fig. 3 Profiles of α as a function of radius, measured from Rd. The spiral stress increases with increasing Ṁinf and with decreasing initial angular momentum ℒout.
In the text

	Fig. 4 α measured at R = Rd/ 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).
In the text