EDP Sciences
Free Access
Issue
A&A
Volume 604, August 2017
Article Number A28
Number of page(s) 9
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201730668
Published online 27 July 2017

© ESO, 2017

1. Introduction

The long-standing problems of identifying the mechanism responsible for transporting angular momentum in accretion discs and of explaining their observed lifetimes have been tackled since the early works of Shakura & Sunyaev (1973) and Lynden-Bell & Pringle (1974) were published. The physical process that has been identified as 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 that well-ionised media trigger the MRI is generally not fulfilled in protoplanetary discs (Gammie 1996).

Several purely hydrodynamical instabilities have been proposed as the drivers of angular momentum transport, 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 more generally applicable mechanisms 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 also applied 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) and in fully radiative discs (Stoll & Kley 2014). Within the context of the planet formation process, these instabilities can strongly impact the dust dynamics, for example by increasing 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 on the planet and force the planet to migrate through the disc (see e.g. Kley & Nelson 2012, and references therein). The direction and magnitude of this migration depend 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 on 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 α-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 Sect. 2 we describe the numerical method used to model the disc, the planet, and the parameter space studied. In Sect. 3 we outline the results of our analysis, and in Sect. 4 we discuss our major findings and draw our conclusions.

2. Set-up

We used the pluto code (Mignone et al. 2007) to model a 3D section of a locally isothermal accretion disc in spherical coordinates (r, θ, φ) 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 α-value adapted to that of the VSI turbulent disc. Radially, the disc section extends from 0.4 rp to 2.5 rp, where rp = 5.2 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, centred 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 Nr = 600,Nθ = 128, and Nφ = 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 Δr = 0.003, and the ratio of gridsizes near the location of the planet is given by Δr:rΔθ:rΔφ = 1.0:1.28:2.012. A comparison model with double spatial resolution is discussed in Sect. 3.5. The main parameters of the simulations are summarised in Table 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,φ). The initial density profile created by force equilibrium is given by (1)where ρ0 is the gas mid-plane density at R = Rp = rp, p = −1.5 the density exponent, and cs the isothermal sound speed. In our simulations we use ρ0 = 2.07 × 10-11 g/cm3, which translates into Σ(R0) = 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 (2)with q = −1 and T0 = 121 K.

The gas moves with an azimuthal velocity given by the Keplerian speed around a 1 M star, corrected by the pressure support (Nelson et al. 2013)(3)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] rp and [2.3, 2.5] rp. For the kinematic viscosity in the α-disc model we use ν = 2/3αcsH, and adopt a constant α parameter of α = 5 × 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 orbits before embedding the planet.

Table 1

Model parameters.

2.2. The planet

We embed a planet that orbits its parent star on a circular orbit with a radius rp = 1 and a mass in the range [5,10,30,100] M. These four planets have the following Hill radii: [0.017,0.022,0.031,0.045] Rp. 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 is not affected by the disc so its orbital parameters remain fixed during the whole simulation. The gravitational potential of the planet is implemented with a cubic expansion in the vicinity of the planet location (Klahr & Kley 2006)(4)where d is the distance of a gas element from the planet location, and the smoothing length drsm = 0.5 RHill 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 orbits of the simulation. We ran the simulations for 200 orbital periods of the planet.

3. Results

3.1. Effective viscosity of the disc with embedded planets

For the VSI models the α-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 α-parameter through the Reynolds stress resulting from the turbulence. We calculate the Reynolds stress in cylindrical coordinates (R,Z,φ) (5)where δuφ is the local difference from the equilibrium azimuthal velocity, which we calculate by time averaging uφ for each grid cell from orbit 60 to orbit 200 over 70 snap shots.

We then calculate the dimensionless α-parameter as a function of the radius, (6)where P is the pressure and ⟨⟩t,θ,φ denotes the average over time (140 orbits) and the whole vertical and azimuthal domain.

In Fig. 1 we compare the α-parameter for our simulations of VSI-active discs with embedded planets of different mass. The inner damping region from 0.4 rp to 0.5 rp 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 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 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.5cs0 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.

thumbnail Fig. 1

Comparison of the α-parameter for our VSI-active simulations with embedded planets of different mass. The planet is located at r = rp. “Res2” indicates a model with double resolution. The dotted line is at α = 5 × 10-4.

Open with DEXTER

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 α-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] rp. For the 30 M case, the variation of the surface density profile is sufficient to reduce the effectiveness of the VSI in close proximity to the planet.

thumbnail Fig. 2

Vertical velocity in the midplane of the VSI disc after 198 orbits with embedded planets of different mass. Upward motion is indicated in red and downwards motion in blue. Larger planets clearly disrupt the usual VSI modes close to the planet. In the last panel (100 M planet) vortices are visible at both sides of the planet. Also visible in the last panel is the more extended VSI activity near the inner boundary compared to the other panels.

Open with DEXTER

Finally, in the outer region, beginning at 2 rp, the VSI can operate unaffected by the planet and leads to an α-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 is the VSI slightly suppressed, due to the formation of a large vortex outside of the planetary orbit.

In these simulations we can see that the overall amplitude and radial dependence of α is very similar to our previous runs (see e.g. Fig. 3 in Stoll & Kley 2014). The fact that α(r) drops towards smaller radii appears to be the result of two effects. First, the inner damping zone reduces the strength of the turbulence and second, the spatial resolution of the models may not be sufficient, as the result for a model with double resolution indicates (shown for the 10 M planet and marked “Res2” in Fig. 1). More discussion on the spatial resolution is given in Sect. 3.5.

3.2. Vortices and vorticity

thumbnail Fig. 3

Density in the midplane relative to the initial density after 198 orbits for the inviscid VSI-active simulations and viscous disc models using α = 5 × 10-4. Upper panels: after 198 orbits for Mp = 30 M. Lower panels: same for Mp = 100 M.

Open with DEXTER

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 α = 5 × 10-4. In fact, in the simulations with VSI and Mp = 100 M we see multiple vortices developing at the outer edge of the gap after the slow introduction of the planet (over 20 orbits). They merge to a single vortex over the next 100 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 simulation1, is shown in Fig. 3 on the right side of the lower panels. Not only do we 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 orbits. This is in contrast to the α-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 α = 5 × 10-4. The upper panel of Fig. 3 shows the simulations with Mp = 30 M. As expected, the α-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: (7)For undisturbed Keplerian rotation this leads to independent of radius, where ωθK is the vorticity at r = rp.

thumbnail Fig. 4

First row: normalised vorticity (Eq. (7) in the disc midplane averaged over the last 100 orbits in the VSI disc. The increased vorticity close to the planet is clearly visible. Second and third rows: normalised vorticity in the midplane for the VSI disc after 198 orbits. Last row: normalised vorticity in the midplane for the viscous α-disc after 198 orbits.

Open with DEXTER

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 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 M this leads, on average, to a 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 see in Sect. 3.4. For the 10 M 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 M 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 M model a vortex forms at the inner edge after 100 orbits at 0.8 rp and starts to migrate inwards. The migration is initially slow, but it increases after leaving the gap edge. After a further 100 orbits it migrates at a speed of approximately 0.001rp per orbit to 0.73 rp, where it can be seen in the figure, which is of 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, which 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 M model also has 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 α-disc simulations, where for the 30 M case no vortices are seen and only weak signatures can be noticed for the 100 M case. The last row shows the vorticity for the α-disc, where small vortices are present in the horseshoe region for the 30 M 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 α-disc in the 100 M 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.

thumbnail Fig. 5

Surface density for the different simulations after 200 orbits. The vertical dotted black lines indicate the position of the planet and the initial surface density profile is marked by the solid black line. For each planet mass the VSI models are compared to the α-disc cases.

Open with DEXTER

In the VSI disc the planet is able to carve slightly deeper gaps with respect to the α-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 discuss in more detail in Sect. 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 orbits before adding the planet and also due to higher α-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 α-disc models. The gaps have the same widths and depths, confirming the estimate for the effective viscosity, α = 5 × 10-4.

3.4. Torques acting on the planet

We calculate the torques using a smoothed force of the planet onto the disc (8)where d is the distance of a grid cell with mass Δm from the planet and ε = 0.5 the smoothing parameter. The smoothing reduces the potential inside the Hill radius, where otherwise we get excessive contributions 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 normalisation for the total torque (9)and for the torque distribution per unit disc mass as a function of radius (10)We display the torques for the VSI and α-disc models over time in Fig. 6, where we quote the total torque averaged over the last 50 orbits in the legend. We smoothed the torques over four orbits with a Gaussian window, to reduce the strong fluctuations. The sources 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 100 M 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 Pp, which translates into a position of the vortex at r ≈ 1.25 Rp, in good agreement with Fig. 3.

thumbnail Fig. 6

Torque acting on the planet over time normalised by Γ0. We smooth over four orbits with a Gaussian window and compare it to the simulation of D’Angelo & Lubow (2010, red dashed line). In the legend we give the total torque averaged over the last 50 orbits. Owing to the turbulence in the VSI disc the torque strongly fluctuates, while for the α-disc the evolution is much smoother. The oscillations in the last panel for the 100 M planet are due to the presence of a vortex.

Open with DEXTER

We compare the torques to the simulations of D’Angelo & Lubow (2010), who obtained for locally isothermal 3D discs the empirical function (11)where β = −dlnΣ/dlnr and ζ = −dlnT/ dlnr, which correspond for our simulations to β = 0.5, ζ = 1, and lead to an expected torque of Γ = −2.15Γ0.

The torques on the planets in the viscous disc with α-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 M 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 the case of the smallest planet with five Earth masses the average of the absolute torques is even five times higher than the torques in the α-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.

thumbnail Fig. 7

Histogram taken over 1000 samples of the occurrence of certain torque values for the two low mass planets in the VSI turbulent disc. The data were collected for the last 100 orbits of the simulations. The torques are normalised to Γ0 in each case and the average values and standard deviation over this time span are quoted in legend.

Open with DEXTER

The time to reach an equilibrium torque is about 50 orbits for the 10 M 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 centred around their mean values, which are quoted in the legend. The standard deviations derived from fitting the Gaussians are 5.0 and 3.7 Γ0 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 orbits and we compare the result to that of 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 Mp = 5 M the torque is five times stronger with active VSI than in the simulation with α-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 can 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 on the outer side in front of the planet. This can be seen in Fig. 9 where we average the density over 100 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.

thumbnail Fig. 8

Torques on the planets normalised by (dΓ/dM)0, see Eq. (10). Upper panel: planet embedded into a disc with VSI. The inner torques are smaller then expected and the outer torques are stronger. Lower panel: planet embedded into a disc with generic α-viscosity. The torques are very similar to the results of D’Angelo & Lubow (2010). In both simulations the torques for the heaviest planet are reduced, due to the lower surface density in the gap.

Open with DEXTER

3.5. Model with double resolution

To check the convergence with resolution, we repeated the simulation for Mp = 10 M and VSI with doubled resolution in every direction. The resolution was then 1200 × 256 × 2048 in (r, θ, ϕ) or (15, 11, 7) cells per RHill 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 α-profile, which we present in Fig. 1. We can see that in the inner region the α-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 rp, 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, where 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 can be expected from simulations with turbulence. No other notable differences have been observed.

thumbnail Fig. 9

Upper panel: density averaged over the last 100 orbits for the VSI simulation with Mp = 5 M in units of initial density. Lower panel: density averaged over the last 100 orbits for the VSI simulation with Mp = 10 M in units of initial density. In both cases an increase in density before the planet is visible, if only barely for the 10 M case. The bump in density does not move over time.

Open with DEXTER

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 as follows:

  • 1.

    Vortices and VSI Larger planetary cores(Mp> 10 M) are easily able to open a gap in the disc, which leads to the generation of vortices at the gap edges. These vortices are observed both in VSI and α-disc models; however, due to the interaction between the vortices and the VSI, their lifetimes are 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), which 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 can be seen by comparing the viscous and turbulent runs in the bottom panel of Fig. 4.

  • 2.

    The VSI active region The region with active VSI is not affected sensibly by the embedded planets for planetary masses less than 30 M. 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 be 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. Another possible 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.

  • Planet migration Concerning the migration of planets we find that the torques acting on small planets for the α-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 M, the time-averaged torques are negative throughout. For the lowest planet mass we studied, 5 M, the inward migration is about five times faster than for the α-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 in our small planet case (5 M), it may nevertheless become important for sub-Earth mass objects. As 3D inviscid simulations of discs will always generate turbulence through the VSI, embedded planets will experience substantial inward migration. This feature cannot be captured by modelling planet migration in 2D inviscid discs that typically show a stalling of migration after a sufficiently wide gap has been formed (Li et al. 2009; Fung & Chiang 2017). In contrast, realistic 3D 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 α-disc models. In addition, small planets (Mp ≤ 10 M) 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.

  • 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 (χ): it slows down upon increasing χ. 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 M planet, the inner vortex has a migration rate after an initial adjustment of the order of 0.001 rp/orbit, in agreement with Richard et al. (2013), while the outer vortex 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 χ< 4). However, its lifetime is much more extended in our runs, meaning that the presence of the VSI is beneficial and extends its life. On the other hand, the outer vortex is much broader and has a larger aspect ratio. Both vortices also have a considerable vertical extent, and only in the outer corona are they effectively dissipated.

  • 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 corresponds to an α = 5 × 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 performing 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 the disc.


1

At 200 orbits the vortex sits at the periodic boundary, thus we show one earlier time slice where it is in the middle of the domain.

Acknowledgments

We thank the anonymous referee for the 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.

References

All Tables

Table 1

Model parameters.

All Figures

thumbnail Fig. 1

Comparison of the α-parameter for our VSI-active simulations with embedded planets of different mass. The planet is located at r = rp. “Res2” indicates a model with double resolution. The dotted line is at α = 5 × 10-4.

Open with DEXTER
In the text
thumbnail Fig. 2

Vertical velocity in the midplane of the VSI disc after 198 orbits with embedded planets of different mass. Upward motion is indicated in red and downwards motion in blue. Larger planets clearly disrupt the usual VSI modes close to the planet. In the last panel (100 M planet) vortices are visible at both sides of the planet. Also visible in the last panel is the more extended VSI activity near the inner boundary compared to the other panels.

Open with DEXTER
In the text
thumbnail Fig. 3

Density in the midplane relative to the initial density after 198 orbits for the inviscid VSI-active simulations and viscous disc models using α = 5 × 10-4. Upper panels: after 198 orbits for Mp = 30 M. Lower panels: same for Mp = 100 M.

Open with DEXTER
In the text
thumbnail Fig. 4

First row: normalised vorticity (Eq. (7) in the disc midplane averaged over the last 100 orbits in the VSI disc. The increased vorticity close to the planet is clearly visible. Second and third rows: normalised vorticity in the midplane for the VSI disc after 198 orbits. Last row: normalised vorticity in the midplane for the viscous α-disc after 198 orbits.

Open with DEXTER
In the text
thumbnail Fig. 5

Surface density for the different simulations after 200 orbits. The vertical dotted black lines indicate the position of the planet and the initial surface density profile is marked by the solid black line. For each planet mass the VSI models are compared to the α-disc cases.

Open with DEXTER
In the text
thumbnail Fig. 6

Torque acting on the planet over time normalised by Γ0. We smooth over four orbits with a Gaussian window and compare it to the simulation of D’Angelo & Lubow (2010, red dashed line). In the legend we give the total torque averaged over the last 50 orbits. Owing to the turbulence in the VSI disc the torque strongly fluctuates, while for the α-disc the evolution is much smoother. The oscillations in the last panel for the 100 M planet are due to the presence of a vortex.

Open with DEXTER
In the text
thumbnail Fig. 7

Histogram taken over 1000 samples of the occurrence of certain torque values for the two low mass planets in the VSI turbulent disc. The data were collected for the last 100 orbits of the simulations. The torques are normalised to Γ0 in each case and the average values and standard deviation over this time span are quoted in legend.

Open with DEXTER
In the text
thumbnail Fig. 8

Torques on the planets normalised by (dΓ/dM)0, see Eq. (10). Upper panel: planet embedded into a disc with VSI. The inner torques are smaller then expected and the outer torques are stronger. Lower panel: planet embedded into a disc with generic α-viscosity. The torques are very similar to the results of D’Angelo & Lubow (2010). In both simulations the torques for the heaviest planet are reduced, due to the lower surface density in the gap.

Open with DEXTER
In the text
thumbnail Fig. 9

Upper panel: density averaged over the last 100 orbits for the VSI simulation with Mp = 5 M in units of initial density. Lower panel: density averaged over the last 100 orbits for the VSI simulation with Mp = 10 M in units of initial density. In both cases an increase in density before the planet is visible, if only barely for the 10 M case. The bump in density does not move over time.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.