Free Access
Issue
A&A
Volume 573, January 2015
Article Number A132
Number of page(s) 12
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201424162
Published online 14 January 2015

© ESO, 2015

1. Introduction

In current planet formation theory, the agglomeration of microscopic dust into kilometre-sized objects remains poorly understood (Chiang & Youdin 2010). The direct gravitational collapse scenario (Goldreich & Ward 1973; Weidenschilling 1980; Nakagawa et al. 1981; Hayashi et al. 1985) requires a midplane dust-to-gas ratio larger than the interstellar medium value by about three orders of magnitude, which is not easily achievable by sedimentation alone (Weidenschilling & Cuzzi 1993; Sekiya 1998; Turner et al. 2010). Growth resulting from ballistic collisions and sticking works efficiently for small particles (Chokshi et al. 1993; Dominik & Tielens 1997; Poppe et al. 2000), but laboratory and numerical experiments suggest two limiting grain sizes (~1 mm at 1 AU), known as the bouncing and fragmentation barriers, beyond which colliding particles do not stick (Weidenschilling & Cuzzi 1993; Benz 2000; Kothe et al. 2010; Wurm et al. 2005; Zsom et al. 2010; Windmark et al. 2012; Seizinger & Kley 2013; Garaud et al. 2013). We note, however, the recent calculations of the evolution of a dust population presented by Garaud et al. (2013), who argue that these barriers are less insurmountable than first thought. For example, the incorporation of a realistic particle velocity distribution ameliorates the problem, as does the realization that collisions between particles of highly different mass can transfer mass to the largest of two impactors.

A separate problem is that the solid content of a protoplanetary (PP) disk rapidly drains out of the disk (Stepinski & Valageas 1997; Kornet et al. 2001; Takeuchi & Lin 2002; Brauer et al. 2007; Hughes & Armitage 2012) because the gas is partly supported by the disk’s pressure gradient and rotates at slightly sub-Keplerian frequencies, while dust grains rotate faster at the Keplerian angular velocity. Consequently, dust grains feel a head-wind and thereby lose their angular momentum, ultimately spiralling into the central star. The inward drift velocity is highest for particles of a few tens of centimetres and such particles only survive for about 100 years (Weidenschilling 1977). This phenomenon constitutes an additional barrier known as the radial-drift barrier.

Among the many scenarios that have been discussed to overcome these problems, some of the most promising involve the pressure bumps formed at planetary gap edges (de Val-Borro et al. 2007; Lin & Papaloizou 2011) or at the interface between a laminar (dead) and a turbulent (active) region (Kretke & Lin 2007; Brauer et al. 2008; Kretke et al. 2009). Friction between gas and dust particles vanishes at the location of the pressure maximum, naturally providing a trap where the disk’s solid content can accumulate. In addition, such pressure extrema are unstable to the vortex or Rossby wave instability (RWI, Lovelace et al. 1999; Li et al. 2000; Meheut et al. 2010; Lin 2012). It is well known that vortices significantly concentrate dust in a short time (Barge & Sommeria 1995; Tanga et al. 1996; Johansen et al. 2004; Klahr & Bodenheimer 2006; Inaba & Barge 2006; Meheut et al. 2012b,a), thus promoting planet formation at these locations. We finally note that the presence of vortices has been suggested by high angular resolution imaging of PP disks that exhibit significant asymmetric features (Brown et al. 2009; Casassus et al. 2012; van der Marel et al. 2013; Isella et al. 2013; Pérez et al. 2014).

In this context, the dead zone inner edge of PP disks is a special location that deserves detailed investigation. We define the inner edge to be the radius where the radially decreasing midplane temperature drops below the critical value at which thermal ionization fails to sustain the magnetorotational instability (MRI, Balbus & Hawley 1991, 1998). As a result, the flow is turbulent inward of that interface, but laminar beyond it (Gammie 1996). Recent 3D magnetohydrodynamics (MHD) simulations have reported the formation of pressure bumps in locally isothermal models of PP disks (Dzyurkevich et al. 2010; Lyra & Mac Low 2012), confirming earlier 2D hydrodynamical viscous simulations. (Varnière & Tagger 2006; Lyra et al. 2008; Regály et al. 2012). This work has established the feasibility of the dead zone inner edge as a trap for solids; however, the robustness of these results is curtailed by the use of simple thermal physics, namely a locally isothermal equation of state. This assumption is problematic because of the pervasive interpenetration of dynamics and thermodynamics in this region, especially at the midplane. Temperature depends on the turbulence via the dissipation of its kinetic and magnetic fluctuations, but the MRI turbulence, in turn, depends on the temperature through the ionization fraction, which is determined thermally via dissipation (Pneuman & Mitchell 1965; Umebayashi & Nakano 1988). In addition, the onset and non-linear evolution of the RWI should depend on the PP disk’s global thermal structure through its radial potential-vorticity profile (Umurhan 2010).

Using a simplified mean field model, Latter & Balbus (2012) found that, if the interplay between thermal and turbulent dynamics is taken into account, the dead-active interface is not static, but rather moves radially before stalling at a well-defined radius. Recently, we confirmed this behaviour using MHD simulations that self-consistently accounted for both turbulent heating and the feedback between temperature and magnetic diffusivity (Faure et al. 2014). However, in order to reduce the computational cost of these simulations, vortex formation was artificially inhibited by using a reduced azimuthal domain. The point of the present paper is to examine the onset and development of vortices in thermally structured models of PP disks. We do so by increasing the azimuthal extent of our previous simulations.

In order to isolate and understand the basic physics we adopt simplified geometry and thermodynamics. Our PP disk is cylindrical, and hence vertically unstratified. Consequently, disk cooling is approximated by a cooling law rather than via a detailed radiative transport model. As expected, we find that a vortex forms at the dead zone inner edge via the RWI, but contrary to expectations (e.g. Paardekooper et al. 2010) the vortex radially migrates inwards, ploughing through the pressure bump and into the active zone where it is ultimately destroyed by turbulent motions. A few orbits later a new vortex forms at the pressure bump and it too follows the same cycle of formation, migration, and disruption. This behaviour fails to appear in isothermal or adiabatic runs, and seems connected to the details of the PP disk’s heating and cooling. It is yet unclear how robust this vortex cycle is, in particular how sensitive it is to the approximate cooling law we have employed. Future vertically structured global simulations using the set-up of Flock et al. (2013) will aid in testing this. Taken at face value, however, these results complicate planet formation theories that appeal to dust trapping at the dead zone inner edge unless a fast formation process is at work within vortices.

The paper is organized as follows. In Sect. 2 we present the results of a 3D MHD simulation which exhibits the basic vortex cycle. Section 3 contains the results of a 2D run of a simple viscous model that reproduces the vortex behaviour observed in the 3D simulation. This simplified model is a useful tool with which to analyse the vortex cycle in more detail. In Sect. 4, we discuss the physical mechanisms at work in both simulations. Finally, we conclude and speculate on the implications of our simulations for planet formation in Sect. 5.

2. The 3D turbulent MHD simulation

We first present results from a 3D simulation subject to MRI-induced turbulence. We focus on an annular region centred around the dead zone inner edge, as in Faure et al. (2014), to which the reader is referred for further details. However, by extending the azimuthal extent of the domain we observe the onset and development of the RWI. The behaviour of the vortices so formed is described in this section and comprises the main result of the paper.

2.1. Set-up

The 3D simulation was performed using a uniform grid version of the code RAMSES (Teyssier 2002; Fromang et al. 2006). We solved the MHD equations in cylindrical coordinates (R,φ,Z) with units vectors (eR,eφ,ez)where ρ is the density, v is the velocity, B is the magnetic field, P is the pressure, and E is the sum of kinetic, magnetic, and thermal energy. The parameter Φ is the gravitational potential. In the cylindrical approximation, it is given by Φ = −GM/R where G is the gravitational constant and M is the stellar mass. We used a perfect gas equation of state to close the above set of equations. The thermal energy is related to the pressure through the relation eth = P/ (γ − 1) in which γ = 1.4. The magnetic diffusivity is denoted by η and is responsible for the resistive flux η that appears in Eq. (3) (see Balbus & Hawley 1998). We captured turbulent heating by solving the equation of total energy conservation and used the same gas cooling function as in Faure et al. (2014), (5)where T is temperature, and Tmin is the temperature associated with radiative equilibrium. We modelled the rapid variation of the resistivity η with the temperature by a step function, (6)where TMRI is the activation temperature for the MRI, typically ~1000 K.

The initial magnetic field configuration is purely toroidal. Its vertical profile is such that the integrated magnetic flux through a vertical slice of the disk vanishes. The computational domain is R ∈ [ R0,8 R0 ], φ ∈ [ 0/ 2 ], and Z ∈ [ −0.3 R0,0.3 R0 ] and has a resolution of [ 320,160,80 ].

The main difference with the set-up of Faure et al. (2014) is the removal of the source term in the continuity equation. We found that the term could interfere with the development of density features and hence with the RWI. In addition, we tuned the radial profile of the viscosity in the buffer zones in order to avoid accretion discontinuities at the buffer edges. For completeness we give here the functional form of the viscosity profile we used, (7)where αbuff is the value of α measured at R = Rbuff and averaged since the beginning of the simulation. The parameter Rbuff is the position of the boundary between the buffer zone and the active domain. Together with free-flow boundary conditions this technique prevents unphysical mass depletion in the domain.

2.2. Notation and units

In the following we denote by X0 the value of any quantity X at the inner edge of the domain. Units are identical to those defined in Faure et al. (2014) and chosen such that where Ω stands for the gas angular velocity at radius R, and ρ is density. Time is measured in units of the inner orbital period.

Density and temperature profiles are initialized with radial power laws, (8)where p = −1.5 and q = −0.75.

In the cooling law, Eq. (5), we chose σ = 1.1 × 10-4. For a typical simulation, this yields a disk aspect ratio of H0/R0 ~ 0.1 and a disk cooling time of about 25 local orbits. Finally, we used η0 = 10-3, which is large enough to prevent the development of the MRI in resistive regions, and TMRI = 0.6.

2.3. Results

We forced the simulation to undergo three consecutive phases, with a different resistivity configuration in each. This approach introduces the requisite physics in a controlled fashion. We describe each step in turn.

2.3.1. The ideal phase

Initially the resistivity is set to zero everywhere except in the buffer zones. The aim is to obtain a fully turbulent disk free of unnecessary transient behaviour associated with the initialization of the MRI. Once this is achieved we then introduce a dead zone in the following phase. The model equations are integrated for a few cooling times and until the temperature reaches a quasi-steady state. The radial profile of the equilibrium temperature that results is plotted in Fig. 1.

2.3.2. The static dead-zone phase

Once the simulation has reached a quasi-steady equilibrium (t = 600) we set η = η0 over the region 3.5 <R< 8 and integrate the equations for another 500 orbits when thermal equilibrium is reached. Magnetohydrodynamics turbulence quickly dies outward of R = 3.5 and a static dead zone forms. Gas cools in this region as a result of the absence of any turbulent heating. Soon, because of the accretion mismatch, we notice the formation of density and pressure maxima around the interface at R = 3.5, and this saturates when the density reaches about twice its initial value. In a few tens of orbits a vortex forms (see the first panel of Fig. 2) at this location because of the RWI.

2.3.3. The self-consistent dead-zone phase

thumbnail Fig. 1

Temperature profile at the end of the ideal phase averaged over 200 inner orbits. The black dashed lines show where the mean temperature is equal to (5 / 4)TMRI. The hashed regions on both sides of the figure correspond to the buffer zones.

Open with DEXTER

thumbnail Fig. 2

From left to right and top to bottom: snapshot of the relative vorticity perturbation at the 1087th, 1124th, 1173th, 1208th, 1216th, 1225th, 1228th, and 1310th orbit. The relative vorticity perturbation has been vertically averaged. We draw with a black line the iso-contour corresponding to 80% of the density maximum ρm. The location of the density maximum is shown by the black cross.

Open with DEXTER

In the final phase of the simulation we restart the simulation at t = 1100 inner orbits and close the feedback loop; the resistivity is now a function of temperature according to Eq. (6). With this configuration, Faure et al. (2014) have shown that the dead zone inner edge moves radially before stalling at a critical radius where the gas temperature is (9)This critical temperature is shown in Fig. 1 and gives Rc = 3.7. In practice, we found that the dead-active interface initially located at R = 3.5 remains at that position over the simulation.

We found that the vortex formed during the static dead-zone step immediately begins to migrate inward. Its behaviour is illustrated by the sequence of successive snapshots shown in Fig. 2 where we plot the vertically integrated vertical relative vorticity δωz of the gas (10)where vK stands for the Keplerian velocity1. The first five panels of Fig. 2 clearly show the migration of the vortex into the active zone. We also note that as it moves toward the star, the vortex becomes smaller and smaller until it disappears after 150 orbits (panel 6). By this stage a new vortex has begun to form at the dead zone inner edge (panel 7) and begins its inward migration about ten inner orbits later (panel 8). The successive snapshots (not shown) indicate that it will also penetrate into the active zone before being similarly disrupted.

In Fig. 3, we plot the evolution of the radial position of the vortex Rvort with time. This is calculated from the position of the density maximum rather than the vorticity minima because the turbulent fluctuations in the active zone can complicate the identification of the vortex. Figure 3 shows that vortices follow a cycle of formation, migration, and disruption with a period of about 200 inner orbits. Vortices migrate inwards with a velocity ~10-3RfΩf (where Rf and Ωf stand for the vortex formation radius and the angular velocity at that location) from R ~ 3.7 to R ~ 2.6, which corresponds to about 3 scale heights at that location.

By definition, the temperature within the vortex at the formation location is below the MRI activation threshold (simply because the vortex forms in the dead zone) and the flow is relatively laminar. The heating rate inside the vortex is consequently low while the cooling rate is slightly increased because it is denser than its surroundings. The combination of such reduced heating and increased cooling rates induces a temperature decrease within the vortex during its migration, as shown in Fig. 4. This causes the temperature in the vortex to remain lower than TMRI during its lifetime. The typical cooling rate is of the order of 4.5 × 10-3 TMRIΩ0, which results in the vortex temperature decreasing by about 20% by the time it disrupts. We also note the secular decrease in the vortex temperature over many cycles. This indicates that the vortices’ properties are still evolving cycle after cycle, which raises the question of the existence of a stationary regime. We come back to this question in Sect. 3 where we perform 2D long-term simulations.

thumbnail Fig. 3

Radial position of the density maximum in the 3D run during the static dead-zone phase (in black) and the self-consistent phase (red curve).

Open with DEXTER

thumbnail Fig. 4

Temperature averaged over the vortex area in the 3D run during the static dead-zone phase (in black) and the self-consistent phase (red curve).

Open with DEXTER

thumbnail Fig. 5

Snapshot of the relative density perturbation at the 1310th orbit. The relative density perturbation has been vertically averaged. The iso-contour of 80% of density maximum ρm is drawn with a black line. The location of the density maximum is shown by the black cross.

Open with DEXTER

Finally, we note that the temperature in the dead zone is strongly correlated with the presence of the vortex. The dead zone’s temperatures increase from T ~ 0.4, when the vortex is weak at the end of its lifespan, to T ~ 0.7 just after a new one has formed. This is due to the presence of a pattern of strong density waves excited by the vortex (see the map of the density perturbations δρ/ ⟨ ρ in Fig. 5)2. The density increase at the crest of these waves is significantly larger (~50% of the local density) than described in Faure et al. (2014) in the absence of a vortex (only ~10% of the local density). These waves provide an additional source of heating and explain the warmer dead zone temperatures we obtain here.

3. The 2D simulation

In the 3D simulation, the flow is complex and the different physics – turbulence, MHD, and thermodynamics – are difficult to disentangle. In addition, the large computational cost associated with such simulations precludes long-term integration. For example, the temperature evolution shows a systematic drift over a few cycles (see Fig. 4), thus raising the question of whether the system is able to settle into a quasi-steady state, or whether the cycle described in the previous section is transient. In order to answer these questions, we present a non-magnetized laminar 2D viscous simulation which also reproduces the vortex cycles. The simpler set-up eases the interpretation of the results and their long-term relevance, and highlights potential limitations to our simulations. In this section, we briefly present the 2D model before describing in detail the formation, migration, and disruption phases of a typical cycle in that situation.

3.1. Set-up and run parameters

We performed a purely hydrodynamic laminar 2D simulation of a PP disk, similar in every other way to the 3D simulations in Sect. 2. The evolution of the disk model described below is calculated using the uniform grid version of the RAMSES code. The equations we solve are where τ is the Navier-Stokes stress tensor and Σ the disk surface density. We use the viscous prescription as a crude model of both turbulent angular momentum transport and heating. The kinematic viscosity radial profile is calculated using (14)where αt is the α profile azimuthally, vertically, and time averaged over the last 200 orbits of the 3D static dead-zone step. Here, cs and Ω are the local sound speed and the angular frequency. In the 3D static dead-zone step the vortex launches waves that significantly increase the angular momentum transport in the dead zone (see Faure et al. 2014). Hence the αt profile has two contributions coming from the turbulence and from the waves. In order to remove the wave contribution in the 2D simulations, which will emerge self-consistently from the vortex, we take αt = 0 outward of R = 3.5. In the disk active zone, αt steadily decreases from about 0.04 at R ~ 2 to 0.01 at R = 3.5. The viscosity jump at R = 3.5 leads to the formation of a pressure and a density bump at the outer edge of the viscous region (or, equivalently, at the inner edge of the inviscid or dead zone). We checked that the bump mass-loading rate is identical to that obtained in the 3D simulations, thus validating the α-prescription. This corresponds to the static dead-zone phase in Sect. 2.

As soon as the density bump has reached five times the initial local density (which occurs at t = 1000), we add a random velocity perturbation at the position of the bump to trigger the RWI and then close the feedback loop between resistivity and temperature by imposing (15)The perturbation amplitude equals 10% of the sound speed so that it mimics the velocity fluctuations induced by the turbulence on the bump. The influence of the initial bump size on the results is discussed in Sect. 3.3. We have entered the analogue of the self-consistent dead-zone phase described in Sect. 2. A vortex forms at R ~ 4 in a few tens of orbits. In the following section we concentrate on a detailed analysis of this phase, for which the idealized conditions of the 2D simulations provide a favourable environment.

3.2. Vortex formation

thumbnail Fig. 6

Radial profiles of at the 1000th orbit (red curve) and after the simulation has reached a quasi-steady state (black curve). is averaged azimuthally for both cases.

Open with DEXTER

thumbnail Fig. 7

Evolution of the logarithm of the vorticity perturbation δω = ωmax(Rvort) − ωmin(Rvort) at the bump location in the 2D viscous simulation. The best fit of the linear part of the amplitude growth is plotted in black and has a slope equal to γ/ Ω(Rf) = 0.2.

Open with DEXTER

The origin of the vortex forming at the dead/active interface is most likely due to the RWI (Lyra & Mac Low 2012; Regály et al. 2012). A 2D adiabatic disk is unstable to the RWI where (16)a conserved quantity similar to the potential vorticity, possesses a local extremum (Lovelace et al. 1999). Here Σ is the equilibrium surface density and γ is the adiabatic index. In Fig. 6 we show the profile of at t = 1000 (i.e. when random velocity fluctuations are added to the flow). Aside from small-scale variations, possesses a deep and long-lasting trough at the vortex location, strong evidence that the RWI is active and generates the vortex. We note that our simulations are diabatic, and hence do not strictly conserve , and that this does not impinge on (and in fact probably aids) instability here. In the 3D simulation, the vortex and the bump form simultaneously which complicates a clear identification of the RWI’s linear phase. Nonetheless, the profile of at the beginning of every cycle shows a similar trough at the position where the next vortex will form (not shown). Taken together, the results of the 2D and 3D simulations consistently point toward the RWI as being responsible for the vortex growth in our simulations.

The growth of vorticity perturbations are presented in Fig. 7. The perturbation grows exponentially with a growth rate of about γ/ Ω(Rf) ~ 0.2. For the density bump properties (height and width), the empirical law of Meheut et al. (2013) gives a growth rate about twice as large, such that γth/ Ω(Rf) ~ 0.5. We note, however, that this relation was obtained for globally isothermal disks with a small density perturbation such as δρ/ρ < 0.3. In our 2D simulation the density bump is 10 times larger and the temperature gradient at the inner edge of the dead zone is significant. The empirical relation still provides a useful rough estimate and sanity check.

When the period of the circular motion around the vortex centre becomes comparable to the growth timescale (at t ~ 1100) the growth ends. The vorticity perturbation is then δω ≃ −0.05Ω0 ≃ 0.4Ω(Rf), corresponding to twice the growth rate, as noticed by Meheut et al. (2013) in their simulations. By this stage the vortex has already begun to migrate while capturing an important fraction of the initial bump’s mass, (17)thus leaving about 60% of the bump mass behind at R = 4.

3.3. Vortex migration and disruption

thumbnail Fig. 8

Radial position of the density maximum in the 2D viscous simulation using the RAMSES code.

Open with DEXTER

thumbnail Fig. 9

From left to right and top to bottom: snapshot of the relative vorticity perturbation at the t = 2502, 2568, 2626, and 2668. The relative vorticity perturbation has been vertically averaged. The iso-contour of 80% of density maximum ρm is drawn with a black line. The location of the density maximum is shown by the black cross.

Open with DEXTER

After the first appearance of a vortex, we observe during the remainder of the 2D simulation a vortex life cycle similar to that obtained in the 3D simulation. As in the 3D simulation, the cycle properties change over the first few cycles, before the disk relaxes into a well-defined quasi-periodic state (see Fig. 8). During that quasi-periodic evolution, the profile of at the beginning of a cycle has converged to the black curve in Fig. 6. The existence of a minimum of at a well-defined radius during the remainder of the disk evolution shows that the conditions for vortex formation persist to late times, presumably driven by the continuing accretion mismatch at this point. The important difference is that the profile of the black curve is self-consistently obtained as a result of the disk evolution. Additional simulations starting with different bump sizes were found to exhibit different relaxation phases, but all converged to the same quasi-periodic state. These results, obtained using the 2D simulation, strongly suggest that the 3D simulation described in Sect. 2 would also evolve toward a quasi-stationary state if evolved for longer. Indeed, during this relaxation phase, the vortex migration range is similar to that observed in the 3D simulation. Forming at R ~ 3.7, vortices migrate inwards to R ~ 2.7. However, such a similarity does not extend to all the simulation diagnostics. For example, the cycle period is two times shorter and vortices migrate two times faster than vortices in the 3D run. We come back to this point in the Sect. 4.

We now investigate the properties of the quasi-periodic phase by focusing on one cycle, starting at t = 2500. At that time, a vortex has formed and is about to start migrating inwards. Its relative density perturbation at this point is ~0.4. Figure 9 presents a series of snapshots showing the vorticity in the disk at four different times. The first panel illustrates the state of the disk at the beginning of the vortex: it is clear that a vortex exists at R = 3.7. In panels 2 and 3 (50 and 100 orbits later) the vortex has drifted closer to the star and has shrunk as it does so. Finally, the last panel shows that the vortex has disappeared, but a new vortex has formed at a larger radii at the interface. Overall, the evolution is very similar to that obtained in the 3D simulation (see Fig. 2). As shown in Fig. 10, the similarity extends to the temperature evolution of the vortex: as in the 3D simulation, the vortex cools while it moves to hotter region. The cooling rate during the stationary cycles is found to be 6 × 10-3 TMRIΩ0, i.e. slightly larger than the value obtained in the 3D simulations. Finally, we show in Fig. 11 the evolution of the mass trapped inside the vortex. As suggested by the snapshots shown in Fig. 9, the vortex loses mass during the cycle.

3.4. A numerical test

In order to check the reliability of our simulations, we have reproduced the simulation presented in this section using the PLUTO code (Mignone et al. 2007). For this run, we used a simplified function for αt, (18)where R1 = 2.65R0, R2 = 3.44R0, αt = 7 × 10-2, and αe = 3.5 × 10-3 but otherwise solve the same set of Eqs. (11)–(13). Figure 12 shows the position of the vortex that forms in this simulation. The vortex cycle is clearly reproduced with a similar vortex migration timescale and amplitude as obtained with RAMSES. The quantitative differences between the PLUTO and the RAMSES simulations probably arise because of the ad hoc radial profile of αt we used in PLUTO that is slightly different from the profile used in RAMSES. The difference is, however, small and the similarity between the two simulations strengthens our main result: a vortex cycle emerges from the interplay between dynamics and thermodynamics at the inner edge of the dead zone.

4. Discussion

In this section we use the results of the previous two sections to illuminate the physical mechanism responsible for the vortex cycle.

4.1. Vortex destruction

thumbnail Fig. 10

Evolution of the mean temperature inside the vortex during the cycle starting at the t = 2500th inner orbit.

Open with DEXTER

thumbnail Fig. 11

Evolution of the vortex total mass over the cycle starting at t = 2500, given in units of the initial disk mass. It is found to decrease steadily during the vortex migration.

Open with DEXTER

thumbnail Fig. 12

Radial position of the density maximum in the 2D viscous simulation using the PLUTO code.

Open with DEXTER

First we address the issue of vortex disruption, which is a little easier to understand. In both the 2D hydro and 3D MHD simulations, the vortex forms in the dead zone (modelled either as a highly resistive zone or as an inviscid region) before migrating into the active region (turbulent in the 3D case and viscous in the 2D simulation) where it gradually disrupts. Throughout its migration the interior of the vortex is either laminar or inviscid and therefore encounters turbulent interference at its outer surface. In both cases this is because the temperature in the core is below the MRI activation threshold TMRI, and so turbulence (or viscosity) switches off. For this reason vortices survive for relatively long times; in a sense, they are cool bubbles of the dead zone moving inside the hot active region.

As is clear, however, from Figs. 2 and 8, vortices shrink in size as they migrate. Quantitative measurements using the 2D simulations showed that the vortex size changed from about 0.6 H to about 0.25 H, at which point they dissipate. This can be understood by realizing that vortices are overwhelmed by MRI turbulence or viscous diffusion once their radius falls below a certain size. An estimate on the critical size can be obtained by equating the turbulent (or viscous) speed external to the vortex (~ν/s) with the typical vortex circulation speed, estimated from simulations (s Ω). We find the critical vortex size scrit is given by (19)This rough estimate is broadly in agreement with the results of the simulation, and so we conclude that once a vortex shrinks by about three times its original size it diffuses away in the active zone.

An estimate on the vortex lifetime is tied to the speed at which the vortex evolves, in particular to the speed at which it loses mass and shrinks. The shrinking of the vortex is possibly due to the gradual breakdown of the vortical flow via turbulent (or viscous) diffusion at the vortex surface. Here there is a strong shear layer (see the first panel in Fig. 13). Thus the vortex is destroyed gradually from the outside in, and as the outer layers disintegrate they release their mass into the surrounding active medium. A lower limit for the vortex lifetime is thus provided by the diffusion timescale of vorticity over the vortex size. Using the α prescription, it can be written as (20)For the parameters of the vortex we obtained in both the 2D and 3D simulations (s ~ 0.5H, Rf ~ 4), this gives an estimated destruction timescale of about 30 orbits. This is much shorter than the typical lifetime of about 300 orbits we obtained in the 3D simulations (see Fig. 3) or 150 orbits found in the 2D simulation (see Fig. 11), because vorticity diffusion only occurs over a thin layer of size Δs at the surface of the vortex. As a result, the viscosity over the vortex section is reduced by a factor of f ~ 2Δs/s. The thickness of the layer Δs is very difficult to estimate. As illustrated by the snapshots in Fig. 9, it is probably different in the 2D and 3D simulations. It might also be affected by numerical diffusion in our simulations. Nevertheless, its smallness significantly increases the vortex lifetime compared to the estimate given above, and might also account for, at least partly, the different vortex lifetimes we find in the 2D and 3D simulations.

thumbnail Fig. 13

Top panels, from left to right: snapshots of the relative vorticity perturbation, the temperature and the resistivity in the vicinity of the vortex at the 1208th orbit, in the 3D simulation. The fields mapped on top panels have been vertically averaged. Bottom panels, from left to right: snapshots of the relative vorticity perturbation, the temperature and the viscosity in the vicinity of the vortex at the 2597th orbit, in the 2D simulation. We draw with a black line the iso-contour of temperature T = TMRI.

Open with DEXTER

4.2. Vortex migration

We now address the question of the vortex migration process. That we find any migration at all in our simulations is surprising in itself. Indeed, all vortices start their journey at a surface density maximum, which according to the isothermal simulations in Paardekooper et al. (2010) should fix the vortex in place. Generally, isothermal vortices migrate toward high pressure regions, because of asymmetric density wave launching. This result appears to be confirmed by Lyra & Mac Low (2012), who report no migration of vortices in their locally isothermal MRI-turbulent simulations of the inner dead-zone interface, and also by Meheut et al. (2012a), who considered the long-term evolution of a RWI vortex in barotropic isentropic disks. However, recently Richard et al. (2013) reported inward migration in their adiabatic runs of RWI vortex formation. Here, however, the vortex destroys the adverse pressure gradient that would prevent its migration by absorbing all the bump material. In our runs, the vortical perturbation is similar, but the bump size at the beginning of a cycle is significantly bigger and the vortices start their migration before they completely absorb the bump mass. As sanity checks, we have successfully reproduced the results of Richard et al. (2013) albeit in a 2D simulation (see model RBD in Table 1). Finally, the migration rate we measure is about 10 times higher than in Richard et al. (2013). Questions therefore remain: why are vortices migrating in our simulations, and why do they migrate so fast?

We focus on the 2D simulations, as the vertical dimension and MHD turbulence are likely to complicate the picture. Several additional numerical experiments were performed, detailed in Table 1. First, at t = 1000 in our standard 2D run (see Sect. 3.1), we have frozen the temperature and calculated the subsequent evolution of the vortex using a locally isothermal equation of state (model STD2D). We found the vortex remained at its formation location. Immediately, it is clear that the gas thermodynamics is crucial to the migration process.

To further test this idea, we performed a set of simulations that reproduce the set-up described by Paardekooper et al. (2010). Within an isothermal and inviscid 2D disk, we initialized a strong vortex by introducing a vortical velocity perturbation. For three different exponents of the background density radial profile (models PLP1, PLP2, and PLP3), we measured three vortex migration velocities that are in agreement with the results reported by Paardekooper et al. (2010). In particular, the zero pressure-gradient case (model PLP3) is almost neutral in the sense that the migration rate of the vortex is vanishingly small.

Next, in order to interrogate the role of the thermodynamics, we relaxed the assumption of isothermality and modified the cooling function in Eq. (13) so that it takes the form (21)The parameter Tb is the gas background temperature, averaged azimuthally3. This is model PLPCOOL, and its piece-wise cooling law, by depending on the strength of vorticity, forces the vortex to possess a different temperature from its surrounding. If we take τv = 10 orbits and τd = 1 orbit, the vortex is cooler than its surroundings which should relax to the initial temperature profile (see top panel in Fig. 14). In this case, the vortex stays at its initial position, showing that a change in the bulk temperature of the vortex is not enough by itself to affect its migration. We then introduce an azimuthal asymmetry in the vortex cooling by further modifying the cooling function: (22)The resulting relative temperature perturbations differ widely as illustrated in the middle and bottom panels of Fig. 14. With this new cooling prescription, we found significant inward or outward migration of the vortex depending on the sign in Eq. (22), despite the vanishing radial pressure gradient. If the vortex is hotter in the region downstream of its core (in the sense of the background rotation) than the surrounding gas while being cooler in the upstream region (case PLPCOOL +), the vortex migrates inward. In the opposite situation (case PLPCOOL -), the vortex migrates outward with the same velocity. The migration speed is significant and is comparable in magnitude (to less than a factor of two) to the isothermal case when p = −1.5 (see Table 1).

Table 1

Numerical experiments reproducing the set-up of Richard et al. (2013) and Paardekooper et al. (2010) with various surface density profiles and thermodynamics.

thumbnail Fig. 14

Snapshots of the temperature (normalized by the mean temperature) in the vortex in model PLPCOOL (top panel), PLPCOOL- (middle panel), and PLPCOOL+ (bottom panel). The black line shows an iso-vorticity contour where the vorticy is equal to the opposite of the vorticity of the background Keplerian shear profile.

Open with DEXTER

We conclude from this series of experiments that diabaticity can play an important role in vortex migration. In particular, an azimuthal temperature gradient may help the vortex to overcome pressure gradients that would otherwise hold it in place. The detailed and quantitative understanding of these effects is largely beyond the scope of this paper and requires additional more specialized simulations. However, it seems likely that a migration mechanism is at work in our simulations that differs from that discussed in earlier works. It is likely that the new effect is associated with the baroclinic term in the vorticity equation, which vanishes in the barotropic flow of model PLP3, but is non-zero in model PLPCOOL ±. This constitutes the only difference between the two experiments. We speculate that the cooler temperatures in our vortex solutions and the sharp temperature gradient at its surface would mean that the vortex circulation is modified, with fluid parcels having to turn abruptly. This could modify the position of the sonic lines and, along with the modified vortex shape (see the black lines in Fig. 14), might strengthen their asymmetry and trigger radial migration.

In addition, the importance of thermodynamics effects is consistent with the differences we find between the 2D and 3D simulations, particularly when it comes to comparing the migration velocities of the vortex. Obviously, the thermal history of the vortices in these simulations differ on account of the two distinct heating processes (turbulent and ohmic dissipation in the 3D case, viscous dissipation in the 2D simulation). The fact that the different baroclinic terms might result in different vortensity evolution is partly supported by the left panels in Fig. 13, which shows the different vorticity field distribution within a single vortex in the two cases.

4.3. The case of the static dead-zone phase

We end the discussion by noting that our findings naturally explain the difference between the static dead-zone step and the self-consistent step. Indeed, a vortex is observed to grow in the former. The question is, why does this vortex not migrate into the dead zone? We have checked this issue and we have found that, in fact, vortices do migrate inward in this case as well. This is illustrated in Fig. 15. The difference in this case is that the vortex migration range is much smaller than in the self-consistent case. The vortex disappears upon reaching R ~ 3.1 (as opposed to ~2.7 in the self-consistent case) because the resistivity is not a function of temperature in this case, but a function of position only. As the vortex penetrates into the active zone, its resistivity drops to zero which renders the flow vulnerable to MRI-induced velocity fluctuations. As a result, vortices formed during the static dead-zone step have a smaller lifetime (of the order of 100 orbits, as opposed to the 300 orbits measured during the self-consistent step) when they enter the active zone.

5. Conclusion

In this paper we have performed a non-ideal MHD simulation of the region centred around the dead zone inner edge of a protoplanetary disk. In accordance with previously published work (Lyra & Mac Low 2012), we find that a vortex forms at the inner edge of the dead zone. Our simulation reveals that the vortex is not fixed to the pressure maximum as one would expect from the results of Paardekooper et al. (2010); instead, it migrates inward and penetrates into the active zone where it is gradually destroyed by turbulent motions. A few orbits later a new vortex forms at the interface and follows the same evolution, thereby creating what we have called a vortex cycle: formation, migration, and disruption follow each other and comprise a quasi-periodic disk evolution.

thumbnail Fig. 15

Radial position of the density maximum in the 3D run during the static dead-zone step.

Open with DEXTER

We find that the vortex cycle does not exist in isothermal simulations of the dead zone inner edge because the vortex stays at its formation location without migrating inward. We have been unable to design any test in which a vortex moves in a uniform pressure background when its temperature is equal to that of its surroundings. Vortex migration seems to occur only if the vortex temperature is different from (and maybe even free to evolve independantly of) the temperature of its environment. A systematic study and a detailed understanding of vortex migration rates in such specific environments is needed. We caution the reader that simulations using a different cooling function may exhibit cycles with very different properties. Moreover, we completely neglected heat diffusion in the present paper, particularly radiative diffusion. Even if vortices forming in an optically thick region such as the dead zone inner edge are likely to be cooler than the turbulent region in which they will penetrate, radiative diffusion will act as a heating source that will help decrease the difference between the vortex temperature and the temperature of its surroundings. If the heat diffusion timescale is shorter than the vortex cycle, the ionization threshold will be crossed at some point. This will quickly activate the magneto-elliptic instability (Lebovitz & Zweibel 2004; Mizerski & Lyra 2012), the growth rate of which (about one local orbit) is much faster than the vortex cycle timescale, and will disrupt the vortex. An accurate measurement of the cycle period and the vortex migration rate should be done using simulations including radiative transfer that properly account for the vortex thermodynamics.

An additional limitation of our work comes from geometry itself. The 3D simulation uses the cylindrical approximation. Taking into account the disk vertical structure (Richard et al. 2013; Meheut et al. 2012a; Lesur & Papaloizou 2009) will affect both the vortex properties and the shape of the dead-active interface. Along the same lines, the magnetic configuration is also known to influence the development of vortices (Lyra & Klahr 2011; Yu & Li 2009; Yu & Lai 2013). For example, the presence of a net vertical flux in the inner parts might change the results presented here in surprising ways.

Before closing this paper, we speculate on the possible consequences of the vortex cycle for the dynamics of dust particles in the disk. As we have already discussed, vortices such as those we see in our simulations are known to concentrate dust grains and help planet formation processes. However, we have shown here that vortices do not stay at their formation location but migrate inward, most likely carrying the dust particles they captured at the density bump. If the particles are still small once they are released by the disrupted vortex, they might continue to migrate inward as a result of the gas friction. The vortex cycle would then help the dust to pass across the pressure bump barrier. From this discussion, a number of questions arise:

  • 1.

    What dust concentration can the vortex achieve?

The 3D bi-fluid simulations of Meheut et al. (2012b) have shown that the dust-to-gas ratio may change from 0.01 to 1 inside a vortex in three local orbits. This is much shorter than the cycle period, suggesting that vortices formed at the dead zone inner edge would trap the entire content of dust initially present in the bump.

  • 2.

    Can dust grains embedded in the vortex become large enough asa result of collisions alone such that friction becomes negligible atthe time of vortex disruption?

According to 1+1D models of dust growth (e.g. Brauer et al. 2008a), the dust growth rate in laminar regions is mainly due to differential settling. Such a growth timescale can be taken as a proxy for the typical time required to grow centimetre-sized particles into metre-sized bodies within the vortex. It amounts to about a thousand years, which is much longer than the cycle period of a few hundred years. This seems to suggest that particles transported by a vortex across the bump barrier would not grow significantly over one cycle. Particles would then quickly drift toward the central star and take no part in the process of planetesimal formation. A fast mechanism is needed to prevent a loss of the disk’s solid content. Given the large dust-to-gas ratios that are likely to be found inside vortices, such fast processes could be the streaming instability (Youdin & Johansen 2007; Johansen & Youdin 2007), gravitational collapse (Lyra et al. 2008), or a combination of both phenomena (Johansen et al. 2007). However, we caution the reader that gas choatic motions in the vortex resulting from sound waves may disrupt embryos formed by these mechanisms. A detailed study of the outcome of these instabilities is needed.

For such a non-linear and complex problem, it is difficult to go beyond these simple qualitative statements. Clearly, self-consistent simulations including the vortex cycle phenomenon along with dust dynamics (including dust growth) are needed if we want to make any quantitative statement regarding the fate of dust particles at the dead zone inner edge. As we have already argued, such simulations will also have to properly include radiative effects since the vortex migration is sensitive to the gas thermodynamics. Such multifluid radiative MHD simulations are very demanding with present day computing resources.


1

Using the background velocity (i.e. taking into account the modifications to Keplerian rotation induced by pressure) makes very little difference to δωz.

2

In the remainder of the paper, the symbol . denotes an azimuthal, vertical average.

3

The second part of the cooling function in Eq. (21) might seem unnecessary. When this part was omitted we found that the disk cooled down as a whole because the density waves launched by the vortex create large areas in the disk where the vorticity is large. Such areas cool down as a result of the vorticity dependent cooling function despite not being associated with the vortex itself, an artefact that is taken care of by the second part of the cooling function.

Acknowledgments

The authors acknowledge a useful report from Dr Wladimir Lyra that significantly strengthened the results presented in this paper. J.F., S.F., and H.M. acknowledge funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement No. 258729. HL acknowledges support via STFC grant ST/G002584/1. The simulations presented in this paper were granted access to the HPC resources of Cines under the allocation x2013042231 and x2014042231 made by GENCI (Grand Equipement National de Calcul Intensif).

References

All Tables

Table 1

Numerical experiments reproducing the set-up of Richard et al. (2013) and Paardekooper et al. (2010) with various surface density profiles and thermodynamics.

All Figures

thumbnail Fig. 1

Temperature profile at the end of the ideal phase averaged over 200 inner orbits. The black dashed lines show where the mean temperature is equal to (5 / 4)TMRI. The hashed regions on both sides of the figure correspond to the buffer zones.

Open with DEXTER
In the text
thumbnail Fig. 2

From left to right and top to bottom: snapshot of the relative vorticity perturbation at the 1087th, 1124th, 1173th, 1208th, 1216th, 1225th, 1228th, and 1310th orbit. The relative vorticity perturbation has been vertically averaged. We draw with a black line the iso-contour corresponding to 80% of the density maximum ρm. The location of the density maximum is shown by the black cross.

Open with DEXTER
In the text
thumbnail Fig. 3

Radial position of the density maximum in the 3D run during the static dead-zone phase (in black) and the self-consistent phase (red curve).

Open with DEXTER
In the text
thumbnail Fig. 4

Temperature averaged over the vortex area in the 3D run during the static dead-zone phase (in black) and the self-consistent phase (red curve).

Open with DEXTER
In the text
thumbnail Fig. 5

Snapshot of the relative density perturbation at the 1310th orbit. The relative density perturbation has been vertically averaged. The iso-contour of 80% of density maximum ρm is drawn with a black line. The location of the density maximum is shown by the black cross.

Open with DEXTER
In the text
thumbnail Fig. 6

Radial profiles of at the 1000th orbit (red curve) and after the simulation has reached a quasi-steady state (black curve). is averaged azimuthally for both cases.

Open with DEXTER
In the text
thumbnail Fig. 7

Evolution of the logarithm of the vorticity perturbation δω = ωmax(Rvort) − ωmin(Rvort) at the bump location in the 2D viscous simulation. The best fit of the linear part of the amplitude growth is plotted in black and has a slope equal to γ/ Ω(Rf) = 0.2.

Open with DEXTER
In the text
thumbnail Fig. 8

Radial position of the density maximum in the 2D viscous simulation using the RAMSES code.

Open with DEXTER
In the text
thumbnail Fig. 9

From left to right and top to bottom: snapshot of the relative vorticity perturbation at the t = 2502, 2568, 2626, and 2668. The relative vorticity perturbation has been vertically averaged. The iso-contour of 80% of density maximum ρm is drawn with a black line. The location of the density maximum is shown by the black cross.

Open with DEXTER
In the text
thumbnail Fig. 10

Evolution of the mean temperature inside the vortex during the cycle starting at the t = 2500th inner orbit.

Open with DEXTER
In the text
thumbnail Fig. 11

Evolution of the vortex total mass over the cycle starting at t = 2500, given in units of the initial disk mass. It is found to decrease steadily during the vortex migration.

Open with DEXTER
In the text
thumbnail Fig. 12

Radial position of the density maximum in the 2D viscous simulation using the PLUTO code.

Open with DEXTER
In the text
thumbnail Fig. 13

Top panels, from left to right: snapshots of the relative vorticity perturbation, the temperature and the resistivity in the vicinity of the vortex at the 1208th orbit, in the 3D simulation. The fields mapped on top panels have been vertically averaged. Bottom panels, from left to right: snapshots of the relative vorticity perturbation, the temperature and the viscosity in the vicinity of the vortex at the 2597th orbit, in the 2D simulation. We draw with a black line the iso-contour of temperature T = TMRI.

Open with DEXTER
In the text
thumbnail Fig. 14

Snapshots of the temperature (normalized by the mean temperature) in the vortex in model PLPCOOL (top panel), PLPCOOL- (middle panel), and PLPCOOL+ (bottom panel). The black line shows an iso-vorticity contour where the vorticy is equal to the opposite of the vorticity of the background Keplerian shear profile.

Open with DEXTER
In the text
thumbnail Fig. 15

Radial position of the density maximum in the 3D run during the static dead-zone step.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.