- Split View
-
Views
-
Cite
Cite
Min-Kai Lin, Testing large-scale vortex formation against viscous layers in three-dimensional discs, Monthly Notices of the Royal Astronomical Society, Volume 437, Issue 1, 01 January 2014, Pages 575–587, https://doi.org/10.1093/mnras/stt1909
- Share Icon Share
Abstract
Vortex formation through the Rossby wave instability (RWI) in protoplanetary discs has been invoked to play a role in planet formation theory and suggested to explain the observation of large dust asymmetries in several transitional discs. However, whether or not vortex formation operates in layered accretion discs, i.e. models of protoplanetary discs including dead zones near the disc mid-plane – regions that are magnetically inactive and the effective viscosity greatly reduced – has not been verified. As a first step towards testing the robustness of vortex formation in layered discs, we present non-linear hydrodynamical simulations of global 3D protoplanetary discs with an imposed kinematic viscosity that increases away from the disc mid-plane. Two sets of numerical experiments are performed: (i) non-axisymmetric instability of artificial radial density bumps in viscous discs and (ii) vortex formation at planetary gap edges in layered discs. Experiment (i) shows that the linear instability is largely unaffected by viscosity and remains dynamical. The disc–planet simulations also show the initial development of vortices at gap edges, but in layered discs, the vortices are transient structures which disappear well into the non-linear regime. We suggest that the long-term survival of columnar vortices, such as those formed via the RWI, requires low effective viscosity throughout the vertical extent of the disc, so such vortices do not persist in layered discs.
INTRODUCTION
Recent observations have revealed a class of transition discs – circumstellar discs which are dust poor in its inner regions – with non-axisymmetric dust distributions in its outer regions (Brown et al. 2009; Mayama et al. 2012; Isella et al. 2013; van der Marel et al. 2013). One interpretation of such a non-axisymmetric structure is the presence of a large-scale disc vortex, which is known to act as a dust trap (Barge & Sommeria 1995; Bracco et al. 1999; Chavanis 2000; Inaba & Barge 2006; Ataiee et al. 2013; Birnstiel, Dullemond & Pinilla 2013; Lyra & Lin 2013). Because of its occurrence adjacent to the inner dust hole, i.e. a cavity edge, it has been suggested that such a vortex is a result of the Rossby wave instability (RWI): a hydrodynamical instability that can develop in radially structured discs.
Modern work on the RWI began with two-dimensional (2D) linear stability analysis (Lovelace et al. 1999; Li et al. 2000). These studies show that a disc with radially localized structure, such as a surface density enhancement of ≳ 10 per cent over a radial length-scale of the order of the local disc scaleheight, is unstable to non-axisymmetric perturbations, which grow on dynamical (orbital) time-scales. Early 2D non-linear hydrodynamic simulations showed that the RWI leads to multivortex formation, followed by vortex merging into a single large vortex in quasi-steady state (Li et al. 2001; Inaba & Barge 2006).
While these studies consider disc models with artificial radial structure, it has recently been established that a natural site for the RWI is the edge of gaps induced by disc–planet interaction (Koller, Li & Lin 2003; Li et al. 2005; de Val-Borro et al. 2007; Li et al. 2009; Lyra et al. 2009; Lin & Papaloizou 2010, 2011). Indeed, this has been the proposed explanation for the lopsided dust distribution observed in the Oph IRS 48 transition disc system (van der Marel et al. 2013).
An important extension to the aforementioned studies is the generalization of the RWI to three-dimensional (3D) discs. Both non-linear 3D hydrodynamic simulations (Meheut et al. 2010, 2012b; Lin 2012a; Lyra & Mac Low 2012; Richard, Barge & Le Dizes 2013) and 3D linear stability calculations (Umurhan 2010; Lin 2012b, 2013b; Meheut, Yu & Lai 2012a) have been carried out. These studies reveal that the RWI is a 2D instability, in that there is negligible difference between growth rates obtained from 2D and 3D linear calculations. The associated density and horizontal velocity perturbations have weak vertical dependence, and vertical velocities are small. In non-linear hydrodynamic simulations, the vortices are columnar and extend throughout the vertical extent of the disc (Richard et al. 2013).
The RWI therefore appears to be a global instability in the direction perpendicular to the disc mid-plane: the vortical perturbation involves the entire fluid column. Thus, conditions away from the disc mid-plane may have important effects on vortex formation via the RWI. For example, Lin (2013a) only found linear instability for certain upper disc boundary conditions. This issue is relevant to protoplanetary disc models including ‘dead zones’.
It is believed that mass accretion in protoplanetary discs is driven by magnetohydrodynamic (MHD) turbulence as a result of the magnetorotational instability (MRI; Balbus & Hawley 1991, 1998). However, it is not clear if the MRI operates throughout the vertical extent of the disc, because the mid-plane of protoplanetary disc is dense and cold (Armitage 2011). As a result, Gammie (1996) proposed the layered disc model: accretion due to MHD turbulence is small near the mid-plane (the dead zone), while MHD turbulence-driven accretion operates near the disc surface (the active zone). The layered accretion disc model has been subject to numerous studies (e.g. Fleming & Stone 2003; Terquem 2008; Oishi & Mac Low 2009; Dzyurkevich et al. 2010; Kretke & Lin 2010; Okuzumi & Hirose 2011; Flaig et al. 2012; Landry et al. 2013). If MRI-driven accretion can be modelled through an effective viscosity (Balbus & Papaloizou 1999), this corresponds to a low-viscosity mid-plane and high-viscosity atmosphere. It is therefore valid to ask how such a vertical disc structure would affect large-scale vortex formation via the RWI.
This problem is partly motivated by viscous disc–planet simulations which show that gap-edge vortex formation only occurs when the viscosity is sufficiently small (de Val-Borro et al. 2006, 2007; Edgar & Quillen 2008). What happens if the effective viscosity near the mid-plane is sufficiently low for the development of Rossby vortices, but is too high away from the mid-plane?
In this work, we examine vortex formation through the RWI in layered discs. As a first study, we take an experimental approach through customized numerical hydrodynamic simulations. We simulate global 3D protoplanetary discs with an imposed kinematic viscosity that varies with height above the disc mid-plane.
The central question is whether or not applying a viscosity only in the upper layers of the disc damps the RWI and subsequent vortex formation. The purpose of this paper is to demonstrate, through selected simulations, the potential importance of layered disc structures on vortex formation. We defer a detailed parameter survey to a future study.
This paper is organized as follows. The accretion disc model is set up in Section 2 and the numerical simulation method described in Section 3. Results are presented in Section 4 for viscous discs initialized with a density bump. These simulations employ a special setup such that the density bump is not subject to axisymmetric viscous diffusion. This allows one to focus on the effect of layered viscosity on the linear non-axisymmetric instability. Section 5 revisits vortex formation at planetary gap edges, but in 3D layered discs, where it will be seen that vortex formation can be suppressed by viscous layers. Section 6 concludes this work with a discussion of important caveats of the present disc models.
DISC MODEL AND GOVERNING EQUATIONS
Disc model and initial conditions
The numerical disc model occupies r ∈ [rin, rout], |$\theta \in [\theta _\mathrm{min}, \pi /2]$| and ϕ ∈ [0, 2π] in spherical coordinates. Only the upper disc is simulated explicitly (ψ > 0), by assuming symmetry across the mid-plane. The maximum angular height is |$\psi _\mathrm{max}\equiv \pi /2 - \theta _\mathrm{min}$|. The extent of the vertical domain is parametrized by nh ≡ tan ψmax/h, i.e. the number of scaleheights at the reference radius.
The initial cylindrical radial velocity vRi and the viscosity profile ν depends on the numerical experiment and will be described along with simulation results. Note that vRi and ν are not independent if one additionally requires a steady state (see Section 4).
Damping
Planet potential
NUMERICAL EXPERIMENTS
The necessary condition for the RWI – a PV extremum (Li et al. 2000) – is either set as an initial condition via a density bump or obtained from a smooth disc by evolving it under disc–planet interaction. The setup of each experiment is detailed in subsequent sections.
We adopt units such that G = M* = 1 and the reference radius r0 = 1. We set σ = 0.5 for the initial surface density profile and apply frictional damping within the shells r < rd, in = 1.25rin and r > rd, out = 0.84rout.
The fluid equations are evolved using the pluto code (Mignone et al. 2007) with the FARGO algorithm enabled (Masset 2000; Mignone et al. 2012). We employ a static spherical grid with (Nr, Nθ, Nϕ) zones uniformly spaced in all directions. For the present simulations, the code was configured with piecewise linear reconstruction, a Roe solver and second-order Runge–Kutta time integration.
Boundary conditions are imposed through ghost zones. Let the flow velocity parallel and normal to a boundary be v∥ and v⊥, respectively. Two types of numerical conditions are considered for the (r, θ) boundaries: (a) reflective: ρ and v∥ are symmetric with respect to the boundary, while v⊥ is antisymmetric and (b) unperturbed: ghost zones retain their initial values. The boundary conditions adopted for all simulations are unperturbed in r, reflective in θ and periodic in ϕ.
Diagnostics
We list several quantities calculated from simulation data for use in results visualization and analysis.
Density perturbations
Vortical structures
Potential vorticity
Perturbed kinetic energy density
We define the perturbed kinetic energy as |$W\equiv \rho |\boldsymbol {v}|^2/[\rho _{\rm i}|\boldsymbol {v}(t =0)|^2] - 1$|, and its Fourier transform |$W_m\equiv \int _0^{2\phi } W\exp {(-\mathrm{i}m\phi )}{\rm d}\phi$|. We will examine |Wm(r, θ)| averaged over subportions of the (r, θ) plane.
NON-AXISYMMETRIC INSTABILITY OF ARTIFICIAL RADIAL DENSITY BUMPS IN LAYERED DISCS
We first consider strictly isothermal discs (q = 0) initialized with a density bump (A > 1). Our aim here is to examine the effect of (layered) viscosity on the RWI through the linear perturbation. In general, a density bump in a viscous disc will undergo viscous spreading (Lynden-Bell & Pringle 1974), but we can circumvent this by choosing the viscosity profile ν and initial cylindrical radial velocity vRi appropriately. Although artificial, this setup avoids the simultaneous evolution of the density bump subject to axisymmetric viscous spreading and growth of non-axisymmetric disturbances; only the latter of which is our focus.
Viscous equilibria for a radially structured disc
In choosing ρi and Ωi, we neglected radial velocities and viscous forces in the steady-state vertical and cylindrical radial momentum equations (equations 7 and 11, respectively). This is a standard practice for accretion disc modelling (e.g. Takeuchi & Lin 2002).
However, vR and ν cannot be ignored in the azimuthal momentum equation. Indeed, if a steady state is desired, then the conservation of angular momentum in a viscous disc implies special relations between the viscosity, cylindrical radial velocity and density field.
Initial cylindrical radial velocity
Viscosity profile for a steady state
If the initial conditions corresponds to a steady state, then RρivRi can only be a function of z. With vRi chosen by equation (20), this implies |$R\rho _{\rm i}\nu \Omega _{\rm i}^\prime /\Omega _{\rm i}$| is only a function of z. We are therefore free to choose the vertical dependence of viscosity.
Equation (21) implies that at the fixed cylindrical radius R = r0, the dimensionless viscosity increases from |$\hat{\nu } = \hat{\nu }_0$| at the mid-plane to |$\hat{\nu } = A_\nu \hat{\nu }_0$| for z > ζνH0. An example of such a layered viscosity profile is depicted in Fig. 1.
Simulations
We consider discs with radial extent [rin, rout] = [0.4, 2.0]r0, vertical extent nh = 2 scaleheights and aspect ratio h = 0.1 at R = r0. We use (Nr, Nθ, Nϕ) = (256, 64, 512) grid points. The resolution at the reference radius is then 16, 32 and 8 cells per scaleheight in (r, θ, ϕ) directions, respectively. The planet potential is disabled for these runs (Mp ≡ 0). We apply a damping rate |$\hat{\gamma }=1$| with the reference velocity |$\boldsymbol {v}_\mathrm{ref}=\boldsymbol {v}(t=0)$|.
The bump parameters are set to A = 1.25 and ΔR = 0.05r0 for all runs in this section. The corresponding PV profile is shown in Fig. 2. The spherical radial velocity is subject to random perturbations of magnitude 10− 4cs a few time-steps after initialization.
Linear growth rates and frequencies
The present setup allows us to define a linear instability in the usual way: exponential growth of perturbations measured with respect to an axisymmetric steady equilibrium. A proper linear stability analysis, including the full viscous stress tensor, is beyond the scope of this paper, but we can nevertheless extract linear mode frequencies from the non-linear simulations.
In a linear stability problem, σm is a constant eigenvalue. However, when extracted numerically from a non-linear simulation, we will generally obtain σm = σm(r, θ, t). Thus, we compute 〈σm〉r = mωm + iqm, where ωm is the mode frequency and qm is the growth rate. We normalize the linear frequencies by Ω0 ≡ Ωi(r0) ≃ Ωk(r0).
Results
Table 1 summarizes the simulations presented in this section. For reference, we simulate an effectively inviscid disc, case B0, with the viscosity parameters |$\hat{\nu }_0=10^{-9}$| and Aν = 1. Thus, viscosity is independent of z at R = r0. Inviscid setups similar to case B0 have previously been simulated both in the linear and non-linear regimes (Meheut et al. 2012a; Lin 2013b).
. | t = 10P0 (linear phase) . | t = 100P0 . | |||||||
---|---|---|---|---|---|---|---|---|---|
Case . | |$\log {\hat{\nu }_0}$| . | Aν . | ζν . | m . | ωm/Ω0 . | qm/Ω0 . | m . | 102am . | min[Ro(z = 0)] . |
B0 | −9 | 1 | n/a | 4 | 0.985 | 0.199 | 1 | 8.5 | −0.15 |
V0 | −6 | 1 | n/a | 4 | 0.985 | 0.199 | 1 | 6.8 | −0.11 |
V1 | −6 | 100 | 1.5 | 4 | 0.986 | 0.191 | 1 | 7.8 | −0.19 |
V2 | −6 | 100 | 1.0 | 4 | 0.986 | 0.182 | 1 | 4.9 | −0.21 |
V3 | −4 | 1 | n/a | 4 | 0.986 | 0.131 | 3 | 3.7 | −0.25 |
. | t = 10P0 (linear phase) . | t = 100P0 . | |||||||
---|---|---|---|---|---|---|---|---|---|
Case . | |$\log {\hat{\nu }_0}$| . | Aν . | ζν . | m . | ωm/Ω0 . | qm/Ω0 . | m . | 102am . | min[Ro(z = 0)] . |
B0 | −9 | 1 | n/a | 4 | 0.985 | 0.199 | 1 | 8.5 | −0.15 |
V0 | −6 | 1 | n/a | 4 | 0.985 | 0.199 | 1 | 6.8 | −0.11 |
V1 | −6 | 100 | 1.5 | 4 | 0.986 | 0.191 | 1 | 7.8 | −0.19 |
V2 | −6 | 100 | 1.0 | 4 | 0.986 | 0.182 | 1 | 4.9 | −0.21 |
V3 | −4 | 1 | n/a | 4 | 0.986 | 0.131 | 3 | 3.7 | −0.25 |
. | t = 10P0 (linear phase) . | t = 100P0 . | |||||||
---|---|---|---|---|---|---|---|---|---|
Case . | |$\log {\hat{\nu }_0}$| . | Aν . | ζν . | m . | ωm/Ω0 . | qm/Ω0 . | m . | 102am . | min[Ro(z = 0)] . |
B0 | −9 | 1 | n/a | 4 | 0.985 | 0.199 | 1 | 8.5 | −0.15 |
V0 | −6 | 1 | n/a | 4 | 0.985 | 0.199 | 1 | 6.8 | −0.11 |
V1 | −6 | 100 | 1.5 | 4 | 0.986 | 0.191 | 1 | 7.8 | −0.19 |
V2 | −6 | 100 | 1.0 | 4 | 0.986 | 0.182 | 1 | 4.9 | −0.21 |
V3 | −4 | 1 | n/a | 4 | 0.986 | 0.131 | 3 | 3.7 | −0.25 |
. | t = 10P0 (linear phase) . | t = 100P0 . | |||||||
---|---|---|---|---|---|---|---|---|---|
Case . | |$\log {\hat{\nu }_0}$| . | Aν . | ζν . | m . | ωm/Ω0 . | qm/Ω0 . | m . | 102am . | min[Ro(z = 0)] . |
B0 | −9 | 1 | n/a | 4 | 0.985 | 0.199 | 1 | 8.5 | −0.15 |
V0 | −6 | 1 | n/a | 4 | 0.985 | 0.199 | 1 | 6.8 | −0.11 |
V1 | −6 | 100 | 1.5 | 4 | 0.986 | 0.191 | 1 | 7.8 | −0.19 |
V2 | −6 | 100 | 1.0 | 4 | 0.986 | 0.182 | 1 | 4.9 | −0.21 |
V3 | −4 | 1 | n/a | 4 | 0.986 | 0.131 | 3 | 3.7 | −0.25 |
We then simulate discs with floor viscosity |$\hat{\nu }_0=10^{-6}$|. The control run case V0 has Aν = 1. Thus, case V0 is the viscous version of case B0. We then consider models where the kinematic viscosity increases by a factor Aν = 100 for z > ζνH0 at the bump radius. We choose ζν = 1.5and 1.0 for cases V1 and V2, respectively. This gives an upper viscous layer of thickness 0.5H and H at R = r0. (See Fig. 1 for a plot of the kinematic viscosity profile for case V2.) For cases V1 and V2, the transition thickness is fixed to Δζν = 0.2. Finally, we consider a high-viscosity run, case V3, with |$\hat{\nu }_0=10^{-4}$| and Aν = 1. This is equivalent to extending the viscous layer in case V1/V2 to the entire vertical domain.
Inviscid reference case
Fig. 3 shows the density fluctuation and Rossby number for case B0. The dominant linear mode is m = 4 with a growth rate 0.2Ω0, consistent with recent 3D linear calculations (Meheut et al. 2012a; Lin 2013b). The non-linear outcome of the RWI is vortex formation (Li et al. 2000). Four vortices develop initially and then merge on a dynamical time-scale into a single vortex. Case B0 evolves similarly to previous simulations of the RWI in an inviscid disc (e.g. Meheut et al. 2010, 2012b, where more detailed analyses are given). This, together with the agreement with linear calculations, demonstrates the ability of the plutocode to capture the RWI.
The effect of a viscous layer
We now examine viscous cases V0–V3. Recall from Table 1 that the viscous layer (with |$\hat{\nu }\sim 10^{-4}$|) occupies the uppermost 0, 25, 50 and 100 per cent of the vertical domain at R = r0 for cases V0, V1, V2 and V3, respectively.
We first compare the viscous case V0 to the inviscid case B0. Table 1 shows that despite increasing the viscosity by a factor of 103, the change to the linear mode frequencies is negligible. The value of am and minimum Rossby number indicate that the final vortex in V0 is only slightly weaker than that in B0. This is also reflected in Fig. 3 (case B0) and the leftmost column in Fig. 4 (case V0). Case V0 develops a more elongated vortex with smaller |Δρ| than that in B0.
As we introduce and thicken the viscous layer from case V0 to V3, the dominant linear mode remains at m = 4 (Table 1), but linear growth rate does appreciably decrease (by ∼34 per cent from case V0 to V3). However, these linear growth time-scales are still ∼P0. We thus obtain the important result that viscosity (layered or not) does not significantly affect the linear instability because the RWI grows dynamically even in the high-viscosity disc.
The effect of layered viscosity in the non-linear regime is more complicated. The bottom row of Fig. 4 compares the Rossby number associated with the overdensities. Thickening the viscous layer decreases the vortex aspect ratio. Since their widths remain at ∼2H0, the vortices become smaller with increasing viscosity. This is partly attributed to fewer vortex merging events having occurred as viscosity is increased, which usually results in larger but weaker vortices (smaller |Ro|). If merging is resisted, then each vortex can grow individually. Strangely, vortices become stronger (more negative Ro) as viscosity is increased.
Fig. 5 compares the perturbed kinetic energy for cases B0, V0 and V1, which are all dominated by a single vortex in quasi-steady state. We compute W1 and compare its average over the disc atmosphere and over the disc bulk. There is only a minor difference between the perturbed kinetic energy densities between the disc bulk and atmosphere, even in case V1 where the kinematic viscosity in the two regions differs by a factor of ∼102. This suggests that the vortex evolves two dimensionally.
The energy perturbation in cases B0 and V0 are both subject to slow decay (a result also observed by Meheut et al. 2012a). By contrast, case V1, which includes a high-viscosity layer, does not show such a decay. We discuss this counter-intuitive result below.
Order of magnitude comparison of time-scales
Meheut, Lovelace & Lai (2013) argued that tRWI is also the vortex turn-over time tturn when the instability saturates and the linear phase terminates. Then, tν ∼ 10tturn, implying that viscous effects are unimportant over one turn-over time. However, if we estimate a vortex turn-over time as |$t_\mathrm{turn} \sim 2\pi /|Ro|\Omega$|, then |$t_\nu \sim (h^2|Ro|/2\pi \hat{\nu })t_\mathrm{turn}$|. Inserting |$h=0.1,\,\hat{\nu }=10^{-4}$| and |Ro| = 0.25 (case V3) gives tν ∼ 4tturn. Therefore, depending on the vortex shape, tν may not be much larger than tturn.
In any case, our high-viscosity simulations span several local viscous time-scales, tsim ∼ 10tν (for |$\hat{\nu }\sim 10^{-4}$|), so viscous damping should have taken place, making the observation that Ro becomes more negative as the viscous layer increases from case V0 to V3, a surprising result. However, recall that we imposed a stationary, radially structured viscosity profile consistent with a steady-state disc containing a density bump. We suggest that for such setups, viscosity attempts to restore the initial disc profile, i.e. the initial PV minimum, thereby acting as a vorticity source.
Potential vorticity evolution
When viscosity is small, the local viscous time-scale tν is long compared to our simulation time-scale tsim. Then, vortex formation weakens the PV minimum with viscosity playing no role. Increasing viscosity eventually makes tν < tsim. This means that over the course of the simulation, our spatially fixed viscosity profile can act to recover the initial PV minimum.
We also notice reduced vortex migration in Fig. 4 with increased viscosity (e.g. the vortex in case V0 has migrated inwards to R ≃ 0.9r0, while that in case V3 remains near R ≃r0). Paardekooper et al. (2010) have shown that vortex migration can be halted by a surface density bump which, in our case, can be sourced by the radially structured viscosity profile.
We conjecture that in the non-linear regime there is competition between destruction of the background PV minimum by the vortices and reformation of the initial radial PV minimum by the imposed viscosity profile. The latter effect should favour the RWI, since the anticyclonic vortices are regions of local vorticity minima. In this way, viscosity acts to source vorticity, and this effect outweighs viscous damping of the linear perturbations. We discuss additional simulations supporting this hypothesis in Appendix A.
VORTEX FORMATION AT PLANETARY GAP EDGES IN LAYERED DISCS
The previous simulations, while necessary to isolate the effect of viscosity on the linear RWI, has the disadvantage that the radially structured viscosity profile can act to source radial disc structure in the non-linear regime. In this section, we employ a radially smooth viscosity profile and use disc–planet interaction to create the disc structure required for instability. Then, we expect viscosity to only act as a damping mechanism.
Vortex formation at gap edges is a standard result in 2D and 3D hydrodynamical simulations of giant planets in low-viscosity discs (de Val-Borro et al. 2007; Lin & Papaloizou 2010, 2011; Lin 2012b; Zhu et al. 2013). The fact that this is due to the RWI has been explicitly verified through linear stability analysis (de Val-Borro et al. 2007; Lin & Papaloizou 2010). Here, we simulate gap-opening giant planets in 3D discs where the kinematic viscosity varies with height above the disc mid-plane. Our numerical setup is similar to those that in Pierens & Nelson (2010), but our interest is gap stability in a layered disc.
Radially smooth viscosity profile for disc–planet interaction
Disc–planet simulations
We simulate locally isothermal discs with constant aspect ratio h = 0.05 (by choosing q = 1), vertical extent nh = 3 scaleheights and radial extent [rin, rout] = [0.4, 2.5]r0. Initially, the surface density is smooth (A = 1) with zero meridional velocity (vr = vθ = 0). The standard resolution is (Nr, Nθ, Nϕ) = (256, 96, 768), corresponding to 6, 32 and 6 cells per H along the r, θ and ϕ directions at the reference radius. We apply a damping rate |$\hat{\gamma }=2$| with the reference velocity field |$\boldsymbol {v}_\mathrm{ref}=(0,0,v_\phi )$| in spherical coordinates.
We insert into the disc a planet of mass Mp = 10− 3M*, which corresponds to a Jupiter mass planet if M* = M⊙. The softening length for the planet potential is ϵp = 0.5rh. The planet potential is switched on smoothly over t ∈ [0, 10]P0. We note that the disc can be considered as two-dimensional for gap-opening giant planets, because the Hill radius rh exceeds the local scaleheight H (rh/H ≃ 1.4 in our cases). Fig. 7 shows a typical PV profile associated with the gap induced by the planet.
We remark that, apart from the viscosity prescription, the above choice of physical and numerical parameter values is typical for global disc–planet simulations (e.g. de Val-Borro et al. 2006; Mignone et al. 2012).
Results
Table 2 summarizes the disc–planet simulations. The main simulations to be discussed are cases P0–P1, with a floor viscosity of |$\hat{\nu }_0=2.5\times 10^{-7}$|. The fiducial run P0 has Aν = 1, i.e. no viscous layer, so that α ∼ 10−4 everywhere. The more typical viscosity value adopted for disc–planet simulations, |$\hat{\nu }\sim 10^{-5}$| or α ∼ 10−3, is known to suppress vortex formation (de Val-Borro et al. 2007; Mudryk & Murray 2009). Thus, vortex formation is expected in case P0. For cases P0.5 and P1, we set Aν = 100 with transition angle ζν = 2.5h and 2h, respectively, so the viscous layer with α ∼ 10−2 occupies the uppermost 0.5H and H of the vertical domain. Case P0R is case P0 restarted from t = 100P0 with the layered viscosity profile of case P1.
Case . | |$10^6\hat{\nu }_0$| . | Aν . | Visc. layer . | |$10^2\overline{a}_1$| . | 102a1(200P0) . | Comment . |
---|---|---|---|---|---|---|
P0 | 0.25 | 1 | 0 | 18.3 | 14.5 | Single vortex by t = 130P0 and persists until end of sim. |
P0.5 | 0.25 | 100 | 0.5H | 17.4 | 8.6 | Single vortex by t = 90P0 and persists until end of sim. |
P0R | 0.25 | 1→100 | 0 → H | 10.1 | 2.5 | Single vortex by t = 130P0, disappears after t = 180P0 |
P1 | 0.25 | 100 | H | 12.5 | 2.1 | Single vortex by t = 80P0, disappears after t = 170P0 |
Pb0 | 1.0 | 1 | 0 | 9.2 | 5.4 | Single vortex by t = 130P0, disappears after t = 160P0 |
Pb0.5 | 1.0 | 10 | 0.5H | 7.3 | 3.0 | Similar to Pb0 |
Pb1 | 1.0 | 10 | H | 5.7 | 2.2 | Single vortex by t = 120P0, disappears after t = 140P0 |
Pc0.5 | 1.0 | 100 | 0.5H | 4.2 | 2.3 | Two vortices at t ∼ 100P0, no vortices after t ∼ 120P0 |
Pc1 | 1.0 | 100 | H | 2.2 | 1.6 | Two weak vortices at t ∼ 60P0, no vortices after t ∼ 80P0 |
Case . | |$10^6\hat{\nu }_0$| . | Aν . | Visc. layer . | |$10^2\overline{a}_1$| . | 102a1(200P0) . | Comment . |
---|---|---|---|---|---|---|
P0 | 0.25 | 1 | 0 | 18.3 | 14.5 | Single vortex by t = 130P0 and persists until end of sim. |
P0.5 | 0.25 | 100 | 0.5H | 17.4 | 8.6 | Single vortex by t = 90P0 and persists until end of sim. |
P0R | 0.25 | 1→100 | 0 → H | 10.1 | 2.5 | Single vortex by t = 130P0, disappears after t = 180P0 |
P1 | 0.25 | 100 | H | 12.5 | 2.1 | Single vortex by t = 80P0, disappears after t = 170P0 |
Pb0 | 1.0 | 1 | 0 | 9.2 | 5.4 | Single vortex by t = 130P0, disappears after t = 160P0 |
Pb0.5 | 1.0 | 10 | 0.5H | 7.3 | 3.0 | Similar to Pb0 |
Pb1 | 1.0 | 10 | H | 5.7 | 2.2 | Single vortex by t = 120P0, disappears after t = 140P0 |
Pc0.5 | 1.0 | 100 | 0.5H | 4.2 | 2.3 | Two vortices at t ∼ 100P0, no vortices after t ∼ 120P0 |
Pc1 | 1.0 | 100 | H | 2.2 | 1.6 | Two weak vortices at t ∼ 60P0, no vortices after t ∼ 80P0 |
Case . | |$10^6\hat{\nu }_0$| . | Aν . | Visc. layer . | |$10^2\overline{a}_1$| . | 102a1(200P0) . | Comment . |
---|---|---|---|---|---|---|
P0 | 0.25 | 1 | 0 | 18.3 | 14.5 | Single vortex by t = 130P0 and persists until end of sim. |
P0.5 | 0.25 | 100 | 0.5H | 17.4 | 8.6 | Single vortex by t = 90P0 and persists until end of sim. |
P0R | 0.25 | 1→100 | 0 → H | 10.1 | 2.5 | Single vortex by t = 130P0, disappears after t = 180P0 |
P1 | 0.25 | 100 | H | 12.5 | 2.1 | Single vortex by t = 80P0, disappears after t = 170P0 |
Pb0 | 1.0 | 1 | 0 | 9.2 | 5.4 | Single vortex by t = 130P0, disappears after t = 160P0 |
Pb0.5 | 1.0 | 10 | 0.5H | 7.3 | 3.0 | Similar to Pb0 |
Pb1 | 1.0 | 10 | H | 5.7 | 2.2 | Single vortex by t = 120P0, disappears after t = 140P0 |
Pc0.5 | 1.0 | 100 | 0.5H | 4.2 | 2.3 | Two vortices at t ∼ 100P0, no vortices after t ∼ 120P0 |
Pc1 | 1.0 | 100 | H | 2.2 | 1.6 | Two weak vortices at t ∼ 60P0, no vortices after t ∼ 80P0 |
Case . | |$10^6\hat{\nu }_0$| . | Aν . | Visc. layer . | |$10^2\overline{a}_1$| . | 102a1(200P0) . | Comment . |
---|---|---|---|---|---|---|
P0 | 0.25 | 1 | 0 | 18.3 | 14.5 | Single vortex by t = 130P0 and persists until end of sim. |
P0.5 | 0.25 | 100 | 0.5H | 17.4 | 8.6 | Single vortex by t = 90P0 and persists until end of sim. |
P0R | 0.25 | 1→100 | 0 → H | 10.1 | 2.5 | Single vortex by t = 130P0, disappears after t = 180P0 |
P1 | 0.25 | 100 | H | 12.5 | 2.1 | Single vortex by t = 80P0, disappears after t = 170P0 |
Pb0 | 1.0 | 1 | 0 | 9.2 | 5.4 | Single vortex by t = 130P0, disappears after t = 160P0 |
Pb0.5 | 1.0 | 10 | 0.5H | 7.3 | 3.0 | Similar to Pb0 |
Pb1 | 1.0 | 10 | H | 5.7 | 2.2 | Single vortex by t = 120P0, disappears after t = 140P0 |
Pc0.5 | 1.0 | 100 | 0.5H | 4.2 | 2.3 | Two vortices at t ∼ 100P0, no vortices after t ∼ 120P0 |
Pc1 | 1.0 | 100 | H | 2.2 | 1.6 | Two weak vortices at t ∼ 60P0, no vortices after t ∼ 80P0 |
Density evolution
Fig. 8 compares the time evolution of the mid-plane density perturbation δρ(z = 0) for cases P0, P0.5 and P1. In all the cases, we observed that the RWI with m = 4 develops early on (t ≃ 20P0), consistent with the limited effect of viscosity on the linear instability, as found above. The no-layer case P0 and layered case P0.5 (viscous layer of 0.5H) behave similarly, showing that a thin viscous layer has little effect on the evolution of the unstable gap edge, at least over the simulation time-scale of 200P0.
Case P1 evolves quite differently from case P0. While a single vortex does form at t ∼ 100P0, it is transient, having disappeared at the end of the simulation for P1. The final m = 1 amplitude is about three times smaller than that in case P0 (Table 2). This result is significant because the upper viscous layer in case P1, of thickness H, only occupies ∼4 per cent of the total column density, but the vortex is still destroyed. This suggests that vortex survival at planetary gap edges requires low effective viscosity throughout the vertical fluid column.
Kinetic energy density
Here, we compare the m = 1 component of the kinetic energy density (W1) between the no-layer case P0, layered case P1 and case P0R which is P0 resumed from t = 100P0 with a viscous layer. Fig. 9 shows W1(t) averaged over the outer gap edge. For each case, we average W1 over the disc bulk and the atmosphere, and plot them separately in the figure.
The m = 1 component does not emerge from the linear instability, but is a result of non-linear vortex merging. Fig. 9 shows that merging is accelerated by a viscous layer: the single vortex appears at t ∼ 70P0 for case P1 but only forms at t ∼ 120P0 for case P0. Also note that for all cases, W1 in the disc bulk (thick lines) is similar to that in the disc atmosphere (thin lines), implying that the m = 1 disturbance (i.e. the vortex) evolves two dimensionally. We checked that this is consistent with the Froude number Fr ≡ |Ro|H/z < 1 away from the mid-plane (Barranco & Marcus 2005; Oishi & Mac Low 2009).
Case P0R shows that introducing a viscous layer eventually destroys the vortex. The local viscous time-scale is tν ≡ H2/ν ≳ 16P0, so on short time-scales after introducing the viscous layer (t = 100P0), vortex merging proceeds in case P0R similarly to case P0 (t ∈ [100, 110]P0). However, W1 decays for t > 110P0 and evolves towards that of case P1. We expect viscosity to damp the m = 1 disturbance in the disc atmosphere between t ∈ [110, 200]P0 because this corresponds to ∼6tν, but the disturbance in the disc bulk is also damped out: the evolution remains two-dimensional.
We emphasize that the kinetic energy is dominated by horizontal motions, with |$\mathrm{max}(|v_z|/|\boldsymbol {v}|) < 0.03$| at the outer gap edge (r ∈ [1.2, 1.6]r0). Vertical motions are well subsonic. When averaged over z ∈ [0, 2]H and z ∈ [2, 3]H, the vertical Mach number Mz ≡ |vz|/cs ≃ 0.05, 0.08 (P0), Mz ≃ 0.04, 0.06 (P0R) and Mz ≃ 0.05, 0.06 (P1), respectively.
Potential vorticity
We examine the PV evolution for case P0R in Fig. 10. To highlight the vortices, which are positive (negative) density (vertical vorticity) perturbations, we show the inverse PV perturbation, |$\delta \eta _z^{-1}\equiv \eta _z(t=0)/\eta _z - 1$|. As noted above, a single vortex still forms despite introducing a viscous layer at t = 100P0. However, it decays rapidly compared to case P0. The region with |$\delta \eta _z^{-1}>0$| (i.e. the vortex) elongates and shifts outwards from R ≃1.38r0 at t = 140P0 to R ≃1.5r0 at t = 200P0, by which the vortex has disappeared. (A similar evolution was observed for case P1.) The vortex is stretched azimuthally much more than radially. This is not surprising since the imposed viscosity profile is axisymmetric. The important point is that viscosity is only large near the disc surface, but still has a significant effect on the vortex.
Resolution check
We repeated simulations P0 and P1 with resolution (Nr, Nθ, Nϕ) = (512, 96, 1536), corresponding to 12 and 32 cells per scaleheight in (r, ϕ) and θ, respectively. We denote these runs as P0HR and P1HR below.
We observe similar evolution in P0HR and P1HR as their standard resolution versions. However, due to lower numerical diffusion, we find stronger vortices in P0HR. Although the vortex in P1HR persisted longer than the standard resolution run, it was still subject to rapid decay in comparison with P0HR. At t = 200P0, we find the m = 1 amplitude to be a1 = 0.29 and 0.10 at the outer gap edge, respectively, for P0HR and P1HR, a similar contrast as that between P0 and P1. A weak overdensity was still observed in P1HR at t = 200P0, but it further decays to a1 = 0.06 at t = 230P0 and there is no vortex. By contrast, P0HR was simulated to t = 250P0 and the vortex survived with little decay (a1 = 0.25).
Interestingly, we observe small-scale (≲ H) disturbances inside the vortex in P0HR. This is shown in the left-hand panel of Fig. 11 in terms of the (inverse) PV perturbation. We checked that the density field remains smooth, so this small-scale structure is due to vorticity variations. This is unlikely the elliptic instability (Lesur & Papaloizou 2009), though, because the numerical resolution is still insufficient for studying such instabilities, especially since the vortex is very elongated with large aspect ratio ∼10 (so the elliptic instability is weak; Lesur & Papaloizou 2009). Despite the disturbances, the vortex overdensity in P0HR remains coherent until the end of the simulation, possibly because the planet maintains the condition for RWI. On the other hand, the vortex in the layered case P1HR does not develop small-scale structure (Fig. 11, right-hand panel), yet it is destroyed by the end of the simulation.
Additional simulations
Locally isothermal, low-viscosity discs are vulnerable to the so-called vertical shear instability because ∂zΩi ≠ 0 (Nelson, Gressel & Umurhan 2013). Nelson et al. employed a radial resolution ≳60 cells per H to resolve this instability because it involves small radial wavelengths (≪H). Our numerical resolution is unlikely to capture this instability. Nevertheless, we have performed additional simulations designed to eliminate the vertical shear instability.
Larger floor viscosity
We performed several simulations with |$\hat{\nu }_0=10^{-6}$|. A viscosity of |$\hat{\nu }\sim 10^{-6}$| is expected to damp the vertical shear instability (Nelson et al. 2013), while still permitting the gap-edge RWI. Table 2 summarizes these cases with Aν = 10 (‘Pb’ runs) and Aν = 100 (‘Pc’ runs).
In these simulations, we find that vortices eventually decay, even in the no-layer case Pb0. For Aν = 10, the layered cases Pb0.5 and Pb1 evolve similarly to Pb0: three vortices formed by t ∼ 30P0, merging into two vortices by t ∼ 40P0, then finally into a single vortex by t ∼ 130P0, which subsequently decays. However, the final vortex decays faster in the presence of a viscous layer. This is shown in Fig. 12, which compares the m = 1 kinetic energy density for cases Pb0 and Pb1. The evolution only begins to differ after the single vortex has formed.
For Aν = 100 (cases Pc0.5 and Pc1), we find that the m = 2 amplitude dominated over m = 1, so a single-vortex configuration never forms. For both Pc0.5 and Pc1, the m = 2 (two-vortex configuration) amplitude decreases from t ∼ 50P0. For case Pc1, the vortices are transient features and are entirely absent for t ≳ 80P0.
Strictly isothermal discs
We repeated simulations Pb0, Pb1 and Pc1 with a strictly isothermal equation of state (q = 0). These are summarized in Table 3. Fig. 13 compares their m = 1 kinetic energy density evolution at the outer gap edge. Consistent with the above simulations, a viscous layer causes a faster decay in this quantity. Most interesting, though, is that we found case Iso2 (with a viscous layer of ∼H) only shows very weak non-axisymmetric perturbations early on (t ≲ 50P0): vortex formation is suppressed.
Case . | |$10^6\hat{\nu }_0$| . | Aν . | Visc. layer . | |$10^2\overline{a}_1$| . | Vortex . |
---|---|---|---|---|---|
Iso0 | 1.0 | 1 | 0 | 19.3 | Yes |
Iso1 | 1.0 | 10 | H0 | 12.8 | Yes |
Iso2 | 1.0 | 100 | H0 | 1.1 | No |
Case . | |$10^6\hat{\nu }_0$| . | Aν . | Visc. layer . | |$10^2\overline{a}_1$| . | Vortex . |
---|---|---|---|---|---|
Iso0 | 1.0 | 1 | 0 | 19.3 | Yes |
Iso1 | 1.0 | 10 | H0 | 12.8 | Yes |
Iso2 | 1.0 | 100 | H0 | 1.1 | No |
Case . | |$10^6\hat{\nu }_0$| . | Aν . | Visc. layer . | |$10^2\overline{a}_1$| . | Vortex . |
---|---|---|---|---|---|
Iso0 | 1.0 | 1 | 0 | 19.3 | Yes |
Iso1 | 1.0 | 10 | H0 | 12.8 | Yes |
Iso2 | 1.0 | 100 | H0 | 1.1 | No |
Case . | |$10^6\hat{\nu }_0$| . | Aν . | Visc. layer . | |$10^2\overline{a}_1$| . | Vortex . |
---|---|---|---|---|---|
Iso0 | 1.0 | 1 | 0 | 19.3 | Yes |
Iso1 | 1.0 | 10 | H0 | 12.8 | Yes |
Iso2 | 1.0 | 100 | H0 | 1.1 | No |
SUMMARY AND DISCUSSION
We have performed customized hydrodynamic simulations of non-axisymmetric instabilities in 3D viscous discs. We adopted height-dependent kinematic viscosity profiles, such that the disc mid-plane is of low viscosity (α ∼ 10−4) and the disc atmosphere is of high viscosity (α ∼ 10−2). We were motivated by the question of whether or not the RWI, and subsequent vortex formation, operates in layered accretion discs.
We first considered viscous disc equilibria with a radial density bump and varied the vertical dependence of viscosity. This setup can isolate the effect of viscosity on the linear RWI. We found that the linear RWI is unaffected by viscosity, layered or not. The viscous RWI remains dynamical and leads to vortex formation on time-scales of a few 10 s of orbits. We continued these simulations into the non-linear regime, but found that vortices became stronger as the viscous layer is increased in thickness. We suggest that this counter-intuitive result is an artefact of the chosen viscosity profile because it is radially structured: viscosity attempts to restore the equilibrium radial density bump, which favours the RWI. This effect outweighs viscosity damping the linear instability.
We also simulated vortex formation at planetary gap edges in layered discs with a radially smooth viscosity profile. Although vortex formation still occurs in layered discs, we found that the vortex can be destroyed even when the viscous layer only occupies the uppermost scaleheight of the vertical domain which is 3 scaleheights. This is significant because most of the disc mass is contained within 2 scaleheights (i.e. the low-viscosity layer), but simulations show that a viscous atmosphere inhibits long-term vortex survival. We found that the non-axisymmetric energy densities have weak vertical dependence, so the disturbance evolves two dimensionally. It appears that applying a large viscosity in the disc atmosphere is sufficient to damp the instability throughout the vertical column of the fluid.
Barranco & Marcus (2005) have described two 3D vortex models: tall columnar vortices and short finite-height vortices. Rossby vortices are columnar, i.e. the associated vortex lines extend vertically throughout the fluid column. One might have expected an upper viscous layer to damp out vortex motion in the disc atmosphere, leading to a shorter vortex. This, however, requires vortex lines to loop around the vortex (the short vortex of Barranco & Marcus). Such vortex loops form the surface of a torus (see, for example, fig. 1 in Barranco & Marcus 2005), instead of ending on vertical boundaries. This implies significant vertical motion near the vertical boundaries of the vortex, which would be difficult in our model because of viscous damping applied there. We suspect that this is why short/tall vortices fail to form/survive in our layered disc–planet models. We conclude that vortex survival at planetary gap edges require low viscosity (α ≲ 10−4) throughout the vertical extent of the disc.
Relation to other works
Pierens & Nelson (2010) simulated the orbital migration of giant planets in layered discs by prescribing a height-dependent viscosity profile. They considered significant reduction in kinematic viscosity in going from the disc atmosphere (the active zone, with α ∼ 10−2) to the disc mid-plane (the dead zone, with α ∼ 10−7). According to previous 2D simulations, such a low kinematic viscosity should lead to the RWI (de Val-Borro et al. 2006, 2007). However, Pierens & Nelson (2010) did not report vortex formation, nor are vortices visible from their plots. Very recent MHD simulations of giant planets in a layered disc also did not yield vortex formation (Gressel et al. 2013). These results are consistent with our simulations.
Oishi & Mac Low (2009) carried out MHD shearing box simulations with a resistivity profile that varied with height to model a layered disc: the disc atmosphere was MHD turbulent, while the disc mid-plane remained stable against the MRI. They envisioned the active zone as a vorticity source for vortex formation in the mid-plane. Although their setup is fundamentally different to ours, they also reported a lack of coherent vortices in the dead zone. They argued that the MHD turbulence in the active layer was not sufficiently strong to induce vortex formation in the dead zone. If MHD turbulence can be represented by a viscosity, the lack of tall columnar vortices in Oishi & Mac Low (2009) is consistent with our results. That is, even when MRI turbulence is only present in the disc atmosphere, it is able to damp out columnar vortices.
Caveats and outlooks
The most important caveat of the current model is the viscous prescription to mimic MRI turbulence. In doing so, an implicit averaging is assumed (Balbus & Papaloizou 1999). The spatial averaging should be taken on length-scales no less than the local disc scaleheight, and the temporal average should be taken on time-scales no less than the local orbital period. These are, however, the relevant scales for vortex formation via the RWI. Furthermore, our viscosity profile varies on length-scales comparable to or even less than H (e.g. the vertical transition between high- and low-viscosity layers). Nevertheless, our simulations demonstrate the importance of disc vertical structure on the RWI. That is, damping, even confined to the disc atmosphere, can destroy Rossby vortices.
Another drawback of a hydrodynamic viscous disc model is the fact that it cannot mimic magnetoelliptic instabilities (MEI), which are known to destroy vortices in magnetic discs (Lyra & Klahr 2011; Mizerski & Lyra 2012). A natural question is how Rossby vortices are affected by the MEI when it only operates in the disc atmosphere. Extension of this work to global non-ideal MHD simulations will be necessary to address RWI vortex formation in layered discs.
However, some improvements can be made within the viscous framework. A static viscosity profile neglects the back-reaction of the density field on the kinematic viscosity. Thus, our simulations only consider how Rossby vortices respond to an externally applied viscous damping. A more physical viscosity prescription should depend on the local column density (Fleming & Stone 2003), with viscosity decreasing with increasing column density. The effective viscosity inside Rossby vortices would be lowered relative to the background disc because disc vortices are overdensities. If the overdensity is large, then it is conceivable that vortex formation itself may render the effective viscosity to be sufficiently low throughout the fluid column to allow long-term vortex survival. Preparation for this study is underway and results will be reported in a follow-up paper.
This work benefited from extensive discussion with O. Umurhan. I also thank R. Nelson for discussion and M. de Val-Borro for a helpful report. Computations were performed on the CITA Sunnyvale cluster, as well as the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, Ontario Research Fund Research Excellence and the University of Toronto.
REFERENCES
APPENDIX A: ARTIFICIAL RADIAL DENSITY BUMPS WITH A RADIALLY SMOOTH VISCOSITY PROFILE
In Section 4, we found that vortices became stronger as the viscous layer thickness is increased, even though linear growth rates were reduced. Here, we present additional simulations to support the hypothesis that this is due to the localized radial structure in the viscosity profile.
This simulation is shown as the dotted line in Fig. A2 in terms of the m = 1 component of the kinetic energy density. We compare it to the corresponding case using the radially structured viscosity profile in Section 4 (i.e. the original case V2 but with lowered floor viscosity). Vortex formation occurs in both runs. With a radially smooth viscosity profile, the vortex decays monotonically after |W1| reaches maximum value of ∼0.05. Using the radially structured viscosity profile (solid line) gives a larger disturbance amplitude at the linear stage (max |W1| ∼ 0.08), and although it subsequently decays, the decay is halted for t ≳ 110P0.
The contrast between these cases show that the radial structure in the viscosity profile helps vortex survival. This experiment indicates that the dominant effect of viscosity is its influence on the evolution of the axisymmetric part of background disc. The radially structured viscosity profile is a source for the radial PV minimum, which is needed for the RWI.
Our result here is qualitatively similar to that in Regály et al. (2012), where a sharp viscosity profile was imposed in a 2D simulation and vortex formation ensues via the RWI. The vortex eventually disappears, but redevelops after the system returns to an axisymmetric state. This is because the imposed viscosity profile causes the disc to develop the required PV minimum for the RWI.