© ESO, 2016

1. Introduction

To drive mass flow in accretion discs, an anomalous source of angular momentum is required (Frank et al. 2002). A strong candidate is the magneto-rotational instability (MRI), which gives rise to turbulent magnetohydrodynamical (MHD) flows that create an outward angular momentum transport disc (Balbus & Hawley 1998). Driven by magnetic fields, the MRI requires a sufficient level of ionization to sustain a turbulent state within the disc. However, protoplanetary discs only have a very low temperature regime and insufficient thermal ionization. Even considering external sources of ionization, there appears to be a region of insufficient ionization level such that the MRI cannot operate, as shown by resistive MHD simulations, including radiative transport (Flaig et al. 2012). Hence, there may exist a dead zone somewhere between 2−20au (Armitage 2011), where the MRI can only produce very weak turbulence. Recent simulations that also included, in addition to Ohmic resistivity, ambipolar diffusion have even shown no signs of turbulence at all in this region (Gressel et al. 2015). The Hall effect creates strong winds on the surface of the disc and may even reintroduce angular momentum transport in the dead zone, but this depends on the sign of the magnetic field (Bai 2014, 2015). Thus another origin of instability inside the dead zones is warranted to drive accretion in protoplanetary discs.

As an alternative to the MRI, different examples of purely hydrodynamic instabilities in discs have been suggested, such as gravitational instability (Lin & Pringle 1987), convective instability (Ruden et al. 1988), or baroclinic instability (Klahr & Bodenheimer 2003), but they do not operate under general conditions. One possibility, that has recently attracted more attention is the vertical shear instability (VSI) suggested for accretion discs by Urpin (2003). The mechanism was first examined in relation to differential rotating stars (Goldreich & Schubert 1967; Fricke 1968), and it is also known as Goldreich-Schubert-Fricke instability. In the context of discs, first simulations have been carried out by Arlt & Urpin (2004). While there is a much stronger radial gradient in the angular velocity, Ω, to feed instabilities, most instabilities cannot overcome the stabilising effect of rotation. In the context of the VSI, it is the vertical shear in Ω, created by a radial temperature gradient, that allows the disc to become unstable. The numerical work of Nelson et al. (2013) shows that a small turbulent α-value in the range of a few 10^-4 was possible for isothermal discs. For non-isothermal discs, Nelson et al. (2013) point out that radiative cooling (diffusion) and viscosity reduce the instability, and they developed a theoretical model describing the initial vertical elongated modes that destabilise the disc. In simulations that included full radiative transport, Stoll & Kley (2014) showed that, for situations typical in protoplanetary discs, a sustained VSI was possible providing an α ~ 10^-4. They found the development of a global wave pattern within the disc whose wavelength was determined partly by viscous effects. Later, Barker & Latter (2015) analysed the VSI through linear analyses of locally isothermal discs and support the modal behaviour seen in the non-linear simulations by Nelson et al. (2013) and Stoll & Kley (2014). They also stress the importance of viscosity to set the smallest length scale.

Recently, Lin & Youdin (2015) shed light on the cooling requirements of the VSI that arise because the VSI has to compete with the stabilising vertical buoyancy. Their theoretical models predict activity in regions with long cooling time only for very large wavenumbers. Thus the VSI is limited by the cooling time on large scales and by viscosity on small scales. With this in mind they predict VSI activity for typical disc models only between 5 and 100 au.

We test this idea by expanding our previous work (Stoll & Kley 2014), where we used a self-consistent radiation transport module and a vertically irradiated disc in the simulations. Here, we treat the irradiation in a more realistic way as originating from the central star, and show that, even under this condition, the VSI can be sustained.

The turbulence in protoplanetary discs is also critical for the initial dust evolution that leads eventually to planet formation. Numerical simulations performed by Johansen & Klahr (2005) and Fromang & Nelson (2005) showed that particles can be caught in the local pressure maxima generated by the MRI turbulence. This clustering of dust particles can then trigger a streaming instability (Youdin & Goodman 2005) that will lead to further clustering and subsequently to the formation of km-planetesimals. In the context of the VSI the large-scale velocity patterns of the corrugation mode promises interesting behaviour for embedded dust grains and larger particles. To investigate the impact of the VSI modes on the dust particles we add particles into our disc model and follow their dynamical evolution.

The paper is organized as follows. In Sect. 2 we present our numerical and physical setup. We present a detailed analysis of an isothermal disc model in Sect. 3, and we discuss in Sect. 4 the results for the particle evolution is this model. In Sect. 5 we describe the results of a viscous model. The simulations with radiation transport and stellar irradiation are presented in Sect. 6 and in Sect. 7 we conclude.

2. The model setup

We use the same equations and physical disc setup, as described in detail in our first paper (Stoll & Kley 2014) and give here only a very brief outline. In summary, for the integration of the hydrodynamical equations we use PLUTO, a publicly available code, based on a Godunov scheme for viscous hydrodynamical flow (Mignone et al. 2007), extended by flux-limited diffusion module for radiation transport and a ray-tracing method for stellar irradiation (Kolb et al. 2013). Having used a two-dimensional, axisymmetric disc setup in our previous work, we now extend the computational domain in the azimuthal direction to three dimensions (3D) and add particles to the flow. For this purpose we added a particle solver based on the method by Bai & Stone (2010), that we describe in the next section.

2.1. Particle Solver

We use Lagrangian particles with drag and gravitation: (1)where f is an acceleration due to an external force, here the gravitation of the star. v_p and u are the particle and gas velocity and t_s is the stopping time.

We treat all particles as if they were in the Epstein regime, where the mean free path of the gas molecules is typically greater than the particle cross section (Epstein 1924). The stopping time is then (2)where r_p is particle radius, ρ_p the particle bulk density, ρ_g the gas density, and c_s is the sound speed. In addition we use τ_s = t_sΩ_K for the dimensionless stopping time.

To solve the equation of motion of the particles (1) we follow Bai & Stone (2010). If the stopping time is longer than the time step, Δt, of the simulation, we solve the semi-implicit equations

where n denotes the timestep level, and x^′ an intermediate position of the particle. The particle acceleration, a, is a function of the particle velocities, v_p, see Eq. (1).

For stopping times smaller than the time step we solve the following implicit equation, where the velocity update does not depend on the old velocity.

This enables the drag force to damp the particle velocity without unphysical oscillations. The test simulations to verify the correct implementation of the particle solver are described in the appendix in Sect. A.

2.2. Physical setup

In order to study the importance of radiative effects we decided to perform first a sequence of isothermal simulations. We describe briefly the setup of our fiducial disc model, which consists of a 3D isothermal model that simulates one eighth (covering 45° in azimuth) of the disc at a resolution of 1024 × 256 × 64 (in r,θ,φ) without viscosity. The full radiative model is described in Sect. 6.

We use the same general disc setup as in Stoll & Kley (2014), where we started with a disc in force equilibrium, (3)with a temperature that is constant on cylinders (4)Here, denotes the isothermal sound speed, and H = Rh = c_s/ Ω_K is the local pressure scale height of the accretion disc. We use (R,Z,φ) for cylindrical coordinates and (r,θ,φ) for the spherical coordinates that we use in our simulations. Typical values for the density and temperature exponents are p = −1.5 and q = −1.

The computational domain of the fiducial model is limited to r = 2−10au in the radial and ± 5 scale heights in the meridional direction. In contrast to Stoll & Kley (2014), where we only simulated 2D, axisymmetric discs, we use here in the azimuthal direction one eighth of the full disc (φ_max = π/ 4) to capture the complete 3D physics of the turbulent disc.

3. The turbulent isothermal disc

Before embedding the particles we first describe the turbulent properties of the VSI for several disc models. We start with the fiducial isothermal and inviscid model for one eighth of the disc and compare this below to models with viscosity and a larger azimuthal domain. In the subsequent section we use the fiducial model for the simulations with particles. Finally, in the last section we also show the results of a radiative, irradiated disc model.

As shown in previous numerical simulations of the VSI the turbulent state is characterized by vertically elongated flow structures as shown in Fig. 1 for the meridional velocity, for more details see Nelson et al. (2013) and Stoll & Kley (2014). To characterize the turbulence of the disc we measure the Reynolds stress that we define here as (5)where is the deviation of the angular velocity from the (time averaged) mean azimuthal velocity. We denote averages taken over certain variables by ⟨⟩ with the appropriate indices. In this definition, T_r,φ is a 2D array in r,θ. The time average, , needed in Eq. (5) is in general not known in advance but one can rewrite the Reynolds stress as (6)where the right hand side can in principle be calculated “on the fly” during the simulation. For convenience, we chose to store φ-averaged 2D data sets at regular time intervals, and then average these over time to obtain T_r,φ(r,θ).

From these we calculate the dimensionless α-parameter as a function of radius (7)where P = ⟨p⟩_φ,t is the azimuthal and time averaged pressure. For the vertically dependent α(z)-parameter at a certain radius r_c, (8)we only integrate over a small radial domain around the desired radius, r_c. This averaging procedure to calculate α has to be used for general discs, for example the radiative discs below.

However, for the isothermal simulations, one can approximate the time averaged in the calculation of the Reynolds stress, by the analytically calculated solution for the equilibrium angular velocity (Nelson et al. 2013), that can be obtained from the initial equilibrium disc setup (9)where is the Keplerian angular velocity. In Fig. 2 (upper panel) we can see that this is indeed a valid approximation. Aside from reducing noise in α, since it no longer depends on the time averaged velocity, this has the further advantage of allowing us to directly calculate the α-parameter at each timestep and we thus only need to store 1D arrays, which we then can later average over arbitrary timespans.

In Fig. 3 we compare the radial α-parameter for different isothermal simulations. The averaging was done from 1000 yr to 1800 yr with 800 snapshots taken at regular intervals. Shown are cases with different resolutions, where the labels indicate the radial number of grid cells, and also viscous models. The dimensionless viscosity is given in units of (u_K,1au·1au), where u_K,1au is the Kepler velocity at 1 au. We can see that in contrast to the 2D simulations in Stoll & Kley (2014), where the α-parameter differed by more than 50% when we doubled the resolution, these 3D simulations show no clear dependence on resolution, with all curves being within 10% of each other at 6au. We ran a further, shorter simulation without viscosity and double resolution of 2048 × 512 × 128 and could not see a difference in wavelength or Reynolds stress in the early equilibrium phase from 500 to 800 yr. Additionally, we see only a weak dependence on viscosity with noticeable differences beginning with the strongest kinematic viscosity, which starts to suppress the VSI in the inner region. A further increase in viscosity would suppress the instability completely, as shown already in Nelson et al. (2013).

Since we are interested in the particle drift, we also show in Fig. 4 the radial velocity in the disc at 5au as a function of height. We averaged from 1000 yr to 1800 yr, which is roughly the quasi-stationary phase. The radial velocity is inwards in the midplane and outwards in the corona. There are only minor differences in height were the direction of the flow changes. Interestingly, this profile is opposite to that of a laminar viscous flow, where it is outwards in the midplane, and inward near the disc surfaces (Urpin 1984; Kley & Lin 1992). Our findings are in agreement with results of isothermal MHD simulations of global turbulent accretion discs without a net magnetic vertical flux that also show gas inflow near the disc midplane and outflow in the disc’s surface layers (Flock et al. 2011). Within the framework of viscous discs the vertical variation of α (shown in the lower panel of Fig. 2) will play a role in determining the u_r(z)-profile, see also Kley & Lin (1992) and Takeuchi & Lin (2002), who studied the u_r(z)-profile for constant α.

3.1. 3D-simulation: full disc

In this subsection we present a full disc, meaning an azimuthal domain from 0 to 2π in contrast to the 0 to π/ 4 of our fiducial model, in order to check the validity of our results obtained from calculations with the reduce domain. Since this full simulation is computationally expensive, we used a lower resolution of 512 × 128 × 512. The full simulation is compared to the one eighth simulation of the disc with a resolution of 512 × 128 × 64, which has the same azimuthal extent as the fiducial model, but not the resolution. We added a small dimensionless physical viscosity, ν = 10^-7, to be independent of the unknown numerical viscosity, allowing better comparison with other simulations.

To give an impression of possible differences in the flow structures, between the full and fiducial model, we display in Fig. 5, from top to bottom, the fluctuations of the density, and the radial and vertical velocity components in the midplane of the disc, where the top inset in each panel refers to the fiducial model with φ_max = π/ 4 and the bottom part to the full disc. The two top panels for the density and radial velocity clearly show non-axisymmetric, wave-like features in the disc. The bottom panel seems to indicate a more axisymmetric structure which is a result of the VSI eigenmode dominating the vertical motion in the disc, as seen above in Fig. 1.

To analyse the turbulent structure in more detail we calculated azimuthal power spectra of the radial and vertical kinetic energy fluctuations in the disc midplane for two different times, spatially averaged from 3 to 7au. We can see in Fig. 6 that during the initial growth phase (top panel) the instability is driven at two length scales. The smaller one, at azimuthal wave numbers m ≈ 200, is most likely due to the initial noise as given by the finite discretization which is enhanced by the growth phase of the instability. The larger one (at m ≈ 10) is on the scale of the wavelength of the strongest growing VSI mode. This feature is also visible in the azimuthal direction even though the VSI should be axisymmetric. We speculate that this is due to a disturbance created in the radial direction by a Kelvin-Helmholtz instability that is sheared into the azimuthal direction. We can also see that in the quasi stationary phase (bottom panel) the turbulence decays faster than Kolmogorov turbulence for wavenumbers around m ≈ 20, which is the scale at which the VSI is driven. There the energy is concentrated in the VSI modes and not in the turbulent kinetic energy and thus the model of Kolmogorov decay is not applicable (Dubrulle 1992). Only at small scales the vertical kinetic energy fluctuations decay with a Kolmogorov spectrum. This can be seen in both simulations. From this we infer that the local properties of the turbulence generated by the VSI are very similar for the restricted and full azimuthal domain and is already fully captured by the smaller disc.

To analyse the locality of the VSI further, we also compare the two point correlation functions for the two disc models. This is defined as (10)where f is the quantity to be correlated, which has a zero mean, ⟨f⟩ = 0. We evaluate the correlation in the disc midplane and take the radial domain from 3 to 7au, which we treat for this calculation as periodic. The correlations are evaluated on a logarithmic grid, to better capture the properties under investigation.

In Fig. 7 we present the results of the two point correlations for the density fluctuations and the radial and meridional velocity after 1700 yr. The fluctuations are clearly non-isotropic and correlated not only on a local scale but also weakly over the whole domain. We can see again that the smaller domain is a reasonable approximation to the larger domain, even though the correlations are enhanced in the smaller domain, but for the vertical velocity we have a global correlation for both cases. This again strengthens the impression that the fluctuations in azimuthal direction are driven by a Kelvin-Helmholtz instability that feeds off the strongest VSI mode. From this we conclude that our reduced domain captures all of the important physics even though it enforces stronger correlations¹.

The models for different levels of background viscosity were shown already in Fig. 3. The full model has a slightly smaller α, but is apart from this very similar. We also compared the results from the particle motions (not shown), but also only found minor differences. Hence, we conclude that we can use the reduced model with φ_max = π/ 4 to analyse the motion of embedded particles.

4. Isothermal discs with dust

The particles are added into the fiducial model after 200 yr when the VSI is reaching the quasi-stationary phase at the inner part of the disc. After further 800 yr both VSI and particles are then in quasi-equilibrium. The particles are inserted in the midplane, randomly distributed over the radius and azimuthal angle, with the velocity of the gas at their current position. If a particle leaves the domain in the inner edge we insert it again randomly positioned over 1au at the outer edge. This ensures that a clump of particles that leaves the domain is sufficiently smoothed out when added again. We add 10 000 particles per size with 20 different sizes, beginning with 1 μm up to 3000 m. The Epstein regime is strictly valid only up to sizes of 10 m, and particles greater than that are added as a numerical experiment. Note, that for the isothermal discs the scale for the density is not fixed and thus the particle size regime may be different for different choices of the disc density. This is not the case for the disc with radiation transport as shown in the next section, where the chosen disc mass (the density) fixes the disc temperature, and hence the disc scale height, through specific values of the opacity.

We begin with the results for the radial drift and diffusion of the dust particles, then we discuss the vertical diffusion and finally the relative velocity distribution for colliding particles. We discuss our results either in terms of the physical size of the particles or the stopping time. For the models in this section the correspondence between these can be read off from Table 1.

4.1. Radial particle drift

The radial drift velocity for the dust in the midplane in the Epstein regime is given by (Nakagawa et al. 1986) (11)This is due to the drag force resulting from the difference between the Kepler velocity of the particles and the gas velocity, which is modified by the pressure support p. We use the midplane pressure for our theoretical curves, since most particles are in the midplane, at least the larger ones. Here τ_s is the dimensionless stopping time (see Eq. (2)) and u_K is the Kepler velocity.

We start by comparing the radial drift of the dust particles in our simulation with this theoretical prediction for the drift velocity in Fig. 8. The results from the simulation are extracted by fitting a linear function to the mean radial position of the particles starting with a distance to the star of 5 ± 0.5au. We fit over the span of 1000 yr beginning with 800 yr after inserting the particles or over the time the particles need to travel 0.5au, whichever is smaller. The inwardly directed drift is plotted in green (dots), and outward drift in blue (squares). Since the error in the measurement of the velocity stems from the random walk of the particles due to the turbulence, we estimate the error from the radial diffusion coefficient (see Fig. 12). The error is then given by the half-width of the distribution divided by the square-root of the number of involved particles.

We can see that the speed of large particles is similar to the predictions. A difference in speed of approximately 20% can be seen for the smaller particles in a size range between 0.1 cm to 10 cm. This deviation can be partly attributed to the spread of the particles in the vertical direction, because particles not in the midplane have a longer stopping time (lower gas density), and the prediction is calculated using the stopping time of particles in the midplane. An additional factor is that they can be caught temporarily in small scale vortices (Johansen & Klahr 2005).

The drift velocity for the smallest particles can even be positive for some time intervals. This is shown for particles at 5au in Fig. 8. At distances closer to the star they clearly drift inwards even though the gas momentum is the same (not shown). We return to the analysis for those particles later in Sect. 5.1, where the sign of the migration direction is better constraint.

An interesting behaviour of the particles is shown in the histogram in Fig. 9. For this plot we divide the distance from the star from 2 to 10au into 800 equal sized bins and count the number of particles per bin. We can clearly see that 1500 yr after we inserted the particles, the distribution clearly deviates from the initial Poisson distribution and instead they clump together. This happens only for particles with a dimensionless stopping time of the order of unity.

To further illustrate this feature, we show the number of particles per radial bin over time for different particles sizes in Fig. 10. We note that the colour scale is logarithmic. For this image and the following analysis we correct the particle density per 1au to remain constant, as it was at the beginning of the simulations, by weighting the number of particles per bin by the number of particles per 1au. Since the particles move faster in the inner region and the particles that leave the inner region are added in the outer region, we would produce overdensities otherwise. This does not change the bunching statistics, however.

After we insert the particles, the VSI is only active in the inner region, but is quickly spreading out to the whole disc, until after 1000 yr the whole domain is active. One can see that the onset of the VSI leads to bunching of the particles. Due to the bunching we can clearly see the different drift velocities of the particles. But we can also see that at certain radii the particles are sometimes caught in the pressure fluctuations for a short time before moving on, which leads to the visible lines in the image.

To make the dynamics involved clearer, we calculate another statistic property. This time we count how often a certain number of particles in a radial bin occurs. We average over 50 snapshots, each 10 yr apart, beginning 1300 yr after we inserted the particles. We then normalise by the total number of particles to find the probability for certain number of particles in a radial bin with width Δr = 0.01. This is shown in Fig. 11.

We can see that for the limit of small stopping time we still follow the initial Poisson distribution, which has its peak at 12.5 particles per bin (the average number) and decays with ∝exp(−n²). These particles are tightly coupled to the gas and can not be caught in pressure fluctuations. In contrast, particles with (dimensionless) stopping time near unity decay with ∝exp(−n). This increases the likelihood to find bins with a high number of particles which can easily lead to overdensities in dust by a factor of around 10. This is caused by short lived pressure fluctuations originating from the VSI, which briefly slow down the crossing particles. The largest particles again revert to the Poisson statistics, since they are not coupled to the gas.

4.2. Radial diffusion of particles

Next we compare the radial diffusion of the embedded particles with theoretical predictions. Since the power spectra of u_r,δu_φ,u_θ are similar, the radial diffusion coefficient for particles (Youdin & Lithwick 2007) is given by (12)For the simulation used for Fig. 12 we measure and near the midplane, indicating isotropic turbulence. With these values and Eq. (11) for the radial drift we can use the radial diffusion to measure t_eddy. Surprisingly, from this we calculate a small dimensionless, τ_eddy = t_eddyΩ_K of 0.1, compared with the large-scale oscillations of the VSI on a timescale of 5 orbits per oscillation.

In Fig. 12 we show the radial distribution of particles with different sizes from which we extracted the radial drift and diffusion during a timespan of 100 yr for different particle sizes, again for particles starting at 5 ± 0.5au. While the smallest particles (blue and green curves) follow directly the prediction, the particles with τ_s = 0.19 (red curve) seem to lag behind the theoretical curve. As described below this difference is caused by the fact that the particles are more spread out vertically. Even larger particles (yellow) again show bunching behaviour and are collected in two distinct peaks, due to being caught in different VSI waves.

From this we can also calculate the radial Schmidt number, which Youdin & Lithwick (2007) determine for homogeneous isotropic turbulence in the xy-plane with , to be (13)where for the gas diffusion coefficient, D_g,r, we use in our case the dust diffusion of the smallest particles.

For smaller particles we find good agreement with Eq. (13) but for τ_s of unity and larger, we measure Schmidt numbers smaller as predicted by Eq. (13) by a factor up to three at a stopping time of τ_s = 200. This is just noticeable in Fig. 12 for the stopping time of τ_s = 1.9. There we can see that the predicted diffusion fits a single peak well, but not the whole curve, which is 50% wider. This can be explained by particles crossing large-scale VSI modes which are not expected by the theoretical model of homogeneous isotropic turbulence used to calculate the predictions for the Schmidt number.

4.3. Vertical diffusion of particles

After the radial diffusion shows a small eddy lifetime, we are now interested in the vertical diffusion. In Fig. 13 we selected particles between 4.5 and 5.5au and calculated a histogram of the vertical position to show the vertical distribution of different sized dust particles. For better statistics we added up the 50 last snapshots, spanning from 1200 yr to 1700 yr. While the small particles with stopping time smaller than unity show the expected Gaussian distribution with a scale height equal to the gas scale height, the particles with stopping time around unity clearly deviate from this. This can be explained by the large-scale velocity pattern of the gas with active vertical shear as displayed in Fig. 1. The corrugation mode of the gas will move the particles upward away from the disc’s midplane, up to the point where the density of the gas is so small that the drag force can no longer overpower gravity. The particle then can swiftly “surf” on the updraft of the large-scale VSI-mode. Since these corrugation modes oscillate, the particles will also oscillate around the midplane.

For isotropic turbulence Dubrulle et al. (1995) and Zhu et al. (2015) calculate a dust scale height of (14)where h = H/R is the relative gas scale height. Together with the measured velocity dispersion we can use this equation to calculate again the eddy timescale by fitting a Gaussian to the data for different stopping times, thus extracting the dust scale height. We note that for a very small stopping time this does not work, since the dust scale height is equal to the gas scale height, independent of the eddy time scale. Also we superimpose two Gaussian with the same scale height for the distributions with two peaks, which then fits well. These peaks are then usually two scale heights apart for particles with dimensionless stopping time around one. Without this scheme, we would get a larger scale height and thus a longer eddy lifetime, which makes sense, since they are caused by the large-scale structures that have a long life time. We present the results in Table 1. Again we find τ_eddy ≈ 0.2, which is similar for other types of turbulence, for example MRI has τ_eddy ≈ 1 (Youdin & Lithwick 2007; Carballido et al. 2011). We note that the longer eddy lifetimes for the larger particles indicate that the eddy lifetime is longer in the midplane, since the smaller particles also see the eddy lifetime of the corona.

4.4. Collision statistics

In this section we evaluate the relative velocity distribution for colliding particles. Since we do not have enough particles to directly measure this distribution, we take each particle between 4 and 7au and check for other particles in a sphere with radius smaller than 0.05au. To make them independent of translation we remove the Kepler velocity (for independence of radius) and rotate all particles into the same r-θ plane, after we calculate the distance (for independence of the azimuthal angle, since they are all in circular orbits). We then calculate histograms of the relative velocity between two particles depending on the size of the particles. We normalise by the number of particles to get a probability.

This can be seen in Fig. 14. Lost by the normalization procedure is the fact that there are around 300 same sized particles in the sphere for small particles, while there are around 11 000 same sized particles with stopping time around unity in the sphere. Since the median relative velocity for larger particles is reduced by an order of magnitude relative to the smaller particles, the number of collisions will also be reduced. The blue dashed line represents the relative velocity between small and large particles, since those are dominated by the velocity difference that is created by the strongly coupled small particles, moving with the gas, and weakly coupled particles moving with Keplerian velocity.

Using an alternative way of estimating the bunching behaviour of particles, we also calculate the pair correlation function g(r). (15)where V is the area of integration, d_ij is the distance between particle i and j and δ(r) is one if | r | < Δr/ 2 and we use Δr = 0.01au. This function returns one for Poisson distributed particles and greater than one if there is an increased surface density in the ring around the particles at this radius and thus picks up 2D clustering instead of the 1D clustering in the earlier analysis. We calculate this property for particles between 4 and 7au projected in the r−φ plane, and show it for different particle sizes in Fig. 15. We can see that small particle positions are uncorrelated, but the particle positions with stopping time near unity display a clear correlation, as we could already infer from Fig. 11. Particles with a stopping time closer to unity have a larger correlation length. This makes the VSI a possible candidate to trigger the streaming instability.

5. 3D-simulations: viscosity

In this section we present the same simulation but using now a dimensionless kinematic viscosity coefficient of ν = 5 × 10^-7, in order to check the influence of viscosity, which in turn influences the α-parameter and the velocity dispersion. We also add the same amount of viscosity to the model with radiation transport, as shown below. This viscosity corresponds to an α-value of 9 × 10^-5 at 5 au or 4 × 10^-5 at 25 au. There it limits the smallest length scale of the VSI, which can not be resolved otherwise. We only show the results if there is a clear difference to the previous simulation. From the hydro-dynamic perspective they are very similar, but the α-parameter is smaller by a factor of 2 and the wavelength of the instability slightly smaller.

5.1. Radial drift

If we include viscosity and repeat the analysis for radial drift, here from 2700 yr to 3700 yr, due to the slower diffusion, we can see in Fig. 16 that particles larger than about 0.1 cm drift inward with approximately the theoretical speed (see Eq. (11)) with slightly higher deviations than in the inviscid case. The smallest particles now clearly drift away from the star at r = 5au. This is true for different radii of the disc.

In Fig. 17 we see that the particles are moving inwards at the midplane and outwards otherwise. The smallest particles follow the gas velocity and larger particles are moving outwards away from the midplane faster than the gas, similar to Takeuchi & Lin (2002), where they move outwards even though the gas is moving inwards, due to the gas being super-Keplerian in the disc’s corona. This leads to a massflow (shown in Fig. 18) that is inward in the midplane and outward farther away from the midplane. As seen in Figs. 17 and 18 small particles (blue and green line) show the same radial velocity profile and mass flux behaviour as the disk’s gas flow (black dashed lines). To estimate the overall mean flow of the particles it is necessary to consider the vertical dependence of the mass flow (Fig. 18) rather than the radial velocity z-profile (Fig. 17). For example, from Fig. 17 we notice that the particles with stopping time τ = 6.1 × 10^-5 (green line, 3.2 × 10^-3 cm radius) move inward slower than the particles with stopping time τ = 6.1 × 10^-4 (red line, 3.2 × 10^-2 cm radius) even though the mean radial dust velocities in Fig. 16 are smaller for the larger particles. This is due to the different vertical distribution of the particles (see Fig. 18 and Table 2). In addition, the velocity and mass flow profiles of the smallest displayed particles (blue lines in Figs. 17 and 18, 3.2 × 10^-4 cm radius) look very similar to the particles with stopping time τ = 6.1 × 10^-5 (green line), while having drift rates in the opposite direction in Fig. 16.

Thus the net particle flow is very sensitive to the gas flow profile and the ratio of particles near the midplane, which is decided by the stopping time. This also means that, independent of the mean flow, there will always be a small fraction of particles drifting away from the star, faster than one would expect from diffusion alone.

5.2. Radial diffusion

The increase in viscosity leads to an decrease of the velocity dispersion to which is smaller by a factor of two. In Fig. 19 we compare our results to theoretical predictions with τ_eddy = 0.1. The diffusion is clearly smaller than in the inviscid case, as predicted by the decreased velocity dispersion. This is in contrast to the vertical eddy lifetime were we measure an increase in eddy lifetime due to vertical diffusion.

We can also see (right wing of green and blue curves) that a small fraction of the small particles diffuse faster outwards than the rest of the particles. These are particles far away from the midplane, where the stopping time is magnitudes longer, typically longer than 10^-2, and the gas flow is in average outwards. These weaker coupled particles can quickly travel a short distance away from the star, before drifting back inwards nearer to the midplane. The difference in radial drift between theory and simulation noticed in Fig. 16 causes the offset for the results of the longer stopping time (red curve).

5.3. Vertical Diffusion

For vertical diffusion we obtain slightly longer eddy lifetimes as can be seen in Table 2 where we measured from 2700 yr to 3700 yr. In these simulations it took four times as long for the dust scale height to converge to the gas scale height for the smallest particles, even though the vertical velocity dispersion is identical to the inviscid case.

5.4. Clustering

Finally we present the particle distribution in the r−φ plane in Fig. 20. We show the particles with dimensionless stopping time of τ_s = 0.6 and τ_s = 1.9 at 5au. One can see for both displayed particles sizes that the VSI modes have produced nearly axisymmetric clusters. The bunching leading to the ring structure is strongest for the simulation with high viscosity. We also show in Fig. 21 the same effect for the full disc with small viscosity of 10^-7 and resolution of 512 × 128 × 512. For this simulation we also increased the number of particles to 500 000. We can see that the ring structure already seen in the histogram of Fig. 9 does indeed persist even in a full disc. Particles with more than a magnitude longer or smaller stopping time do not show this features.

6. 3D-simulations: radiative model

In this section we present the results for a radiative disc with radiation transport and irradiation from the central star. In contrast to our first paper (Stoll & Kley 2014), we model the stellar irradiation in a more realistic fashion as coming from the central star, similar to the treatment in Bitsch et al. (2013). This star has a temperature of 4000 K and a radius of 4 R_⊙.

6.1. Setup

For this simulation we have to take into account that the cooling time has to be sufficiently small for the VSI to be active (Nelson et al. 2013; Lin & Youdin 2015), which made changes in the domain necessary. We moved the radial extent of the disc for the inner boundary from 2 to 8au and for the outer boundary from 10 to 80au. This radial range is expected to be the active region of the VSI (Lin & Youdin 2015). We simulate again one eighth of the disc in the azimuthal direction and this domain is resolved by 1200 × 260 × 60 grid cells. We also changed the density profile exponent from p = −1.5 to a value that is more in line with the observations p = −1.8 (Williams & Cieza 2011). Initially the temperature drops with T = T₀·r^-1/r₀, thus we have ρ = 10^-9 g / cm³·r^-1.8/r₀. This translates to Σ = 1700 g / cm² at 1au, which corresponds to is the MMSN-model with a shallower decay of the density. The radiation transport then quickly leads to a new equilibrium with T = 900 K·r^-0.6/r₀ in the corona and T = 700 K·r^q/r₀ in the midplane, where the temperature gradient exponent q varies slightly around the mean of q = −0.9, from q = −1.1 in the inner region to q = −0.6 in the outer region.

During the evolution to the new equilibrium we damp the velocities in the whole disc.

We add a small viscosity of ν = 5 × 10^-7. This suppresses the VSI in the inner region, where it would otherwise be weakly active, but due to a wavelength on grid scale clearly not resolved, which in turn would lead to unphysical numerical artifacts. As shown above, we observe only a small change of the VSI activity in the active domain with viscosity enabled compared to the inviscid case. Thus we see no harm in adding it.

In our first paper on the behaviour of the VSI in radiative discs we only considered vertical irradiation onto the disc surfaces (Stoll & Kley 2014). Here, we make the simulations more realistic and irradiate the disc from a central stellar source from the origin along the radial direction, see Bitsch et al. (2013). In this procedure the inner rim of the disc in our simulation is directly exposed to the stellar irradiation. To prevent unphysical heating of the midplane at the inner boundary, we absorb the irradiation flux coming from the star in a fictitious ghost cells with a width 0.25 au  using the gas properties of the adjacent innermost active cells of the domain.

For the irradiation opacity we choose a value 10 times higher than the gas opacity, to compensate for the fact, that this radiation is emitted by a hot star and not the surrounding gas. This leads to a heated corona with a cooler midplane as can be seen in Fig. 22, instead of the cooler corona in Stoll & Kley (2014). At the boundary of the corona we can also see a change in the VSI mode. They have a larger wavelength in the hotter corona region and split where the temperature changes to a smaller wavelength in the midplane, see lower panel Fig. 24.

To make a direct comparison between the isothermal and radiative case, we ran an additional isothermal model with p = −1.8 and q = −0.9 and ν = 5 × 10^-5 and damping in the vertical and radial velocity in the region between 8 and 10 au to avoid boundary effects. In principle one could also compare the radiative case directly to the isothermal models from Sect. 4, because in isothermal simulations the unit of length is not fixed and can be scaled to a different regime. However, the gradients in density and temperature are not the same. We thus included a new isothermal model that can be directly compared to.

6.2. Hydrodynamic properties

We begin by presenting the α-parameter in Fig. 23, here calculated by time averaging the azimuthal velocity, since the equilibrium velocity cannot be computed analytically for radiative discs with a vertically varying temperature. In the inner region the VSI is suppressed by the viscosity of ν = 5 × 10^-7 on small wavelengths and by the high cooling time on large wavelengths. Compare this to the isothermal simulation, where the same viscosity is not able to suppress even at 4au (see Fig. 3). This is followed by an active region beginning at 15au, where we reach α = (1−4) × 10^-4, which is still smaller than the isothermal simulations. The drop off in the outer region may be linked to the reduced activity in this region, see also Fig. 24, but is also visible in the isothermal model.

As the VSI is critically dependent on small cooling times, we analyse the cooling times due to radiative diffusion in the irradiated disc models. The radiative diffusion coefficient is given by: (16)where λ is the flux limiter, a the radiation constant, c the speed of light and κ_R the Rosseland mean opacity. To calculate the cooling time we also need the appropriate length scale. For the optically thick region we simply take the length scale of the perturbation, which we approximate as a fourth of the scale height l_thick = H/ 4. In the optically thin region we use the optical mean free path l_thin = 1 /κ_Rρ for the length scale. This leads to a combined dimensionless cooling time of (17)In Fig. 24 we compare the cooling time, τ_cool, as calculated from our numerical irradiated disc models with the critical cooling time, τ_crit, as estimated by Lin & Youdin (2015), who compared the destabilising vertical shear rates with the stabilising vertical buoyancy frequency. They obtained (18)We see a good agreement in the inner region between the active regions as predicted by the critical cooling time and the active regions in our simulation. The inner midplane region up to 10au is completely inactive and the following region which is also predicted to be inactive is only active with a higher order mode. We note that without viscosity one expects modes with higher wavenumber in this region.

In the outer region beyond 60au the VSI is inactive despite a small enough cooling time. This may be due the dynamics of the VSI that shows larger wavelengths in the outer region, thus requiring a smaller cooling time, or to boundary effects.

One can also see that the jump in temperature and cooling time, that also defines the boundary between disc and corona, creates a boundary for the VSI, where the surface modes can attach to (Barker & Latter 2015).

6.3. Dust properties

For this simulation we add 20 000 particles per size after 1000 yr. In Fig. 25 we present a histogram of the distribution of particles with a size of 31 cm, after 13 500 yr. We see that in the outer region with the inactive VSI the particles are still Poisson distributed, but in the active region they are caught in the eddies. The particles in the outer region are only collected weakly, since the VSI is reduced, due to the long cooling time. The isothermal case we compare to has clustering throughout the whole disc and the overdensities are stronger by a factor of around two, even though the velocity dispersion is higher by a factor of 5 to 10.

These results show that, even with realistic cooling times, the VSI can create small axisymmetric regions with overdensities in the dust by a factor of three. This is the right range of metallicity and size of particles which is needed for the streaming instability to set in (Youdin & Goodman 2005). This instability can further enhance the clumping until gravity is strong enough to directly form planetesimals out of the cluster of particles.

In Fig. 26 we repeat the statistical analysis for the distribution of particles. We take into account all particles in the region from 15 to 40au, in the timespan between 11 000 and 13 500 yr over 50 snapshots. We note that the average number of particles per bin is the same as in the isothermal case in the earlier section, since we increased the size of the bins to compensate for the lower density of particles. Again we see a clear deviation from the initial Poisson distribution for the particles with stopping time around unity, even though it is weaker. Interestingly the effect is now most powerful for τ_s = 6.3 × 10^-2 (particles with 10 cm radius), where the inward drift velocity for the particles and the inward drift of the vertical motion of the VSI mode is the same, thus the particles move with the bunching gas mode instead of through the mode. Those are only bunched at around 30au, and are diffused again after they have passed this region. This resonance does not exist in the isothermal case, because the wavelength is greater than in the radiative case. The isothermal case in general behaves very similar to the isothermal case in the earlier section. Both show the strongest bunching for particles with stopping time close to unity.

The radial drift shown in Fig. 27 measured at 20 ± 2au from 13 500 yr to 18 500 yr is similar to the isothermal case with the same viscosity. While the outward migration is no longer as clear as in the isothermal case with viscosity, there is still a trend to outward migration. That the effect is weaker can be explained by the weaker effect the viscosity has at 20 au.

More important for the radial motion of a single particle is the diffusion. For the radial and vertical velocity dispersion in the region at 20 ± 5au we measure for the radiative case and and for the isothermal case we measure and . Both values lead to a dimensionless eddy time of τ_eddy = 1.0 even though the velocity dispersion differs by a factor of ten. The larger difference between prediction and simulation in Fig. 28 results from the error in the measurement of q, the exponent in the radial temperature distribution. Here, T(r) is determined through the radiation transport and q varies now with radius. For the plot we use an average value of q = −1.1.

In Table 3 we can see that in this simulation the dust scale height is smaller than the gas scale height even for the smallest particles. We averaged from 13 500 yr to 18 500 yr. In this simulation the radial and vertical calculated eddy lifetimes are again very similar, despite the turbulence not being isotropic.

For the collision statistics we increased the cutoff distance within which we compare particle velocities to 0.2au to compensate the decreased density of particles. As the distribution of particles already indicated the clustering is indeed weaker. This is reflected additionally in Fig. 30, where we can see that the correlation is slightly weaker. For this radiative case the effect appears to be strongest for a dimensionless stopping time τ_s ≈ 0.1 instead of 1 for the isothermal case. The correlation length is larger for the particles with stopping time around one, reflecting the larger wavelength of the VSI in the isothermal case.

The histogram of the relative velocities between particles as displayed in Fig. 29 illustrates this situation. For τ_s ≈ 0.1 the particles have about an order of magnitude smaller relative velocities than for large and small values. We can also see that the wider velocity dispersion in the isothermal model leads to greater relative velocities.

7. Summary and conclusions

In the paper we analysed the dynamics of particles embedded in hydrodynamic discs that show fully developed turbulence as induced by the VSI.

In a first step we calculated isothermal disc models in full three dimensions and analysed the properties of the turbulence generated by the VSI. Our standard model consisted of an eighth of a full circle (φ_max = π/ 4) and showed in the fully developed turbulent state α-values around 6 × 10^-4, which is of the same order of magnitude or even slightly larger than the corresponding 2D models (Stoll & Kley 2014). The 3D models shows variations in the azimuthal direction and these fluctuations follow a Kolmogorov-type spectrum. The mean radial velocity of the gas in a VSI turbulent disc turned out to be directed inward in the disc midplane and outward in the upper layers, in agreement with global MHD simulations using zero net vertical magnetic flux (Flock et al. 2011). This flow is opposite to viscous laminar discs (Urpin 1984; Kley & Lin 1992) or MHD discs with non-zero vertical magnetic field (Suzuki & Inutsuka 2014). For 3D discs covering the full circle (φ_max = 2π) we found very similar results, which allowed us to treat particle evolution in the reduced domain.

In addition to the isothermal case we studied fully radiative models including heating from the central star. To allow for regimes where the VSI instability can operate we extended to radial domain from 8−80 au. The temperature structure in the disc displayed a central disc region with a nearly constant temperature in the vertical direction and hotter surface layers produced by the stellar irradiation. The vertically varying opacity in the disc resulted in different cooling times and the turbulence turned out to be slightly weaker in comparison to the purely isothermal situation. For the effective α-parameter values of around 10^-4 were reached in the active state that extended from about 10 to 60 au.

After having reached the equilibrium state we inserted particles of different sizes to study their motion in the disc, where the drag force between gaseous disc and particles was treated in the Epstein regime. Overall we found for both, isothermal and radiative discs, comparable results. On average the particles drift inwards with the expected speed. For all disc models we found that the smallest particles show an outwardly directed radial drift. This comes about because the small particles are coupled more to the gas flow and are lifted upward by the vertical motions of the VSI induced large-scale flows. Since the average flow direction in the upper layers is positive small dust particles that are elevated above the disk’s midplane are dragged along and move outwards. Particles below about 1 mm in size experience this fate. This outward drift might be beneficial in transporting strongly heated solid material to larger radii as required to explain for example the presence of chondrules at larger radii in the Solar System (Bockelée-Morvan et al. 2002). The upward drift of small particles in the disc by the VSI modes will also help to explain the observed presence of a population small particles in the later stages of the disc evolution that were produced by a fragmentation process (Dullemond & Dominik 2005).

Using the information of histograms, probability functions and pair correlation functions we analysed the spatial re-distribution of particles in the disc that were initially homogeneously distributed. We found that the particles are strongly “bunched” together by the large-scale motions of the VSI turbulence. The bunching effect is strongest for particles with a stopping time of the order unity and the maximum overdensities reached were about 5 times the average initial density of the particles. The relative velocity between particles of the same size is smallest (about a few m/s) for those particles that show the strongest bunching. This combination of high density and low relative speed is highly beneficial for the early formation process of planetary precursors. First, at these relative speeds collisions between two particles can lead to sticking collisions (Blum & Wurm 2008; Meru et al. 2013). The higher relative velocities between particles of different sizes does not necessarily lead to fragmentation. The experiments of Teiser & Wurm (2009) have shown that particles with different size can stick to each other even for collisions up to 50 m/s and possibly more. Secondly, through the concentration of particles it is possible to trigger streaming instabilities in the disc which can further increase the particle concentration and growth (Youdin & Goodman 2005).

The two dimensional distribution of particles in the disc shows axisymmetric ring-like concentration zones of the particles resembling very roughly the features observed recently in the disc around HL Tau (Brogan et al. 2015). Even though the strongest effect is seen here in our simulations for particles about one meter in size, it is possible that through collisions of nearly equal sized bodies much smaller particles that could generate the observed emission can be produced and which follow a similar spatial distribution. Obviously the observed spacing of the “bright” rings in our simulations is smaller than those observed in HL Tau but the inclusion of variations in opacity or chemical abundances may create larger coherent structures.

---

Acknowledgments

Moritz Stoll received financial support from the Landesgraduiertenförderung of the state of Baden-Württemberg and through the German Research Foundation (DFG) grant KL 650/16. Simulations were performed on the bwGRiD cluster in Tübingen or the ForHLR cluster in Karlsruhe (both funded by the DFG and the Ministry for Science, Research and Arts of the state Baden-Württemberg), and the cluster of the Forschergruppe FOR 759 “The Formation of Planets: The Critical First Growth Phase” funded by the DFG.

References

Appendix A: Particle solver

To verify the correct implementation of our particle solver we repeat some of the tests of Zhu et al. (2014).

Appendix A.1: Orbit test

We release the particle at r = 1, at the midplane with a velocity of u_φ = 0.7 and integrate for 20 orbits. The presented timesteps are Δt = 0.1 and Δt = 0.01 for an orbital time of 2π. Even though the orbit precesses for the longer timestep, the geometric property is conserved. There is no visible precession in Fig. A.1 for the timestep of Δt = 0.01, which we use in our simulations.

Appendix A.2: Settling test

We release particles with different stopping times at one scale height from the midplane. For particles with τ_s< 1 we can see in Fig. A.2 the exponential decay of the vertical position. Particles with τ_s> 1 oscillate around the midplane and instead the amplitude decays exponentially.

Appendix A.3: Drift test

For the drift test we use the disc in hydrostatic equilibrium and release particles with different stopping times at r = 5 au on Keplerian orbits in Fig. A.3. We compare to the theoretical expected drift velocity of Eq. (11) (black line).

All Tables

All Figures