# Vortices in stratified protoplanetary disks: from baroclinic instability to vortex layers

Pierre Barge, Samuel Richard, Stéphane Le Dizès

► To cite this version:

Pierre Barge, Samuel Richard, Stéphane Le Dizès. Vortices in stratified protoplanetary disks: from baroclinic instability to vortex layers. Astronomy & Astrophysics - A&A, 2016, 592, pp.A136. 10.1051/0004-6361/201628381. hal-01387029

HAL Id: hal-01387029

https://hal.science/hal-01387029v1

Submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Vortices in stratified protoplanetary disks

## From baroclinic instability to vortex layers

P. Barge<fnref n="1" />, S. Richard<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, and S. Le Dizès<fnref n="3" />

Received 25 February 2016 / Accepted 12 June 2016

## ABSTRACT

*Context.* Large-scale vortices could play a key role in the evolution of protoplanetary disks, particularly in the dead-zone where no turbulence associated with magnetic field is expected. Their possible formation by the subcritical baroclinic instability is a complex issue because of the vertical structure of the disk and the elliptical instability.

*Aims.* In 2D disks the baroclinic instability is studied as a function of the thermal transfer efficiency. In 3D disks we explore the importance of radial and vertical stratification on the processes of vortex formation and amplification.

*Methods.* Numerical simulations are performed using a fully compressible hydrodynamical code based on a second-order finite volume method. We assume a perfect gas law in inviscid disk models in which heat transfer is due to either relaxation or diffusion.

*Results.* In 2D, the baroclinic instability with thermal relaxation leads to the formation of large-scale vortices, which are unstable with respect to the elliptic instability. In the presence of heat diffusion, hollow vortices are formed which evolve into vortical structures with a turbulent core. In 3D, the disk stratification is found to be unstable in a finite layer which can include the mid-plane or not. When the unstable layer contains the mid-plane, the 3D baroclinic instability with thermal relaxation is found to develop first in the unstable layer as in 2D, producing large-scale vortices. These vortices are then stretched out in the stable layer, creating long-lived columnar vortical structures extending through the width of the disk. They are also found to be the source of internal vortex layers that develop across the whole disk along baroclinic critical layer surfaces, and form new vortices in the upper region of the disk.

*Conclusions.* In 3D disks, vortices can survive for a very long time if the production of vorticity by the baroclinic amplification balances the destruction of vorticity by the elliptical instability. However, this possibility is strongly dependent on the disk properties. Such baroclinic vortices could play a significant role in the global disk evolution and in participating to the decoupling of solids from the gas component. They could also contribute to the formation of new out-of-plane vortices by a critical layer excitation mechanism.

**Key words.** hydrodynamics – instabilities – protoplanetary disks

## 1. Introduction

Large-scale vortices in protoplanetary disks could play an impor-
tant role in the evolution of gas and dust at the decoupling stage
of planet formation. For example, they could participate in the
transport of angular momentum, providing a natural mechanism
for accretion onto the star; they could also stop the systematic
drift of the solid particles thanks to the trapping mechanism iden-
tified by Barge & Sommeria (1995). This is the reason why the
formation, stability, and evolution of such vortices in protoplan-
etary disks deserve special attention. Recent observations with
high resolution in the millimeter and centimeter range, have re-
cently raised a strong interest on dust trapping and vortices. Vor-
tex formation in protoplanetary disks has been studied a lot in
the last ten years. Their formation by the Rossby-wave instabil-
ity was addressed first (Lovelace et al. 1999; Li et al. 2000) and
studied in detail in 2D and 3D disks (Li et al. 2001; Meheut et al.
2012b,a; Lin 2012, 2013; Richard et al. 2013). Other mecha-
nisms (like the baroclinic instability) were also proposed for
the formation of vortices (Klahr 2004; Petersen et al. 2007a,b;
Lesur & Papaloizou 2010).

The baroclinic instability is well known in geophysical fluid
dynamics where it is responsible for the formation of cyclones
in the atmosphere of the Earth (Pedlosky 1987). It originates
in the systematic bending of the iso-pressure surfaces with re-
spect to the iso-density surfaces, which leads to a source term
in the evolution equation of the vorticity field, the so-called
baroclinic source term. It is only recently that this instability
was actually addressed in the context of protoplanetary disks.
Klahr & Bodenheimer (2003) explored the possibility of finding
a baroclinic instability in 2D protoplanetary disks with global
numerical simulations. They reported the formation of large-
scale vortices, but the actual baroclinic origin of these vor-
tices is unclear since thermal transfer is absent from their sim-
ulations, contrary to what is normally required (Petersen et al.
2007a,b). The non-linear nature of this instability was under-
lined by Lesur & Papaloizou (2010, thereafter LP10) using in-
compressible shearing box simulations, clearly explaining why
the works of Klahr (2004) and Johnson & Gammie (2005) failed
to describe the baroclinic amplification of the vortices. In the
2D context, Raettig et al. (2013) have tried to estimate the gen-
eration of baroclinic turbulence in protoplanetary disks. They

<footnote n="1">¹ Aix Marseille Université, CNRS, LAM UMR 7326, 13013 Marseille, France
e-mail: pierre.barge@lam.fr</footnote>

<footnote n="2">² Aix Marseille Université, CNRS, Centrale Marseille, IRPHE UMR 7342, 13013 Marseille, France</footnote>

<footnote n="3">³ Astronomy Unit, Queen Mary University of London, Mile End Road, London E1 4NS, UK</footnote>
studied the dependence of vortex amplification as a function of the entropy gradient, the thermal transfer efficiency and the numerical resolution of the computations.

Lyra & Klahr (2011) also investigated numerically the growth of the subcritical baroclinic instability in a magnetic disk; they concluded that baroclinic vortices are destroyed by the magneto-rotational instability and can only survive in the dead zone of the protoplanetary disks.

This instability was then studied by LP10 who noted that the baroclinic instability in protoplanetary disks is basically non-linear and differs from the linear instability that forms cyclones in geophysical flows. They also found that, in appropriate conditions, 3D baroclinic vortices may survive the elliptical instability. Klahr & Hubbard (2014) recently proposed that a convective overstability could develop following a linear amplification of the epicyclic oscillations in axisymmetric and vertically un-stratified disks and Lyra (2014) showed using 3D simulations, it could form vortices. In this paper, this mechanism is not considered. Here, we focus on the subcritical baroclinic instability as described by Lesur & Papaloizou (2010).

LP10 performed 3D simulations to study the stability of baroclinic vortices without vertical stratification and with a radial stratification only. However, to our knowledge, the baroclinic instability has never been studied in a fully stratified 3D disk. In this paper, we revisit the problem using fully compressible and long term numerical simulations. First, we start from numerical simulations in a non-homotropic disk and show that vortices can be created by the baroclinic instability in a 2D disk. Then we check the stability of the 2D vortices with respect to the 3D elliptic instability. Finally, we study the formation and growth of 3D vortices in a 3D stratified disk.

The paper is organized in the following way: the basic equations, the disk model and the numerical method are presented in Sect. 2. Section 3 is devoted to the 2D simulations. The effect of the nature and strength of the thermal transfer on the baroclinic instability is analyzed in detail. 3D simulations are performed in Sect. 4. The structure of the vertical stratification is first examined. Then, the formation of long-lived vortices in a fully 3D simulation with an unstable mid-plane layer is demonstrated. The results are briefly summarized and discussed in Sect. 5.

## 2. Fluid equations and numerical method

### 2.1. Standard disk assumptions

The gas of the nebula is assumed to be a mixture of molecular hydrogen and helium with a mean molecular weight $\mu = 2.34$ g/mol. The low pressure conditions inside the disk justify the use of the perfect gas law as the appropriate equation of state.

In the steady state, the gas is stratified with an hydrostatic equilibrium in the vertical direction $z$ and a pressure modified centrifugal equilibrium in the radial direction $r$. The gas flows around the star at nearly the Keplerian velocity and fills a flared disk in which the pressure scale height is $H(r) = c_s(r)/\Omega_K(r)$ where $c_s = \sqrt{\gamma P/\rho}$ is the sound speed and $\Omega_K = 1/r^{3/2}$ is the Kepler angular velocity.

The surface density and the temperature of the disk are assumed to decrease as simple power laws of the distance to the star: $r^{-p}$ and $r^{-q}$, respectively; numerical values follow the model of the minimum mass solar nebula (e.g. Hayashi 1981). The self-gravity of the gas is neglected and we focus on optically thick regions of the nebula where the coupling of the gas with the magnetic field is negligible due to weak ionization. In these regions the presence of dust particles embedded in the gas also justifies

the presence of thermal transfer which is a crucial ingredient for the baroclinic amplification of the vortices. The problem is addressed by solving numerically the full set of the compressible hydrodynamical equations in cylindrical coordinates with an energy equation that includes thermal transfer.

### 2.2. The governing equations

The standard equations for an inviscid gas flowing around a central gravitational potential read

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{V}) = 0 \qquad (1) $$

$$ \frac{\partial \rho \mathbf{V}}{\partial t} + \nabla \cdot (\rho \mathbf{V} \mathbf{V}) + \nabla \cdot \mathbf{P} = \rho \nabla \phi \qquad (2) $$

$$ \frac{\partial \rho e}{\partial t} + \nabla \cdot (\mathbf{V}(\rho e + P)) = \rho \mathbf{V} \cdot \nabla \phi + Q \qquad (3) $$

$$ \rho e = \frac{P}{\gamma - 1} + \frac{1}{2} \rho \mathbf{V}^2 \qquad (4) $$

where $\rho$ and $P$ are the density and pressure of the gas; $\mathbf{V}$ is the gas velocity and $e$ its total specific energy; $\phi$ is the gravitational potential, $Q$ the heat transfer and $\gamma = 1.4$ the adiabatic index.

The equations will be used in non-dimensional form thanks to the following normalization

$$ \mathbf{r} = r_o \tilde{\mathbf{r}} \qquad \mathbf{V} = v_o \tilde{\mathbf{V}} \qquad (5a) $$

$$ \rho = \rho_o \tilde{\rho} \qquad P = \rho_o v_o^2 \tilde{P} \qquad T = \frac{v_o^2}{R} \tilde{T} \qquad (5b) $$

where the $o$ index refers to the values at the reference radius $r_o$ and $R$ is the ideal gas constant; $v_o = (GM/r_o)^{1/2}$ and $\rho_o = \rho(r_o, z=0)$ stand for the values of the velocity and the density in the mid-plane and at the radius $r_o$.

In the following, the tilde over the new dimensionless variables is dropped. The form of the Eqs. (1)–(4) is unchanged, except the gravitational potential which now reads $\phi = 1/(r^2 + z^2)^{1/2}$. The important parameter is the Mach number $M_{A_o} = v_o/\sqrt{\gamma R T_{o}}$ where $T_o$ is the gas temperature at radius $r_o$.

### 2.3. Equilibrium state of the disk

In 2D the steady state is assumed to be simple power law in temperature and surface density

$$ T_D = \frac{r^{-q}}{\gamma M_{A_o}^2} \qquad (6) $$

$$ \Sigma_D = r^{-p} \qquad (7) $$

where $p = 1.5$ and $q = 0.5$ in the case of the standard MMSN hypothesis. The perfect gas equation of state gives

$$ P_D = \Sigma_D T_D = \frac{r^{-p-q}}{\gamma M_{A_o}^2}. \qquad (8) $$

The radial velocity is assumed to be null and the azimuthal velocity is given by the radial hydrodynamic equilibrium

$$ v_D = \sqrt{r \frac{\partial \phi}{\partial r} + \frac{r}{\Sigma_D} \frac{\partial P_D}{\partial r}} \qquad (9) $$

For large Mach number, the disk angular velocity $\Omega_D = v_D/r$ is close to the Kepler profile $\Omega_K$ in the mid-plane. This equilibrium
solution only depends on a small number of parameters: the indices $p$ and $q$, and the Mach number $M_{A_o}$. It is also interesting to note that the structure of the disk is independent of the disk mass.

The pressure scale height is given by

$$H(r) = \frac{c_s}{\Omega_K} = \frac{r^{(3-q)/2}}{M_{A_o}}. \qquad (10)$$

## 2.4. Radial stratification

In the radial direction, disk stratification is measured by the Brunt-Väisälä frequency

$$N_r^2 = -\frac{1}{\gamma \rho} \frac{\partial P}{\partial r} \frac{\partial}{\partial r} \ln \left( \frac{P}{\rho^\gamma} \right), \qquad (11)$$

and the Schwarzschild criterion states that this stratification is stable if $N_r^2 > 0$ and unstable if $N_r^2 < 0$.

In a 2D disk at equilibrium, pressure and density can be expressed as simple power-laws and the Brunt-Väisälä frequency writes

$$N_r^2 = \frac{(p+q)(p(\gamma-1)-q)}{\gamma r^2} T. \qquad (12)$$

Then, if pressure and density are decreasing outward (i.e $p$ and $q$ are positive), we find that stratification is unstable when $p < q/(\gamma-1)$. In term of the Richardson number $Ri = N_r^2/S^2$, where $S = -3\Omega_K/2$ is the Kepler shear, the condition of instability is then equivalent to

$$Ri > 0. \qquad (13)$$

## 2.5. Thermal transfer

Thermal transfer in optically thick regions of the nebula results from an energy dissipation either in the form of a black body radiation escaping from the disk or a diffusive process through the scattering of IR photons in the disk. These complex phenomena are not considered here. In the present paper, we mainly focus on the physical consequences of thermal transfer. As commonly done, we assume that thermal transfer only takes place in the form of thermal relaxation or heat diffusion. Moreover, we clearly distinguish between the two cases by assuming either pure relaxation or pure diffusion.

- Thermal relaxation. In this case the gas temperature variations are assumed to relax toward an equilibrium temperature $T_0$ at a rate

$$Q = \frac{\gamma}{\gamma-1} \rho \frac{T-T_0}{\tau}, \qquad (14)$$

where $\tau$ is the cooling time normalized to the inverse of the orbital frequency. Thermal relaxation then affects equally all spatial scales. However, it does not conserve the total energy which is now a function of time.

- Heat diffusion. The thermal energy of the gas is assumed to diffuse in the disk at a rate

$$Q = \frac{\gamma}{\gamma-1} \frac{\rho}{Pe} \Delta(T-T_0), \qquad (15)$$

where $\Delta$ is the Laplacian, and $Pe$ is the Péclet number, a dimensionless number that quantifies the importance of the transfer in the disk and defined as

$$Pe = \frac{H(r_*) r_* \Omega(r_*)}{\kappa}, \qquad (16)$$

where $\kappa$ is a diffusion coefficient and $r_*$ the mean orbital radius of the computational box (in what follows we will use $Pe/\rho(r_*)$ instead of $Pe$). Heat diffusion operates at all spatial scales but its effect is stronger on smaller scales. We shall see the consequences of this property on the evolution and structure of large size vortices. We also note that, in contrast to thermal relaxation, heat diffusion conserves total energy.

## 2.6. Numerical method

The system of non-linear Eqs. (1)–(4) is solved using a second order finite volume scheme, the MUSCL Hancock scheme, and an exact Riemann solver. The 2D version of the code is described in Inaba et al. (2005) (see also Surville 2013). The 3D version is the same as used in Richard et al. (2013). The two versions were used in this paper to study the baroclinic instability in a 2D or a 3D context. We do not use the standard shearing sheet approximations and our simulations correspond to local integrations of the full set of the fluid equations.

No explicit dissipation mechanism is present in our simulations. Viscous effects (numerical diffusion) originate either at a global scale owing to the numerical interpolation of the field on a finite size grid, or at small scales, owing to the existence of a cut-off length in numerical integration. They can be described with effective viscosities: $v_a$ at large scales and $v_n$ at small scales. The first one has been estimated by analyzing the time evolution of known viscous solutions e.g. Pringle (1981); numerical benchmarks with the 2D version of our code (Surville, Ph.D. Thesis) have shown that $\Delta\Sigma/(\Delta t\Sigma\Omega) \approx 3.2 \times 10^{-6}$ so that $v_a \approx 7.2 \times 10^{-6} r_*^2 \Omega_* \approx 10^9 m^2 s^{-1}$ (or $Re \approx 5.2 \times 10^5$). The second one scales as $\Omega_* \delta_*^2$ where $\delta_*$ is the mesh size; in the radial direction, this estimate gives $v_n \approx 10^{10} m^2 s^{-1}$ (or $Re \approx 5.2 \times 10^4$).

# 3. Two dimensional simulations of the baroclinic instability

In this section, we perform 2D simulations to recover and extend a few results already obtained in the literature. As explained in LP10 the necessary conditions for the subcritical baroclinic instability (SBI) are: (i) a radially unstable entropy gradient or negative Richardson number; (ii) a thermal transfer at a rate that can sustain the baroclinic generation of vorticity; (iii) a seeding of the disk with finite amplitude perturbations to trigger the nonlinear response of the disk.

To satisfy the above requirements our 2D simulations are performed with the following initial conditions:

(a) a background temperature with an index $q = 0.5$;

(b) a fixed value of the Mach number $M_{A_o} = 8$ (at $r_o$) or correspondingly a scale height $H = 1.55$ at $r = 7.5$;

(c) a perturbation of the equilibrium state obtained by randomly superimposing on the background density small Gaussian bumps with amplitude ~10%.

The computational domain of the simulations was limited in the radial direction to $7 < r < 8$ and in the azimuthal direction to $0 < \theta < \pi/4$. The numerical resolution was set up to 300 × 300, that is 465 cells by scale height in the radial direction and
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Vorticity field after 20, 200, 400 and 600 rotations (from left to right, respectively) in a 2D disk with thermal relaxation with a cooling time $τ = 160$. Top row: stable stratification ($p = 1.5, q = 0.5$). Bottom row: unstable stratification ($p = 1, q = 0.5$). Vortices are rapidly decaying in the case of a stable stratification while they are slowly but systematically amplified in the case of an unstable stratification.</figcaption></figure>

80 cells by scale height in the azimuthal direction. The results of the simulations are presented below.

## 3.1. A two stage development

In order to illustrate the development of the baroclinic instability in 2D disks, we have first carried out simulations with thermal relaxation. The crucial importance of radial stratification was tested by comparing simulations performed for background densities with significantly different radial gradients. We have chosen two power-laws with different values of the index $p$.

For $p = 1.5$, stratification is stable since the Richardson number is positive ($Ri \approx 2 \times 10^{-3}$ at $r = 7.5$) whereas for $p = 1$, stratification is unstable since the Richardson number is negative ($Ri \approx -1.5 \times 10^{-3}$ at $r = 7.5$). Baroclinic feedback is expected only in the case of unstable stratification that is when $p < q/(\gamma - 1)$ (here $p < 1.25$).

In Fig. 1, the vorticity is plotted at different instants of the simulation. Two different stages can be distinguished in the evolution of the vorticity: a first one during which small vortices are formed from the initial density perturbations, nearly irrespective of stratification; a second one, after approximately 50 rotations, during which the evolution strongly depends on stratification. If stratification is stable, the small vortices formed during the first stage gradually decay in a few hundred rotations whereas, if stratification is unstable, they are amplified and grow as a function of time. These observation confirms known results.

Figure 2 shows the variation of the enstrophy $\iint \omega_z^2 dS/2$ as a function of time. This enstrophy evolution is consistent with the evolution of the vorticity field. The peak at the very beginning corresponds to the formation of vortices from the initial density seeds. The rest of the plot corresponds either to a gradual decay of the vortices or to their baroclinic amplification depending on whether stratification is stable or not, respectively.

The baroclinic amplification of the vortices (thereafter BVA) was first identified by Petersen et al. (2007b). LP10 also explained why unstable disk stratification is required for the vortex amplification. The mechanism operates along the cyclic path of a gas parcel around the vortex center; it consists in a thermalization of the parcel during its azimuthal excursion followed by an acceleration of the parcel during its radial excursion. It must be noted that, in the case of a stable stratification, the mechanism

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Evolution of the enstrophy in a 2D stratified disk with thermal relaxation. The full line corresponds to an unstable stratification ($N_r^2 < 0$) whereas the dotted line corresponds to a stable stratification ($N_r^2 > 0$). Enstrophy is computed for $τ = 160$.</figcaption></figure>

becomes ineffective since the radial motion of the parcel is not accelerated but, on the contrary, decelerated; this leads to a gradual weakening of the vortex instead of an amplification. This explains why vortices weaken when $N_r^2 > 0$. We want to stress that this gradual decay of the vortices is very different from the breakdown of the vortex which may occur in 3D disks due to the elliptic instability. In the rest of the section, we only consider an unstable stratification by fixing the parameter $p$ and $q$ to the values $p = 1.5$ and $q = 0.5$. Our goal is to analyze the importance of the thermal transfer in the baroclinic instability. We consider first thermal relaxation then heat diffusion.

## 3.2. Thermal relaxation

The results are presented in Fig. 3 which shows the evolution of the enstrophy as a function of time and for different values of the cooling time $τ$. The cooling time is given in term of the orbital period at $r = 7.5$. In the beginning of the simulation, the production of vorticity is very effective with the formation of a peak of enstrophy which is associated with the birth of small vortices
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Enstrophy as a function of time, in a 2D stratified disk ($p = 1.5, q = 0.5$) with thermal relaxation for different values of the cooling time. Only intermediate values of $\tau$ lead to an efficient vortex amplification.</figcaption></figure>

from the density perturbations initially seeded in the disk. We also observe that the longer the cooling time, the higher the initial peak of enstrophy. The generation of vortices is therefore more efficient for longer cooling times.

The subsequent evolution of the enstrophy depends on the cooling time. For both short and large cooling times, enstrophy decreases monotonically while for intermediate values it increases after the rapid phase of decay. For intermediate cooling time, the growth of enstrophy corresponds to a baroclinic amplification of the vortices. This amplification is effective on a time scale much longer (some hundred rotations) than the time necessary for the creation of the vortices. During this amplification phase, vortices are gradually merging one another to give a single vortex that continues to grow until it reaches the maximum size compatible with the limits of the computational box.

When $\tau = \infty$ or $\tau = 0.1$ no baroclinic amplification is observed. In the first case the flow is adiabatic (no thermal transfer) and we find that the small vortices formed during the first phase cannot grow. This confirms that thermal transfer is required to get a baroclinic amplification. It is, however, interesting to stress that the adiabatic assumption does not inhibit vortex formation but only the amplification mechanism. In the second case ($\tau = 0.1$) enstrophy does not increase because no vortex has emerged during the first phase. No baroclinic amplification is possible since no vortex is present. To summarize, the formation of large and strong vortices is possible for intermediate values of $\tau$ ($\tau \approx \Omega_K^{-1}$) only. After studying the global production of vorticity with the enstrophy, it is also interesting to look at the growth of a single vortex. For this purpose we follow the evolution of the strongest vortex formed at the end of the birth period. The strength of this vortex is estimated using the Rossby number

$$Ro = \frac{\omega_z}{2\Omega_D(r_v)} \qquad (17)$$

where $\omega_z$ is maximum relative vorticity of the vortex and $r_v$ is the radial position of the vortex. Figure 4 shows the evolution of the Rossby number as a function of time and for various cooling times. We clearly see the strength increase of the vortex when the baroclinic amplification is active. The strongest vortex also increases in size. This is illustrated in Fig. 5 where is shown the evolution of the radial extent of the strongest vortex during its amplification. The aspect ratio $\chi$ of the vortices decreases as a function of time (see Fig. 6), as expected from the systematic increase of the vorticity. The strong fluctuations observed in the

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Rossby number of the strongest vortex as a function of time in a 2D stratified disk ($p = 1.5, q = 0.5$) with thermal relaxation for different values of the cooling time.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Radial extent of the strongest vortex as a function of time in a 2D stratified disk ($p = 1.5, q = 0.5$) with thermal relaxation for different values of the cooling time.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Aspect ratio of the strongest vortex as a function of time in a 2D stratified disk ($p = 1.5, q = 0.5$) with thermal relaxation for different values of the cooling time.</figcaption></figure>

beginning of the evolution are due to the difficulty to evaluate the aspect ratio of the vortices as they are merging with one another. After some hundred rotations, a single vortex remains in the disk and the dispersion of the measures strongly reduces, uncovering a slow decrease of $\chi$ as a function of time. Figure 7 shows the Rossby number as a function of the aspect ratio and for various cooling times. A clear relation exists between the two
<figure><img src="image_8.png" /><figcaption>Fig. 7. Rossby number of the strongest vortex as a function of its aspect ratio in a 2D stratified disk ($p = 1.5$, $q = 0.5$) with thermal relaxation for different values of the cooling time and different instants. The black curve represents the relation (18) obtained for Kida's vortex model.</figcaption></figure>

parameters, nearly independent of the cooling time. It is well fit- ted by the analytical expression derived for Kida's vortex model (Kida 1981)

$$Ro = -\frac{3}{4} \frac{\chi + 1}{\chi(\chi - 1)} \qquad (18)$$

### 3.3. Heat diffusion

In this subsection, thermal transfer is due to heat diffusion. The strength of heat diffusion is quantified by the Péclet number. The smaller this number, the stronger the diffusion. The non-diffusive case corresponds to an infinite Péclet number. The computa- tional domain has been chosen to be a box with $7 < r < 8$ and $0 < \theta < \pi/2$ except in the simulation presented in Fig. 9 in which the box is extended to $6.5 < r < 8.5$ and $0 < \theta < \pi$. The re- sults we get for heat diffusion are quite similar to those obtained for thermal relaxation. Two stages can also be distinguished in the evolution of the enstrophy plotted in Fig. 8. The peak at the very beginning corresponds to a first stage with the formation of small vortices from the density bumps initially seeded in the disk. The subsequent evolution corresponds to a second stage during which vortices are growing under the baroclinic feedback mechanism. Our simulations show that baroclinic amplification is effective only for large enough values of the Péclet number ($Pe > 2 \times 10^4$). However, as for thermal relaxation, the amplifi- cation is weakened for very large $Pe$ and disappears at $Pe = \infty$. There is therefore an optimal heat diffusion for the amplification which corresponds to $Pe \approx 10^5$.

The main differences between an evolution under thermal re- laxation and under heat diffusion clearly appear during the sec- ond stage, after the vortices have grown and merged into a sin- gle large vortex. Vorticity during this time period is presented in Fig. 9 in the form of 6 different panels corresponding to 300, 400, 500, 700, 860 and 960 rotations, respectively. In the three first panels a strong peak develops in the center, surrounded by an annular bump. It must be noted that this peculiar structure looks like the one observed by Raettig et al. (2013) in their Fig. 8. In the fourth panel, after 700 rotations, a hollow vortex has formed: the negative vorticity in the center becomes lower than the neg- ative vorticity on the periphery. In the two last panels the hollow vortex is strongly perturbed in the center but not in the outer rims that preserves a coherent shape. At this stage the hollow vortex is transformed in a vortex with a turbulent core that we observed

<figure><img src="image_9.png" /><figcaption>Fig. 8. Enstrophy as a function of time in a 2D stratified disk, in the case of heat diffusion and for different values of the Péclet number.</figcaption></figure>

till the end of the simulations. This is a significantly different evolution than in the case of thermal relaxation.

The differences observed in the evolution of the vortices lie in the details of the amplification mechanism. Indeed, BVA can work if the gas warms up when moving azimuthally in the in- ner side (the side of the star) and cools down when moving az- imuthally in the outer side (opposite side of the star). Inside a large size vortex, this condition is always satisfied in the case of thermal relaxation but not in the case of heat diffusion. In the presence of heat diffusion, energy exchanges are not uniform inside in the vortex. This point is illustrated in Fig. 10 which displays the distribution of heat transfert (Q) inside a growing vortex at two different instants, after 250 and 500 orbital peri- ods. This figure shows that when the vortex reaches a large size, the distribution of thermal transfer changes sign in the vortex core while it remains unchanged in the vortex border.

So, the baroclinic amplification mechanism is still active in the vortex periphery, but no more in the vortex core where an opposite mechanism now leads to the damping of the vorticity. This could explain the formation of hollow vortices. These hol- low vortices turn out to be unstable as observed in the last two panels of Fig. 9.

This simulation shows that, in the case of heat diffusion, the structure of a baroclinic vortex depends on its spatial extent and evolves in a complex way. At the end of the simulations the vor- tex fills a large part of the available computational domain and stabilizes in a vortical structure with a big turbulent core.

## 4. Three dimensional simulations

In this section we leave the 2D setup and consider 3D ef- fects. We fix the nature of the thermal transfer to thermal re- laxation with a cooling time $\tau = 1$. Previous 3D simulations (Lesur & Papaloizou 2010; Lyra & Klahr 2011) were carried out using Boussinesq approximation and constant radial stratifica- tion. Here we use a fully compressible code first in an vertically unstratified configuration and then in a non-uniformly stratified configuration. In the z direction boundary conditions are set up by imposing the values derived from the equilibrium state of the disk.

### 4.1. Vertically unstratified disk

The continuous growth of the vortices by baroclinic amplifica- tion obviously raises the question of their evolution when their aspect-ratio reaches a value small enough to be affected by the
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Vorticity map in a 2D stratified disk, in the case of heat diffusion after 300, 400, 500, 700, 860 et 960 rotations, from left to right and from top to bottom. A hollow vortex forms and evolves in a vortex with a turbulent core.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Map of the heat transfert (Q) in a 2D stratified disk in the case of heat diffusion after 250 and 500 rotations, from left to right, respectively.</figcaption></figure>

elliptic instability (Lesur & Papaloizou 2009). In the case of un-stratified disks, LP10 found that vortices can survive this 3D instability if the baroclinic production of vorticity balances the mixing induced by the elliptical instability. Here, the problem is revisited in the case of compressible flow. In order to avoid the long (time consuming) integration necessary to produce a single and strong enough vortex we started the simulations from an approximate vortex solution. This solution is build up in a two steps procedure: first, a 2D baroclinic vortex is formed using the 2D version of the code and we let it grow until it reaches an aspect ratio $\chi = 4$; then, the 2D vortex is considered in a 3D setting by piling up the vortex profile in the third direction. The resulting vortex is a columnar vortex in an unstratified 3D disk. Figure 11 shows the evolution of a 3D baroclinic vortex constructed in such a way. It is clearly affected by the elliptical instability with some vertical wavy motions but the instability is never strong enough to completely destroy the vortex. We observe that the vortex re-organizes in a vortical structure with a turbulent core that can survive till the end of the simulation. The results we found are consistent with those of Lesur & Papaloizou (2010) and Lyra & Klahr (2011) who found the emergence of vertical motions inside the vortices when their aspect ratio satisfies $\chi < 4$.

## 4.2. Vertically stratified disk

We now consider a fully stratified 3D disk.

### 4.2.1. Radial stratification and layering

The disk structure is governed by a vertical equilibrium between gravity and pressure gradient

$$ \frac{\partial P}{\partial z} = -\rho \frac{z}{(r^2 + z^2)^{3/2}}. \qquad (19) $$

As in the 2D case, the radial profiles of the temperature and the surface density are assumed to be simple power laws. It is then easy to derive the gas density

$$ \rho_e = r^{-\beta} \exp\left(\frac{1}{T}\left(\frac{1}{\sqrt{r^2 + z^2}} - \frac{1}{r}\right)\right) \approx r^{-\beta} e^{-\frac{1}{2T}(z/r)^2} \qquad (20) $$

in which $\beta = p + (3 - q)/2$ depends on the two indices p and q. Radial stratification is characterized by the Brunt-Väisälä frequency

$$ N_r^2 = -\frac{RT}{\gamma} \left( \frac{\partial \ln \rho_e}{\partial r} - \frac{q}{r} \right) \left( (1-\gamma) \frac{\partial \ln \rho_e}{\partial r} - \frac{q}{r} \right) \qquad (21) $$

which in the thin disk approximation ($z/r \ll 1$) and using (10) leads to

$$ N_r^2 \approx \gamma(\gamma-1) \left(\frac{3-q}{2}\right)^2 \frac{RT}{r} (z^2 - z_2^2) (z^2 - z_1^2) \qquad (22) $$

where

$$ z_1 = H(r) \sqrt{\frac{2 \beta(\gamma-1) - q}{\gamma(\gamma-1) (3-q)}}, \qquad (23) $$

$$ z_2 = H(r) \sqrt{\frac{2 \beta + q}{\gamma (3-q)}}. \qquad (24) $$

So, Schuchardt criterion leads to four possible roots ($\pm z_1, \pm z_2$) that indicate the altitudes in the disk where gas stability is changing.

In the case of protoplanetary disks, $q < 3$ is a reasonable assumption; therefore, $z_2$ is always a real root but $z_1$ may be imaginary or real, depending on the sign of $\beta(\gamma - 1) - q$.

- If $\beta(\gamma - 1) - q < 0$, $z_1$ is imaginary: stratification is unstable in a mid-plane layer for $z < |z_2|$ and stable in the upper layer for $z > |z_2|$.

- If $\beta(\gamma - 1) - q > 0$, $z_1$ is real: stratification is stable, in both the mid-plane layer and the upper layer, but it is unstable in the intermediate layer for $|z_1| < z < |z_2|$.

This layering of the stability regions in a 3D disk contrasts with the uniform distribution of stability in a 2D disk in which $N_r^2$ keeps the same sign at any distance from the star. In the following, we restrict the analysis to the first case ($\beta(\gamma - 1) - q < 0$) for which the disk is unstable in a mid-plane layer.

We set the power law indices for the density and the temperature to: $\beta = 1$ and $q = 2$, respectively. For these values, we get $|z_2| = 2H$. The disk is then unstable for $z < 2H$ and stable for $z > 2H$. The computational domain is chosen to include both layers. It is defined in the non-dimensional variables as $7 < r < 8$, $0 < \theta < \pi/4$ and $0 < z < 1$ (that is $0 < z < 2.90H$ at $r = 7.5$). The numerical resolution we used in this case is $450 \times 300 \times 300$, that is $155 \times 18 \times 100$ cells by scale height in the radial, the azimuthal and the vertical directions, respectively.

As for the 2D simulations, the baroclinic instability is initiated by adding perturbations of sufficiently large amplitude to the equilibrium state.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Vorticity of a 3D baroclinic vortex with aspect ratio $\chi = 4$: at $t = 0$ and after 20, 22, 24 and 40 rotations; in a 3D perspective (top), in the $r-\theta$ plane (middle) and in the $\theta - z$ plane (bottom).</figcaption></figure>

### 4.2.2. Development of the instability

In a 3D disk the baroclinic instability is found to develop in nearly the same way as in a 2D disk. The first stage corresponding to the formation of small vortices from the density perturbations initially seeded in the background is recovered. These vortices mainly appear in the mid-plane of the disk like in the 2D case but they also possess a small vertical extent. Then, once formed, they tend to grow slowly in size and strength.

Figure 12 shows the evolution of the Rossby number of the strongest vortex formed at the end of the first stage. It clearly illustrates the vorticity amplification process occurring during the formation and growth of the vortex.

Figure 13 shows how vortices grow and also how they extend in the vertical direction. We observe the formation of slightly tilted columns of vorticity that progressively fill the whole disk thickness, including the stable upper layer. This spread of the vorticity in the stable layers of the disk is an interesting and new result. The aspect ratio of the vortices decreases as the vorticity increases, and reaches 10 at the end of the simulation. This large aspect ratio makes the vortex relatively stable with respect to the elliptical instability.

### 4.2.3. Formation of vortex layers

In this section, we focus on the weaker features observed on Fig. 13 after a few hundred rotations. We clearly see on the vorticity contours in ($r-z$) plane (middle row of Fig. 13) the presence of linear rays whose amplitude increases in time, especially in the upper-disk region. These rays exhibit an azimuthal structure which is visible on the 3D plots shown in Fig. 14. According to the inclination of the ray, we observe 2, 3 or 4 vortices, corresponding to the azimuthal wavenumbers $m = 16, 24$ or $32$ in our $\pi/4$ angular computation domain. The amplitude of the vorticity perturbation increases along the rays, and is found to be maximum at the boundary of the computation box, where it forms well-defined anticyclonic vortices (see Fig. 14 at 400 rotations) that tend to merge two by two (see Fig. 14 at 450 rotations).

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Rossby number as a function of time of the strongest vortex formed by baroclinic instability in a 3D disk.</figcaption></figure>

We think that these observations could correspond to the first step of what has been christened the Zombie Vortex Instability (ZVI) by Marcus et al. (2015). Barranco & Marcus (2005) were the first to observe the formation of off-midplane vortices during the destabilization of a mid-plane vortex in a stratified disk. Marcus et al. (2013b) attributed the formation of these vortices to the presence of baroclinic critical layers. They further showed that once formed these vortices could force new critical layers that create new vortices and so on by a self-replication process. It is this instability process which has been called ZVI. Although Marcus et al. (2015) demonstrated a certain robustness of the process, it has never been observed in a fully compressible simulation of 3D Keplerian disks.

We do not see any replication process in the present simulation, but claim that the rays do correspond to baroclinic critical layers. They clearly resemble the structures shown in Marcus et al. (2013a). The position of the baroclinic critical
<figure><img src="image_14.png" /><figcaption>Fig. 13. Vorticity map after 50, 250, 400 et 450 rotations, in the mid-plane (top), and in cuts through the strongest vortex in the $r-z$ plane (middle) and in the $\theta-z$ plane (bottom). The dashed line shows the height where $N_r^2$ become positive.</figcaption></figure>

layers is theoretically given by the relation

$$ \omega - m\Omega_D(r, z) = \pm N_z(r, z) \quad (25) $$

where $\omega$ is the frequency of the perturbation, $m$ its azimuthal wavenumber and $N_z$ the Brunt-Väisälä frequency in the vertical direction (see also the recent work of Umurhan 2016). With the present vertical structure of the disk and assuming that the mid-plane vortex is the source of the perturbations, we get

$$ m[\Omega_D(r_v, 0) - \Omega_D(r, z)] = \pm \sqrt{\gamma - 1} \Omega_K(r) \frac{z}{H(r)} \quad (26) $$

where $r_v$ is the radial position of the mid-plane vortex. For each value of $m$, we obtain two (almost linear) curves in the $(r-z)$ plane. In Fig. 15, we have plotted these curves for $m = 16, 24, 32$ on the vorticity contours of the perturbation. We obtain a quite good agreement between the inclinations of the vortex layers and the slopes of the critical layers derived analytically from Eq. (26). The azimuthal wavenumber of the theoretical critical layer also matches the observed azimuthal structure of the vortex layer: we do observe the formation of two, three and four vortices in the vortex layers associated with the critical layer of azimuthal wavenumber $m = 16, 24$, and 32 respectively (see Fig. 14). There is however an uncertainty on the position of the source which does not seem to be very localized. Several sources between 7.4 to 7.6 seem to contribute. This is perhaps related to the other vortex which is present in the simulation (see Fig. 14) but located at a different radius, outside the cross section plotted in Fig. 15. It then generates a different vortex-layer pattern, which should leave a trace in Fig. 15. This is clearly visible in Fig. 14: we do observe the formation of a double rows of two vortices and of three vortices corresponding to vortex layers of

each mid-plane vortex. This confirms that both the vortex layers and the vortices that form at the boundary of the computation box are not spurious numerical artifacts. They are real physical phenomena which deserve further attention.

Before going into deeper studies, it is straightforward to look at other changes in the disk associated with the presence of vortex layers. For example, we have plotted in Fig. 16 the distribution of density, pressure and velocity divergence in the same section as in Fig. 15. The most striking feature of this figure is the difference between the plots of density and pressure with those of velocity divergence and vorticity. The density/pressure contours are mainly limited to the unstable layer of the disk while the vorticity/divergence contours exhibit linear rays that extend through the whole thickness of the disk. The rays are also visible on the density contours in the mid-plane layer but they rapidly fade away with increasing amplitude. No rays are by contrast visible on the pressure contours.

## 5. Discussion and conclusion

In this paper, we have revisited the subcritical baroclinic instability. We have studied a number of unexplored situations, including 2D simulations with heat diffusion and 3D simulations in a completely stratified disk.

In the case of 2D disks, our simulations have confirmed the existence of two main stages in the evolution of the vortices: (i) a first stage during which small vortices form from density bumps seeded in the background; (ii) a second stage during which the small vortices grow and expand by the baroclinic amplification mechanism. We have considered two situations in which thermal transfer is due to either relaxation or heat diffusion. In the case of thermal relaxation, the evolution is found to be the same as
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** 3D distribution of the vorticity in a fully stratified disk, after 250, 400 and 445 rotations (from top to bottom, respectively). The colored arrows correspond to the radial (red), the azimuthal (green) and to the vertical (blue) directions.</figcaption></figure>

described by Petersen et al. (2007a,b). In the case of heat diffu-
sion, the evolution of the vorticity is more complex and leads to
the formation of hollow vortices. These vortices were found to
be unstable, but they were not destroyed. Instead, they gave rise
to coherent vortical structures with a turbulent core that persist
for a large number of rotation periods.

In the case of 3D disks, we have first analyzed the influ-
ence of the radial distributions of density and temperature on the
stratification properties of the disk. We have shown that the disk
could be either unstable in a mid-plane layer or in a layer dis-
tant from the mid-plane. We have examined the first situation
in which an unstable mid-plane layer is sandwiched between
two stable layers. We have found that, at the beginning of the
evolution, the baroclinic instability follows the 2D scenario, i.e.
small vortices are formed in the mid-plane and are then baro-
clinically amplified. However, during their growth, the vortices
were found to be stretched out in the vertical direction, leading
to columnar and long-lived vortex structures extending through
the whole width of the disk (including the stable layers).

It would be interesting to examine the second situation where
the unstable layer is far from the mid-plane. This situation is a
priori possible in the framework of the standard power law model
where *q* ≈ 0.5 and *p* > 0 (Andrews et al. 2009). It is possible
to imagine that vortices would first develop in the intermediate
unstable layer before being advected toward the mid-plane and
the upper layer.

The baroclinic instability that we have studied is non-linear
in nature. It can be triggered only in the presence of finite ampli-

tude perturbations. This requires strong enough perturbations of
the disk whose origin can only be speculated. It is unlikely that
they could be created by the magneto-rotational instability in the
vicinity of the dead-zone since this instability tends to inhibit
vortex formation (Lyra & Klahr 2011). However, they could be
due to another hydrodynamical instability such as the linear ver-
sion of the baroclinic instability proposed by Klahr & Hubbard
(2014). The vertical shear instability was proposed by Nelson
et al. (2013) as a possible way to trigger the SBI, but recently,
Richard et al. (2016) showed that the vertical shear instability
leads to the formation of vortices via the Rossby wave instabil-
ity rather than the SBI.

The recent discovery by ALMA of an asymmetry in the
dust distribution around Oph IRS 48 (van der Marel et al. 2013;
Armitage 2013) has been interpreted as an assembly of particles
captured by a large-scale Rossby vortex. The authors suggested
that the vortex had formed by a Rossby wave instability occur-
ring at the edge of the gap opened by an unseen companion or-
biting the inner disk. However, it cannot be excluded that such a
vortex is the result of a baroclinic instability.

Finally, we have also shown that 3D baroclinic vortices are
sources of internal vortex layers. These layers develop on baro-
clinic critical layer surfaces and form vortices in the upper-disk
region. Whether these vortices could force other vortex lay-
ers and follow the self-replication mechanism of Marcus et al.
(2013b) remains among the open questions that would be inter-
esting to answer.
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Vorticity map after 445 rotations in a $r-z$ cross section through the strongest vortex (same plot as in Fig. 13 but the contrast has been enhanced to make more salient weak features. The horizontal line (black dashed) indicates the height where $N_r^2$ becomes positive. The colored lines (dark-blue, cyan and green) indicate the position of the baroclinic critical layers, as computed from Eq. (26) with a source at the center of the mid-plane vortex for $m = 16, 24$ and $32$, respectively. The dashed-colored lines are other critical layer curves with different positions of the source.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Maps of the pressure, the density and the divergence of the flow after 445 rotations in a $r-z$ cross section through the strongest vortex observed in the simulation. $P - \rho_e T(r)$ is plotted on the left, $\rho - \rho_e$ is plotted on the middle and $0.5\nabla \cdot V - 1/4$ on the right.</figcaption></figure>

*Acknowledgements.* Computations were performed on a Bull multicore machine funded by AMU (Aix-Marseille University) and LAM (Laboratoire d'Astrophysique de Marseille); the most demanding computations were performed using HPC resources from GENCI [TGCC - (grant2013-2014) - t2013046832 - t2014046832 - t2015047407]. 3D visualizations were performed with *glnemo2* developed at CeSAM - http://projets.lam.fr/projects/glnemo2. Support by LABEX MEC (ANR-11-LABX-0092) and ANR LIPSTIC (ANR-13-JS05-0004-01) is also acknowledged. We thank an anonymous referee for the review of our manuscript.

## References

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502

Armitage, P. J. 2013, Science, 340, 1179

Barge, P., & Sommeria, J. 1995, A&A, 295, L1

Barranco, J. A., & Marcus, P. S. 2005, ApJ, 623, 1157

Inaba, S., Barge, P., Daniel, E., & Guillard, H. 2005, A&A, 431, 365

Johnson, B. M., & Gammie, C. F. 2005, ApJ, 626, 978

Kida, S. 1981, J. Phys. Soc. Jap., 50, 3517

Klahr, H. 2004, ApJ, 606, 1070

Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21

Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869

Lesur, G., & Papaloizou, J. C. B. 2010, A&A, 513, A60

Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023

Li, H., Colgate, S. A., Wendroff, B., & Liska, R. 2001, ApJ, 551, 874

Lin, M.-K. 2012, ApJ, 754, 21

Lin, M.-K. 2013, ApJ, 765, 84

Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805

Lyra, W. 2014, ApJ, 789, 77

Lyra, W., & Klahr, H. 2011, A&A, 527, A138

Marcus, P. S., Jiang, C.-H., Pei, S., & Hassanzadeh, P. 2013a, in Eur. Phys. J. Web Conf., 46, 3006

Marcus, P. S., Pei, S., Jiang, C.-H., & Hassanzadeh, P. 2013b, Phys. Rev. Lett., 111, 084501

Marcus, P. S., Pei, S., Jiang, C.-H., et al. 2015, ApJ, 808, 87

Meheut, H., Keppens, R., Casse, F., & Benz, W. 2012a, A&A, 542, A9

Meheut, H., Yu, C., & Lai, D. 2012b, MNRAS, 422, 2399

Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610

Pedlosky, J. 1987, Geophys. Fluid Dynamics (Springer-Verlag)

Petersen, M. R., Julien, K., & Stewart, G. R. 2007a, ApJ, 658, 1236

Petersen, M. R., Stewart, G. R., & Julien, K. 2007b, ApJ, 658, 1252

Pringle, J. 1981, ARA&A, 19, 137

Raettig, N., Lyra, W., & Klahr, H. 2013, ApJ, 765, 115

Richard, S., Barge, P., & Le Dizès, S. 2013, A&A, 559, A30

Umurhan, O. M., Shariff, K., & Cuzzi, J. N. 2016, ApJ, submitted [arXiv:1601.00382]

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199