© ESO, 2016

1. Introduction

Large-scale vortices in protoplanetary disks could play an important 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 identified by Barge & Sommeria (1995). This is the reason why the formation, stability, and evolution of such vortices in protoplanetary disks deserve special attention. Recent observations with high resolution in the millimeter and centimeter range, have recently raised a strong interest on dust trapping and vortices. Vortex formation in protoplanetary disks has been studied a lot in the last ten years. Their formation by the Rossby-wave instability 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 mechanisms (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 respect 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 vortices is unclear since thermal transfer is absent from their simulations, contrary to what is normally required (Petersen et al. 2007a,b). The non-linear nature of this instability was underlined by Lesur & Papaloizou (2010, thereafter LP10) using incompressible 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 generation of baroclinic turbulence in protoplanetary disks. They 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 unstratified 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-homentropic 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 μ = 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)/Ω_K(r) where is the sound speed and Ω_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

where ρ and P are the density and pressure of the gas; V is the gas velocity and e its total specific energy; φ is the gravitational potential, Q the heat transfer and γ = 1.4 the adiabatic index.

The equations will be used in non-dimensional form thanks to the following normalization 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 ρ_o = ρ(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 φ = 1/(r² + z²)^1/2. The important parameter is the Mach number 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

where p = 1.5 and q = 0.5 in the case of the standard MMSN hypothesis. The perfect gas equation of state gives (8)The radial velocity is assumed to be null and the azimuthal velocity is given by the radial hydrodynamic equilibrium (9)For large Mach number, the disk angular velocity Ω_D = v_D/r is close to the Kepler profile Ω_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_Ao. 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 (10)

2.4. Radial stratification

In the radial direction, disk stratification is measured by the Brunt-Väisälä frequency (11)and the Schwarzschild criterion states that this stratification is stable if and unstable if .

In a 2D disk at equilibrium, pressure and density can be expressed as simple power-laws and the Brunt-Väisälä frequency writes (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/ (γ − 1). In term of the Richardson number , where S = −3Ω_K/ 2 is the Kepler shear, the condition of instability is then equivalent to (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₀ at a rate (14)where τ 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 (15)where Δ 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 (16)where κ is a diffusion coefficient and r_∗ the mean orbital radius of the computational box (in what follows we will use Pe/ρ(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: ν_a at large scales and ν_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 ΔΣ/(ΔtΣΩ) ≃ 3.2 × 10^-6 so that (or Re ≃ 5.2 × 10⁵). The second one scales as where δ_∗ is the mesh size; in the radial direction, this estimate gives ν_n ≃ 10¹⁰m²s^-1 (or Re ≃ 5.2 × 10⁴).

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 non-linear 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_Ao = 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 <θ<π/ 4. The numerical resolution was set up to 300 × 300, that is 465 cells by scale height in the radial direction and 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 ≈ 2 × 10^-3 at r = 7.5) whereas for p = 1, stratification is unstable since the Richardson number is negative (Ri ≈ −1.5 × 10^-3 at r = 7.5) . Baroclinic feedback is expected only in the case of unstable stratification that is when p<q/ (γ − 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 approximatively 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 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 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 . 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 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 τ = ∞ or τ = 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 (τ = 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 τ () 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 (17)where ω_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 χ 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 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 χ 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 parameters, nearly independent of the cooling time. It is well fitted by the analytical expression derived for Kida’s vortex model (Kida 1981) (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 computational domain has been chosen to be a box with 7 <r< 8 and 0 <θ<π/ 2 except in the simulation presented in Fig. 9 in which the box is extended to 6.5 <r< 8.5 and 0 <θ<π . The results 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 × 10⁴). However, as for thermal relaxation, the amplification is weakened for very large Pe and disappears at Pe = ∞. There is therefore an optimal heat diffusion for the amplification which corresponds to Pe ≈ 10⁵.

The main differences between an evolution under thermal relaxation and under heat diffusion clearly appear during the second stage, after the vortices have grown and merged into a single 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 negative 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 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 inner side (the side of the star) and cools down when moving azimuthally 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 periods. 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 hollow 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 vortex 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 effects. We fixe the nature of the thermal transfer to thermal relaxation with a cooling time τ = 1. Previous 3D simulations (Lesur & Papaloizou 2010; Lyra & Klahr 2011) were carried out using Boussinesq approximation and constant radial stratification. 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 amplification obviously raises the question of their evolution when their aspect-ratio reaches a value small enough to be affected by the elliptic instability (Lesur & Papaloizou 2009). In the case of unstratified 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 χ = 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 χ< 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 (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 (20)in which β = p + (3 − q)/2 depends on the two indices p and q. Radial stratification is characterized by the Brunt-Väisälä frequency (21)which in the thin disk approximation (z/r ≪ 1) and using (10) leads to (22)where

So, Schwarchild criterion leads to four possible roots (± z₁, ± z₂) 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₂ is always a real root but z₁ may be imaginary or real, depending on the sign of β(γ − 1) − q.

• If β(γ − 1) − q< 0 , z₁ is imaginary: stratification is unstable in a mid-plane layer for z< | z₂ |and stable in the upper layer for z> | z₂ |.

• If β(γ − 1) − q> 0 , z₁ 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₁ | <z< | z₂ |.

This layering of the stability regions in a 3D disk contrasts with the uniform distribution of stability in a 2D disk in which keeps the same sign at any distance from the star. In the following, we restrict the analysis to the first case (β(γ − 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: β = 1 and q = 2, respectively. For these values, we get | z₂ | = 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 <θ<π/ 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 × 300 × 300, that is 155 × 18 × 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.

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 expend 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 π/ 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).

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 layers is theoretically given by the relation (25)where ω 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 (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 described by Petersen et al. (2007a,b). In the case of heat diffusion, 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 influence 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 distant 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 baroclinically 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 amplitude 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 version 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 instability 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 occurring at the edge of the gap opened by an unseen companion orbiting 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 baroclinic critical layer surfaces and form vortices in the upper-disk region. Whether these vortices could force other vortex layers and follow the self-replication mechanism of Marcus et al. (2013b) remains among the open questions that would be interesting to answer.

Acknowledgments

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 ressources 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

All Figures