The role of gap edge instabilities in setting the depth of planet gaps in protoplanetary discs

TL;DR

This study investigates how gap edge instabilities, especially Rossby wave instability, influence the depth of planet-induced gaps in protoplanetary discs, revealing the importance of two-dimensional effects and torque distribution modeling.

Contribution

It demonstrates that gap edge instabilities critically affect gap depth and highlights the importance of two-dimensional modeling and improved torque approximations for accurate predictions.

Findings

Two-dimensional gaps are shallower than one-dimensional gaps for the same initial conditions.

Rossby wave instability influences the equilibrium gap depth in two-dimensional discs.

Improved torque density models from three-dimensional simulations reduce discrepancies between 1D and 2D models.

Abstract

It is known that an embedded massive planet will open a gap in a protoplanetary disc via angular momentum exchange with the disc material. The resulting surface density profile of the disc is investigated for one dimensional and two dimensional disc models and, in agreement with previous work, it is found that one dimensional gaps are significantly deeper than their two dimensional counterparts for the same initial conditions. We find, by applying one dimensional torque density distributions to two dimensional discs containing no planet, that the excitement of the Rossby wave instability and the formation of Rossby vortices play a critical role in setting the equilibrium depth of the gap. Being a two dimensional instability, this is absent from one dimensional simulations and does not limit the equilibrium gap depth there. We find similar gap depths between two dimensional gaps formed by torque density distributions, in which the Rossby wave instability is present, and two dimensional planet gaps, in which no Rossby wave instability is present. This can be understood if the planet gap is maintained at marginal stability, even when there is no obvious Rossby wave instability present. Further investigation shows the final equilibrium gap depth is very sensitive to the form of the applied torque density distribution, and using improved one dimensional approximations from three dimensional simulations can go even further to reducing the discrepancy between one and two dimensional models, especially for lower mass planets. This behaviour is found to be consistent across discs with varying parameters.

Tables

Table 1. Table 1: Values of the parameter p in Equation 12 .

$p_n$	Value
$p_1$	$0.0297597$
$p_2$	$1.09770$
$p_3$	$0.938567$
$p_4$	$0.0421186$
$p_5$	$0.902328$
$p_6$	$1.03579$
$p_7$	$0.0981183$
$p_8$	$4.68108$

Equations

$$\frac{\partial \Sigma }{\partial t}+\nabla \cdot (\Sigma \mathbf{v})=0,$$

$$\Sigma \left(\frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}\cdot \nabla \mathbf{v}\right)=-\nabla P-\nabla \cdot \mathbf{T}-\Sigma \nabla \Phi ,$$

$$\begin{gathered} \frac{\partial v_R}{\partial t}+v_R\frac{\partial v_R}{\partial R}+\Omega \frac{\partial v_R}{\partial \varphi }-R{\Omega }^2=-\frac{1}{\Sigma }\frac{\partial P}{\partial R} \\ -\frac{2}{\Sigma R}\frac{\partial }{\partial R}\left(\nu \Sigma R\frac{\partial v_R}{\partial R}\right)-\frac{1}{\Sigma R}\frac{\partial }{\partial \varphi }\left[\nu \Sigma \left(R\frac{\partial \Omega }{\partial R}+\frac{1}{R}\frac{\partial v_R}{\partial \varphi }\right)\right] \\ -\frac{\partial \Phi }{\partial R} \end{gathered}$$

$$\begin{gathered} \frac{\partial (R\Omega )}{\partial t}+v_{R}R\frac{\partial \Omega }{\partial R}+\Omega \frac{\partial (R\Omega )}{\partial \varphi }+2v_R\Omega =-\frac{1}{R\Sigma }\frac{\partial P}{\partial \varphi } \\ -\frac{2}{\Sigma R}\frac{\partial }{\partial \varphi }\left(\nu \Sigma \frac{\partial \Omega }{\partial \varphi }\right)+\frac{1}{\Sigma R^2}\frac{\partial }{\partial R}\left[\nu \Sigma R^2\left(R\frac{\partial \Omega }{\partial R}+\frac{1}{R}\frac{\partial v_R}{\partial \varphi }\right)\right] \\ -\frac{1}{R}\frac{\partial \Phi }{\partial \varphi }, \end{gathered}$$

$$\Sigma R\frac{\partial (R\Omega )}{\partial t}+\Sigma v_R\frac{\partial (R^2\Omega )}{\partial R}=\frac{1}{R}\frac{\partial }{\partial R}\left(\nu \Sigma R^3\frac{\partial \Omega }{\partial R}\right)+\Sigma \Lambda .$$

$$\Sigma v_{R}R=-3R^{\frac{1}{2}}\frac{\partial }{\partial R}(\nu \Sigma R^{\frac{1}{2}})+\frac{2\Sigma \Lambda R^{\frac{3}{2}}}{(GM{)}^{\frac{1}{2}}}.$$

$$R\frac{\partial \Sigma }{\partial t}+\frac{\partial }{\partial R}(\Sigma v_{R}R)=0,$$

$$\frac{\partial \Sigma }{\partial t}=\frac{1}{R}\frac{\partial }{\partial R}\left[3R^{\frac{1}{2}}\frac{\partial }{\partial R}\left(\nu \Sigma R^{\frac{1}{2}}\right)-\frac{2\Lambda \Sigma R^{\frac{3}{2}}}{(GM{)}^{\frac{1}{2}}}\right].$$

$$\Lambda =\text{sign}(R-R_0)\frac{f{q}^{2}GM}{2R}{\left(\frac{R}{\left|{\Delta }_0\right|}\right)}^4,$$

$$\text{sign}(R-R_0)=\{\begin{array}{ll}-1,& \text{if}R-R_0<-H\\ \frac{R-R_0}{H},& \text{if}\left|R-R_0\right|\le H\\ 1,& \text{if}R-R_0>H\end{array}$$

$$\Lambda =\mathcal{F}\left(x,\beta ,\zeta \right){\Omega }^2{R}_0^2{q}^2{\left(\frac{R_0}{H}\right)}^4,$$

$$\begin{gathered} \mathcal{F}(x,\beta ,\zeta )=\left(p_1{e}^{\left(-\frac{(x+p_2{)}^2}{p_3^2}\right)}+p_4{e}^{\left(-\frac{(x-p_5{)}^2}{p_6^2}\right)}\right) \\ \cdot \text{tanh}(p_7-p_{8}x). \end{gathered}$$

$$\Lambda =\frac{{\int }_{-\pi }^{\pi }\Sigma (R,\varphi )\frac{\partial {\Phi }_p}{\partial \varphi }d\varphi }{2\pi {\Sigma }_{Av}(R)},$$

$${\Phi }_p=-\frac{G{M}_p}{\sqrt{R^2+R_0^2-2R{R}_0\cos (\varphi -{\varphi }_0)+s^2}},$$

$$\frac{{\Sigma }_{\text{Gap}}}{{\Sigma }_0}=0.14{\left(\frac{q}{10^{-3}}\right)}^{-2.16}{\left(\frac{\alpha }{10^{-2}}\right)}^{1.41}{\left(\frac{h}{0.05}\right)}^{6.61},$$

Full text

††pubyear: 2016††pagerange: The role of gap edge instabilities in setting the depth of planet gaps in protoplanetary discs.–The role of gap edge instabilities in setting the depth of planet gaps in protoplanetary discs.

1 Introduction

A planet situated in a protoplanetary disc will excite density waves in the disc material, which transport angular momentum away from the planet. The angular momentum is then deposited in the disc as these waves dissipate. Therefore, the planet exerts a torque on the surrounding disc material, forming a gap when the torque is strong enough.

For a planet to open a gap it must fulfil two conditions. Firstly, it must be massive enough (Rafikov, 2002; Crida et al., 2006) such that this torque overcomes the viscous torque from the disc material attempting to diffuse back into the low density gap region, e.g. Kley & Nelson (2012). This is known as the viscous condition and it is the balance of these torques that sets the equilibrium profile of the gap (Lin & Papaloizou, 1979; Goldreich & Tremaine, 1980; Takeuchi et al., 1996; Crida et al., 2006). Secondly, the planet’s Hill radius must be larger than the disc scale height, such that the excited density waves deposit angular momentum close to the planet (Lin & Papaloizou, 1993; Bryden et al., 1999). This is known as the thermal criterion for gap opening. It has also been shown that for small viscosities low mass planets can open gaps further away in the disc, where the excited density wave shocks (Goodman & Rafikov, 2001). Gap opening has an important impact on both theoretical and observational studies of protoplanetary discs.

Recent results from the Atacama Large Millimetre Array (ALMA) have provided our first high resolution images of protoplanetary discs and the features present across the course of their lifetime (ALMA Partnership et al., 2015; Andrews et al., 2016). Of particular note is the presence of a number of gaps visible in the disc. We cannot currently say with certainty whether or not these gaps are planetary in origin, due to lack of resolution and understanding of possible gap opening processes. Despite this it is possible that a number of the gaps we see could be the result of planetary interaction and this furthers the desire to understand this phenomena (Dipierro et al., 2015; Gonzalez et al., 2015; Zhang et al., 2015; Rapson et al., 2015).

In the context of planet formation theory, gap formation plays a vital role in the transition between two regimes of planet migration, Type I and Type II. Type I migration occurs when the planet is not massive enough to open a gap in the disc, hence the disc is almost unperturbed by the presence of a planet. Type II migration occurs when a more massive planet opens a gap in the disc. The migration time-scale for the planet will then depend on the viscous evolution time-scale of the disc (Lin & Papaloizou, 1986). Type I migration is considerably faster than Type II migration (Ward, 1997). These types of migration are important to theoretical models used to explain planetary system formation as they dictate the final position of a planet and can help to explain the features of the numerous exoplanetary systems discovered. Gap formation is also a very important factor when considering how massive a planet can grow. When an accreting planet becomes massive enough to open a gap in the disc its accretion rate will fall significantly, as the planet has emptied its accretion zone. Hence opening a gap can act as a limit on the final mass of a planet. These factors play an important role in population synthesis models (Ida & Lin, 2008; Mordasini et al., 2009) and also the more recent pebble accretion models (Bitsch et al., 2015; Ali-Dib, 2016).

Gap formation as a result of planet-disc interactions is an extensively studied subject area, and has been explored in one, two and three dimensions. The results of one dimensional models show thin deep gaps (Lin & Papaloizou, 1986; Kanagawa et al., 2015), comparative to their higher dimensional analogues. In two dimensions gap formation criteria and shape have been studied (Crida et al., 2006; Duffell & MacFadyen, 2013) and a power law scaling relation for gap depth has been deduced (Fung et al., 2014; Kanagawa et al., 2015). It is well known that the depth and shape of two dimensional gaps from numerical simulations are inconsistent with expected results from one dimensional analytical and numerical models (Crida et al., 2006; Kanagawa et al., 2015). Crida et al. (2006), however, show that including the pressure torque in the balance of torques on the discs returns good agreement between semi-analytical and numerical gap profiles. The differences between two and three dimensional simulations are a less explored area, but nonetheless has received some investigation (Fung & Chiang, 2016) and improved one dimensional models of torque distributions have been determined from three dimensional hydrodynamic simulations (D’Angelo & Lubow, 2010).

The evolution of gaps in two and three dimensional discs can be heavily influenced by the presence of instabilities and the resultant formation of an unstable gap edge. Two relevant instabilities are the Rayleigh instability and the Rossby wave instability. The Rayleigh instability occurs due to violation of the Rayleigh stable condition, $d(R^{2}\Omega)/dR\geq 0$ (Chandrasekhar, 1961), which can be due to the deviation from Keplerian velocity of material within a deep gap formed by a massive planet. If this condition is violated the gap would become unstable and promote angular momentum transfer from the gap edge, lowering the surface density gradient. The effect of this on the aforementioned discrepancy has been investigated by Kanagawa et al. (2015). The Rossby wave instability occurs due to the formation of Rossby vortices by the shear in velocity of the material at the edge of the gap. These vortices only form when the velocity shear is significant, which corresponds to a gap edge with a steep surface density gradient (Lovelace et al., 1999; Li et al., 2000).

We focus our efforts on investigating the discrepancy between the results of one dimensional and two dimensional simulations, specifically the significant difference in gap depth. Here we explore a new approach to explaining this discrepancy, we remove the gap forming planet and instead apply a gap forming one dimensional torque density distribution radially across the two dimensional disc. This now becomes the mechanism for gap formation, and is the same distribution that forms the gap in one dimensional simulations. We proceed to investigate the resultant gap depth for a number of different forms of this one dimensional torque density distribution. We also extend our investigation to both lower viscosity and higher aspect ratio discs.

This paper is arranged as follows. In Section 2 we derive the relevant equations solved to simulate disc evolution. In Section 3 we discuss the code used and the numerical setup of our simulations. In Section 4 we present the results of our one and two dimensional simulations. In Section 5 we present the results of varying the disc parameters and compare these to our previous results. In Section 6 we discuss our findings in the context of prior work, while justifying assumptions made and postulating any impact they may have. Finally, in Section 7 we present our conclusions.

2 Basic Equations

2.1 Two dimensional

The continuity equation for the evolution of a protoplanetary discs surface density, $\Sigma$, is

[TABLE]

where $\mathbf{v}$ is the velocity field. We simulate the evolution of a protoplanetary discs surface density due to the presence of a planet by solving the two dimensional Navier-Stokes equation for the motion of the disc planet system,

[TABLE]

where $\mathbf{T}$ is the Newtonian viscous stress tensor, $P$ is the pressure and $\Phi$ is the gravitational potential of the planet and star system. We assume the disc is locally isothermal, with equation of state $P=c_{s}^{2}\Sigma$, where $c_{s}(R)$ is the sound speed at a radius $R$. We impose a constant aspect ratio $h=H/R=0.05$ by selecting a profile for the sound speed, where $H$ is the disc scale height. A cylindrical coordinate system is used, such that $\mathbf{v}=(v_{R},R\Omega)$ where $v_{R}(R)$ and $\Omega(R)$ are the radial and angular velocities respectively, at a given radius. Now splitting Equation 2 into its component forms we find,

[TABLE]

and

[TABLE]

where $\nu$ is the kinematic viscosity of the disc.

2.2 One dimensional

We reduce the two dimensional equations to standard one dimensional equations for disc evolution outlined in Pringle (1981), modified to account for the presence of a planet. We begin by averaging azimuthally Equation 4 so that only the radial dimension remains. This removes the potential term due to the planet, so we reintroduce an approximation to this as a torque density distribution, $\Lambda$,

[TABLE]

Assuming Keplerian rotation, $\Omega=(GM/R^{3})^{1/2}$, where $M$ is the stellar mass, Equation 5 can be written as,

[TABLE]

Using the azimuthally averaged form of Equation 1,

[TABLE]

we eliminate the $\Sigma v_{R}R$ term and find,

[TABLE]

We then solve this form of the continuity equation using standard finite difference techniques (Richtmyer & Morton, 1967).

2.3 Torque density distributions

We investigate three major forms of the torque density distribution, $\Lambda$. These are an impulse approximation (Equation 9), an improved model of the one dimensional torque density distribution from D’Angelo & Lubow (2010) (Equation 11) and an azimuthally averaged torque density distribution calculated from the planet itself in two dimensional simulations (Equation 13).

The impulse approximation for $\Lambda$ is in the manner of Lin & Papaloizou (1986),

[TABLE]

where $q=M_{p}/M$ is the ratio of planet mass to stellar mass and $R_{0}$ is the location of the planet. We take the dimensionless constant $f=0.23$, and $\Delta_{0}$ to be the greater of the disk scale height $H$ or $\left|R-R_{0}\right|$.

We split the disk into three sections, defined as inner disk, outer disk and inside gap. Discontinuities are removed by ensuring equality at the boundaries of these regions. To this end the $\textrm{sign}(R-R_{0})$ term is approximated to

[TABLE]

providing a continuous form for $\Lambda$. In Equation 10 and the definition of $\Delta_{0}$ we replace $H$ with the Hill Radius, $R_{\textrm{Hill}}$, for planet masses of which $R_{\textrm{Hill}}>H$.

The improved model of the one dimensional torque density distribution from D’Angelo & Lubow (2010) is of the form,

[TABLE]

where $x=(R-R_{0})/H$, $\mathcal{F}$ is a dimensionless function and $\beta=0.5$ and $\zeta=1.0$ are the surface density and temperature radial gradients respectively. $\mathcal{F}$ is determined by D’Angelo & Lubow (2010) to be,

[TABLE]

$p_{1}-p_{7}$ are constant for a given $\beta$, $\zeta$ and can be found in Table 1 of this paper and table 1 of D’Angelo & Lubow (2010). The function $\mathcal{F}$ is found by fitting the results of three dimensional simulations.

We then consider the torque density distributions presented in the upper panel figure 15 of D’Angelo & Lubow (2010). These distributions are determined from three dimensional discs in which the tidal torques are stronger than the viscous torques, such that a gap is formed. In this regime Equation 11 requires some modification to match the distributions shown in D’Angelo & Lubow (2010). From their figure 15 we can see a non-linear progression in the magnitude of the local maxima and minima of the torque density distributions, despite the linear mass progression. To account for this deviation from linearity we introduce a scaling factor such that our applied torque follows the trend present in their results. This factor is a cubic polynomial in planet mass and is applied in the region immediately before transition from using the disc aspect ratio to the hill radius of the planet, smoothing the transition between the low mass region (with no scaling factor) and the high mass region (in which a scaling factor of $1/2$ is applied as discussed in D’Angelo & Lubow (2010)).

The last torque density distribution we investigate is calculated from the torque applied to the disc from the planet in the two dimensional planet simulations. We average this torque azimuthally to form a one dimensional torque density distribution. The torque density distribution is given by

[TABLE]

where $\Sigma_{Av}(R)$ is the average surface density at a given radius and the gravitational potential of the planet $\Phi_{p}$ is given by

[TABLE]

where $s=0.6H$ is the smoothing length and $\phi_{0}$ is the azimuthal position of the planet. This provides the average torque the disc feels at a given radius due to the presence of the planet. A comparison of these two distributions can be seen in Figure 1 for a $q=10^{-3}$ planet.

3 Numerical Setup

3.1 One dimensional

One dimensional simulations run over a radial domain of $0.4\leq R/R_{0}\leq 2.5$ using $2100$ radial elements, providing a high resolution of $\Delta R_{\textrm{Cell}}=10^{-3}$ necessary for resolving the extremely steep gap edges, while maintaining computational efficiency. This is possibly an unnecessarily large domain for a one dimensional simulation however it allows for more accurate and direct comparisons with our two dimensional simulations discussed in Section 3.2. The initial surface density distribution was such that $\Sigma_{\textrm{int}}\propto R^{-1/2}$ and a constant viscosity $\nu=10^{-5}$ was used. This corresponds to a dimensionless viscosity $\alpha=4\cdot 10^{-3}$ at $R=R_{0}$, using the alpha prescription $\nu=\alpha c_{s}H$ of Shakura & Sunyaev (1973). We impose boundary conditions such that the surface density at the edge of the simulation is constant in time and hence feels no perturbation due to gap evolution. Simulations are run until equilibrium is attained.

3.2 Two dimensional

Two dimensional simulations were run using the FARGO3D code. This is a magnetohydrodynamic code developed with specific emphasis on simulating disc evolution for one dimensional to three dimensional models. The hydrodynamics equations are solved using a finite differences method. Importantly, FARGO3D uses a C to CUDA translator to allow simulations to run on GPU’s. GPU’s have limited memory, however decrease computational time significantly, making them a good choice for running moderate resolution simulations for extended periods of time. This makes FARGO3D a very good choice for investigating gap forming planets, as these are often computationally expensive and need to be run until equilibrium is attained. For further detail see Benítez-Llambay & Masset (2016).

The two dimensional simulation uses, wherever possible, the same parameters as the one dimensional simulation to draw accurate comparisons between them. The two dimensional simulation extends from $0.4\leq R/R_{0}\leq 2.5$ radially and $-\pi\leq\phi\leq\pi$ azimuthally with $256$ by $768$ cells respectively, with an initial surface density distribution of $\Sigma_{\textrm{int}}\propto R^{-1/2}$. This corresponds to a radial resolution of $\Delta R_{\textrm{Cell}}=8.203\cdot 10^{-3}$. The radial range was chosen to allow density waves excited by the embedded planet to propagate and potentially impact gap evolution, while the azimuthal range provides a view of the entire disk. This resolution was found to give consistent results with higher resolution runs, whilst still providing computational efficiency. Reflecting boundary conditions were used with parabolic damping wave killing zones in the regions $0.4\leq R/R_{0}\leq 0.505$ and $2.29\leq R/R_{0}\leq 2.5$, as per de Val-Borro et al. (2006). The combination of this closed boundary and wave killing zones acts as a good approximation to setting the boundaries at infinity, as excited density waves are significantly damped before reaching the boundaries and no wave reflections are observed. The planet is held on a fixed circular orbit ($e=0$) with no migration and no disc self gravity. We use this setup as the basis for simulating disc evolution in both planet and applied torque density distribution simulations.

We now consider the further modifications to this setup for our new approach to investigating the prior observed discrepancy in gap depth between one and two dimensional simulations. We apply a gap forming one dimensional axisymmetric torque density distribution radially across a two dimensional disc devoid of embedded planet. As a result the formation of the gap will be due to a torque applied to every cell, with magnitude dependent on the radial position of the cell. We attempt to closely mimic the torque due to the presence of the planet and measure resultant gap depth. To this end we use the parameters of the simulation as described here, with the addition that random noise of magnitude $10^{-5}$ was applied to the disc surface density, to remove axisymmetry. For the simulations containing planets the noise would be swiftly damped. We investigate the formation of a gap due to the one dimensional torque density distributions listed in Section 2.3.

To calculate gap depth from both planet and torque density distribution simulations two dimensional results were averaged azimuthally to form a one dimensional profile and the gap depth was taken to be the minimum of this. A small region of $10$ cells either side of the planet’s position was removed when taking the azimuthal average so the presence of the planet does not artificially decrease gap depth, similarly to Fung et al. (2014). See Figure 2 and Section 4.2. Gap depth was always calculated at equilibrium, defined after $1000$ orbits had elapsed for a $\nu=10^{-5}$ simulation, and always calculated from the average of the gap depth per orbit over $100$ orbits. This eliminates any short term variance in equilibrium gap depth due to gap turbulence. Figure 3 shows that beyond $1000$ orbits gap depth is, on average, approximately constant and so provides good reasoning for this definition of equilibrium.

4 Results

4.1 One dimensional results

The equilibrium surface density profiles for a number of sample runs of the one dimensional simulation can be seen in Figure 4. This shows the variation in shape of the gap profile with planet mass, namely the wider and deeper gaps are resultant of higher mass planets. The $q=10^{-3}$ planet case shows how sensitive the gap depth is to planet mass, forming a gap of depth $\Sigma/\Sigma_{0}\approx 10^{-24}$, considerably deeper than the other profiles shown in Figure 4. It was also found that higher mass planets approach equilibrium faster.

For confirmation of accurate definition of equilibrium, Equation 8 was solved analytically for $\partial\Sigma/\partial t=0$. These results can be seen in Figure 5, showing good agreement over a large range of planet masses. This prompts the conclusion that the gap depths from the one dimensional simulations are at an acceptable degree of accuracy.

4.2 Two dimensional planet results

Gap depths in two dimensions show the same expected mass-depth trend as they do in one dimension. The most interesting result of the two dimensional simulations arises from the comparison between one and two dimensional results. A direct comparison of one and two dimensional gap profiles for the same planet and disc can be seen in Figure 6. The gap profiles are distinctly different as the one dimensional case is thinner and significantly deeper, with $\Sigma/\Sigma_{0}\approx 10^{-24}$. This discrepancy is also illustrated in Figure 7, which shows us it is present for any mass planet that can open a significant gap in the disc.

We check the consistency of our two dimensional gap depth results against the gap depth scaling law of Fung et al. (2014), valid for $10^{-4}\leq q\leq 5\cdot 10^{-3}$,

[TABLE]

where $\Sigma_{\textrm{Gap}}$ is the gap depth. This scaling was calculated from two dimensional disc-planet models. From Figure 7 we can see that, for masses of roughly $q\geq 3\cdot 10^{-4}$, the results of our two dimensional simulations follow closely the scaling law.

4.3 Two dimensional impulse approximation

We now apply the impulse approximation, described by Equation 9, directly from the one dimensional simulations radially across a two dimensional disc devoid of any planet. We do this for a range of masses comparable to both prior simulations. The results can be seen in Figure 8, which compares the results of this simulation to the prior two simulations. A sample simulated output can be seen in Figures 9 and 10. Figure 9 is a parallel to Figure 2, however using the one dimensional torque density distribution rather than the gap forming planet. Particularly noticeable is the presence of a turbulent gap edge, absent in Figure 2 and unaccounted for in Section 4.1. The ramifications of this can be seen in Figure 10, which shows the significantly reduced depth of the two dimensional impulse approximation gap compared to the one dimensional simulation.

Figure 8 shows the significant difference in gap depth resulting from using the impulse approximation in two dimensions. Applying this approximation to the torque of the planet accounts for the majority of the discrepancy, with only roughly an order of magnitude in difference between the two dimensional planet simulations and the two dimensional impulse approximation simulations.

It is not expected that the applied torque density distribution would return a gap depth of equal magnitude, as it is only an approximation to the planet, however it is important to understand why the formed gap is significantly shallower in two dimensions than in one. The disc material is forced by the torque density distribution to attempt to form a gap of similar shape to those in Figure 4, namely deep, thin and axisymmetric, characteristic of one dimensional simulations. This results in a steep surface density gradient at the gap edge, which causes Rossby vortices to form and the gap edge to become unstable (Lovelace et al., 1999; Li et al., 2000). This results in the excitation of density waves which redistribute angular momentum. This lessens the surface density gradient at the gap edge forming a more stable gap and reducing the gap depth, providing the shallower gap we see. The Rossby vortices outside the gap edge will often merge and weaken before equilibrium is reached. They are however, present throughout the discs lifetime after excitation to ensure the disc does not return to an axisymmetric state. The profile of this shallower gap can be seen in Figure 10, not dissimilar to that of the two dimensional planet case.

Despite the angular momentum transfer creating a shallower gap, a quick comparison of Figures 2 and 9 show that the one dimensional torque density distribution forms a significantly more turbulent gap. The presence of this turbulence gives rise to short term variability in equilibrium gap depth, which makes averaging gap depth over $100$ orbits an important factor in determining true equilibrium gap depth.

While this rectifies a large portion of the discrepancy between these simulations, the gap depths are still close to an order of magnitude lower than in the two dimensional planet simulations. We proceed to make further amendments to the applied torque density distribution in order to investigate the remaining difference.

4.4 Improved one dimensional torque density distribution

Using the results of D’Angelo & Lubow (2010) we proceed to modify the torque density distribution to the form given in Equation 11 and measure the impact on gap depth. The variation of the minimum surface density with mass under the influence of this torque density distribution can be seen in Figure 11. Again we can see better agreement between the model and the two dimensional planet model to a larger mass, but still a large discrepancy in the higher mass regime.

Returning to the torque density distributions presented by D’Angelo & Lubow (2010) in their figure 15, we now consider the lower panel. The torque density distributions here correspond to high mass planets, the regime in which our simulations are inconsistent with the two dimensional planet results. Comparing the shapes of these distributions we see interesting phenomena present in the larger mass distributions that are otherwise absent in the lower mass distributions. As the mass of the planet increases, the peak/trough of the torque density distribution becomes narrower, it approaches zero at a greater distance from the planet, and remains approximately zero for longer. This is a significant deviation from the torque density distribution currently being applied at these high masses. Therefore, we crudely reproduce this phenomena in our model by setting our torque density distribution to zero for $\lvert R-R_{0}\rvert/R_{\textrm{Hill}}<0.95$ with a sharp transition region from $0.95\leq\lvert R-R_{0}\rvert/R_{\textrm{Hill}}\leq 1.05$. The difference between these two torque models is clearly visible in Figure 12. This new torque model is applied to the high mass simulations, where the previous model begins to diverge from the two dimensional planet simulation. The results of this can be seen in Figure 13. This shows greater consistency with the two dimensional planet simulation to even larger masses, compared to previous simulations, however we still observe a trend away from the two dimensional planet simulation as the mass increases. The discontinuous region present at roughly $q=3\cdot 10^{-4}$ is a result of the regime change from the low mass to the high mass torque density distribution shown in Figure 12.

A conclusion we can draw from this is that the shape of the applied torque density distribution has a significant impact on the equilibrium depth of the resultant gap. Additionally it appears that the shape of the distribution at close proximity to the planet plays an important role in determining gap depth. Further investigation shows that while the shape of the distribution in this region does impact the gap depth, it is currently unknown how important and what other factors are contributing to the gap depth. This is discussed in greater detail in Section 6. Nevertheless we can clearly see the form of the torque density distribution is important, and worth further investigation.

4.5 One dimensional form of the planet’s torque density distribution

Concluding that the torque density distribution is an important factor, we now apply the azimuthally averaged shape of this distribution calculated from our two dimensional planet simulations, as presented in Equation 13. For a sample case of a $q=10^{-3}$ mass planet, this torque distribution can be seen in Figure 1. This was calculated after equilibrium had been attained, resulting in a constant torque distribution orbit to orbit. This distribution is comparable to the previous distributions used, although shows interesting detail at very close proximity to the planet. We now apply this torque to a two dimensional disc in the same manner as our previous distributions. In Figure 14 we compare the azimuthally averaged gap profile of this with the corresponding profiles from prior simulations. In Figure 15 we compare the calculated one dimensionally averaged torque distributions for a range of planet masses with our prior simulations.

This clearly shows that using the exact torque distribution is an improvement over the results of the one dimensional model, but does not match the previous accuracy of the modified improved torque density distribution. As a result, we have to conclude here that the torque density distribution is significant to the depth of the gap and using the correct form in two dimensions can greatly reduce the discrepancy between one and two dimensional models. However, there appears to be a limit to which the torque density distribution can improve over the one dimensional model, as the torque density distribution is only an approximation to the planet. We now investigate the behaviour of the gap under applied torque density distributions in discs with different parameters.

5 Results of varying Disc Parameters

5.1 Low viscosity disc

We now proceed to investigate the behaviour of a disc with $\nu=10^{-6}$, corresponding to an $\alpha=4\cdot 10^{-4}$ at $R=R_{0}$. We initialise the disc as discussed in Section 3.2 with the only difference being the lower viscosity and follow a procedure largely the same as that described in Section 4, however we do not investigate the improved torque density distribution that we discuss in Section 4.4 as unfortunately this disc is outside the range of discs studied in D’Angelo & Lubow (2010) and therefore they do not have fit parameters pertaining to it. We investigate this disc for one planet mass, which we choose to return similar gap depths to our previous $q=10^{-3}$ planet simulations. This corresponds to $q=2.77\cdot 10^{-4}$ in this case. In two dimensions, simulations now run to $10000$ orbits before reaching equilibrium due to this low viscosity.

The cases investigated were the one dimensional, two dimensional planet, two dimensional impulse approximation and two dimensional azimuthally averaged calculated torque distribution. The azimuthally averaged gap profiles for these can be seen in Figure 16. This shows a similar result to those presented in Figure 14, the one dimensional gap is significantly deeper than its two dimensional counterparts, and the application of one dimensional torques to the two dimensional disc make up a sizeable amount of the discrepancy, but still remain roughly an order of magnitude lower than the two dimensional planet case. A deviation from the results of Figure 14 is that here the azimuthally averaged calculated torque distribution provides better agreement than the impulse approximation, whereas in the higher viscosity disc this is not the case. Despite this we see good agreement between the low viscosity disc and our standard disc model, in that the gaps formed by one dimensional torque density distributions in two dimensions are significantly shallower than in one dimension, however some discrepancy still remains between these gap depths and the two dimensional planet case.

5.2 High aspect ratio disc

We now proceed to investigate the behaviour of a disc with $h=H/R=0.063$ and $\nu=10^{-5}$, corresponding to an $\alpha=2.52\cdot 10^{-3}$ at $R=R_{0}$. We initialise the disc as discussed in Section 3.2 with the only difference being the higher aspect ratio. We investigate this disc for the same cases as the low viscosity disc discussed in Section 5.1. Again we investigate for one planet mass, $q=2\cdot 10^{-3}$, selected with the same reasoning as in Section 5.1.

The azimuthally averaged gap profiles for the cases investigated can be seen in Figure 17. Again we see a similar result to that presented in Figure 14, with the one dimensional gap still significantly deeper than its two dimensional counterparts, and the application of one dimensional torque density distributions to the two dimensional disc make up a sizeable amount of the discrepancy, but still remain roughly an order of magnitude lower than the two dimensional planet case. However in contrast to Section 5.1, here we see the impulse approximation providing a better agreement to the two dimensional planet case than the azimuthally averaged calculated torque distribution. This is the exact result found in Figure 14 and hence we see the two dimensional gaps formed by one dimensional torque density distributions are significantly closer in depth to the two dimensional planet gap than their one dimensional counterparts.

6 Discussion of Results

In this paper we find that the discrepancy between equilibrium gap depths from one and two dimensional simulations can be explained by the absence of the Rossby wave instability in one dimension. We determine this using one dimensional torque density distributions as the gap forming mechanism in two dimensional simulations. We find good agreement with prior work regarding our one and two dimensional planet simulations (Lin & Papaloizou, 1986; Crida et al., 2006; Kanagawa et al., 2015), specifically our two dimensional planet gap depths closely follow the established gap depth scaling law (Fung et al., 2014).

Our results show the excitement of the Rossby wave instability in two dimensional simulations when a one dimensional torque density distribution is used as a gap forming mechanism. Prior work (de Val-Borro et al., 2006; de Val-Borro et al., 2007) shows that the Rossby wave instability is present in two dimensional planet simulations for low viscosity discs or high mass planets. It is also shown that these instabilities die off before equilibrium is reached (Fu et al., 2014; Hammer et al., 2016), suggesting that the gap is maintained at marginal stability set by the point at which the excitement of the Rossby wave instability occurs.

We also find that when the planet is not massive enough to open a gap in the disc there is good agreement between one and two dimensional simulations as shown in Figure 15. This is because there is no Rossby wave instability present in the two dimensional simulations that do not open a gap. This result could be important for modelling the behaviour of low mass planets in one and two dimensions.

For population synthesis models, e.g. Mordasini et al. (2009), especially those regarding the formation of giant planets (Ida & Lin, 2008), the depth of the gap formed and the feedback of the planet onto the disc structure is very important. If the depth of the gap is set by the instabilities observed, such as the Rayleigh instability investigated in Kanagawa et al. (2015) and the Rossby wave instability presented here, this could lead to an improved recipe for gap formation, which would be very important in future population synthesis models. The fact that our results hold for differing disc parameters (Figure 16, Figure 17) reinforces the relevance of our results to such a recipe for gap formation.

We can find some similarity to our investigation of the gap depth discrepancy in Kanagawa et al. (2015). They examine gap depth in viscous, one dimensional discs accounting for deviation from Keplerian disc rotation and include the Rayleigh stable condition. They find a shallowing of the gap due to angular momentum transfer from the gap edge, arising as a result of the deviation from Keplerian disc rotation, or violation of the Rayleigh stable condition. We see a similar situation in our two dimensional simulations, however the shallowing of the gap is a result of the excitement of the Rossby wave instability, which is unaccounted for in one dimension. Kanagawa et al. (2015) also show that a significantly large wave propagation length before angular momentum deposition occurs can make the one dimensional gap shallower and wider. They find a wave propagation length that gives a comparable depth to the two dimensional model for high mass planets, however the increase in mass corresponds to a long propagation length and a wider gap. The correct width of the two dimensional gap is unknown, with prior work both agreeing (Varnière et al., 2004) and disagreeing (Duffell & MacFadyen, 2013) with their result.

Closely related to this is the work of Crida et al. (2006). They find that the equilibrium gap profile of two dimensional simulations can be accurately reproduced semi-analytically by balancing the viscous, gravitational and pressure torques on the disc. The pressure torque here arises from the fraction of the gravitational torque that is carried away by density waves instead of being locally damped and is important to increasing the width of the analytical gaps to match their numerical counterparts. An ansatz for the pressure torque is determined from reference numerical simulations, hence the result is not fully analytic.

The results of both Crida et al. (2006) and Kanagawa et al. (2015) imply that wave propagation can account for the discrepancy between one and two dimensional simulations, however both approaches introduce parameters that must be fitted in order to determine the wave propagation length, or the magnitude of the pressure torque. From our results we find that two dimensional gaps formed by one dimensional torque density distributions have depths significantly closer to those of two dimensional planet gaps than their one dimensional counterparts. This is a consequence of the excitement of the Rossby wave instability. Considering this together with the result that high mass planets appear to have gaps maintained at marginal stability (Fu et al., 2014; Hammer et al., 2016), we suggest that the pressure torque required to form the characteristic two dimensional gaps is such that the gap edge is maintained at marginal stability.

Our results from exploring different discs provide good agreement between themselves. We find in each case that using a one dimensional torque density distribution as an approximation to the torque on the disc due to the planet provides significantly better agreement with the planet case in two dimensions than in one. We also find that in the $\nu=10^{-6}$ case, the averaged one dimensional torque array proves to be a better approximation than the impulse approximation, however in both the $\nu=10^{-5}$ discs we find the opposite to be true.

The results of Section 4.4 give the impression that the treatment of the torque at close proximity to the location of the planet has a major impact on the depth of the resultant gap. This was investigated using the azimuthally averaged calculated torque density distribution, largely due to the detailed shape of the distribution at close proximity to the planet (see Figure 1). It was found that setting this detailed region to zero, similarly to the modification made in the second part of Section 4.4, had very minimal impact on the gap depth. Widening this area, such that the transition to this region is steeper, moves the gap depth in the correct direction, but not significantly. This was unexpected, as the shape of the azimuthally averaged calculated torque density distribution was crudely forced until it was very close to the modified improved torque density distribution from Section 4.4, yet it still does not return a similar gap depth. In order to determine where the difference between these two torque density distributions arises, we split these torque density distributions into two regions, the inner disc ranging from the inner edge of the simulation to the beginning of the zero torque region and the outer disc ranging from the end of the zero torque region to the outer edge of the simulation. We then create two more torque density distributions, consisting of the inner part of one torque density distribution and the outer part of another. In each of these distributions one region is modelled by the improved torque density distribution while the other is the azimuthally averaged calculated torque density distribution. We find that the Inner Disc Modelled simulation has minimal effect in changing the gap depth, whereas the Outer Disc Modelled returns a gap depth extremely close to the original, modified improved torque density distribution. At this point we can conclude that it is not only the torque close to the planet that has a large impact on the gap depth, there are a variety of other factors. Currently, however, we are unable to deduce exactly what all of these factors are or specifically determine how these impact the depth of the formed gap.

One such factor that could be important to consider is the presence of material on horseshoe orbits and its impact on gap depth. Horseshoe orbits are present in simulations containing planets, as the planet perturbs the orbit of material at similar orbital radius to itself as usual. However, this is absent in the case of a torque density distribution forming the gap. The interaction between the planet and material on a horseshoe orbit gives rise to a corotation torque, due to the vortensity (Ward, 1991; Masset, 2001) and entropy (Baruteau & Masset, 2008; Paardekooper & Papaloizou, 2008; Paardekooper et al., 2010) gradients across this region. The presence of this torque in the planet case may provide an explanation as to why there is consistently a discrepancy between gap depths from planet formed gaps and torque density distribution formed gaps.

Recent studies into the excitement of the Rossby wave instability have shown that using more realistic (longer) planetary growth times induce shorter lived and lower amplitude vortices (Hammer et al., 2016). We do not increase the intensity of our applied torque distributions over any number of orbits, instead we apply our torque density distributions from the beginning of each simulation. We do not believe the results of Hammer et al. (2016) will effect our findings, as in the torque density distribution simulations the instabilities are present throughout the whole simulation. If the instabilities did die off the disc would return to an axisymmetric state and the gap would again try to become deep and narrow, which would in turn re-excite the Rossby wave instability. Hence the Rossby wave instability will be present throughout the discs lifetime.

We make a number of simplifying assumptions during this investigation. Planets were held on a constant, circular orbit ($e=0$) and we do not include planetary migration. We use a locally isothermal equation of state, as radiative transfer and cooling is difficult to accurately model and is very computationally expensive. Hence, disc thickness $H$ is linear in $R$ and gives us a constant aspect ratio $h=H/R$, for $H\ll R$. We use a constant kinematic viscosity across the disc and neglect disk self gravity. The effect of disc self gravity on vortex instabilities has been investigated (Lin & Papaloizou, 2011) and it was found that larger disc masses form more vortices. With sufficiently large self gravity, however, vortex formation was suppressed. It was also found that self gravity acts to delay vortex merging. The fact that vortex formation is still present with self gravity gives confidence that our results would still be applicable in that regime.

We also neglect magnetohydrodynamic (MHD) effects on the disc, which would otherwise cause the disc to become turbulent. Gap formation in MHD turbulent discs has been investigated (Papaloizou et al., 2004) and have been found to return wider and shallower gaps than their purely hydrodynamic counterparts (Winters et al., 2003). The impact of turbulence would be difficult to account for in the applied torque distribution cases and would be computationally expensive. Papaloizou et al. (2004) show that in the vicinity of the planet, local shearing box simulations are a good approximation to global disc models and are less computationally expensive, which would be useful if accounting for MHD was pursued. These factors all constitute a very idealised case of gap formation and removing some of these simplifying assumptions could also impact our results.

We compare results of two dimensional and one dimensional simulations for one dimensional gap forming mechanisms such as the impulse approximation. Currently it is unknown as to how a gap formed by a one dimensional torque density distribution such as the impulse approximation would behave in three dimensions and how this would compare to both its one and two dimensional counterparts. Recent studies into three dimensional gap formation has found three dimensional gaps have similar depth and shape to their two dimensional counterparts (Fung & Chiang, 2016). With this in mind, and considering the gap depth is set by the excitement of the Rossby wave instability, a two dimensional instability, we can predict with some confidence that there would be minimal change in our findings if it was taken to three dimensions. Of course this area is far from fully explored and a number of assumptions are still involved, so further investigation would be required before any conclusions could be drawn with certainty.

7 Conclusions

We have investigated gap formation in viscous protoplanetary discs in both one and two dimensions, with intent to measure the depth in surface density of the gap formed. We encountered a well known discrepancy (Crida et al., 2006; Kanagawa et al., 2015) in depth between one and two dimensional gaps. We proceeded to investigate this discrepancy using a new method, the application of a one dimensional gap forming torque density distributions to a two dimensional disc devoid of a planet. We studied this for a number of forms of the torque density distribution;

• •

An impulse approximation (Equation 9), as of Lin & Papaloizou (1986).

• •

An improved one dimensional torque density distribution using results from three dimensional models (Equation 11), as of D’Angelo & Lubow (2010).

• •

A modified version of the above torque density distribution from D’Angelo & Lubow (2010), in which $\Lambda=0$ at close proximity to the planet.

• •

An azimuthally averaged one dimensional torque distribution calculated due to the presence of the planet in two dimensional simulations (Equation 13).

We found that applying a one dimensional torque density distribution across the disc as a gap forming mechanism results in the formation of a gap significantly shallower than the one dimensional gap, however still roughly an order of magnitude deeper than the equivalent two dimensional planet simulation. This occurs as a deep thin gap is attempted to be formed, much like in the one dimensional case, however as the two dimensional gap edge becomes too steep Rossby vortices begin to form and the gap edge becomes unstable. This excites density waves that redistribute angular momentum across the disc, which acts to reduce the steepness of the gap edge and as a result the gap becomes shallower. de Val-Borro et al. (2006); de Val-Borro et al. (2007); Fu et al. (2014); Hammer et al. (2016) show that the Rossby wave instability is present in low viscosity discs or high mass planet simulations, however the instability dies off before equilibrium is reached. This suggests that the gap is maintained at marginal stability, set by the point at which the excitement of the Rossby wave instability occurs. Our results show a similar gap depth between two dimensional gaps formed by torque density distributions, in which the Rossby wave instability is present, and two dimensional planet gaps, in which no Rossby wave instability is present. This can be understood if the planet gap is maintained at marginal stability, even when there is no obvious Rossby wave instability present.

We observe very similar behaviour across all applied torque density distributions. We see that while the choice of torque density distribution does effect the resultant gap depth, the change is often minimal and cannot make up for the order of magnitude discrepancy remaining. The only torque density distribution that makes any significant progression towards the two dimensional planet model is the modified version of the improved torque density distribution from D’Angelo & Lubow (2010), in which the torque at close proximity to the planet is set to zero. This would propose that the torque close to the planet is the dominant factor in the resultant gap depth, however while it does impact the gap depth, this does not appear to be the case. Therefore, while we can say the gap depth is sensitive to the torque density distribution, the exact dependencies are still unknown. We are considering here only the gap depths in the higher mass regime, where the discrepancy still exists. For lower mass planets, roughly $q\leq 3\cdot 10^{-4}$, the improved torque density distribution and the azimuthally averaged calculated torque density distribution are significant improvements over the impulse approximation, and almost entirely eliminate the discrepancy. Hence in this mass range, these are very good approximations to the effect of a two dimensional planet.

When we extend this study into discs with different parameters we find very good agreement with our prior results. While we only study two additional discs, changing the aspect ratio and viscosity, and study each case for a single planet mass, they return the expected results. With this we can confidently say our results are not constrained solely to the disc we have been investigating.

Acknowledgements

Sijme-Jan Paardekooper is supported by a Royal Society University Research Fellowship.