Parametric Study of the Rossby Wave Instability in a Two-dimensional Barotropic Disk. II. Nonlinear Calculations

Tomohiro Ono; Takayuki Muto; Kengo Tomida; Zhaohuan Zhu

doi:10.3847/1538-4357/aad54d

1. Introduction

Recent observations have revealed protoplanetary disks with lopsided structures, especially in transitional disks (e.g., Fukagawa et al. 2013; van der Marel et al. 2013; Casassus et al. 2015). It is important to know how the lopsided structures are formed and how they are related to the disk evolution and planet formation. One possible mechanism for producing such lopsided structures is capturing of dust particles (∼mm in size) by a large-scale gas vortex. Theoretically, it has been known that a large-scale vortex of gas can survive for a long time (Godon & Livio 2000) and can efficiently trap dust particles due to gas drag (e.g., Barge & Sommeria 1995; Chavanis 2000). With a vortex induced by an embedded planet, Zhu & Stone (2014) showed that some ALMA observations can be reproduced with three-dimensional (3D) MHD simulations, including dust particles.

Several vortex formation mechanisms have been proposed; the Rossby wave instability (RWI; e.g., Lovelace & Hohlfeld 1978; Lovelace et al. 1999), the baroclinic instability (e.g., Klahr & Bodenheimer 2003), the vertical shear instability (e.g., Goldreich & Schubert 1967; Nelson et al. 2013), the zombie vortex instability (e.g., Marcus et al. 2013), and pebble accretion onto planets (Owen & Kollmeier 2017). As far as the observations show, all the protoplanetary disks with lopsided structures are transitional disks, which have an inner cavity. In the case of a transitional disk, the existence of sharp variations of physical quantities (e.g., surface density) is naturally expected. When a protoplanetary disk has sharp radial variations, a large-scale vortex of gas is expected to be formed by the RWI. Therefore, the RWI is one of the most promising mechanisms for explaining the observed lopsided structures.

The RWI has been studied with linear stability analyses (Lovelace et al. 1999; Li et al. 2000; Umurhan 2010; Meheut et al. 2012b; Lin 2013). The existence of a vortensity local minimum is necessary for the RWI to take place (Lovelace et al. 1999). However, the necessary and sufficient conditions for the onset of the RWI have been unknown until recent years. In Ono et al. (2016, hereafter Paper I), we performed linear stability analyses of the RWI within the framework of two-dimensional (2D), barotropic, and purely hydrodynamic disks. We revealed the parameter sets where the disks are unstable against the RWI and derived the critical condition for the RWI in a semi-analytic form. The RWI has been also studied with numerical simulations (Li et al. 2001; Varnière & Tagger 2006; Lyra et al. 2008, 2009; Meheut et al. 2010, 2012a; Richard et al. 2013). However, our knowledge of the evolution and final outcomes of the RWI is still limited, partially due to a lack of systematic parameter surveys in 2D cylindrical coordinates. In this paper, we perform numerical simulations of the RWI in 2D, barotropic, and purely hydrodynamic disks. We explore a wide parameter space of initial surface density profiles and the disk temperature. We investigate the RWI evolution and the properties and evolution of the vortices formed by the RWI.

This paper is organized as follows. We describe our disk model and numerical setup in Section 2. We present results and discussions of the RWI evolution in Section 3 and the vortices formed by the RWI in Section 4. Section 5 is for the summary.

2. Disk Model and Numerical Method

2.1. Disk Models and Initial Condition

We consider barotropic and purely hydrodynamic disks orbiting a central star of mass M in global 2D cylindrical coordinates, which are the same as the model used in our linear stability analyses presented in Paper I. We assume that the disks are geometrically thin and neglect the effects of magnetic fields, viscosity, and self-gravity. We employ these rather simple assumptions in order to compare the numerical simulations with the results of Paper I in detail and to perform a systematically controlled parameter survey. Previous works showed that viscosity, self-gravity, an indirect term of gravity force, vertical stratification, and baroclinicity have some effects on the RWI or the vortices formed by the RWI (Lin 2012a, 2012b, 2013, 2014; Lovelace & Hohlfeld 2013; Miranda et al. 2016, 2017; Zhu & Baruteau 2016). Numerical calculations of the RWI with dust particles, planets, and magnetic fields have also been performed (Li et al. 2005; Inaba & Barge 2006; Lin & Papaloizou 2011a, 2011b; Lyra & Mac Low 2012; Regály et al. 2013; Fu et al. 2014a, 2014b; Hammer et al. 2017). However, we can capture the essential physics of the RWI even within the 2D, barotropic, and purely hydrodynamic framework.

Our numerical simulations employ a nonrotating frame centered on a star and a 2D cylindrical coordinate with (r, φ). The gravitational potential of the central star is given by ${\rm{\Phi }}(r)=-{GM}/r$ , where G is the gravitational constant. We denote the surface density by Σ and the (vertically integrated) pressure by P. We assume that the disk is barotropic, i.e., $P=P({\rm{\Sigma }})\ \propto {{\rm{\Sigma }}}^{{\rm{\Gamma }}}$ , where Γ is the effective adiabatic index of the gas. In our simulations, we consider only Γ = 5/3. The continuity equation is

$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}(r{\rm{\Sigma }}{v}_{r})+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial \varphi }({\rm{\Sigma }}{v}_{\varphi })=0,\end{eqnarray} \tag{ 1 }$

where t is the time, ${\boldsymbol{v}}(r,\varphi ,t)\equiv {v}_{r}(r,\varphi ,t)\ \hat{{\boldsymbol{r}}}\ +\ {v}_{\varphi }(r,\varphi ,t)\ \hat{{\boldsymbol{\varphi }}}$ is the velocity field, $\hat{{\boldsymbol{r}}}$ is the unit vector in the r direction, and $\hat{{\boldsymbol{\varphi }}}$ is the unit vector in the φ direction. The equations of motion are

$\begin{eqnarray}&&\displaystyle \frac{\partial {v}_{r}}{\partial t}+{v}_{r}\displaystyle \frac{\partial {v}_{r}}{\partial r}+\displaystyle \frac{{v}_{\varphi }}{r}\displaystyle \frac{\partial {v}_{r}}{\partial \varphi }-\displaystyle \frac{{v}_{\varphi }^{2}}{r}=-\displaystyle \frac{{GM}}{{r}^{2}}-\displaystyle \frac{\partial {\rm{\Pi }}}{\partial r},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&\displaystyle \frac{\partial {v}_{\varphi }}{\partial t}+{v}_{r}\displaystyle \frac{\partial {v}_{\varphi }}{\partial r}+\displaystyle \frac{{v}_{\varphi }}{r}\displaystyle \frac{\partial {v}_{\varphi }}{\partial \varphi }+\displaystyle \frac{{v}_{r}{v}_{\varphi }}{r}=-\displaystyle \frac{1}{r}\displaystyle \frac{\partial {\rm{\Pi }}}{\partial \varphi },\end{eqnarray} \tag{ 3 }$

where Π is the pressure function. For the barotropic flow with ${\rm{\Gamma }}\ne 1$ , Π is written as

$\begin{eqnarray}&&{\rm{\Pi }}\equiv \displaystyle \frac{{\rm{\Gamma }}}{{\rm{\Gamma }}-1}\displaystyle \frac{P}{{\rm{\Sigma }}}.\end{eqnarray} \tag{ 4 }$

From Equations (1)–(3), the equation of the vortensity conservation is obtained as

$\begin{eqnarray}&&\displaystyle \frac{\partial q}{\partial t}+{v}_{r}\displaystyle \frac{\partial q}{\partial r}+\displaystyle \frac{{v}_{\varphi }}{r}\displaystyle \frac{\partial q}{\partial \varphi }=0,\end{eqnarray} \tag{ 5 }$

where $q{(r,\varphi ,t)\equiv (\mathrm{rot}{\boldsymbol{v}})}_{z}/{\rm{\Sigma }}$ is the vortensity.

We perform numerical calculations with various initial conditions and investigate the RWI and vortices formed by the RWI. We adopt stationary (∂/∂t = 0), axisymmetric (∂/∂φ = 0), and circular (v_r = 0) flow as the initial conditions, which are denoted by the subscripts "0," e.g., Σ₀(r), P₀(r), and ${{\boldsymbol{v}}}_{0}(r)={v}_{\varphi 0}(r)\hat{{\boldsymbol{\varphi }}}$ . The initial surface density Σ₀(r) is given by a Gaussian bump on a uniform profile,

$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{{{\rm{\Sigma }}}_{{\rm{n}}}}=1+{{ \mathcal A }}_{0}\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{r-{r}_{{\rm{n}}}}{{\rm{\Delta }}{w}_{0}}\right)}^{2}\right],\end{eqnarray} \tag{ 6 }$

where Σ_n is the surface density of the uniform profile and r_n is the representative radius of the initial bump. This initial profile is the same as the "GB"-type profile in Paper I. There are two parameters to characterize the initial bump profile Σ₀: the contrast ${{ \mathcal A }}_{0}$ and the radial half-width Δw₀.

Since we consider the barotropic flow, ${P}_{0}(r)$ follows

$\begin{eqnarray}&&{P}_{0}(r)={S}_{0}{{\rm{\Sigma }}}_{0}^{{\rm{\Gamma }}},\end{eqnarray} \tag{ 7 }$

where S₀ is the entropy and constant. We define a dimensionless parameter h by

$\begin{eqnarray}&&h\equiv \displaystyle \frac{\sqrt{{\rm{\Gamma }}{S}_{0}{{\rm{\Sigma }}}_{{\rm{n}}}^{({\rm{\Gamma }}-1)}}}{{r}_{{\rm{n}}}{{\rm{\Omega }}}_{{\rm{n}}}},\end{eqnarray} \tag{ 8 }$

where ${{\rm{\Omega }}}_{{\rm{K}}}(r)\equiv \sqrt{{GM}/{r}^{3}}$ is the Kepler angular velocity and ${{\rm{\Omega }}}_{{\rm{n}}}\equiv {{\rm{\Omega }}}_{{\rm{K}}}({r}_{{\rm{n}}})$ . In this case, the initial entropy is written as

$\begin{eqnarray}&&{S}_{0}=\displaystyle \frac{{h}^{2}}{{\rm{\Gamma }}}\displaystyle \frac{{({r}_{{\rm{n}}}{{\rm{\Omega }}}_{{\rm{n}}})}^{2}}{{{\rm{\Sigma }}}_{{\rm{n}}}^{({\rm{\Gamma }}-1)}}=\mathrm{constant}.\end{eqnarray} \tag{ 9 }$

It is noted that h can be regarded as the dimensionless disk scale height, or, equivalently, the dimensionless sound speed. The value of h also represents the disk temperature. From Equation (2), the initial velocity field in the azimuthal direction ${v}_{\varphi 0}(r)$ is obtained as

$\begin{eqnarray}&&{v}_{\varphi 0}(r)=\sqrt{{v}_{{\rm{K}}}^{2}+r\displaystyle \frac{d{{\rm{\Pi }}}_{0}}{{dr}}},\end{eqnarray} \tag{ 10 }$

where ${{\rm{\Pi }}}_{0}(r)\equiv {\rm{\Gamma }}{S}_{0}{{\rm{\Sigma }}}_{0}^{{\rm{\Gamma }}-1}/({\rm{\Gamma }}-1)$ is the initial pressure function.

The initial conditions are characterized by three parameters: h, Δw₀, and ${{ \mathcal A }}_{0}$ . First, we fix h and Δw₀ and vary ${{ \mathcal A }}_{0}$ . The larger ${{ \mathcal A }}_{0}$ is, the more unstable against the RWI the system is. For an unstable configuration against the RWI, the largest linear growth rate of the RWI, ${\gamma }_{* }({{ \mathcal A }}_{0}:\ h,{\rm{\Delta }}{w}_{0})$ , monotonically increases with ${{ \mathcal A }}_{0}$ (see Paper I). If, however, ${{ \mathcal A }}_{0}$ exceeds a certain value, ${{ \mathcal A }}_{0,\max }(h,{\rm{\Delta }}{w}_{0})$ , the system violates Rayleigh's condition and is prone to rotational instability (see Appendix B.1). Since the linear growth rate of the rotational instability is typically larger than that of the RWI, we expect that the system that is unstable against the rotational instability immediately transfers to the marginally stable configuration of the rotational instability ( ${{ \mathcal A }}_{0}={{ \mathcal A }}_{0,\max }(h,{\rm{\Delta }}{w}_{0})$ ). We therefore consider the cases where the system does not violate Rayleigh's condition. In other words, we consider the cases with ${{ \mathcal A }}_{0}\lt {{ \mathcal A }}_{0,\max }(h,{\rm{\Delta }}{w}_{0})$ as the initial conditions. The maximum of the largest linear growth rate of the RWI is limited below the value of that with ${{ \mathcal A }}_{0}={{ \mathcal A }}_{0,\max }(h,{\rm{\Delta }}{w}_{0})$ , which we denote by ${\gamma }_{* ,\max }({{ \mathcal A }}_{0,\max }:\ h,{\rm{\Delta }}{w}_{0})$ .

When the three parameters $(h,{\rm{\Delta }}{w}_{0},{\gamma }_{* })$ are given, ${{ \mathcal A }}_{0}$ is uniquely determined. We vary h and Δw₀ in the ranges of h = [0.05, 0.1, 0.15, 0.2] and Δw₀/r_n = [0.02, 0.0356, 0.0632, 0.112, 0.2]. We also vary γ_* in the ranges of γ_*/Ω_n = [γ_*,max/Ω_n, 0.2, 0.15, 0.1, 0.05] for h = 0.1 and γ_*/Ω_n = [γ_*,max/Ω_n, 0.1] for $h\ \ne 0.1$ . We run 54 models in total, whose parameter sets are shown in Table 1. Note that we do not have the "h10w5g2" model because γ_*,max/Ω_n is smaller than 0.2 for h = 0.1 and Δw₀ = 0.2r_n. We show the most unstable azimuthal mode m_*, as well as the largest linear growth rate γ_* in Table 1. In addition, we calculate the linear growth rate of the RWI for each azimuthal mode m, γ_m, in the same manner as described in Paper I. The setup of the linear stability analyses and linear growth rates for 1 ≤ m ≤ 10 are shown in Appendix C.1. In this paper, we regard the "h10w3g1" model (h = 0.1, Δw₀ = 0.0632r_n, γ_*/Ω_n = 0.246, ${{ \mathcal A }}_{0}=0.439$ , and m_* = 4) as a fiducial case. When we investigate the overall properties of the RWI and vortices formed by the RWI, we always refer to the outcome of the "h10w3g1" model.

Table 1. The Parameter Sets of the Models

Name	h	Δw₀/r_n	γ_*/Ω_n	${{ \mathcal A }}_{0}$	m_*
h10w1g1	0.1	2.00E−2	0.227	4.04E−2	9
h10w1g2	0.1	2.00E−2	0.200	3.40E−2	8
h10w1g3	0.1	2.00E−2	0.150	2.41E−2	7
h10w1g4	0.1	2.00E−2	0.100	1.57E−2	6
h10w1g5	0.1	2.00E−2	0.050	8.70E−3	4
h10w2g1	0.1	3.56E−2	0.242	1.31E−1	6
h10w2g2	0.1	3.56E−2	0.201	1.06E−1	5
h10w2g3	0.1	3.56E−2	0.150	7.81E−2	5
h10w2g4	0.1	3.56E−2	0.100	5.46E−2	4
h10w2g5	0.1	3.56E−2	0.050	3.42E−2	3
h10w3g1	0.1	6.32E−2	0.246	4.39E−1	4
h10w3g2	0.1	6.32E−2	0.200	3.58E−1	4
h10w3g3	0.1	6.32E−2	0.150	2.77E−1	3
h10w3g4	0.1	6.32E−2	0.100	2.05E−1	3
h10w3g5	0.1	6.32E−2	0.050	1.42E−1	2
h10w4g1	0.1	1.12E−1	0.227	1.57E+0	3
h10w4g2	0.1	1.12E−1	0.200	1.38E+0	2
h10w4g3	0.1	1.12E−1	0.150	1.06E+0	2
h10w4g4	0.1	1.12E−1	0.100	8.02E−1	2
h10w4g5	0.1	1.12E−1	0.050	6.00E+1	2
h10w5g1	0.1	2.00E−1	0.191	5.66E+0	2
h10w5g3	0.1	2.00E−1	0.150	4.69E+0	2
h10w5g4	0.1	2.00E−1	0.100	3.39E+0	1
h10w5g5	0.1	2.00E−1	0.050	2.33E+0	1

h20w1g1	0.2	2.00E−2	0.209	1.00E−2	8
h20w1g4	0.2	2.00E−2	0.100	3.97E−3	4
h20w2g1	0.2	3.56E−2	0.225	3.17E−2	5
h20w2g4	0.2	3.56E−2	0.100	1.23E−2	3
h20w3g1	0.2	6.32E−2	0.237	1.00E−1	3
h20w3g4	0.2	6.32E−2	0.100	4.15E−2	3
h20w4g1	0.2	1.12E−1	0.237	3.17E−1	2
h20w4g4	0.2	1.12E−1	0.100	1.53E−1	2
h20w5g1	0.2	2.00E−1	0.192	9.79E−1	2
h20w5g4	0.2	2.00E−1	0.100	5.38E−1	1

h15w1g1	0.15	2.00E−2	0.216	1.78E−2	8
h15w1g4	0.15	2.00E−2	0.100	6.89E−3	5
h15w2g1	0.15	3.56E−2	0.233	5.69E−2	5
h15w2g4	0.15	3.56E−2	0.100	3.56E−2	3
h15w3g1	0.15	6.32E−2	0.240	1.83E−1	4
h15w3g4	0.15	6.32E−2	0.100	8.13E−1	2
h15w4g1	0.15	1.12E−1	0.233	6.00E−1	2
h15w4g4	0.15	1.12E−1	0.100	2.96E−1	2
h15w5g1	0.15	2.00E−1	0.193	1.96E+0	2
h15w5g4	0.15	2.00E−1	0.100	1.13E+0	1

h05w1g1	0.05	2.00E−2	0.244	1.68E−1	11
h05w1g4	0.05	2.00E−2	0100	7.14E−2	7
h05w2g1	0.05	3.56E−2	0.248	5.86E−1	7
h05w2g4	0.05	3.56E−2	0.100	2.77E−1	5
h05w3g1	0.05	6.32E−2	0.240	2.30E+0	4
h05w3g4	0.05	6.32E−2	0.100	1.18E+0	4
h05w4g1	0.05	1.12E−1	0.225	1.01E+1	3
h05w4g4	0.05	1.12E−1	0.100	5.28E+0	2
h05w5g1	0.05	2.00E−1	0.188	4.12E+1	2
h05w5g4	0.05	2.00E−1	0.100	2.52E+1	1

Note. Name: the name of the model. h: the dimensionless disk aspect ratio. Δw₀/r_n: the radial half-width of the initial bump normalized by r_n. γ_*/Ω_n: the largest linear growth rate against the RWI normalized by Ω_n. ${{ \mathcal A }}_{0}$ : the radial surface density contrast of the initial bump. m_*: the most unstable azimuthal mode.

Download table as: ASCII Typeset image

2.2. Numerical Method

We use the Athena++ code (J. M. Stone et al. 2018, in preparation) with the HLLC approximate Riemann solver, the second-order piecewise linear reconstruction, and the second-order van Leer time integrator. We assume barotropic flows for simplicity, and therefore we overwrite the pressure after every time step to satisfy Equation (7). The computational domain extends r_in < r < r_out, where we set r_in = 0.3 r_n and r_out = 2.5 r_n, in the radial direction and covers a full 2π in the azimuthal direction in the 2D cylindrical coordinates. We choose the radial extension of the numerical domain so that all the effective Lindblad resonances from the corotation point with the vortex center reside within the computational domain unless the vortices become too close to the boundaries.

We set the mesh structure so that the size of a cell is at least smaller than 0.04h in the vortex-forming region (r ∼ r_n). For h = 0.1, 0.15, and 0.2, the mesh has 576 cells in the radial direction and 1596 grids in the azimuthal direction. For h = 0.05, the mesh has 1296 cells in the radial direction and 3744 cells in the azimuthal direction. While the azimuthal spacing of cells is uniform, we make the radial spacing logarithmically constant and keep an aspect ratio of cells about unity. From the resolution study, the calculations with the mesh structure are high-resolution enough to discuss the results of this paper (see Appendix D).

We adopt the nonreflective boundary conditions (Godon 1996) in the radial direction and the periodic boundary conditions in the azimuthal direction. The nonreflective boundary conditions are designed to be nonreflective only for one-dimensional simple waves. Even for 2D nonlinear simulations, however, we have observed the strong reduction of the wave reflection at the radial boundaries. These nonreflective boundary conditions are also used in previous works (e.g., Paardekooper et al. 2010). Note that these nonreflective boundary conditions cannot vanish the wave reflection perfectly. However, the inner boundary does not have significant effects on the vortices formed by the RWI (see Appendix D).

We have further modified the original Athena++ code by introducing fast Fourier transform (FFT) filters. We perform the Fourier transform of ${\rm{\Sigma }},{{ \mathcal M }}_{r}$ , and ${{ \mathcal M }}_{\varphi }$ in the azimuthal direction at every radius and time step. We denote the Fourier components for the azimuthal mode m by ${ \mathcal F }{({\rm{\Sigma }})}_{m},{ \mathcal F }{({{ \mathcal M }}_{r})}_{m},$ and ${ \mathcal F }{({{ \mathcal M }}_{\varphi })}_{m}$ . We have developed two kinds of FFT filters, namely, "axisymmetric filter" and "single-mode filter." For the numerical relaxation of the initial conditions before the main calculations, we use the axisymmetric filter, where all the nonaxisymmetric ( $m\ne 0$ ) modes are filtered out. We use the single-mode filter, where all the nonaxisymmetric components except for m = k are filtered out, when we investigate the linear and weakly nonlinear regimes and the saturation of a specific azimuthal mode k (see Sections 3.1 and 3.2). After these filters modify the Fourier components, we recalculate ${\rm{\Sigma }},{{ \mathcal M }}_{r}$ , and ${{ \mathcal M }}_{\varphi }$ by the inverse Fourier transform and update the quantities in the calculations.

Before starting the main calculation of each run, we evolve the disk numerically using the axisymmetric filter for 10 orbits at r = r_n in order to relax the initial profile to a numerical equilibrium. We impose a small initial perturbation on the radial momentum, ${{ \mathcal M }}_{r}$ , to trigger the RWI and start the main calculation. The Fourier component of the radial momentum for an azimuthal mode m is defined by ${ \mathcal F }{({{ \mathcal M }}_{r})}_{m}$ . We perform two types of numerical calculations. The first is a single-mode calculation in which we focus on one specific azimuthal mode k. In the single-mode calculations, the initial perturbations satisfy $| { \mathcal F }{({{ \mathcal M }}_{r})}_{m}| ={10}^{-6}\,| { \mathcal F }{({{ \mathcal M }}_{\varphi })}_{m=0}|$ $\exp [{\{(r/{r}_{{\rm{n}}}-1)/0.2\}}^{2}/2]$ for m = k and $| { \mathcal F }{({{ \mathcal M }}_{r})}_{m=0}| =0$ otherwise, where ${{ \mathcal M }}_{\varphi }$ is the azimuthal momentum and ${ \mathcal F }{({{ \mathcal M }}_{\varphi })}_{m=0}$ is the axisymmetric Fourier component of the azimuthal momentum. At that point, the k mode is not restricted to the most unstable azimuthal mode. We also make the single-mode filter of the k mode work to filter out the other nonaxisymmetric components ( $m\ \ne \ k$ ) during calculations. We use the results of the single-mode calculations for the purpose of investigating the initial evolution and saturation of the RWI. The other is a white-noise calculation. In the white-noise calculations, the initial perturbations satisfy $| { \mathcal F }{({{ \mathcal M }}_{r})}_{m}| \,={10}^{-6}\ | { \mathcal F }{({{ \mathcal M }}_{\varphi })}_{m=0}|$ $\exp [{\{(r/{r}_{{\rm{n}}}-1)/0.2\}}^{2}/2]$ for 1 ≤ m ≤ 128. Note that we set a maximum azimuthal mode of the white-noise perturbation to m = 128 in order to avoid the effects from the numerical resolution. The phase of ${ \mathcal F }{({{ \mathcal M }}_{r})}_{m}$ is randomly varied for each m. In the white-noise calculations, we do not use the single-mode filter. If not stated otherwise, we refer to the white-noise calculations.

3. Evolution of the RWI

First of all, we give an overview of the RWI evolution. Figure 1 shows 2D snapshots of the surface density at τ = 0, 8, 11, 18, 20, 30, 50, 100, and 150 in the fiducial calculation, where $\tau \equiv t{{\rm{\Omega }}}_{{\rm{n}}}/2\pi$ is the time measured in the unit of the orbital period at r = r_n.

After the onset of the RWI, the perturbation shows linear and weakly nonlinear evolution. The saturation occurs when the perturbation becomes comparable to the initial axisymmetric bump. Then, the system enters the fully nonlinear regime at τ ∼ 8.3. At that time, four vortices are formed by fragmentation of the initial axisymmetric bump. The number of vortices formed initially is in accordance with the most unstable azimuthal mode of the RWI, m_*. The vortices coalesce one after another ( $4\to 3$ at τ ∼ 11, $3\to 2$ at τ ∼ 18, $2\to 1$ at τ ∼ 20). In the end, one quasi-stationary vortex remains after the final merger.

In this section, we consider each stage of the RWI evolution individually: the linear and weakly nonlinear regimes in Section 3.1, the saturation in Section 3.2, and the vortex merger in Section 3.3.

3.1. Linear and Weakly Nonlinear Regime of the RWI

Here we direct our attention to the linear and weakly nonlinear regimes of the RWI. We compare the results of the numerical calculations with those of the linear stability analyses in order to confirm the validity of our numerical calculations. We also take a step into the weakly nonlinear regime and study how applicable the linear stability analyses are to understand the RWI evolution.

For the sake of the comparison with the linear stability analyses, we separate the surface density in the numerical calculations into axisymmetric components and nonaxisymmetric components. The axisymmetric component corresponds to the azimuthally averaged surface density $\langle {\rm{\Sigma }}\rangle (r)$ . We define the nonaxisymmetric component of the m mode by ${{\rm{\Sigma }}}_{m}(r,\varphi )\equiv \mathrm{Real}[{ \mathcal F }{({\rm{\Sigma }})}_{m}\exp ({im}\varphi )]$ . The single-mode calculation for the m mode has only the m-mode component and the axisymmetric component. Therefore, the m-mode component is calculated by subtracting $\langle {\rm{\Sigma }}\rangle (r)$ from Σ(r, φ). On the other hand, the white-noise calculation requires using the single-mode filter in a post-process to obtain the nonaxisymmetric components.

We perform single-mode calculations for the m modes (1 ≤ m ≤ 10) and a white-noise calculation of the "h10w3g1" model. Figure 2 shows the time evolution of ${{\rm{\Sigma }}}_{m,\max }$ in the calculations, where ${{\rm{\Sigma }}}_{m,\max }$ is the maximum of Σ_m around r = r_n. The linear growth rates against the RWI are independently derived from the linear stability analyses. All the numerical calculations show excellent agreement with the linear analyses in the linear regime. The white-noise calculation also shows the weakly nonlinear growth of the mode with a small linear growth rate due to the coupling between the two modes with a large linear growth rate. For example, we can observe the mode-coupling regime between the m = 3 mode and the m = 4 mode to produce the m = 1 (=4–3) mode component in 4.4 ≤ τ ≤ 6.9. We find that the linear stability analyses predict the weakly nonlinear evolution precisely.

**Figure 2.** Time evolution of Σ_m,max in the single-mode calculations (red lines) and the white-noise calculation (blue lines) of the "h10w3g1" model for each azimuthal mode m (1 ≤ m ≤ 10). The black dashed lines show the linear growth, and the green dashed lines show the growth estimated from the mode–mode coupling based on the linear stability analyses. The last panel shows the time evolution of the maximum of ${\sum }_{m\geqslant 1}{{\rm{\Sigma }}}_{m}$ in the white-noise calculation (blue line).
Download figure:
Standard image High-resolution image

Figure 3 compares the distribution of Σ_m(r, φ) normalized by Σ_m,max in the single-mode calculations and white-noise calculation at τ = 5.6 and the surface density perturbation normalized by the maximum value derived in the linear stability analyses for 1 ≤ m ≤ 4. The azimuthal phase is shifted so that the point of Σ_m = Σ_m,max is at φ = 0. Note that we also denote the surface density perturbation of the linear stability analyses for the m mode by Σ_m(r, φ). Except for the m = 1 mode in the white-noise calculation, the profiles of Σ_m in the numerical calculations match those in the linear stability analyses. The discrepancy for the m = 1 mode occurs because the coupling between the m = 3 mode and the m = 4 mode becomes significant, and the m = 1 mode already enters the weakly nonlinear regime at τ = 5.6 in the white-noise calculation.

From Figures 2 and 3, our numerical calculations agree with the linear stability analyses in the linear and weakly nonlinear regimes. Therefore, our numerical calculations and linear stability analyses are reliable in these regimes.

3.2. Saturation Mechanism of the RWI

As shown in Figure 2, the RWI saturation occurs when the amplitude of the nonaxisymmetric components becomes comparable to ${{ \mathcal A }}_{0}$ in both the single-mode calculations and the white-noise calculation. During the growth of the nonaxisymmetric components, the axisymmetric components, or the m = 0 mode components, also evolve due to the couplings of the nonaxisymmetric components. For example, a self-coupling of the m = k mode can produce the m = 0 (=k − k) mode component. As another example, the couplings between three or more modes can also produce the m = 0 mode components. Here we attempt to explain the saturation mechanism of the RWI by investigating the time evolution of the axisymmetric components.

We analyze the axisymmetric components in the single-mode calculation for the m = 4 mode and the white-noise calculation of the "h10w3g1" model. Since the radial profiles of the azimuthally averaged surface density $\langle {\rm{\Sigma }}\rangle (r)$ resemble a Gaussian bump during the RWI evolution, as seen in panel (a) of Figure 4, we measure the location of the peak r_p, the contrast ${ \mathcal A }$ , and the half-width Δw of the bump by fitting with $\langle {\rm{\Sigma }}\rangle /{{\rm{\Sigma }}}_{{\rm{n}}}={ \mathcal A }\exp [-{\{(r-{r}_{{\rm{p}}})/{\rm{\Delta }}w\}}^{2}/2]+1$ . From panels (b)–(d) of Figure 4, r_p and ${ \mathcal A }$ start to decrease and Δw starts to increase a few orbits before the saturation in both calculations. While the change of r_p is gradual, the changes of ${ \mathcal A }$ and Δw are rapid. These mean that the axisymmetric components approach the stable configurations against the RWI during the RWI evolution.

**Figure 4.** Time evolution of the azimuthally averaged profiles in the white-noise calculation (solid lines) and single-mode calculation for the m = 4 mode (dashed lines) of the "h10w3g1" model. Panel (a) shows the initial surface density (τ = 0; black line) and the azimuthally averaged surface density at τ = 7.5 (red), 8.3 (blue), and 20 (green). The time evolution of r_p, ${ \mathcal A }$ , and Δw is shown in panels (b), (c), and (d), respectively.
Download figure:
Standard image High-resolution image

**Figure 4.** Time evolution of the azimuthally averaged profiles in the white-noise calculation (solid lines) and single-mode calculation for the m = 4 mode (dashed lines) of the "h10w3g1" model. Panel (a) shows the initial surface density (τ = 0; black line) and the azimuthally averaged surface density at τ = 7.5 (red), 8.3 (blue), and 20 (green). The time evolution of r_p, ${ \mathcal A }$ , and Δw is shown in panels (b), (c), and (d), respectively.
Download figure:
Standard image High-resolution image

In order to investigate quantitatively the time evolution of the axisymmetric components, we use the semi-analytic condition for the onset of the RWI derived in Paper I,

$\begin{eqnarray}&&{\eta }_{m}\equiv {\int }_{{r}_{\mathrm{IR}}}^{{r}_{\mathrm{OR}}}\sqrt{-{D}_{\mathrm{MS},m}(r)}{dr}\gtrsim {\eta }_{{\rm{c}}},\end{eqnarray} \tag{ 11 }$

where D_MS,m is the effective potential of the m mode if the system is assumed to be marginally stable against the RWI of the m mode and r_IR and r_OR are the radii where D_MS,m vanishes. The threshold of the condition, η_c, is roughly equal to $\pi /(2\sqrt{2})$ when the profile of D_MS,m for r_IR < r < r_OR is approximated by a parabolic function. We show the detailed expression for D_MS,m in Appendix B.2. Since η_m depends on the azimuthal mode m and the axisymmetric components, η_m evolves with the axisymmetric components if m is fixed. Calculating η_m every one-tenth orbit, we study the time evolution of the stability of the axisymmetric components against the RWI for the m mode.

First, we look at the time evolution of η₄ in the single-mode calculation for the m = 4 mode because the m = 4 mode is the most unstable azimuthal mode of the "h10w3g1" model. As shown in panel (a) of Figure 5, η₄ is initially larger than $\pi /(2\sqrt{2})$ so that the system is unstable against the RWI of the m = 4 mode. As the RWI evolves, η₄ decreases and becomes smaller than $\pi /(2\sqrt{2})$ at τ ≈ 7.5. This means that the axisymmetric component approaches the stable configuration during the RWI evolution and reaches the marginally stable configuration at the RWI saturation. The same thing occurs in the white-noise calculation, where the RWI is saturated at τ ≈ 8.3, even though the calculation contains all the nonaxisymmetric components, as well as the axisymmetric component, as shown in panel (b) of Figure 5. Therefore, we consider that the RWI saturation occurs when axisymmetric components become marginally stable against the RWI for the most unstable azimuthal mode of the initial conditions. This indicates that the evolution of the axisymmetric components is mainly due to the self-coupling of the most unstable azimuthal mode. We also find that η_m > 4 is smaller and η_m < 4 is larger than $\pi /(2\sqrt{2})$ at the RWI saturation in the white-noise calculation. In other words, the axisymmetric components are stable for the higher modes but still unstable for the lower modes at the RWI saturation.

**Figure 5.** Time evolution of the stability of the axisymmetric components against the RWI in the "h10w3g1" runs. Panel (a) shows the time evolution of η₄ in the single-mode calculation for the m = 4 mode. The time evolution of η_m for the m = 1 mode (black solid lines), the m = 4 mode (red solid lines), and other azimuthal modes m ≤ 10 (dashed lines) is shown for a short period (τ ≤ 10) in panel (b) and a long period (τ ≤ 30) in panel (c). Here ${\eta }_{{\rm{c}}}=\pi /(2\sqrt{2})$ is assumed.
Download figure:
Standard image High-resolution image

**Figure 5.** Time evolution of the stability of the axisymmetric components against the RWI in the "h10w3g1" runs. Panel (a) shows the time evolution of η₄ in the single-mode calculation for the m = 4 mode. The time evolution of η_m for the m = 1 mode (black solid lines), the m = 4 mode (red solid lines), and other azimuthal modes m ≤ 10 (dashed lines) is shown for a short period (τ ≤ 10) in panel (b) and a long period (τ ≤ 30) in panel (c). Here ${\eta }_{{\rm{c}}}=\pi /(2\sqrt{2})$ is assumed.
Download figure:
Standard image High-resolution image

We also observe the similar evolution of η_m in other calculations. However, the time when ${\eta }_{m={m}_{* }}$ becomes smaller than $\pi /(2\sqrt{2})$ deviates from that of the RWI saturation in the models with a small linear growth rate. We consider that the time deviation is due to ${\eta }_{{\rm{c}}}\ne \pi /(2\sqrt{2})$ because the profile of D_MS,m for r_IR < r < r_OR is not approximated very well by a parabolic function when the initial Gaussian bump is weak; i.e., the linear growth rate is small.

Meheut et al. (2013) interpreted the saturation mechanism of the RWI in an analogy of the wave–particle interaction in plasma physics. Our explanation is based on the linear stability of the axisymmetric components and is complementary to that by Meheut et al. (2013). We expect that combining these explanations helps us understand the physical mechanisms of the RWI evolution.

3.3. Vortex Merger

After the RWI saturation, multiple vortices formed as a result of the RWI coalesce one after another. In this section, we investigate the vortex merger regime.

In all the runs, the regimes with more than two vortices continue at most for a few orbits. On the other hand, the lifetime of the two-vortices regime shows some variations. By visual inspection of the surface density distribution, we identify the orbits when the vortex mergers occur. We define an orbit when the number of the vortices becomes two by τ₂ and an orbit when the final vortex merger occurs by τ₁. By that account, τ₂ − τ₁ represents the duration of the two-vortices regime. We show the values of τ₂ and τ₁ in Appendix C.2. Note that these orbits have errors of a few tenths due to the uncertainties of our visual inspection. For all the models with m_* = 1 and some models with m_* = 2, it is difficult to measure τ₂ and τ₁, so we set τ₂ to no data and τ₁ to the orbit number at the RWI saturation.

The values of (τ₁ − τ₂) seem to be random. As shown in Figure 6, however, there is an upper limit in (τ₁ − τ₂),

$\begin{eqnarray}{\tau }_{1}-{\tau }_{2} & \lesssim & 300\exp (-11.5h),\\ & \approx & 95\exp [1-h/0.1],\end{eqnarray} \tag{ 12 }$

within our parameter space. For r_n = 100 au, one orbit corresponds to about 10³ yr. From Equation (12), the lifetime of the two-vortices regime is up to about a few ×0.1 Myr. The duration of the two-vortices regime is one to two orders of magnitude shorter than the disk lifetime, which is 1–10 Myr (e.g., Haisch et al. 2001). It is difficult to observe protoplanetary disks with multiple vortices formed by the same RWI event except at outer disks.

**Figure 6.** Duration of the two-vortices regime. The values of τ₁ − τ₂ are shown with by green squares (h = 0.05), black circles (h = 0.1), red triangles (h = 0.15), and blue diamonds (h = 0.2). The black dashed line shows ${\tau }_{1}-{\tau }_{2}=300\exp (-11.5h)$ .
Download figure:
Standard image High-resolution image

**Figure 6.** Duration of the two-vortices regime. The values of τ₁ − τ₂ are shown with by green squares (h = 0.05), black circles (h = 0.1), red triangles (h = 0.15), and blue diamonds (h = 0.2). The black dashed line shows ${\tau }_{1}-{\tau }_{2}=300\exp (-11.5h)$ .
Download figure:
Standard image High-resolution image

The vortex mergers strongly depend on the perturbations imposed on the initial conditions. In our calculations, the white-noise perturbations always have the same Fourier phase and power spectrum because we use the same random seed and set the maximum azimuthal mode to m = 128. When the Fourier phase or power spectrum of the perturbations is different, the time when the vortex mergers occur varies. Even in those cases, however, τ₂ − τ₁ always satisfies Equation (12).

We turn our attention to the stability of the axisymmetric components during the vortex merger regime. As discussed in Section 3.2, the axisymmetric components are still unstable at the RWI saturation for the lower modes than the most unstable azimuthal mode. After the saturation, r_p and ${ \mathcal A }$ continue to decrease and Δw increases. The rate of change of r_p is similar to that before saturation, but the rates of change of ${ \mathcal A }$ and Δw are slower. From panel (c) of Figure 5, the values of η_m for 1 ≤ m ≤ 3 continue to decrease during the vortex mergers and finally become below the threshold. Particularly, η₁ reaches the threshold just after the final vortex merger at τ ≈ 20. Therefore, the axisymmetric components evolve toward the stable configurations during the vortex merger regime and become marginally stable against the RWI for the m = 1 mode at the final vortex merger.

4. Quasi-stationary Vortex Formed by the RWI

We call the quasi-stationary vortex formed after the final vortex merger the "RWI vortex." Hereafter, we focus on the structure and evolution of the RWI vortex. We provide a method to analyze the RWI vortex in Section 4.1 and show results in Section 4.2. Section 4.3 is for discussion.

4.1. Analyses

4.1.1. Structure of the RWI Vortex

We show definitions of some physical quantities that characterize the vortex structure (vortex center, velocity gradient, vortex size, vortex aspect ratio, and turnover time).

We define the center of the RWI vortex (r_v, φ_v) by

$\begin{eqnarray}&&\delta {v}_{\varphi }({r}_{{\rm{v}}},{\varphi }_{{\rm{v}}})=0,\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&{v}_{r}({r}_{{\rm{v}}},{\varphi }_{{\rm{v}}})=0,\end{eqnarray} \tag{ 14 }$

in the vicinity of the surface density peak, where $\delta {v}_{\varphi }(r,\varphi )\equiv {v}_{\varphi }(r,\varphi )-{v}_{{\rm{K}}}(r)$ . From panels (a) and (b) of Figure 7, the surface density at the vortex center, ${{\rm{\Sigma }}}_{{\rm{v}}}\equiv {\rm{\Sigma }}({r}_{{\rm{v}}},{\varphi }_{{\rm{v}}})$ , is almost the same as the peak value of the surface density.

**Figure 7.** Structure of the RWI vortex in the "h10w3g1" run at τ = 100. Panels (a)–(d) show the profiles of ${\rm{\Sigma }}(r,{\varphi }_{{\rm{v}}})$ , ${\rm{\Sigma }}({r}_{{\rm{v}}},\varphi )$ , $\delta {v}_{\varphi }(r,{\varphi }_{{\rm{v}}})$ , and ${v}_{r}({r}_{{\rm{v}}},\varphi )$ , respectively. The red dashed line in panels (c) and (d) shows the radial (azimuthal) gradient of δv_φ (v_r) at the vortex center. The black points show the values at the extrema expected from the velocity gradients.
Download figure:
Standard image High-resolution image

Next, we turn our attention to the velocity field in the vortex and the vortex size. Panels (c) and (d) of Figure 7 show the radial profile of δv_φ and the azimuthal profile of v_r, respectively. In the vicinity of the vortex center, these velocity profiles are almost on straight lines. We define the radial and azimuthal velocity gradients at the vortex center by

$\begin{eqnarray*}\delta {v}_{\varphi {\rm{v}},r} & = & {\left[\displaystyle \frac{\partial \delta {v}_{\varphi }(r,{\varphi }_{{\rm{v}}})}{\partial r}\right]}_{r={r}_{{\rm{v}}}},\\ {v}_{r{\rm{v}},\varphi } & = & \displaystyle \frac{1}{{r}_{{\rm{v}}}}{\left[\displaystyle \frac{\partial {v}_{r}({r}_{{\rm{v}}},\varphi )}{\partial \varphi }\right]}_{\varphi ={\varphi }_{{\rm{v}}}}.\end{eqnarray*}$

In addition, we define the radial and azimuthal convexities of the pressure function at the vortex center by

$\begin{eqnarray*}{{\rm{\Pi }}}_{{\rm{v}},{rr}} & = & {\left[\displaystyle \frac{{\partial }^{2}{\rm{\Pi }}(r,{\varphi }_{{\rm{v}}})}{\partial {r}^{2}}\right]}_{r={r}_{{\rm{v}}}},\\ {{\rm{\Pi }}}_{{\rm{v}},\varphi \varphi } & = & \displaystyle \frac{1}{{r}_{{\rm{v}}}^{2}}{\left[\displaystyle \frac{{\partial }^{2}{\rm{\Pi }}({r}_{{\rm{v}}},\varphi )}{\partial {\varphi }^{2}}\right]}_{\varphi ={\varphi }_{{\rm{v}}}}.\end{eqnarray*}$

The velocity gradients and the convexities of the pressure function are used to compare the RWI vortices with the analytic solutions of steady vortices (see Section 4.2.1). The velocity profiles gradually deviate from the straight lines with distance from the vortex center and finally have two extrema. At these extrema, the values of $| \delta {v}_{\varphi }|$ and $| {v}_{r}|$ are about two-thirds times as large as $| (r-{r}_{{\rm{v}}})\delta {v}_{\varphi {\rm{v}},r}|$ and $| {r}_{{\rm{v}}}(\varphi -{\varphi }_{{\rm{v}}}){v}_{r{\rm{v}},\varphi }|$ , respectively. We define the radial and azimuthal half-widths Δr and r_vΔφ by the half of the distance between the two extrema of δv_φ and v_r.

The vortex aspect ratio and turnover time in the vortex are important physical quantities of the vortex (e.g., Kida 1981). In this paper, we measure these quantities using streamlines. As shown in panel (a) of Figure 8, the streamlines around the vortex center look like closed loops, indicating that the flow is in a quasi-stationary state. The streamlines are almost elliptic in the (r – r_v) – r_v(φ – φ_v) plane, and the semiminor axes of them are aligned to the radial direction. We measure the semiminor axis b (the radial direction) and the semimajor axis a (the azimuthal direction) for each streamline. We define the aspect ratio of each streamline by χ ≡ a/b. We also measure the turnover time of each streamline normalized by 2π/Ω_v, ψ,

$\begin{eqnarray}&&\psi \equiv \displaystyle \frac{{{\rm{\Omega }}}_{{\rm{v}}}}{2\pi }{\oint }_{{ \mathcal L }}\displaystyle \frac{d{\ell }}{\tilde{v}},\end{eqnarray} \tag{ 15 }$

where ${{\rm{\Omega }}}_{{\rm{v}}}\equiv {\rm{\Omega }}({r}_{{\rm{v}}})={{\rm{\Omega }}}_{{\rm{K}}}({r}_{{\rm{v}}})$ , ${ \mathcal L }$ denotes the integration along the streamline, $d{\ell }$ is the line element along the streamline, and $\tilde{v}=\sqrt{{v}_{r}^{2}+{({v}_{\varphi }-r{{\rm{\Omega }}}_{{\rm{v}}})}^{2}}$ is the magnitude of the velocity field in the rotating frame with the vortex center. Panels (b) and (c) of Figure 8 show the profiles of χ and ψ as functions of the normalized distance from the vortex center, b/r_v. Both χ and ψ are almost constant around the vortex center. At b ≥ 0.1 r_v, these quantities are no longer constant and increase rapidly. Here, 0.1 r_v is close to one disk scale height at the vortex center. We define χ₂ and ψ₂ by χ and ψ at b/r_v = 0.02 as representative vortex aspect ratio and normalized turnover time in the vortex, respectively.

**Figure 8.** Panel (a) shows streamlines around the vortex center in the "h10w3g1" run at τ = 100. Panels (b) and (c) show the profiles of χ and ψ as functions of b/r_v, respectively.
Download figure:
Standard image High-resolution image

4.1.2. Measurements about the RWI Vortex

We terminate all our main calculations after the disks have had a few hundred orbits at r = r_n since the RWI vortex formation. Here we provide a way to measure the properties of the RWI vortices in our calculations.

From panels (a) and (b) of Figure 9, the RWI vortices migrate toward the central star, and their surface densities increase. At that time, the vortices keep ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ almost constant during the vortex migration, as seen in panel (c) of Figure 9. From panels (d)–(g) of Figure 9, the vortex aspect ratio χ₂, vortex turnover time ψ₂, vortensity at the vortex center ${q}_{{\rm{v}}}\equiv [\mathrm{rot}{\boldsymbol{v}}]({r}_{{\rm{v}}},{\varphi }_{{\rm{v}}})/{{\rm{\Sigma }}}_{{\rm{v}}}$ , radial half-width of the vortex Δr, and azimuthal half-width of the vortex r_vΔφ are approximately constant. To investigate the migration speed of the RWI vortex, we define a physical quantity ξ by

$\begin{eqnarray}&&\xi \equiv -\displaystyle \frac{1}{{{\rm{\Omega }}}_{{\rm{v}}}}\displaystyle \frac{{{dr}}_{{\rm{v}}}}{{dt}}=\displaystyle \frac{-1}{2.5\sqrt{{GM}}}\displaystyle \frac{{{dr}}_{{\rm{v}}}^{2.5}}{{dt}},\end{eqnarray} \tag{ 16 }$

where ${{\rm{\Omega }}}_{{\rm{v}}}=\sqrt{{GM}/{r}_{{\rm{v}}}^{3}}$ . The value of ξ shows the distance the vortex moves in the −r direction for a unit time. We define the timescale of the vortex migration τ_mig by ${\tau }_{\mathrm{mig}}\equiv {{\rm{\Omega }}}_{{\rm{n}}}{r}_{{\rm{v}}}/(2\pi \xi {{\rm{\Omega }}}_{{\rm{v}}})$ . Panel (h) of Figure 9 shows that ξ is also almost constant. In this paper, we use the time-averaged values of these physical quantities over 40 orbits in each run as the measurements of the RWI vortex. We summarize the physical quantities and measurements about the RWI vortex in Table 2.

**Figure 9.** Time evolution of the RWI vortex in the "h10w3g1" run. Each panel shows the time evolution of (a) r_v, (b) Σ_v, (c) Σ_vr_v^1.5, (d) χ₂ (red) and ψ₂ (blue), (e) q_v, (f) Δr, (g) r_vΔφ, and (h) ξ. The crosses represent the time-averaged values over 40 orbits at τ = 120.
Download figure:
Standard image High-resolution image

Table 2. List of Physical Quantities and Measurements about the RWI Vortex

Measurement	Meaning
r_v	The distance of the vortex center from the central star.
φ_v	The azimuthal angle at the vortex center.
${{\rm{\Omega }}}_{{\rm{v}}}\equiv \sqrt{{GM}/{r}_{{\rm{v}}}^{3}}$	The angular velocity at the vortex center.
Σ_v	The surface density at the vortex center.
δv_φ	The relative rotation velocity from the Keplerian velocity at the vortex center.
δv_φv,r	The radial gradient of δv_φ at the vortex center.
v_rv,φ	The azimuthal gradient of v_r at the vortex center.
Π_v,rr	The radial convexity of the pressure function at the vortex center.
Π_v,φφ	The azimuthal convexity of the pressure function at the vortex center.
Δr	The radial half-width of the vortex.
r_vΔφ	The azimuthal half-width of the vortex.
χ₂	The vortex aspect ratio in the vicinity of the vortex center.
ψ₂	The turnover time normalized in the vicinity of the vortex center.
q_v	The vortensity at the vortex center.
ξ	The distance the vortex moves in the −r direction for a unit time.
τ_mig	The vortex migration timescale.
χ_2,ini	The value of χ₂ just after the RWI vortex formation.
τ_χ	The decrease timescale of the vortex aspect ratio.

Download table as: ASCII Typeset image

The RWI vortex is quasi-stationary but not completely stationary. In a longer timescale than 1000 orbits, the values of ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ , q_v, and Δr are still almost constant, but the values of χ₂, ψ₂, r_vΔφ, and ξ are not. We will discuss the applicability of our results to the long-term evolution in Section 4.3.2. However, we speculate that the long-term evolution is due to numerical viscosity in our calculations. Since viscosity seems to be generally important for the long-term evolution of the RWI vortices, a detailed investigation of the long-term evolution falls outside the scope of this paper.

4.1.3. Tracer Particle Analysis

In order to obtain more detailed information on the RWI vortex, we perform a tracer particle analysis. We calculate path lines of fluid particles for 30 orbits in the "h10w3g1" run. Some 240 × 360 tracer particles are initially distributed uniformly within the range of 0.6 ≤ r ≤ 1.2 and −π ≤ φ ≤ π at τ = 55, 60, 65, 70, 75, 80, and 85.

4.2. Results

4.2.1. Comparison with the Analytic Steady Vortices

We compare the RWI vortices with known analytic solutions of steady vortices (the Kida and GNG solutions; see Appendix A for details). The steady vortices satisfy

$\begin{eqnarray}&&\delta {v}_{\varphi {\rm{v}},r}=(1.5-{\chi }_{2}/{\psi }_{2}){{\rm{\Omega }}}_{{\rm{v}}},\end{eqnarray} \tag{ 17 }$

$\begin{eqnarray}&&{v}_{r{\rm{v}},\varphi }={{\rm{\Omega }}}_{{\rm{v}}}/({\chi }_{2}{\psi }_{2}),\end{eqnarray} \tag{ 18 }$

$\begin{eqnarray}&&{{\rm{\Pi }}}_{{\rm{v}},{rr}}=(1/{\psi }_{2}^{2}-2{\chi }_{2}/{\psi }_{2}+3.0){{\rm{\Omega }}}_{{\rm{v}}}^{2},\end{eqnarray} \tag{ 19 }$

and

$\begin{eqnarray}&&{{\rm{\Pi }}}_{{\rm{v}},\varphi \varphi }=[1/{\psi }_{2}^{2}-2/({\psi }_{2}{\chi }_{2})]{{\rm{\Omega }}}_{{\rm{v}}}^{2}.\end{eqnarray} \tag{ 20 }$

Equations (17)–(20) correspond to Equations (38), (39), (41), and (42). In the case of the Kida solution, the vortex aspect ratio χ and the vortex turnover time ψ for the steady vortices are related as

$\begin{eqnarray}&&(\chi -1)={\mathfrak{S}}\psi .\end{eqnarray} \tag{ 21 }$

The GNG solution (Goodman et al. 1987) gives another relation:

$\begin{eqnarray}&&({\chi }^{2}-1)=2{\mathfrak{S}}{\psi }^{2}.\end{eqnarray} \tag{ 22 }$

As shown in panel (a) of Figure 10, the χ₂ and ψ₂ of the RWI vortices reasonably satisfy Equations (21) and (22). Panels (b)–(e) of Figure 10 show that the RWI vortices satisfy Equations (17)–(20) within a factor of 1.5. Therefore, the RWI vortices resemble the steady vortices in the Keplerian shear, as shown in previous works (e.g., Surville & Barge 2015).

**Figure 10.** Comparing the RWI vortices with the analytic steady solutions of vortices. The time-averaged values of (a) ψ₂, (b) $\delta {v}_{\varphi {\rm{v}},r}{{\rm{\Omega }}}^{-1}/(1.5-{\chi }_{2}/{\psi }_{2})$ , (c) ${v}_{r{\rm{v}},\varphi }{{\rm{\Omega }}}_{{\rm{v}}}^{-1}/{({\chi }_{2}{\psi }_{2})}^{-1}$ , (d) ${{\rm{\Pi }}}_{{\rm{v}},{rr}}{{\rm{\Omega }}}_{{\rm{v}}}^{-2}/(1/{\psi }_{2}^{2}-2{\chi }_{2}/{\psi }_{2}+3.0)$ , and (e) ${{\rm{\Pi }}}_{{\rm{v}},\varphi \varphi }{{\rm{\Omega }}}_{{\rm{v}}}^{-2}/(1/{\psi }_{2}^{2}-2/({\psi }_{2}{\chi }_{2}))$ against the time-averaged χ₂ are shown with green squares (h = 0.05), black plus signs (h = 0.1), red circles (h = 0.15), and blue triangles (h = 0.2). The Kida vortex lies on the solid line and the GNG solution on the dashed line in panel (a).
Download figure:
Standard image High-resolution image

4.2.2. Empirical Relations of the RWI Vortex

We obtain some empirical relations between the RWI vortices and the initial bump structures so that they can help us predict the nonlinear outcomes from initial structures.

First, we consider the surface density at the vortex center, Σ_v. Compiling all the calculations with the different initial conditions, we find the empirical relations between ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ and the parameters of the initial bumps as (see panel (a) of Figure 11)

$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}\approx \left[1+\displaystyle \frac{{({{ \mathcal A }}_{0}/h)}^{0.8}}{4}\right]{{\rm{\Sigma }}}_{{\rm{n}}}{r}_{{\rm{n}}}^{1.5}.\end{eqnarray} \tag{ 23 }$

**Figure 11.** The values of (a) ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ against ${{ \mathcal A }}_{0}/h$ , (b) χ_2,ini and (c) Δr against γ_*/γ_*,max, and (d) ξ/(hΔr) against χ₂ are shown with green squares (h = 0.05), black plus signs (h = 0.1), red circles (h = 0.15), and blue triangles (h = 0.2). In panel (a), the solid line shows ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}=[1+{({{ \mathcal A }}_{0}/h)}^{0.8}/4]{{\rm{\Sigma }}}_{{\rm{n}}}{r}_{{\rm{n}}}^{1.5}$ . In panel (b), the solid line shows ${\chi }_{2}=5{({\gamma }_{* }/{\gamma }_{* ,\max })}^{-1}$ , and the dashed lines show ${\chi }_{2}=4{({\gamma }_{* }/{\gamma }_{* ,\max })}^{-1}$ and ${\chi }_{2}=6{({\gamma }_{* }/{\gamma }_{* ,\max })}^{-1}$ . In panel (c), the solid line shows ${\rm{\Delta }}r=1.46{h}^{1/5}{({\gamma }_{* }/{\gamma }_{* ,\max })}^{1/3}{\rm{\Delta }}{w}_{0}^{2/3}{r}_{{\rm{n}}}$ . In panel (d), the solid line shows $\xi =1.6h{\chi }_{2}^{-3}{\rm{\Delta }}r$ , and the dashed lines are twice and half of the solid line. The pentagons show the values derived from the results of Paardekooper et al. (2010; cyan) and Richard et al. (2013; orange).
Download figure:
Standard image High-resolution image

Next, we consider the vortex aspect ratio χ₂. Qualitatively, the vortex aspect ratio χ₂ is small for the small growth rate cases, and vice versa. As shown in panel (b) of Figure 11, we find that χ₂ and γ_*/γ_*,max are related by

$\begin{eqnarray}&&{\chi }_{2}\approx {5}^{\pm 1}{({\gamma }_{* }/{\gamma }_{* ,\max })}^{-1}.\end{eqnarray} \tag{ 24 }$

Using Equation (24), we can estimate the aspect ratio of the RWI vortex from the initial bump structures. Once we know the aspect ratio, it is possible to estimate the vortex turnover time and the velocity gradients and convexities of the pressure function at the vortex center using the analytic vortex solutions.

We turn our attention to the radial vortex size. We find that Δr is related to the initial bump parameters by (see panel (b) of Figure 11)

$\begin{eqnarray}&&{\rm{\Delta }}r\approx 1.46{h}^{1/5}{({\gamma }_{* }/{\gamma }_{* ,\max })}^{1/3}{\rm{\Delta }}{w}_{0}^{2/3}{r}_{{\rm{n}}}^{1/3}.\end{eqnarray} \tag{ 25 }$

We also compare Δr with the disk scale height at the vortex center ${H}_{{\rm{v}}}\equiv h{({{\rm{\Sigma }}}_{{\rm{v}}}/{{\rm{\Sigma }}}_{{\rm{n}}})}^{({\rm{\Gamma }}-1)/2}{({r}_{{\rm{v}}}/{r}_{{\rm{n}}})}^{1.5}{r}_{{\rm{n}}}$ . Figure 12 indicates that a maximum value of Δr/H_v is about 2. This maximum value of the radial vortex size is consistent with the values reported by previous works (e.g., Li et al. 2001; Surville & Barge 2015).

Finally, we consider the vortex migration speed. From panel (d) of Figure 11, we find that ξ satisfies

$\begin{eqnarray}&&\xi =1.6h{\chi }_{2}^{-3}{\rm{\Delta }}r\end{eqnarray} \tag{ 26 }$

within a factor of about 2. Using Equation (26), the timescale of the vortex migration ${\tau }_{\mathrm{mig}}$ is obtained as

$\begin{eqnarray}&&{\tau }_{\mathrm{mig}}\approx 2.1\times {10}^{3}{\left(\displaystyle \frac{h}{0.1}\right)}^{-1}{\left(\displaystyle \frac{{\chi }_{2}}{6}\right)}^{3}{\left(\displaystyle \frac{{\rm{\Delta }}r}{0.1{r}_{{\rm{n}}}}\right)}^{-1}{\left(\displaystyle \frac{{r}_{{\rm{v}}}}{{r}_{{\rm{n}}}}\right)}^{2.5}.\end{eqnarray} \tag{ 27 }$

If r_n = 100 au and M = M_⊙, an orbital time at r = r_n is about 1000 yr, so that the timescale of the vortex migration is

$\begin{eqnarray}&&2.1\ \mathrm{Myr}\times {\left(\displaystyle \frac{h}{0.1}\right)}^{-1}{\left(\displaystyle \frac{{\chi }_{2}}{6}\right)}^{3}{\left(\displaystyle \frac{{\rm{\Delta }}r}{10\mathrm{au}}\right)}^{-1}{\left(\displaystyle \frac{{r}_{{\rm{v}}}}{100\mathrm{au}}\right)}^{2.5}.\end{eqnarray} \tag{ 28 }$

Unless the vortex aspect ratio is small ( ${\chi }_{2}\sim 4$ ) and the vortex size is large (Δr ∼ 2H_v), the vortex migration timescale is comparable to or longer than the lifetime of the protoplanetary disks (∼1–10 Myr). In other words, the RWI vortex resulting from a narrow and/or weak initial surface density bump stays within the disk without suffering from the radial migration.

Using Equations (23)–(26), we can estimate the properties of the RWI vortex (the surface density at the vortex center, the aspect ratio, the radial size, and the migration speed) from the initial conditions (h, Δw₀, and γ_*/γ_*,max). We note that ${{ \mathcal A }}_{0}$ and γ_*,max can be calculated from h, Δw₀, and γ_* performing the linear stability analyses.

4.2.3. The RWI Vortex from the Point of View of Tracer Particles

We investigate the structure and evolution of the RWI vortex from the point of view of the tracer particles. We focus on the tracer particles initially distributed at τ = 70 in the "h10w3g1" run. According to the r evolution of those fluid particles (see panel (a) of Figure 13), we categorize fluid particles into four groups (I, II, III, and IV). Group I represents the particles that compose the vortex at τ = 70. All the group I particles remain in the vortex part during the 30 orbits. The fluid particles categorized into group II are in the inner part of the disk initially and move to the outer part of the disk within the 30 orbits. The group III and IV particles remain in the inner and outer parts for the 30 orbits, respectively. The vortex captures a few particles that reside in the inner region, but such particles escape from the vortex to the outer part after only a few turnover motions and are categorized into group II. We find that no particle moves inward (the outer part $\to$ the inner part, the outer part $\to$ the vortex part, and the vortex part $\to$ the inner part) during the 30 orbits. We show the distribution of the particle groups at τ = 70 in panel (b) of Figure 13. The tracer particles in the other calculations are also categorized in the same way.

**Figure 13.** Categorization of the fluid particles into four groups and distribution of the particle groups. The r evolution of (a) the group I particle (red solid line), (b) the group II particle (blue dotted line), (c) the group III particle (the green dot-dashed line), and (d) the group IV particle (orange double-dot-dashed line) that depart from the starting points at τ = 70 in the "h10w3g1" run is shown in panel (a). The black dashed line shows the evolution of the vortex center. Panel (b) shows the distribution of the particle groups at τ = 70 (red: group I; blue: group II; green: group III; and orange: group IV). The black lines show the streamlines in the vortex at τ = 70.
Download figure:
Standard image High-resolution image

We calculate the total mass of the group I particles, M_v, which represents the vortex mass. From Figure 14, M_v is approximately constant within 1% over 30 orbits in the "h10w3g1" run. We therefore expect that, in the course of the inward migration, the RWI vortex carries most of the fluid particles originating from the initial location where the vortex is formed, down to the inner radii. In this sense, the RWI vortex can be considered as a physical entity, like a large fluid particle.

4.3. Discussions

4.3.1. Definition of the Vortex Center

Our definition of the vortex center is based on the velocity field. The vortex center defined from the velocity field corresponds to the center of the streamlines in the vortex. Therefore, we call the vortex center obtained from our definition the "streamline vortex center." However, there is another way to obtain the vortex center using tracer particles. To focus on the time evolution of the tracer particles close to the vortex center, we use three particles that are located near the vortex center at τ = 60, 70, and 80 in the "h10w3g1" run. From panel (a) of Figure 15, those particles migrate inward with small oscillations. Here we can define the center of the oscillations by the "particle vortex center."

**Figure 15.** The r and Σ evolution of the fluid particles that depart from the starting points near the vortex center at τ = 60 (red line), 70 (blue line), and 80 (green line) are shown in panels (a) and (b), respectively. In panel (a), the black solid line shows the evolution of r_v, and the black dotted line shows the evolution of $1.022{r}_{{\rm{v}}}$ . In panel (b), the black solid line shows the evolution of Σ_v, and the black dotted line shows the evolution of 0.974Σ_v.
Download figure:
Standard image High-resolution image

**Figure 15.** The r and Σ evolution of the fluid particles that depart from the starting points near the vortex center at τ = 60 (red line), 70 (blue line), and 80 (green line) are shown in panels (a) and (b), respectively. In panel (a), the black solid line shows the evolution of r_v, and the black dotted line shows the evolution of $1.022{r}_{{\rm{v}}}$ . In panel (b), the black solid line shows the evolution of Σ_v, and the black dotted line shows the evolution of 0.974Σ_v.
Download figure:
Standard image High-resolution image

The distance of the particle vortex center from the central star is larger than r_v by about 2.2%. We consider that this discrepancy originates from the difference between the position of the surface density peak and that of the vortensity minimum. Since the dynamical equilibrium is achieved in the vicinity of the surface density peak, the flow should be stagnant there. Therefore, the streamline vortex center is very close to the surface density peak. On the other hand, the particle vortex center is very close to the vortensity minimum due to the vortensity conservation law.

However, the migration speed of both vortex centers almost corresponds, as shown in panel (a) of Figure 15. We can also confirm this fact in panel (b) of Figure 15. The time-averaged surface densities of the three particles are smaller than Σ_v by about 2.6%, but the growth rates of the surface densities of these tracer particles are almost equivalent to Σ_v. Therefore, the choice of the vortex center has no effect on our results.

4.3.2. Applicability of Empirical Formulae to the Long-term Evolution

In cases with a large linear growth rate, the RWI vortices migrate too fast to survive for 1000 orbits. On the other hand, the vortex migration is slow enough to survive for 1000 orbits in cases with a small linear growth rate. In such a long timescale, not all the physical quantities that are approximately constant in a short timescale are almost constant. In this section, we check the applicability of the empirical formulae obtained in Section 4.2.2 to the long-term evolution.

Here we regard the "h10w3g5" run ( $h=0.1,{\rm{\Delta }}{w}_{0}\,=0.0632{r}_{{\rm{n}}}$ , γ_*/Ω_n = 0.1, ${{ \mathcal A }}_{0}=0.142$ , and m_* = 2) as a representative case with a small linear growth rate. From Figure 16, ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ , q_v, and Δr are still almost constant, but χ₂, ψ₂, and r_vΔφ decrease and ξ increases in the long-term calculation of the "h10w3g5" model. Panel (d) of Figure 16 shows that the decrease of χ₂ is exponential-like following

$\begin{eqnarray}&&{\chi }_{2}={\chi }_{2,\mathrm{ini}}\exp (-\tau /{\tau }_{\chi }),\end{eqnarray} \tag{ 29 }$

where χ_2,ini ≈ 27.6 is the vortex aspect ratio just after the RWI vortex formation and τ_χ ≈ 5040 is the decrease timescale of the vortex aspect ratio. Panels (d) and (g) of Figure 16 show that ψ₂ and r_vΔφ also exponentially decrease on a similar timescale to τ_χ and a slightly longer timescale, respectively. From panel (h) of Figure 16, ξ exponentially increases on a similar timescale to 3τ_χ. We consider that τ_χ shows the timescale of the long-term evolution.

**Figure 16.** Similar to Figure 9 but the long-term calculation of the "h10w3g5" model. The crosses show the time-averaged values over 40 orbits every 100 orbits from τ = 360. The green dashed line in panel (d) shows the best fit of the exponential decreases of χ₂ with χ_2,ini = 27.6 and τ_χ = 5036. The orange dotted lines show ${\psi }_{2}=16.7\ \exp (-\tau /{\tau }_{\chi })$ in panel (d), ${r}_{{\rm{v}}}{\rm{\Delta }}\varphi =1.62{r}_{{\rm{n}}}\ \exp (-\tau /{\tau }_{\chi })$ in panel (g), and $\xi =4.85\times {10}^{-7}{r}_{{\rm{n}}}\ \exp (3\tau /{\tau }_{\chi })$ in panel (h).
Download figure:
Standard image High-resolution image

**Figure 16.** Similar to Figure 9 but the long-term calculation of the "h10w3g5" model. The crosses show the time-averaged values over 40 orbits every 100 orbits from τ = 360. The green dashed line in panel (d) shows the best fit of the exponential decreases of χ₂ with χ_2,ini = 27.6 and τ_χ = 5036. The orange dotted lines show ${\psi }_{2}=16.7\ \exp (-\tau /{\tau }_{\chi })$ in panel (d), ${r}_{{\rm{v}}}{\rm{\Delta }}\varphi =1.62{r}_{{\rm{n}}}\ \exp (-\tau /{\tau }_{\chi })$ in panel (g), and $\xi =4.85\times {10}^{-7}{r}_{{\rm{n}}}\ \exp (3\tau /{\tau }_{\chi })$ in panel (h).
Download figure:
Standard image High-resolution image

We calculate the time-averaged measurements of the RWI vortex over 40 orbits every 100 orbits. Since ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ and Δr are still approximately constant, Equations (23) and (25) are satisfied. We also confirm that the RWI vortex always resembles the steady vortices and satisfies Equation (26) even though χ₂, ψ₂, and ξ vary. However, the value of χ₂ goes away from the value obtained by Equation (24). We observe these trends in the long-term calculations of the other small growth rate cases. Therefore, we conclude that the empirical formulae except for Equation (24) are applicable even to the long-term evolution. On the other hand, Equation (24) is applicable only for the first few hundred orbits after the RWI vortex formation.

Since Δr is approximately constant, the shrink of the RWI vortex in the azimuthal direction can explain the long-term evolution of the vortex aspect ratio. In all the models without a large linear growth rate, τ_χ is longer than 1000. For the cases with a large linear growth rate, it is impossible to measure τ_χ due to the fast vortex migration. However, τ_χ is expected to also be large compared to the migration timescale because the structure of the RWI vortex is almost stationary in the migration timescale. The values of τ_χ are shown in Appendix C.2. We also perform a long-term calculation of the "h10w3g5" model with a twice-coarser resolution than that of the fiducial setup. We find that τ_χ is several times smaller in the coarse calculation than in the fiducial calculation. This indicates that the numerical viscosity has a significant effect on the long-term evolution and the viscosity is important for the long-term evolution in general. Therefore, a detailed investigation of the long-term evolution is outside the scope of this paper.

4.3.3. Reason Why ${{\rm{\Sigma }}}_{v}{{\rm{r}}}_{v}^{1.5}$ Is Approximately Constant

As shown in Sections 4.1.2 and 4.3.2, ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ is approximately constant, as well as q_v in short and long timescales. In our calculations, the vortensity conservation law should be satisfied without taking the numerical viscosity into account. Almost the same tracer particles constitute the vortex center (see Section 4.2.3) so that the vortensity at the vortex center is approximately constant. In this section, we discuss the reason for the invariance of ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ using the vortensity conservation.

Here we assume that the RWI vortex coincides with the Kida vortex. From the Kida solution, we obtain

$\begin{eqnarray}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{v}}}}{{{\rm{\Omega }}}_{{\rm{v}}}}{q}_{{\rm{v}}} & = & 2.0-\displaystyle \frac{{\chi }_{2}}{{\psi }_{2}}-\displaystyle \frac{1}{{\chi }_{2}{\psi }_{2}},\\ & = & \displaystyle \frac{1}{2}-\displaystyle \frac{3}{2}\left(\displaystyle \frac{2}{{\chi }_{2}-1}-\displaystyle \frac{1}{{\chi }_{2}}\right).\end{eqnarray} \tag{ 30 }$

For the large growth rate cases, the second term of Equation (30) is comparable to the first term because χ₂ is small. In these cases, the vortex migration is fast and τ_mig < τ_χ. Therefore, we can safely assume that χ₂ is constant over the vortex lifetime and Σ_vq_v/Ω_v is constant. On the other hand, if linear growth rates are small, the vortex aspect ratio χ₂ is large so that ${{\rm{\Sigma }}}_{{\rm{v}}}{q}_{{\rm{v}}}/{{\rm{\Omega }}}_{{\rm{v}}}$ is approximately constant at 1/2. Therefore, the values of ${{\rm{\Sigma }}}_{{\rm{v}}}{q}_{{\rm{v}}}/{{\rm{\Omega }}}_{{\rm{v}}}$ can be regarded as constant in all cases. Due to the vortensity conservation law and the definition of the vortex center, Σ_vq_v/Ω_v is proportional to ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ . This is the reason why ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ is approximately constant in short and long timescales.

4.3.4. Comparison with Another Formula of Vortex Aspect Ratio

In Section 4.2.2, we obtain the empirical formula of the vortex aspect ratio by Equation (24). Richard et al. (2013) derived another formula assuming that the vorticity of the non-Keplerian motion normalized by that of the background shear flow, ${\tilde{\omega }}_{z}$ , is steady at the peak of the initial bump and the vortex center. Here we compare the two formulae.

When the profile of the initial surface density is given as a Gaussian bump, the value of ${\tilde{\omega }}_{z}$ at r = r_n for the initial conditions is calculated as

$\begin{eqnarray}{\tilde{\omega }}_{z} & = & {\left[\displaystyle \frac{d\{r({v}_{\varphi 0}-{v}_{{\rm{K}}})\}}{{dr}}/\displaystyle \frac{\partial \{r({v}_{{\rm{K}}}-r{{\rm{\Omega }}}_{{\rm{n}}})\}}{\partial r}\right]}_{r={r}_{{\rm{n}}}}\\ & = & \displaystyle \frac{{h}^{2}{{ \mathcal A }}_{0}{(1+{{ \mathcal A }}_{0})}^{{\rm{\Gamma }}-2}}{3{({\rm{\Delta }}{w}_{0})}^{2}}{r}_{{\rm{n}}}^{2}.\end{eqnarray} \tag{ 31 }$

From the Kida solution, the RWI vortex should satisfy

$\begin{eqnarray}&&{\tilde{\omega }}_{z}=\displaystyle \frac{1+\chi }{\chi (\chi -1)}.\end{eqnarray} \tag{ 32 }$

Richard et al. (2013) estimated the vortex aspect ratio, which we denote by χ_r, from the balance between Equations (31) and (32).

From Figure 17, χ₂ matches with χ_r in cases with a large linear growth rate. On the other hand, χ₂ is larger than χ_r by a factor of a few in cases with a small linear growth rate. Therefore, we consider that the formula by Richard et al. (2013) is applicable and tested for the large growth rate cases and that our formula extends their work to the small growth rate cases. We note again that Equation (24) is applicable only for the first few hundred orbits after the RWI vortex formation (see Section 4.3.2).

4.3.5. Generality of the Empirical Formula of Vortex Migration Speed

Now we can estimate the migration speed of the RWI vortex using Equation (26). This empirical formula does not explicitly depend on the initial surface density profile. Here we investigate whether the formula is generally applicable to vortices on disks or not.

Paardekooper et al. (2010) reported the migration speed of the vortex that is formed by imposing a vorticity perturbation in their 2D disk for h = 0.1 at r ≈ r_n. The parameters of the vortex are (χ, Δr, ξ) = (2.5, 0.025 r_n, 5.0 × 10⁻⁴ r_n). Richard et al. (2013) also measured the migration speed of three RWI vortices in their 3D calculations for h ≈ 0.0662 at r ≈ r_n. The parameters of the three vortices are (χ, Δr, ξ) = (7.0, 0.0543 r_n, 3.1 × 10⁻⁵ r_n), (8.5, 0.0531 r_n, 2.0 × 10⁻⁵ r_n), and (14.0, 0.0413 r_n, 3.0 × 10⁻⁶ r_n). We note that the numerical setups in Paardekooper et al. (2010) and Richard et al. (2013) are somewhat different from our setup, where the initial surface density profiles have global radial gradients and the disks are assumed to be locally isothermal. In panel (d) of Figure 11, we plot the result of Paardekooper et al. (2010) with a cyan pentagon and the results of Richard et al. (2013) with orange pentagons. We find that these points are almost on the line of $\xi =3.2h{\chi }_{2}^{-3}{\rm{\Delta }}r$ .

Paardekooper et al. (2010) showed that the inward migration of vortices is faster if the global radial gradient of the surface density is negative and steeper. In our calculations, we assume that the radial profile of the surface density is globally flat. On the other hand, both Paardekooper et al. (2010) and Richard et al. (2013) assumed that the radial profile is globally proportional to r^−1.5 so that ξ is twice larger. In that sense, our results are qualitatively consistent with those of Paardekooper et al. (2010). In addition, Equation (26) is satisfied within a factor of two even in the calculations with such a surface density slope. This indicates that the dependence of the migration speed on the vortex structure and disk aspect ratio seems to be universal.

According to Paardekooper et al. (2010), a pressure bump can prevent the vortex from migrating inward. In their calculations, the pressure bump is stronger and wider than the vortex. The axisymmetric components also have a pressure bump, as shown in Figure 4, in our calculations. However, it seems that the migration of the RWI vortices does occur even in the presence of the pressure bump. We consider that the pressure bump structures seen in our calculations are too weak and narrow to stop the vortex migration. One notable difference between Paardekooper et al. (2010) and our work is that the structure of the pressure bump is determined consistently with the development of the RWI and the formation of the vortex. On the other hand, Paardekooper et al. (2010) used the parameterized model for the pressure bump without calculating its formation process.

In short, we consider that Equation (26) is broadly applicable to estimate the vortex migration speed regardless of the formation mechanism of the vortex in both 2D and 3D disks unless the disks have very steep global slopes of the surface density or strong and wide pressure maxima. Further investigations are necessary to quantitatively study the effects of the global gradients of the surface density and the pressure maxima on the vortex migration.

4.3.6. Mechanism of Vortex Migration

In this section, we discuss the mechanism of the vortex migration. During the vortex migration, the vortex loses angular momentum via density waves (Paardekooper et al. 2010). The velocity perturbations induced by the vortex motion excite density waves in a disk, which carry away negative (inner spirals) or positive (outer spirals) angular momentum, causing the vortex to migrate. The positions of the Lindblad resonances for the m mode are located where the epicyclic frequency, κ, is the m times of the angular velocity of the fluid element in the frame corotating with the vortex (see panel (a) of Figure 18):

$\begin{eqnarray}&&\kappa (r,{\varphi }_{{\rm{v}}})=\pm m[{{\rm{\Omega }}}_{{\rm{v}}}-{v}_{\varphi }(r,{\varphi }_{{\rm{v}}})/r].\end{eqnarray} \tag{ 33 }$

As can be seen in panels (b)–(e) of Figure 18, the density waves are indeed excited around the Lindblad resonances for each mode.

**Figure 18.** Lindblad resonances and surface density distribution for each m mode. Panel (a) shows the epicyclic frequency κ in the "h10w3g1" run at τ = 100. In panel (a), $m({{\rm{\Omega }}}_{{\rm{v}}}-{v}_{\varphi }/r)$ is shown with the red solid line (m = 1), blue dotted line (m = 2), green dot-dashed line (m = 3), and orange double-dot-dashed line (m = 4), and Ω_K is shown with the black dashed line. Panels (b)–(e) show the 2D distribution of Σ_m for m = 1–4, respectively. The points in panel (a) and the vertical dashed lines in panels (b)–(e) indicate the Lindblad resonances of each m mode.
Download figure:
Standard image High-resolution image

The fluid particles of group II (see Section 4.2.3) also contribute to the angular momentum exchange between the vortex and the disk. When they move from the inner part to the outer part, they gain the angular momentum from the vortex. It is analogous to the corotation torque exerted on a planet in the planet–disk interaction (Balmforth & Korycansky 2001). We calculate the variation of the total angular momentum of the group I particles, ΔJ₁, and that of the group II particles, ΔJ₂, over the 30 orbits. We find $| {\rm{\Delta }}{J}_{2}/{\rm{\Delta }}{J}_{1}| \ \sim \ 0.2$ , indicating that about 20% of the total torque exerted on the vortex comes from the contribution of the fluid elements passing through the vortex region. The remaining about 80% of the total torque is interpreted to originate from the density waves. However, we find it difficult to directly measure the total torque that comes from the density waves, because it is difficult to distinguish the density waves from the vortex motion and to precisely measure the angular momentum flux of the density waves due to numerical viscosity. In order to quantify the significance of each mechanism in more detail, we need to precisely measure both types of torque.

4.3.7. Other Physical Effects

In this paper, we consider the simplest disk model to keep the broad parameter search tractable. Various potentially important physical effects on the RWI, such as viscosity, 3D, self-gravity, and dust drag, are not included. Here we briefly discuss how these effects can affect the RWI vortex.

Viscosity has a significant effect on the lifetime of the RWI vortex (Lin 2014). For long-term survival of the RWI vortex, very low viscosity (α ≲ 10⁻⁴) is required, where α is the kinematic viscosity normalized by the sound speed and Kepler time (the α-parameter; Shakura & Sunyaev 1973). We expect that the circumstance with such low α is achieved within the magneto-rotational instability (MRI) dead zone (Gammie 1996). As discussed in Section 4.3.2, we speculate that viscosity is still important for the long-term evolution of the RWI vortex even for very low viscosity. The long-term behavior of the vortex needs further investigation with an explicit viscosity prescription.

Lesur & Papaloizou (2009) reported that the effect of 3D can destroy vortices by the ellipsoidal instability. The analytic steady vortices are strongly unstable for χ ≤ 4 and weakly unstable for χ ≳ 6 from the 3D linear stability analyses and the local numerical simulations in incompressible flow with Keplerian shear. For χ ≲ 4, the velocity field of the vortex violates Rayleigh's condition, and the vortex is destroyed. On the other hand, the vortex with χ ≳ 6 is destroyed due to the resonance between the turnover motion of the vortex and the epicyclic motion of the disk. Richard et al. (2013) performed 3D compressible simulations of vortices formed by the RWI. In those simulations, the destruction of the vortices with χ ≲ 4 is verified. However, the vortices with χ ≈ 7 are not destroyed and survive for a sufficiently long time. Since our calculations are within the 2D framework, the RWI vortex does not suffer from the ellipsoidal instability. In all our calculations except for the "h05w4g1" and "h05w5g1" runs, χ₂ is always larger than 4. Even in the "h05w4g1" and "h05w5g1" runs, the vortex aspect ratio is approximately equal to 4. This is because we only consider the initial conditions that do not violate Rayleigh's criterion. In fact, we confirm that χ₂ is smaller than 4 if the initial conditions violate Rayleigh's condition, but we consider that such initial conditions are not realistic. Therefore, we expect that the evolution of the RWI vortices formed in our calculations almost never changes, even in the 3D frameworks.

The effect of self-gravity prevents the onset of the RWI (Lovelace & Hohlfeld 2013; Yellin-Bergovoy et al. 2016). However, once the RWI vortex is formed, there are possibilities that self-gravity can help the vortex survive for a long time (Lin & Pierens 2018). At that time, self-gravity can also have effects on the RWI evolution.

The effects of the dust particles on the gas flow are negligible in typical protoplanetary disks due to the low dust-to-gas mass ratio (∼10⁻²). When a protoplanetary disk has a gas vortex, the vortex captures the dust particles and concentrates them at the vortex center (Barge & Sommeria 1995). If a sufficiently high dust-to-gas mass ratio is realized at the vortex center, the vortex is destroyed due to the gas–dust interaction (Fu et al. 2014b; Crnkovic-Rubsamen et al. 2015). The motion of dust particles, as well as the hydrodynamics of the gas, should be explored further to understand how the RWI vortices act as the location where dust particles are accumulated.

We have fixed the value of the effective adiabatic index and have not taken into account the baroclinicity and global radial gradient of the initial surface density profile. In addition, the efficiency of disk cooling is also important on the RWI vortex (Pierens & Lin 2018). In order to investigate the applicability of the empirical relations obtained in this paper, a further parameter survey taking into account these other physical effects is needed.

5. Summary

We perform numerical simulations of the RWI in 2D, barotropic, and hydrodynamic disks using the Athena++ code. As initial conditions, we consider axisymmetric disks with a Gaussian surface density bump. We have three parameters to characterize the initial bump: the dimensionless disk aspect ratio h, the radial half-widths Δw₀, and the largest linear growth rate of the RWI γ_*. We vary these parameters in a wide parameter space and explore the nonlinear evolution for 54 models.

First, we investigate the RWI evolution. Perturbations grow as expected by the linear stability analyses not only in the linear regime but also in the weakly nonlinear regime. The axisymmetric component evolves as the RWI develops. When the axisymmetric component becomes marginally stable against the RWI for the most unstable azimuthal mode of the initial condition, the RWI saturation occurs and multiple vortices are formed in accordance with the most unstable azimuthal mode. After the RWI saturation, the vortices coalesce one after another. The axisymmetric component also approaches a stable configuration against the RWI during the vortex mergers. In the end, one quasi-stationary vortex (RWI vortex) remains when the axisymmetric component reaches the marginally stable configuration for the m = 1 mode. The regime with more than two vortices continues at most for a few orbits, and the two-vortices regime continues up to about 100 orbits. We conclude that it is difficult to observe the disks with multiple vortices originating from the RWI of one initial surface density bump except at outer disks.

Next, we turn our attention to the RWI vortex. We confirm that the RWI vortex almost corresponds to the analytic steady vortices on the Keplerian shear (the Kida and GNG solutions), as shown in previous works. Comparing the measurements of the RWI vortex with the initial conditions, we obtain empirical relations between the properties of the RWI vortices (the surface density at the vortex center: Equation (23); aspect ratio: Equation (24); radial size: Equation (25); and migration speed: Equation (26)) and the initial conditions. The radial half-width of the RWI vortex is no larger than twice the disk scale height at the vortex center. Finally, we find that the RWI vortex can be considered as a physical entity like a large fluid particle from the tracer particle analysis.

Our results are not affected by the definition of the vortex center. Even if we take into account the long-term evolution of the RWI vortex, the empirical formulae, except for Equation (24), are still applicable. On the other hand, Equation (24) is applicable only for the first few hundred orbits after the RWI vortex formation. We consider that viscosity is responsible for the long-term evolution. In order to obtain the estimation formula of the surface density at the vortex center, we use the fact that ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ remains almost constant. The vortencity conservation law explains the invariance of ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ . It is likely that the estimation formula of the vortex migration speed is broadly applicable regardless of the formation mechanism of the vortex not only in 2D disks but also in 3D disks, unless the disks have very steep global slopes of the surface density or strong and wide pressure maxima. We also find that the fluid particles passing through the vortex region contribute to about 20% of the total torque exerted on the vortex. In our interpretation, the remaining about 80% of the total torque comes from the density waves.

Our calculations have been performed under a number of simplifying assumptions, but we consider that we have captured some physical aspects of the RWI evolution and vortex. Our results provide solid theoretical ground for quantitative interpretation of the observed lopsided structures in protoplanetary disks. Future studies considering other physical effects would allow us to make models for the vortices that can be compared with observations.

Numerical computations were carried out using the Athena++ code on the Cray XC40 at the Yukawa Institute Computer Facility and the Cray XC30 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. We would like to thank Samuel Richard for showing the snapshots of his calculations. We gratefully acknowledge Hideko Nomura, James Stone, Eugene Chiang, Jeffrey Fung, and Steve Lubow for their comments. We are also grateful to the referee, who helped us improve the quality of the manuscript. This work was partially supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Numbers 15J01554 (TO); 26800106, 23103004, 15H02074, and 17H01103 (TM); and 16H05998 and 16K13786 (KT). This research was also supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) as "Exploratory Challenge on Post-K computer" (Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planets in the Solar System).

Software: Athena++ code (Stone et al. 2018, in preparation; https://princetonuniversity.github.io/athena/).

Appendix A: Analytic Steady Vortex Model

The steady solutions of vortices in shearing flow (Kida 1981; Goodman et al. 1987) are useful to understand the structure of the RWI vortices. In this section, we introduce the analytic steady vortex models.

We consider a 2D vortex orbiting a central star at angular velocity Ω_c = Ω(r_c) under the shearing-box approximation, where r_c is the distance between the central star and the center of the vortex. The shearing box is a rotating Cartesian box centered at r = r_c with an angular velocity of Ω_c. We define $x=r-{r}_{{\rm{c}}}$ and y = r_cφ and neglect the terms arising from the cylindrical geometry. In this rotating frame, the equations of motion are

$\begin{eqnarray}&&\displaystyle \frac{d{\bar{v}}_{x}}{{dt}}=\displaystyle \frac{\partial {\bar{v}}_{x}}{\partial t}+{\bar{v}}_{x}\displaystyle \frac{\partial {\bar{v}}_{x}}{\partial x}+{\bar{v}}_{y}\displaystyle \frac{\partial {\bar{v}}_{x}}{\partial y}=2{\mathfrak{S}}{{\rm{\Omega }}}_{{\rm{c}}}^{2}x+2{{\rm{\Omega }}}_{{\rm{c}}}{\bar{v}}_{y}-\displaystyle \frac{\partial {\rm{\Pi }}}{\partial x},\end{eqnarray} \tag{ 34 }$

$\begin{eqnarray}&&\displaystyle \frac{d{\bar{v}}_{y}}{{dt}}=\displaystyle \frac{\partial {\bar{v}}_{y}}{\partial t}+{\bar{v}}_{x}\displaystyle \frac{\partial {\bar{v}}_{y}}{\partial x}+{v}_{y}\displaystyle \frac{\partial {\bar{v}}_{y}}{\partial y}=-2{{\rm{\Omega }}}_{{\rm{c}}}{\bar{v}}_{x}-\displaystyle \frac{\partial {\rm{\Pi }}}{\partial y},\end{eqnarray} \tag{ 35 }$

where $\bar{{\boldsymbol{v}}}=({\bar{v}}_{x},{\bar{v}}_{y})$ is the velocity field in the local shearing box. In Equation (34), we have defined the mean shear without a vortex by ${\mathfrak{S}}\equiv -r{{\rm{\Omega }}}_{{\rm{c}}}^{-1}{[\partial {\rm{\Omega }}/\partial r]}_{r={r}_{{\rm{c}}}}$ , where ${\mathfrak{S}}=1.5$ for a Keplerian disk.

In the steady vortex, the vorticity ${\omega }_{z}={(\mathrm{rot}\bar{{\boldsymbol{v}}})}_{z}$ is assumed to be uniform. The fluid particles orbit the origin at constant angular velocity Ω_c/ψ, and the shapes of the trajectories are elliptic, where ψ is the turnover time of the vortex normalized by 2π/Ω_c. We denote the semiminor axis of the vortex by b, the semimajor axis of the vortex by a, and the vortex aspect ratio by χ ≡ a/b. According to Kida (1981; Equations (3.2) and (3.3)), the elliptic trajectories have to satisfy the following two conditions to be steady. First, the semiminor axis should be aligned with the x- or y-axis. Otherwise, the elliptic trajectory would precess. When the semiminor axis is aligned with the x-axis (radial direction), the velocity field is written as

$\begin{eqnarray}&&{\bar{v}}_{y}=-\displaystyle \frac{\chi }{\psi }{{\rm{\Omega }}}_{{\rm{c}}}x,\end{eqnarray} \tag{ 36 }$

$\begin{eqnarray}&&{\bar{v}}_{x}=\displaystyle \frac{1}{\chi \psi }{{\rm{\Omega }}}_{{\rm{c}}}y.\end{eqnarray} \tag{ 37 }$

Defining the velocity field in the inertial frame by ${\boldsymbol{v}}=({v}_{x}$ , v_y), we obtain from Equations (36) and (37)

$\begin{eqnarray}&&\displaystyle \frac{\partial ({v}_{y}-{v}_{{\rm{K}}})}{\partial x}=({\mathfrak{S}}-\displaystyle \frac{\chi }{\psi }){{\rm{\Omega }}}_{{\rm{c}}},\end{eqnarray} \tag{ 38 }$

$\begin{eqnarray}&&\displaystyle \frac{\partial {v}_{x}}{\partial y}=\displaystyle \frac{1}{\chi \psi }{{\rm{\Omega }}}_{{\rm{c}}},\end{eqnarray} \tag{ 39 }$

where ${v}_{{\rm{K}}}\,=\,{r}_{{\rm{c}}}{{\rm{\Omega }}}_{{\rm{c}}}+(1-{\mathfrak{S}}){{\rm{\Omega }}}_{{\rm{c}}}x$ is the Keplerian rotation velocity. Second, the vortex aspect ratio χ and the turnover time ψ should satisfy

$\begin{eqnarray}&&(\chi -1)={\mathfrak{S}}\psi \end{eqnarray} \tag{ 40 }$

for the invariance of the vortex aspect ratio. The analytic vortex solution that satisfies these conditions is called the Kida solution.

Goodman et al. (1987) derived another relation between χ and ψ. They assumed the velocity field satisfying Equations (36) and (37) and stationary compressible flow in the shearing box. Substituting Equations (36) and (37) into Equations (34) and (35) and assuming ∂/∂t = 0 (steady state), the pressure function Π should satisfy

$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Pi }}}{\partial {x}^{2}}=\left(\displaystyle \frac{1}{{\psi }^{2}}-2\displaystyle \frac{\chi }{\psi }+2{\mathfrak{S}}\right){{\rm{\Omega }}}_{{\rm{c}}}^{2},\end{eqnarray} \tag{ 41 }$

$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Pi }}}{\partial {y}^{2}}=\left(\displaystyle \frac{1}{{\psi }^{2}}-2\displaystyle \frac{1}{\chi \psi }\right){{\rm{\Omega }}}_{{\rm{c}}}^{2}.\end{eqnarray} \tag{ 42 }$

Substituting Equations (41) and (42) into the continuity equation, another relation between χ and ψ,

$\begin{eqnarray}&&({\chi }^{2}-1)=2{\mathfrak{S}}{\psi }^{2},\end{eqnarray} \tag{ 43 }$

is obtained (Goodman et al. 1987). The analytic vortex solution that satisfies Equation (43) rather than Equation (40) is called the GNG solution.

The Kida solution (Equation (40)) is not compatible with the GNG solution (Equation (43)), except for

$\begin{eqnarray}&&\chi =\displaystyle \frac{2+{\mathfrak{S}}}{2-{\mathfrak{S}}}=7,\end{eqnarray} \tag{ 44 }$

$\begin{eqnarray}&&\psi =\displaystyle \frac{2}{2-{\mathfrak{S}}}=4.\end{eqnarray} \tag{ 45 }$

Here we have assumed ${\mathfrak{S}}=1.5$ (Keplerian shear) in the last equalities. Note that the steady solution gives us the gradients of the velocity field and the convexities of the pressure function around the vortex center but does not give any information about the size or surface density in the vortices.

Appendix B: Diagnostics of the Axisymmetric Disk Profiles

We introduce two criteria to assess axisymmetric disks. In this section, we assume the axisymmetric disks.

B.1. Rayleigh's Condition

When there is a radius at which

$\begin{eqnarray}&&{\kappa }^{2}(r)\equiv \displaystyle \frac{1}{{r}^{3}}\displaystyle \frac{d[{rv}_{\varphi }(r)]}{dr}\lt 0\end{eqnarray} \tag{ 46 }$

is satisfied, where κ is the epicyclic frequency, the gas distribution is unstable to the rotational instability, which is an axisymmetric hydrodynamical instability in differentially rotating disks (Chandrasekhar 1961). This is known as Rayleigh's criterion. We use the term "Rayleigh's condition" when there is a radius where Rayleigh's stability criterion is violated.

In a case where the initial surface density distribution is given by a Gaussian bump, Rayleigh's condition can be regarded as a sufficient condition for the onset of the RWI. We note that Rayleigh's condition is not always a sufficient condition for the RWI in general (see Paper I). When ${{ \mathcal A }}_{0}$ is large or Δw₀ is small, the disk is unstable to the rotational instability. In this paper, we set parameters so that Rayleigh's condition is not violated. In other words, Rayleigh's condition gives the upper limit of ${{ \mathcal A }}_{0}$ for each h and Δw₀ in our parameter space.

B.2. The Critical Condition for the RWI

For a barotropic flow, the RWI requires that the vortensity, $q(r)\equiv {\kappa }^{2}/(2{\rm{\Sigma }}{\rm{\Omega }})$ , has a local minimum (Lovelace et al. 1999; see also Paper I). This is a necessary condition but is not sufficient. We found the way to explore parameters where the disk is marginally stable to the RWI (see Section 5.2 in Paper I). In this paper, we use this condition to determine the lower limit of ${{ \mathcal A }}_{0}$ for each h and Δw₀ in our parameter space.

In Paper I, we also derived the necessary and sufficient conditions for the onset of the RWI with the azimuthal mode m in a semi-analytic form as

$\begin{eqnarray}&&{\eta }_{m}\equiv {\int }_{{r}_{\mathrm{IR}}}^{{r}_{\mathrm{OR}}}\sqrt{-{D}_{\mathrm{MS},m}}{dr}\gtrsim {\eta }_{{\rm{c}}}\approx \displaystyle \frac{\pi }{2\sqrt{2}},\end{eqnarray} \tag{ 47 }$

where

$\begin{eqnarray}&&{D}_{\mathrm{MS},m}(r)\equiv \displaystyle \frac{1}{2}\displaystyle \frac{{{dB}}_{\mathrm{MS},m}}{{dr}}+\displaystyle \frac{{B}_{\mathrm{MS},m}^{2}}{4}-{C}_{\mathrm{MS},m},\end{eqnarray} \tag{ 48 }$

$\begin{eqnarray}&&{B}_{\mathrm{MS},m}(r)\equiv \displaystyle \frac{d\,\mathrm{ln}\,({{rF}}_{\mathrm{MS},m}/{\rm{\Omega }})}{{dr}},\end{eqnarray} \tag{ 49 }$

$\begin{eqnarray}{C}_{\mathrm{MS},m}(r) & \equiv & -\displaystyle \frac{{m}^{2}}{{r}^{2}}-\displaystyle \frac{{\kappa }^{2}-{m}^{2}{({{\rm{\Omega }}}_{q}-{\rm{\Omega }})}^{2}}{{c}_{{\rm{s}}}^{2}}\\ & & -\,\displaystyle \frac{2}{r}\displaystyle \frac{{\rm{\Omega }}}{({{\rm{\Omega }}}_{q}-{\rm{\Omega }})}\displaystyle \frac{d\mathrm{ln}({F}_{\mathrm{MS},m})}{{dr}},\end{eqnarray} \tag{ 50 }$

$\begin{eqnarray}&&{F}_{\mathrm{MS},m}(r)\equiv \displaystyle \frac{{\rm{\Sigma }}{\rm{\Omega }}}{{\kappa }^{2}-{m}^{2}{({{\rm{\Omega }}}_{q}-{\rm{\Omega }})}^{2}},\end{eqnarray} \tag{ 51 }$

${{\rm{\Omega }}}_{q}$ is the angular velocity at the local minimum of the vortensity q, ${c}_{{\rm{s}}}\equiv \sqrt{{\rm{\Gamma }}P/{\rm{\Sigma }}}$ is the adiabatic sound speed, and r_IR and r_OR are the radii where D_MS,m vanishes. Since Equation (47) is derived using the Sommerfeld–Wilson quantization condition (Wilson 1915; Sommerfeld 1916), ${\eta }_{{\rm{c}}}$ is equal to $\pi /(2\sqrt{2})$ only when the azimuthal mode m is large or the shape of the D_MS,m potential well is perfectly parabolic. We expect ${\eta }_{{\rm{c}}}\ \gt \ \pi /(2\sqrt{2})$ for a small m mode and shallower D_MS,m potential well and ${\eta }_{{\rm{c}}}\ \lt \ \pi /(2\sqrt{2})$ for a small m mode and steeper D_MS potential well from the knowledge of quantum mechanics.

We also use Equation (47) to investigate the stability of the axisymmetric components against the RWI during the RWI evolution in Sections 3.2 and 3.3. At that time, we calculate the values defined by Equations (48)–(51) using the radial profiles of the azimuthally averaged surface density and rotation velocity instead of those in the initial profiles.

Appendix C: Supplementary Data

C.1. The Linear Growth Rate of the Initial Conditions of Each Azimuthal Mode

In Table 3, we show the most unstable azimuthal mode, m_*, and the linear growth rates for 1 ≤ m ≤ 10 modes, γ_m, calculated by the same method as in Paper I. For all the linear stability analyses, we adopt the inner radius r_in = 0.3 and outer radius r_out = 2.5. We set the radial grid number N as N = 1024 for h = 0.1, 0.15, and 0.2 and N = 2048 for h = 0.05.

Table 3. The Most Unstable Azimuthal Mode and the Linear Growth Rates of Each Mode

Name	m_*	Linear Growth Rate for Each Azimuthal Mode (γ_m/Ω_n)
		1	2	3	4	5	6	7	8	9	10
h10w1g1	9	0.053	0.100	0.141	0.172	0.195	0.210	0.221	0.226	0.227	0.225
h10w1g2	8	0.053	0.095	0.131	0.158	0.177	0.190	0.197	0.200	0.199	0.194
h10w1g3	7	0.042	0.079	0.108	0.128	0.141	0.148	0.150	0.148	0.141	0.131
h10w1g4	6	0.033	0.061	0.081	0.094	0.0996	0.1000	0.096	0.087	0.074	0.054
h10w1g5	4	0.021	0.038	0.047	0.050	0.047	0.039	0.026	0.008	⋯	⋯
h10w2g1	6	0.074	0.140	0.189	0.220	0.238	0.242	0.237	0.222	0.199	0.168
h10w2g2	5	0.067	0.123	0.165	0.190	0.201	0.200	0.190	0.170	0.142	0.106
h10w2g3	5	0.054	0.100	0.131	0.147	0.150	0.142	0.124	0.098	0.062	0.012
h10w2g4	4	0.042	0.075	0.095	0.100	0.094	0.077	0.052	0.017	⋯	⋯
h10w2g5	3	0.026	0.044	0.050	0.043	0.025	0.002	⋯	⋯	⋯	⋯
h10w3g1	4	0.103	0.184	0.232	0.246	0.230	0.188	0.122	0.040	⋯	⋯
h10w3g2	4	0.089	0.158	0.195	0.200	0.177	0.130	0.059	0.001	⋯	⋯
h10w3g3	3	0.072	0.126	0.150	0.145	0.113	0.059	0.002	⋯	⋯	⋯
h10w3g4	3	0.054	0.090	0.100	0.089	0.041	⋯	⋯	⋯	⋯	⋯
h10w3g5	2	0.033	0.050	0.042	0.009	⋯	⋯	⋯	⋯	⋯	⋯
h10w4g1	3	0.140	0.224	0.227	0.149	0.005	⋯	⋯	⋯	⋯	⋯
h10w4g2	2	0.126	0.200	0.195	0.115	0.002	⋯	⋯	⋯	⋯	⋯
h10w4g3	2	0.098	0.150	0.131	0.044	⋯	⋯	⋯	⋯	⋯	⋯
h10w4g4	2	0.070	0.100	0.065	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h10w4g5	2	0.043	0.050	0.002	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h10w5g1	2	0.174	0.191	0.002	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h10w5g3	2	0.146	0.150	0.001	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h10w5g4	1	0.100	0.080	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h10w5g5	1	0.050	0.00	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯

h20w1g1	8	0.064	0.118	0.153	0.177	0.191	0.201	0.206	0.209	0.208	0.206
h20w1g4	4	0.044	0.076	0.093	0.100	0.099	0.096	0.086	0.076	0.061	0.048
h20w2g1	5	0.093	0.162	0.201	0.220	0.224	0.220	0.205	0.186	0.159	0.132
h20w2g4	3	0.056	0.091	0.100	0.093	0.072	0.046	0.005	⋯	⋯	⋯
h20w3g1	3	0.126	0.207	0.237	0.229	0.195	0.140	0.063	0.001	⋯	⋯
h20w3g4	2	0.069	0.100	0.087	0.045	⋯	⋯	⋯	⋯	⋯	⋯
h20w4g1	2	0.163	0.237	0.213	0.118	⋯	⋯	⋯	⋯	⋯	⋯
h20w4g4	2	0.087	0.100	0.030	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h20w5g1	2	0.191	0.192	0.020	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h20w5g4	1	0.100	0.049	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯

h15w1g1	8	0.059	0.112	0.150	0.177	0.195	0.207	0.213	0.216	0.215	0.212
h15w1g4	5	0.039	0.070	0.089	0.098	0.100	0.097	0.089	0.079	0.065	0.049
h15w2g1	5	0.086	0.155	0.199	0.224	0.233	0.231	0.219	0.200	0.175	0.144
h15w2g4	3	0.050	0.085	0.100	0.098	0.084	0.060	0.028	⋯	⋯	⋯
h15w3g1	4	0.117	0.200	0.238	0.240	0.213	0.162	0.092	0.007	⋯	⋯
h15w3g4	2	0.064	0.100	0.099	0.084	0.060	0.028	⋯	⋯	⋯	⋯
h15w4g1	2	0.153	0.233	0.222	0.136	0.003	⋯	⋯	⋯	⋯	⋯
h15w4g4	2	0.079	0.100	0.048	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h15w5g1	2	0.183	0.193	0.025	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h15w5g4	1	0.100	0.066	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯

h05w1g1	11^a	0.048	0.078	0.116	0.149	0.176	0.199	0.216	0.229	0.238	0.243
h05w1g4	7	0.023	0.044	0.063	0.079	0.090	0.097	0.100	0.099	0.095	0.088
h05w2g1	7	0.057	0.112	0.160	0.198	0.225	0.242	0.248	0.244	0.232	0.211
h05w2g4	5	0.030	0.058	0.080	0.094	0.100	0.098	0.087	0067	0.039	⋯
h05w3g1	4	0.085	0.159	0.213	0.240	0.237	0.206	0.147	0.064	0.011	⋯
h05w3g4	4	0.044	0.079	0.100	0.240	0.237	0.206	0.147	0.064	0.010	⋯
h05w4g1	3	0.128	0.214	0.225	0.151	0.001	⋯	⋯	⋯	⋯	⋯
h05w4g4	2	0.064	0.100	0.082	0.001	⋯	⋯	⋯	⋯	⋯	⋯
h05w5g1	2	0.168	0.188	0.001	⋯	⋯	⋯	⋯	⋯	⋯	⋯
h05w5g4	1	0.100	0088	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯

Notes. Name: the name of the model. m_*: the most unstable azimuthal mode. γ_m/Ω_n: the linear growth rate against the RWI for each azimuthal mode m. The bold values show the largest linear growth rates against the RWI for each model.

^aThe largest linear growth rate (γ_*/Ω_n) is 0.244 in the "h05w1g1" model.

Download table as: ASCII Typeset image

C.2. Table of the Results

We show the values of τ₂, τ₁, and τ₂ − τ₁ in Table 4. For the "h10w5g3"–"h10w5g5," "h20w5g4," "h15w5g4," and "h05w5g4" models, we are not able to measure τ₂ and τ₁ with visual inspection due to low m_*. In those cases, we set τ₂ to no data and τ₁ to the orbit number at the saturation. From the long-term calculations, we also show the values of τ_χ in Table 4. For the "h10w2g1," "h10w3g1," "h10w4g1," "h10w4g2," "h10w5g1," "h10w5g2," "h20w1g1," "h20w2g1," "h20w3g1," "h20w4g1," "h20w5g1," "h15w1g1," "h15w2g1," "h15w3g1," "h15w4g1," "h15w5g1," "h05w1g1," "h05w2g1," "h05w4g1," and "h05w5g1" models, we are not able to measure τ_χ due to fast vortex migration. In those cases, we set τ_χ to no data.

Table 4. Values of τ₂, τ₁, τ₁ − τ₂, and τ_χ

Name	τ₂	τ₁	τ₁ − τ₂	τ_χ
h10w1g1	23	45	22	5.2E3
h10w1g2	33	57	24	3.6E3
h10w1g3	36	96	60	3.7E3
h10w1g4	52	66	14	3.0E3
h10w1g5	54	98	44	1.3E3
h10w2g1	24	84	60	⋯
h10w2g2	21	46	25	2.4E4
h10w2g3	21	111	90	3.6E3
h10w2g4	28	49	21	3.1E3
h10w2g5	41	109	68	2.8E3
h10w3g1	18	20	2	⋯
h10w3g2	20	45	28	5.8E3
h10w3g3	22	42	20	8.9E3
h10w3g4	25	35	10	8.0E3
h10w3g5	47	72	25	5.0E3
h10w4g1	10	17	7	⋯
h10w4g2	12	20	8	⋯
h10w4g3	14	28	14	7.3E3
h10w4g4	20	36	16	1.4E4
h10w4g5	36	45	9	1.7E4
h10w5g1	12	13	1	⋯
h10w5g3	⋯	15	⋯	⋯
h10w5g4	⋯	22	⋯	3.3E3
h10w5g5	⋯	52	⋯	5.0E3

h20w1g1	24	25	1	⋯
h20w1g4	30	59	29	2.2E3
h20w2g1	14	27	13	⋯
h20w2g4	21	37	16	1.3E3
h20w3g1	12	15	3	⋯
h20w3g4	29	42	13	1.4E3
h20w4g1	8	16	8	⋯
h20w4g4	19	26	7	2.7E3
h20w5g1	11	12	1	⋯
h20w5g4	⋯	23	⋯	2.0E3

h15w1g1	37	62	25	⋯
h15w1g4	30	82	52	3.3E3
h15w2g1	16	27	11	⋯
h15w2g4	25	32	7	1.1E3
h15w3g1	14	21	7	⋯
h15w3g4	19	35	16	2.0E3
h15w4g1	8	17	9	⋯
h15w4g4	19	29	10	3.0E3
h15w5g1	11	13	2	⋯
h15w5g4	⋯	23	⋯	2.0E3

h05w1g1	69	177	108	⋯
h05w1g4	57	195	138	1.4E4
h05w2g1	29	199	170	⋯
h05w2g4	36	182	146	2.1E4
h05w3g1	16	45	29	9.2E3
h05w3g4	26	49	23	6.8E3
h05w4g1	10	21	11	⋯
h05w4g4	19	46	27	5.2E3
h05w5g1	13	15	2	⋯
h05w5g4	⋯	24	⋯	2.1E4

Note. Name: the name of the model. τ₂: the orbit when the number of vortices becomes two. τ₁: the orbit when the final vortex merger occurs. τ₁ − τ₂: the duration of the two-vortices regime. τ_χ: the decreasing time of χ₂ in the long-term calculations.

Download table as: ASCII Typeset image

Appendix D: Numerical Test

In order to check the convergence of our numerical calculations, we have additionally performed a high-resolution calculation and a low-resolution calculation of the "h10w3g1" model. The high-resolution calculation has twice as many cells as the fiducial calculation in each direction. In contrast, the low-resolution calculation has half the cells compared to the fiducial calculation in each direction. We also perform wide-domain calculations in order to check the effects of the inner boundary on the RWI vortices. In the wide-domain calculations, the resolution is almost the same as that of the fiducial calculation, but the radius of the inner boundary is set at r_in = 0.2r_n, 0.1r_n, and 0.03r_n instead of 0.3r_n.

In all the calculations, the RWI vortices are formed at almost the same time (τ ≈ 20). We measure five parameters ( ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ , χ₂, ψ₂, Δr, and ξ) that characterize the vortex and its migration and are almost constant in each additional calculation due to the fast vortex migration. Figure 19 shows the deviation of the parameters from those of the high-resolution calculation. Since the RWI vortex in the fiducial calculation shows the same values of the parameters as those in the high-resolution calculation within 5% (see panel (a) of Figure 19), we conclude that the fiducial calculation has sufficiently high resolution and the results are converged. Panel (b) of Figure 19 shows that the values in the wide-domain calculations are also equivalent to those in the fiducial calculation within 5%. This indicates that the inner boundary does not have significant effects on the RWI vortex.

**Figure 19.** Panel (a) shows the deviations of the parameters ( ${{\rm{\Sigma }}}_{{\rm{v}}}{r}_{{\rm{v}}}^{1.5}$ , χ₂, ψ₂, δr, and ξ) from the high-resolution calculation in the high-resolution calculation (red circles), fiducial calculation (blue triangles), and low-resolution calculation (green diamonds). Panel (b) shows the deviations of the parameters from the fiducial calculations in the wide-range calculations ( ${r}_{\mathrm{in}}=0.3{r}_{{\rm{n}}}$ : black crosses; ${r}_{\mathrm{in}}=0.2{r}_{{\rm{n}}}$ : red circles; ${r}_{\mathrm{in}}=0.1{r}_{{\rm{n}}}$ : blue triangles; and ${r}_{\mathrm{in}}=0.03{r}_{{\rm{n}}}$ : green diamonds).
Download figure:
Standard image High-resolution image

Parametric Study of the Rossby Wave Instability in a Two-dimensional Barotropic Disk. II. Nonlinear Calculations

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Disk Model and Numerical Method

2.1. Disk Models and Initial Condition

2.2. Numerical Method

3. Evolution of the RWI

3.1. Linear and Weakly Nonlinear Regime of the RWI

3.2. Saturation Mechanism of the RWI

3.3. Vortex Merger

4. Quasi-stationary Vortex Formed by the RWI

4.1. Analyses

4.1.1. Structure of the RWI Vortex

4.1.2. Measurements about the RWI Vortex

4.1.3. Tracer Particle Analysis

4.2. Results

4.2.1. Comparison with the Analytic Steady Vortices

4.2.2. Empirical Relations of the RWI Vortex

4.2.3. The RWI Vortex from the Point of View of Tracer Particles

4.3. Discussions

4.3.1. Definition of the Vortex Center

4.3.2. Applicability of Empirical Formulae to the Long-term Evolution

4.3.3. Reason Why {{\rm{\Sigma }}}_{v}{{\rm{r}}}_{v}^{1.5} Is Approximately Constant

4.3.4. Comparison with Another Formula of Vortex Aspect Ratio

4.3.5. Generality of the Empirical Formula of Vortex Migration Speed

4.3.6. Mechanism of Vortex Migration

4.3.7. Other Physical Effects

5. Summary

Appendix A: Analytic Steady Vortex Model

Appendix B: Diagnostics of the Axisymmetric Disk Profiles

B.1. Rayleigh's Condition

B.2. The Critical Condition for the RWI

Appendix C: Supplementary Data

C.1. The Linear Growth Rate of the Initial Conditions of Each Azimuthal Mode

C.2. Table of the Results

Appendix D: Numerical Test

4.3.3. Reason Why ${{\rm{\Sigma }}}_{v}{{\rm{r}}}_{v}^{1.5}$ Is Approximately Constant