Planet-disc interaction in laminar and turbulent discs

TL;DR

This study investigates how embedded planets interact with turbulence caused by the vertical shear instability in protoplanetary discs, revealing that turbulence has limited impact on gap formation but can influence planetary migration rates.

Contribution

It provides the first detailed comparison of planet-disc interactions in VSI-turbulent and viscous discs, highlighting the effects on gap structure and migration.

Findings

VSI turbulence has limited impact on gap formation for massive planets.

Low-mass planets experience enhanced inward migration due to over-density behind them.

Vortices at gap edges are stronger and longer-lived in VSI turbulent discs.

Abstract

In weakly ionized discs turbulence can be generated through the vertical shear instability (VSI). Embedded planets feel a stochastic component in the torques acting on them which can impact their migration. In this work we study the interplay between a growing planet embedded in a protoplanetary disc and the VSI-turbulence. We performed a series of three-dimensional hydrodynamical simulations for locally isothermal discs with embedded planets in the mass range from 5 to 100 Earth masses. We study planets embedded in an inviscid disc that is VSI unstable, becomes turbulent and generates angular momentum transport with an effective $\alpha = 5 \cdot 10^{-4}$. This is compared to the corresponding viscous disc using exactly this $\alpha$-value. In general we find that the planets have only a weak impact on the disc turbulence. Only for the largest planet ($100 M_\oplus$) the turbulent activity becomes enhanced inside of the planet. The depth and width of a gap created by the more massive planets ($30, 100 M_\oplus$) in the turbulent disc equal exactly that of the corresponding viscous case, leading to very similar torque strengths acting on the planet, with small stochastic fluctuations for the VSI disc. At the gap edges vortices are generated that are stronger and longer lived in the VSI disc. Low mass planets (with $M_p \leq 10 M_\oplus$) do not open gaps in the disc in both cases but generate for the turbulent disc an over-density behind the planet that exerts a significant negative torque. This can boost the inward migration in VSI turbulent discs well above the Type I rate. Due to the finite turbulence level in realistic three-dimensional discs the gap depth will always be limited and migration will not stall in inviscid discs.

Tables

Table 1. Table 1: Model parameter

Parameter	model
Radial range [$5.2\text{au}$]	$0.4$ - $2.5$
Vertical range [H]	$\pm 5$
Phi range [rad]	$0$ - $2\pi$
Radial grid size	$600$
Theta grid size	$128$
Phi grid size	$1024$
Planet mass [$M_{\oplus }$]	$5$, $10$, $30$, $100$

Equations

$$\rho (R,Z)={\rho }_0{\left(\frac{R}{R_{\mathrm{p}}}\right)}^p\exp \left[\frac{G{M}_{\mathrm{s}}}{c_{\mathrm{s}}^2}\left(\frac{1}{r}-\frac{1}{R}\right)\right],$$

$$T(R)=T_0{\left(\frac{R}{R_{\mathrm{p}}}\right)}^q,$$

$$\Omega (R,Z)={\Omega }_{\mathrm{K}}{\left[(p+q){\left(\frac{H}{R}\right)}^2+(1+q)-\frac{qR}{\sqrt{R^2+Z^2}}\right]}^{\frac{1}{2}}.$$

$${\Phi }_{\mathrm{p}}=\{\begin{array}{ll}-\frac{G{M}_{\mathrm{p}}}{d}\left[{\left(\frac{d}{d_{\mathrm{rsm}}}\right)}^4-2{\left(\frac{d}{d_{\mathrm{rsm}}}\right)}^3+2\left(\frac{d}{d_{\mathrm{rsm}}}\right)\right],& \text{for }d<d_{\mathrm{rsm}}\\ -\frac{G{M}_{\mathrm{p}}}{d},& \text{for }d\ge d_{\mathrm{rsm}},\end{array}$$

$$R_{R,\varphi }=\rho u_R\delta u_{\varphi },$$

$$\alpha (r)=\frac{{⟨R_{R,\varphi }⟩}_{t,\theta ,\varphi }}{{⟨P⟩}_{t,\theta ,\varphi }},$$

$${\overline{\text{@vec}\omega }}_{\theta }={\left(\frac{r}{r_{\mathrm{p}}}\right)}^{3/2}{\left(\nabla \times \text{@vec}u\right)}_{\theta }.$$

$$\text{@vec}{F}_{\mathrm{p}}=\frac{G{M}_{\mathrm{p}}\Delta m}{d^3+(\epsilon R_{\mathrm{Hill}}{)}^3}\text{@vec}d,$$

$${\Gamma }_0=\Sigma (r_{\mathrm{p}}){\Omega }^2(r_{\mathrm{p}})r_{\mathrm{p}}^4{\left(\frac{M_{\mathrm{p}}}{M_s}\right)}^2{\left(\frac{r_{\mathrm{p}}}{H}\right)}^2,$$

$${\left(\frac{d\Gamma }{dm}\right)}_0={\Omega }^2(r_{\mathrm{p}})r_{\mathrm{p}}^2{\left(\frac{M_{\mathrm{p}}}{M_s}\right)}^2{\left(\frac{r_{\mathrm{p}}}{H}\right)}^4.$$

$$\Gamma =-(1.36+0.62\beta +0.43\zeta ){\Gamma }_0,$$

Full text

1 Introduction

The longstanding problem of identifying the mechanism responsible for transporting angular momentum in accretion discs and explain their observed lifetimes has been tackled since the early works of Shakura & Sunyaev (1973) and Lynden-Bell & Pringle (1974). The physical process that has been identified having an important contribution is the magnetorotational instability (MRI) (Balbus & Hawley 1991; Hawley & Balbus 1991). However, within the context of planet formation, the basic requirement of well ionised media to trigger the MRI is generally not fulfilled in protoplanetary discs (Gammie 1996).

Several purely hydrodynamical instabilities have been proposed to drive the angular momentum, even in the non-ionised regions of the disc (also referred to as ”dead zones”), that can develop under special conditions (see e.g. Nelson et al. 2013, and references therein). One of the mechanisms that are more generally applicable is the vertical shear instability (VSI) which was first discovered in the context of differentially rotating stars (Goldreich & Schubert 1967; Fricke 1968) and then later applied also to accretion discs (Urpin & Brandenburg 1998). More recently, the VSI has been modelled by direct numerical simulations and proven effective in locally isothermal discs (Nelson et al. 2013) as well as in fully radiative discs (Stoll & Kley 2014). Within the context of the planet formation process, these instabilities can strongly impact the dust dynamics, by increasing for example the dust concentration inside vortices (Barge & Sommeria 1995; Klahr & Bodenheimer 2006; Baruteau & Zhu 2016) and regulating its settling and growth (Dullemond & Dominik 2005). In particular, Stoll & Kley (2016) studied the effect of the VSI on the dust dynamics and find a strong clustering of dust particles by the large-scale vertical motion induced by the instability. At the same time, the relative speed between individual particles is significantly reduced which can be highly beneficial for the planet formation process. Furthermore, a strong outward migration of sub-mm size dust is observed in the upper layers of the disc that can replenish the outer disc regions of solid material.

An embedded protoplanet can strongly affect the disc structure in its vicinity by creating density waves departing from its location and by depleting the co-orbital region. These asymmetries in the disc’s density generate a gravitational torque acting onto the planet and force it to migrate through the disc (see e.g. Kley & Nelson 2012, and references therein). The direction and magnitude of this migration depends strongly on the local properties (viscosity and thermodynamics) of the disc. The VSI generates strong vertical motions that can perturb the disc equilibrium, potentially affecting the resulting torque acting onto the planet. Previously, direct simulations of planets embedded in turbulent discs have been performed only for the magnetic case of MRI unstable discs (Nelson & Papaloizou 2004; Nelson 2005; Uribe et al. 2011).

In this work we focus on the signatures that are induced by a planet in a protoplanetary disc in which the VSI is operating, and calculate the torque acting on the planet. We compare our results directly to the classical viscous $\alpha$-disc models in order to study the main differences that this simpler prescription can or cannot reproduce. Additionally, we compare the strength of the torques acting on the planet in hydrodynamic VSI-turbulence to that generated in MRI-turbulence. The paper is organised as follows. In Sec. 2 we describe the numerical method used to model the disc, the planet and the parameter space studied. In Sec. 3 we outline the results of our analysis, and in Sec. 4 we discuss our major findings and draw the conclusions.

2 Setup

We used the pluto code (Mignone et al. 2007) to model a three-dimensional section of a locally isothermal accretion disc in spherical coordinates (r, $\theta$, $\phi$) with embedded planets of different masses. To study the planet-disc interaction we compared two sets of simulations. In the first series we ran inviscid models that generated direct turbulence through the VSI and in the second series we ran viscous models using an effective $\alpha$-value adapted to that of the VSI turbulent disc. Radially, the disc section extends from $0.4\,\mathrm{r_{\mathrm{p}}}$ to $2.5\,\mathrm{r_{\mathrm{p}}}$, where $r_{\mathrm{p}}=5.2\,\mathrm{au}$ is the planet position. Lin & Youdin (2015) estimated that the VSI for a minimum mass solar nebula would operate in the range between $5-50$ au, so we are just at the inner boundary of their suggested range. In our situation planets have already formed and considering a later phase of the planet evolution process the surface density in solids is reduced which would make the VSI more efficient due to the enhanced cooling. At the same time, the configuration we study can be considered an exemplary case, centered at 5.2 au, that could be scaled easily to other regimes due to the locally isothermal assumption. In the vertical direction the models cover 5 scale heights above and below the equator, and a full annulus in the azimuthal direction. This domain is covered for the standard model by $N_{r}=600,N_{\theta}=128$, and $N_{\phi}=1024$ gridcells that are spaced logarithmically in the radial and equidistant in the other directions. At the location of the planet ($r=1$ in our units) the radial size of a gridcell is $\Delta r=0.003$, and the ratio of gridsizes near the location of the planet is given by $\Delta r:r\Delta\theta:r\Delta\phi=1.0:1.28:2.012$. A comparison model with double spatial resolution is discussed in section 3.5. The main parameters of the simulations are summarised in Tab. 1.

2.1 Gas component

The initial disc is axisymmetric and, even though we use spherical coordinates in the simulations, we state the initial conditions in cylindrical coordinates $(R,Z,\phi$). The initial density profile created by force equilibrium is given by

[TABLE]

where $\rho_{0}$ is the gas mid-plane density at $R=R_{\mathrm{p}}=r_{\mathrm{p}}$, $p=-1.5$ the density exponent, and $c_{\mathrm{s}}$ the isothermal sound speed. In our simulations we use $\rho_{0}=2.07\cdot 10^{-11}$ g/cm3 which translates into $\Sigma(R_{0})=200$ g/cm2 for our chosen disc thickness. The disc is modelled with a locally isothermal equation of state, which keeps the initial temperature stratification fixed throughout the whole simulation. We assume a constant aspect ratio $H/R=0.05$, which corresponds to a temperature profile

[TABLE]

with $q=-1$ and $T_{0}=121\,\mathrm{K}$.

The gas moves with an azimuthal velocity given by the Keplerian speed around a $1\,\mathrm{M_{\sun}}$ star, corrected by the pressure support (Nelson et al. 2013)

[TABLE]

At the inner and outer boundary we damp the density, radial and vertical velocity to the initial values with a timescale of half a local orbit (see e.g. de Val-Borro et al. 2006), in order to reduce numerical issues due to the interaction of the VSI with the boundary and to prevent reflection of the spiral wave at the boundaries. The damping is applied in the intervals [0.4,0.5 ]$r_{\mathrm{p}}$ and [2.3,2.5 ]$r_{\mathrm{p}}$. For the kinematic viscosity in the $\alpha$-disc model we use $\nu=2/3\alpha c_{\mathrm{s}}H$, and adopt a constant $\alpha$ parameter of $\alpha=5\cdot 10^{-4}$, which is close to the value obtained from the VSI model (see Fig. 1). To bring the VSI model into equilibrium, we evolve it for $200\,\mathrm{orbits}$ before embedding the planet.

2.2 The planet

We embed a planet that orbits its parent star on a circular orbit with a radius $r_{\mathrm{p}}=1$ and a mass in the range $[5,10,30,100\,\mathrm{]}{M_{\oplus}}$. These 4 planets have the following Hill radii, $[0.017,0.022,0.031,0.045\,\mathrm{]}{R_{\mathrm{p}}}$. Using the grid resolution for the standard model in Table 1 this implies that there are about 100 gridcells within the Hill radius for the smallest planet mass and 740 cells for the largest planet mass. Gas accretion is not allowed, and the planet does not feel the disc, so its orbital parameters remain fixed during the whole simulation. The gravitational potential of planet is implemented with a cubic expansion in the vicinity of the planet location (Klahr & Kley 2006)

[TABLE]

where $d$ is the distance of a gas element from the planet location, and the smoothing length $d_{\mathrm{rsm}}=0.5\,\mathrm{R_{\mathrm{Hill}}}$ is adopted in order to avoid singularities. To ensure a smooth start during the initial phase of the simulations, the planetary mass is increased slowly over the first $20\,\mathrm{orbits}$ of the simulation. We run the simulations for $200$ orbital periods of the planet.

3 Results

3.1 The effective viscosity of the disc with embedded planets

For the VSI models the $\alpha$-parameter is not set by the viscosity prescription. Instead, the VSI generates eddies that transport angular momentum self-consistently. We can then calculate the efficiency of this process by evaluating an effective $\alpha$-parameter through the Reynolds stress resulting from the turbulence. We calculate the Reynolds stress in cylindrical coordinates $(R,Z,\phi$)

[TABLE]

where $\delta u_{\phi}$ is the local difference from the equilibrium azimuthal velocity, which we calculate by time averaging $u_{\phi}$ for each grid cell from orbit $60$ to orbit $200$ over $70$ snap shots.

We then calculate the dimensionless $\alpha$-parameter as a function of the radius,

[TABLE]

where $P$ is the pressure and $\left<\right>_{t,\theta,\phi}$ denotes the average over time ($140\,\mathrm{orbits}$) and the whole vertical and azimuthal domain.

In Fig. 1 we compare the $\alpha$-parameter for our simulations of VSI-active discs with embedded planets of different mass. The inner damping region from 0.4 $r_{\mathrm{p}}$ to 0.5 $r_{\mathrm{p}}$ is very effective in suppressing the VSI, which is intended to reduce the interaction with the boundary. The wave killing zone is slightly less effective in suppressing the VSI for an increased grid resolution close to the inner boundary (model ’Res2’ for the 10 $M_{\oplus}$ model), while there is no difference in the outer disc. Interestingly, the VSI is more resilient to the wave killing zone for the 100 $M_{\oplus}$ case. There the VSI is stronger and can access a region much closer to the inner boundary. This behaviour is also highlighted in Fig. 2, where the vertical velocity in the midplane of the disc (induced by the VSI) reaches $0.5c_{\mathrm{s0}}$ for the more massive case. One possibility for this extended turbulent activity might by the action of the spiral wave instability (SWI) where the presence of a massive planet can excite unstable inertial-gravity waves (Bae et al. 2016b, a) which could enhance the turbulent level above that of the VSI alone.

In the top panels of Fig. 2 we see that small planetary masses are not able to influence the strength of the VSI, even close to their location. The strong deviation in the $\alpha$-profile seen in Fig. 1 for the large mass planets can be directly linked to the bottom panels of Fig. 2 with the formation of strong vortices at the gap edges carved by the planets. The vortex at the outer edge is strong enough to suppress the VSI in the region [1.2, 2.0 ]$r_{\mathrm{p}}$. For the 30 $M_{\oplus}$ case, the variation of the surface density profile is sufficient to reduce the effectiveness of the VSI in close proximity to the planet.

Finally, in the outer region, beginning at 2 $r_{\mathrm{p}}$, the VSI can operate unaffected by the planet and leads to an $\alpha$-parameter close to $10^{-3}$, which is approximately the value inferred from observations (see e.g. Andrews et al. 2009). Only for the largest planet the VSI is slightly suppressed due to the formation of a large vortex outside of the planetary orbit.

Altogether, we can see that in these simulations the overall amplitude and radial dependence of $\alpha$ is very similar to our previous runs, see for example Fig. 3 in Stoll & Kley (2014). The fact that $\alpha(r)$ drops towards smaller radii appears to be the result of two effects. First, the inner damping zone reduces the strength of the turbulence, second the spatial resolution of the models may not be sufficient as the result for a model with double resolution (shown for the 10 $M_{\oplus}$ planet and marked ’Res2’ in Fig. 1) indicates. More discussion on the spatial resolution will be given in section 3.5 below.

3.2 Vortices and vorticity

Vortex creation by embedded planets is a phenomenon occurring in nearly inviscid discs (Koller et al. 2003; de Val-Borro et al. 2007), while they are usually suppressed by viscosity over time (de Val-Borro et al. 2006). It is thus unclear what happens in VSI simulations that are able to generate a viscosity on the order of $\alpha=5\cdot 10^{-4}$. In fact, in the simulations with VSI and $M_{\mathrm{p}}=100\,\mathrm{M_{\oplus}}$ we see multiple vortices developing at the outer edge of the gap after the slow introduction of the planet (over $20\,\mathrm{orbits}$). They merge to a single vortex over the next $100\,\mathrm{orbits}$, which then remains stable, and even increases slightly in strength up to the end of the simulation. This vortex, near the end of the simulation111At $200\,\mathrm{orbits}$ the vortex sits at the periodic boundary, thus we show one earlier time slice, where it is in the middle of the domain., is shown in Fig. 3 on the right side of the lower panels. There we not only can see a vortex at the outer gap edge, but also a smaller vortex that has formed at the inner edge of the gap. Both vortices exist over the whole simulated timescale of $200\,\mathrm{orbits}$. This is in contrast with the $\alpha$-disc on the left side of the lower panel, where we can still see two weak vortices that are dissipating in the long run due to the effect of viscosity.

Now we compare directly the inviscid VSI-active simulations with the results of a viscous disc using $\alpha=5\cdot 10^{-4}$. The upper panel of Fig. 3 shows the simulations with $M_{\mathrm{p}}=30\,\mathrm{M_{\oplus}}$. As expected, the $\alpha$-viscosity disc does not show large scale vortices, but the VSI simulations again introduces a smaller vortex, which is also stable over the runtime of the simulation. The runs with smaller planets do not have a gap and thus do not show large scale vortices in this region, but instead smaller vortices are visible.

To further analyse the vortices we calculate a normalised vorticity

[TABLE]

For undisturbed keplerian rotation this leads to $\overline{\@vec{\omega}}_{\theta}=\@vec{\omega}_{\theta K}$ independent of radius, where $\@vec{\omega}_{\theta K}$ is the vorticity at $r=r_{\mathrm{p}}$.

The normalised vorticity in the disc’s midplane is displayed in Fig. 4 for various models. The second and third row show the normalised vorticity for the different planets in the VSI discs after $200\,\mathrm{orbits}$. The creation of small vortices in regions with active VSI turbulence was already noticed and explained by Richard et al. (2016). They are only short lived in locally isothermal discs (Richard et al. 2016) and are only weakly influenced by the lower mass planets in the main part of the disc. However, in its close vicinity the planet impacts the dynamics of the vortices, as shown in the first row where the vorticity of the two low mass planets is plotted averaged in time over 50 snapshots, beginning with orbit 100 and ending with orbit 200. We can see that the characteristic horseshoe motion of the flow transports the vorticity created close to the planet, which leads to an increase in the surface density along the inner horseshoe orbit. For the smallest planet with $5\,\mathrm{M_{\oplus}}$ this leads, on average, to an region of enhanced vorticity in close proximity behind the planet (top left panel in Fig. 4), which affects the torque on the planet as we will see in Sec. 3.4. For the $10\,\mathrm{M_{\oplus}}$ planet the effect becomes somewhat weaker but there is still an enhanced vorticity visible just inside of the planetary orbit spread out in azimuth.

For the larger planet masses with $30$ and $100\,\mathrm{M_{\oplus}}$ the presence of the planet leads to the formation of spiral arms which reduce the strength of the small vortices in the main part of the disc. On the other hand, the vortices created by the more massive planets, which were already visible in Fig. 3, are clearly identifiable as an increase in the vorticity. For the $100\,\mathrm{M_{\oplus}}$ model a vortex forms at the inner edge after $100\,\mathrm{orbits}$ at 0.8 $r_{\mathrm{p}}$ and starts to migrate inwards. The migration is initially slow, but it increases after leaving the gap edge. After further $100\,\mathrm{orbits}$ it migrates with a speed of approximately $0.001r_{\mathrm{p}}$ per orbit to 0.73 $r_{\mathrm{p}}$ where it can be seen in the figure, which is on the same order of magnitude as found by Richard et al. (2013). Vortex migration has already been studied by Paardekooper et al. (2010) and is due to the angular momentum transport by the excited density waves, that are asymmetric due to the density and vorticity gradient, leading to a net torque and migration. The vortex at the outer edge does not migrate notably, because the presence of the planet prevents inward migration.

The $30\,\mathrm{M_{\oplus}}$ model has also an increased vorticity at the edges of the gap. While no large vortices have formed yet, two weak vortices are visible at the outer edge of the gap (see also Fig. 3). This is in contrast to the $\alpha$-disc simulations, where for the $30\,\mathrm{M_{\oplus}}$ case no vortices are seen and only weak signatures can be noticed for the $100\,\mathrm{M_{\oplus}}$ case. The last row shows the vorticity for the $\alpha$-disc, where small vortices are present in the horseshoe region for the $30\,\mathrm{M_{\oplus}}$ model, since the relatively small viscosity cannot dissipate them on a short timescale. The variation in vorticity in the outer region is due to VSI activity, which is not suppressed completely by the viscosity leaving mainly the modes with small wavenumber. Similar to the VSI disc, we see stronger VSI activity in the $\alpha$-disc in the $100\,\mathrm{M_{\oplus}}$ model both in the inner and outer disc.

3.3 Surface density profile

An important observable that can be studied is the surface density profile which, modified by the presence of the planet and local instabilities, can show specific recognisable patterns. We present in Fig. 5 the disc surface density distribution, with respect to the initial model, for the different models and planetary masses, averaged over the azimuthal direction.

In the VSI-disc the planet is able to carve slightly deeper gaps with respect to the $\alpha$-disc models. Moreover, the density is increased at the inner gap edge. These effects are caused by an enhanced vorticity near the gap edges, due to the interplay of the VSI with the Rossby wave instability (see Richard et al. 2016). Surprisingly, even the smaller planets create a perturbation in the surface density profile close to the planet (especially in the VSI models). This is not a typical gap, but instead an increase in the surface density profile close the planet at the inner side of the disc and a decrease at the outer side. This comes from an increased vorticity in this region as we discussed in more detail in the Sec. 3.2. Finally, an overall slight decrease in the surface density in the outer region is visible for all VSI models, partly because they have been evolved for $200\,\mathrm{orbits}$ before adding the planet, and also due to higher $\alpha$-parameter in the outer region (see Fig. 1). For the larger planets the surface density profiles show very good agreement between the VSI models and the corresponding viscous $\alpha$-disc models. The gaps have the same widths and depths, confirming the estimate for the effective viscosity, $\alpha=5\cdot 10^{-4}$.

3.4 Torques acting on the planet

We calculate the torques using a smoothed force of the planet onto the disc

[TABLE]

where $d$ is the distance of a grid cell with mass $\Delta m$ from the planet and $\varepsilon=0.5$ the smoothing parameter. The smoothing reduces the potential inside the Hill radius, where we get too large contributions otherwise from the high gas density due to the missing accretion of the gas onto the planet. Moreover the density distribution close to the planet should be symmetrical with respect to the planet and thus it should not contribute to the total torque, but the finite resolution can make this numerically challenging.

We compare our results to the simulations of D’Angelo & Lubow (2010) and use the same normalization for the total torque

[TABLE]

and for the torque distribution per unit disc mass as a function of radius

[TABLE]

We display the torques for the VSI and $\alpha$-disc models over time in Fig. 6, where we quote the total torque averaged over the last $50\,\mathrm{orbits}$ in the legend. We smoothed the torques over $4\,\mathrm{orbits}$ with a Gaussian window, to reduce the strong fluctuations. The source of the perturbations are different depending on the size of the planet. For the smaller planets they are introduced through the turbulence injected by the VSI. In the case of the $100M_{\oplus}$ planet the more regular oscillations are introduced by a vortex near the outer edge of the gap generated by the planet. From the last panel in Fig. 6 we can estimate an oscillation period in this case of about 3.5 $P_{\mathrm{p}}$ which translates into a position of the vortex at $r\approx 1.25R_{\mathrm{p}}$, in good agreement with Fig. 3.

We compare the torques to the simulations of D’Angelo & Lubow (2010), who obtained for locally isothermal 3D discs the following empirical function

[TABLE]

where $\beta=-d\ln{\Sigma}/d\ln{r}$ and $\zeta=-d\ln{T}/d\ln{r}$, which correspond for our simulations to $\beta=0.5$, $\zeta=1$, and lead to an expected torque of $\Gamma=-2.15\Gamma_{0}$.

The torques on the planets in the viscous disc with $\alpha$-viscosity agree closely with the empirical function, with small differences caused possibly by the different smoothing methods of the planetary potential adopted. Only the $100\,\mathrm{M_{\oplus}}$ planet has a reduced torque due to gap formation.

In contrast, the simulations with VSI strongly disagree with the predictions. While it is not very surprising that the torques can be positive for a few orbits in the context of turbulent discs, it is clearly visible that the average of the torque is far from the predicted value. In case of the smallest planet with five Earth masses the average of the absolute torques is even five times larger than the torques in the $\alpha$-disc, leading to a very rapid inward migration. This effect becomes smaller with increasing planet mass because the larger planets begin to open gaps and generate larger vortices that are able to quench the VSI modes.

The time to reach an equilibrium torque is about 50 orbits for the $10M_{\oplus}$ planet and somewhat longer for the smallest one. To further analyse the stochastic nature of the perturbation on the planets and the resulting torque we display in Fig. 7 the occurrence rate of torque values for the two low mass planets. As shown, they are approximately shaped like a Gaussian and centered around their mean values that are quoted in the legend. The standard deviations derived from fitting the Gaussians are 5.0 and 3.7 for the two low mass planets, such that the strength of the torque fluctuations is comparable to the mean torque for both cases. This is different to simulations of embedded planets in MHD turbulent disc. In the simulations performed for example by Nelson (2005) and Uribe et al. (2011) it is shown that the stochastic part of the torque can exceed the mean value significantly and the timescale for reaching equilibrium is also longer in MHD turbulence. The reason for this difference lies in the fact that the density fluctuations are stronger in the MHD turbulence compared to the VSI case.

To shed more light on this torque increase for the low mass planets we also show the torque distribution per unit disc mass as a function of radius in Fig. 8. Here we average again over the last $50\,\mathrm{orbits}$ and we compare to D’Angelo & Lubow (2008) who give an analytical fit to the radial torque distribution. In the lower panel of Fig. 8 we can again see a good agreement with their results for the viscous disc, except for the heaviest planet. In the simulation with $M_{\mathrm{p}}=5\,\mathrm{M_{\oplus}}$ the torque is five times stronger with active VSI compared to the simulation with $\alpha$-viscosity. This is due to a combination of reduced outer torque and stronger inner torque. This cannot be explained by a variation in the radial surface density, because we could see in Fig. 5 that the density in the outer region is reduced and the density in the inner region is enhanced, leading to the opposite effect. Instead the surface density depends on the azimuthal direction, and there we find an increase in density on the inner side behind the planet and a decrease in density at the outer side in front of the planet. This can be seen in Fig. 9 where we average the density over $100\,\mathrm{orbits}$. The density enhancement is due to the increased vorticity behind the planet, resulting from the interplay between the planet and the VSI as can be seen in the upper panels of Fig. 4.

3.5 A model with double resolution

To check the convergence with resolution, we repeat the simulation for $M_{\mathrm{p}}=10\,\mathrm{M_{\oplus}}$ and VSI with doubled resolution in every direction. The resolution is then $1200\times 256\times 2048$ in (r, $\theta$, $\varphi$) or (15,11,7) cells per $R_{\mathrm{Hill}}$ for this planet. The overall results show no important difference to the results of the simulations with lower resolution. A notable difference is in the $\alpha$-profile, which we present in Fig. 1. We can see that in the inner region the $\alpha$-parameter is larger for the better resolved simulation. This indicates that the simulations using the standard resolution are not properly resolved in the innermost region between 0.5-0.64 $r_{\mathrm{p}}$, which could be a consequence of the fact that the radial extent of the VSI cells increases stronger than linear with radius, as shown in Stoll & Kley (2014). This also indicates that the rest of the domain is sufficiently resolved and not dependent on resolution, in contrast to the 2D simulations in Stoll & Kley (2014). This does not pose a problem for the calculations of the torques, which are only important close to the planet. This can be seen in Fig. 8, were the dotted black line represents the simulation with doubled resolution which can be compared to the green line. Both lines are as close together as one can expect from simulations with turbulence. No other notable differences have been observed.

4 Discussion and Conclusions

We performed a series of hydrodynamical simulations to study the interplay between a growing planet embedded in a protoplanetary disc and the VSI, which is a strong candidate for the main angular momentum transport mechanism in the region of active planet formation. The main results that can be drawn from our study are:

Vortices and VSI

Larger planetary cores ($M_{\textrm{p}}>10\,\mathrm{M_{\oplus}}$) are able to open easily a gap in the disc, which leads to the generation of vortices at the gap edges. These vortices are observed both in VSI and $\alpha$-disc models, however, due to the interaction between the vortices and the VSI, their lifetime is much longer in the VSI models. In our simulations we notice the presence of several weak spiral arms, for example in Fig. 3 (upper row), that can be associated with vortices generated by the Rossby-Wave Instability (RWI) near the edges of the gap created by the planet. They are related to a maximum in the normalised vorticity profile, shown in Fig. 4. Furthermore, the vorticity perturbations generated by the VSI, as suggested by Richard et al. (2016), are effectively increasing the lifetime of the RWI vortices, as one can see comparing the viscous and turbulent runs in the bottom panel of Fig. 4. 2. 2.

The VSI active region

The region with active VSI is not affected sensibly by the embedded planets for planetary masses less than $30\,\mathrm{M_{\oplus}}$. For higher masses the presence of vortices is able to suppress the VSI near their position. On the other hand the active region is extended inwards, at least for the most massive planet studied (see Fig. 4). One possibility for this extended turbulent activity might by the action of the spiral wave instability (SWI) where the presence of a massive planet can excite unstable inertial-gravity waves (Bae et al. 2016b, a) which might enhance the turbulent level above that of the VSI alone. We have analysed the time evolution of the vertical kinetic energy and 2D slices of the vertical velocity but could not find a clear indications for the operation of the SWI. A possible other explanation for this effect is the change in surface density profile due to the creation of a gap. This will in turn change the angular momentum profile and thus influence the VSI. 3. 3.

Planet migration

Concerning the migration of planets we find that the torques acting onto small planets for the $\alpha$-disc models are in good agreement with the prediction of D’Angelo & Lubow (2010). This picture changes drastically for the inviscid models with the VSI operating. For small planets with masses smaller than about $10\,\mathrm{M_{\oplus}}$ the time averaged torques are negative throughout. For the lowest planet mass we studied, $5\,\mathrm{M_{\oplus}}$, the inward migration is about 5 times faster than for the $\alpha$-disc model which is due to a density enhancement directly behind the planet. For increasing planet masses the migration rates approach those of the viscous disc. Due to this special disturbance in the density we do not find any prolonged phases of outward migration as seen for example for planets embedded in MHD-turbulent discs (Nelson 2005; Uribe et al. 2011). In general, because the density fluctuations in VSI turbulent discs are smaller than in MHD turbulence, the stochastic component in the migration process is weaker in VSI discs. Even though the stochastic component does not play a major role for our small planet case with $5\,\mathrm{M_{\oplus}}$, it may nevertheless become important for sub-Earth mass objects. As three-dimensional inviscid simulations of disc will always generate turbulence through the VSI, embedded planets will experience substantial inward migration. This feature cannot be captured by modelling planet migration in two-dimensional inviscid discs that show typically a stalling of migration after a sufficiently wide gap has been formed (Li et al. 2009; Fung & Chiang 2017). In contrast, realistic three-dimensional discs will always show a finite turbulence level that will limit the gap depth. In fact, from our simulations we see that the width and depths of gaps opened by the planets are very similar in the VSI models and the corresponding $\alpha$-disc models. In addition, small planets ($M_{\textrm{p}}\leq 10\,\mathrm{M_{\oplus}}$) that do not open significant gaps are able to modify the surface density profile close to their location which leads to the enhanced migration for the VSI turbulent discs. 4. 4.

Vortex migration and its implications

Vortex migration was first observed by Paardekooper et al. (2010) in isothermal discs, and then further studied by Faure et al. (2015) for radiative discs. Vortices are able to generate spiral arms by compressing the flow around them, and the migration occurs because of an asymmetry in the position of the sonic lines that generates the density waves. According to Richard et al. (2013) the migration speed is directly linked with the vortex aspect ratio ($\chi$), it slows down upon increasing $\chi$. In particular, Faure & Nelson (2016) studied the interaction of a planet with a vortex generated through RWI. Interestingly, they found that a massive vortex can drag the planet with it during the migration process, and it can also cross the planetary gap periodically crossing its location.

In our simulations for the $100\,\mathrm{M_{\oplus}}$ planet, the inner vortex has a migration rate after an initial adjustment on the order of $0.001\,\mathrm{r_{\mathrm{p}}/orbit}$, in accordance with Richard et al. (2013), while the outer one is kept in its orbit by the steep surface density profile carved by the planet. The aspect ratio of the inner vortex is around $3$ which, according to Richard et al. (2013) should be destroyed after a few orbital periods due to elliptical instability (for $\chi<4$). However, its lifetime is much more extended in our runs, meaning that the presence of the VSI is beneficial extending its life. On the other hand the outer vortex is much broader and has a larger aspect ratio. Both vortices have also a considerable vertical extent, and only in the outer corona they are effectively dissipated. 5. 5.

VSI as an angular momentum driving process

In agreement with previous simulations of VSI unstable discs Nelson et al. (2013), we find that an angular momentum transport driven by the VSI that corresponds to an $\alpha=5\cdot 10^{-4}$ which is not strongly affected by embedded planets. Hence, the VSI constitutes a viable candidate for the generation of turbulence in discs where the MRI may be inactive. Recently, we have shown (Stoll et al. 2017) that the turbulence generated by the VSI is anisotropic and can be described by a viscous ansatz for the stress tensor using two different coefficients for the radial and vertical angular momentum transport. It remains to be seen how an embedded planet behaves in such discs with anisotropic viscosities.

Thus, the VSI can on the one hand increase the strength of the vortices forming close to the gap edges for high planetary masses, and on the other hand boost the Type I migration for small planetary cores. The present study can be further improved by relaxing the locally isothermal assumption, and perform radiative simulations over a similar radial range to check the strength of the VSI turbulence and its impact on embedded planets under more realistic conditions. Additionally, in further improvements the planet should be allowed to accrete mass and migrate through disc.

Acknowledgements.

We thank an anonymous referee for useful comments and suggestions. G. Picogna acknowledges the support through the German Research Foundation (”Deutsche Forschungsgemeinschaft”, DFG) grant KL 650/21 within the collaborative research program The first 10 Million Years of the Solar System. M.H.R. Stoll acknowledges the support through the (DFG) grant KL 650/16. This work was performed on the computational resource ForHLR I funded by the Ministry of Science, Research and the Arts of Baden-Württemberg and the DFG.