EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 644 (December 2020) A&A, 644 (2020) A50 Full HTML
Open Access
Issue
A&A
Volume 644, December 2020
Article Number A50
Number of page(s) 8
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202039294
Published online 30 November 2020

© L. Flores-Rivera et al. 2020

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

1 Introduction

Hydrodynamical (HD) simulations in protoplanetary disks have gained substantial attention during the last decades. Due to the low coupling between the magnetic field and the gas in a large fraction of the disk (Turner et al. 2014; Dzyurkevich et al. 2013) we expect that HD instabilities play a more important role in the gas kinematics and the overall gas evolution in protoplanetary disks. In theory, the presence of vertical shear instability (VSI; Goldreich & Schubert 1967; Fricke 1968) in purely HD accretion disks has been well studied by careful perturbational analysis and numerical simulations in several works (e.g., Nelson et al. 2013; Stoll & Kley 2014; Barker & Latter 2015; Lin & Youdin 2015). Several works have shown the importance of this instability to explain the turbulence and angular momentum transport in the disk (Arlt & Urpin 2004; Nelson et al. 2013; Stoll & Kley 2014). The main ingredients of the VSI are a short thermal relaxation timescale and a vertical gradient of the rotational velocity in the disk. However, the effect of the VSI together with realistic thermal profiles in the disk are not well known. The evolution of the VSI has been tested in the context of planet formation and disk evolution. Recent simulations conducted by Lin (2019) have found that the VSI can be suppressed by dust settling and grain growth; however, the streaming instability (SI) may be present when VSI is active in the turbulent disk (Schäfer et al. 2020). In the context of disk evolution its been proposed that photoevaporation by energetic photons is a potential agent of disk dispersal (e.g., Shu et al. 1994; Hollenbach et al. 1994; Clarke et al. 2001; Owen et al. 2012; Alexander et al. 2006a,b, 2014; Ercolano & Pascucci 2017). Investigating the effect of high-energy radiation on protoplanetary disks can be quite complex, but recent efforts have attempted to explain the effect of high-energy radiation onto the disk through observations of winds (e.g., Ercolano et al. 2008; Gorti et al. 2009; Owen et al. 2010, 2011). A prominent signature of photoevaporation processes is the detection of forbidden lines in the protoplanetary disk spectrum. The detection of the NeII line from the surface of TW Hya at around 10 AU (Pascucci et al. 2011), shows the important influence of high-energy photons on the surface layers of the disk; however, the vertical extension of the wind is still quite inconclusive. Moreover, Ballabio et al. (2020) found the synthetic line profile of NeII to be consistent with thermally driven winds, while other forbidden molecules such as OI, OII, and SII require a different scenario, possibly magneto driven winds (Fang et al. 2018; Banzatti et al. 2019). A photoevaporative process in the inner disk (<5 AU) has been tested by Wang et al. (2019) and Gressel et al. (2020) that considered hydromagnetic diffusive effects, thermochemistry, and ray-tracing radiative transfer, but the wind description in terms of the Keplerian rotation and mass loss rate is different. Gressel et al. (2020) determine that the deviation of the winds caused by the magneto-centrifugal mechanism are super-Keplerian, while Wang et al. (2019) found it to be sub-Keplerian. The mass loss in Gressel et al. (2020) is ~ 10−7 M⊙ yr−1, whereas Wang et al. (2019) is ~10−8 M⊙ yr−1.

The coexistence between the VSI and photoevaporation processes by stellar photons is of particular interest, especially in regions where we expect a steep temperature profile and fast thermal relaxation timescales, which few have found that the VSI might operate down to ~1 AU (e.g., Klahr et al. 2018; Lyra & Umurhan 2019; Pfeil & Klahr 2019, 2020). The effect of viscous heating in the disk remains relatively small, making it plausible to use passive irradiated disk models to accurately represent the thermal structure (Flock et al. 2019). As a first step towards including photoevaporation processes in our HD simulations, our aim in this first paper is to study the convergence, the strength, and the region where the VSI operate in the disk, using the full meridional domain.

The structure of this paper is as follows. In Sect. 2 we present the methodology, and we describe the disk model and boundary conditions to resolve the VSI using a global 2D isothermal accretion disk configuration. In Sect. 3 we present the results with a careful analysis of the kinetic energy and vertical velocity from a global and local perspective. We also provide predictions of the influence that extreme ultraviolet (EUV), far-ultraviolet (FUV), and X-ray photons canhave in the inner parts of the disk. In Sect. 4 we discuss and present a perspective of future work regarding the implications for planet formation and winds in the disk. Finally, in Sect. 5 we present our conclusions.

2 Methods

We performed our HD simulations using PLUTO 4.31 by Mignone (2007). The dynamics of the fluids is described by the conservation laws accounting for the divergence of mass density and momentum density. The HD fundamental equations in our setup are (1) (2)

where ρ is the mass density, v is the velocity vector, ρv is the momentum density vector, and P is the gas pressure. We select the isothermal equation of state with , where cs is the isothermal sound speed. The code solves consistently the HD equations in a 2D geometry and considers the system in spherical coordinates (r, θ, ϕ) with axisymmetry in the azimuth. The grid cells are set up with a logarithmic increase in the radial domain and a uniform spacing in cos(θ) for the meridional domain (e.g., Ormel et al. 2015) in order to have a better resolution around the midplane layers. The simulations are scale free. We set the radial domain extending from 0.5 to 5.0 and the vertical domain covering approximately 180°.

Table 1

Grid setup and initial parameters for the different runs.

2.1 Disk model and boundary conditions

Protoplanetary disks, rotating with Keplerian angular velocity and vertical stratification are dynamically stable according to the Solberg-Høiland criteria (e.g., Rüdiger et al. 2002). The accretion disk setup in equilibrium in cylindrical coordinates (R, Z) (Nelson et al. 2013) are defined as (3) (4)

where = 4.5 × 10−10g cm−3 = 2.7 × 1014 cm−3 is the initial density at the midplane, and Σ0 is the initial surface density (see Table 1), R0 is the reference radius, H0 is the reference scale height, G is the gravitational constant, M is the mass of the star, p is the power-law fitting exponent of the density profile, and is the Keplerian frequency. The disk scale height in terms of the radius is (5)

where the reference scale height ratio is H0 /R0 = 0.1 (see Table 1) and cs is related to the temperature through the power-law fitting exponent q. The radial profile of cs is given based on Eq. (5) so that cs = HR/Ωk. The scale length of the vertical mode is proportional to H0.

In order to compute the fluxes in each cell, we employ second-order piece-wise linear spatial reconstruction. This reconstruction method is appropriate for both uniform and non-uniform grid spacing based on the approach presented inMignone (2014). We use the Harten, Lax, Van Leer (HLL) Riemann solver. For the time integration we adopt a second-order Runge-Kutta algorithm and we set the Courant-Friedrichs-Lewy (CFL) number to 0.25. Table 1 summarizes the model parameter and the different runs for each resolution. For 101 and 203 cells per HR, the θ inner limit is shifted slightly inwards due to the high computational precision issue when increasing the resolution at the θ boundary. We investigated the third-order piecewise parabolic method (PPM3) spatial reconstruction using third-order Runge-Kutta (RK3) for the time integration for all resolution cases. We found that the PPM3 could resolve the VSI for even lower resolutions; however, our PPM3 configuration is currently is not yet stable. Therefore, we only show our analysis and discussion for the VSI using the linear reconstruction method for this paper.

2.2 EUV, FUV, and X-ray column density

The inner parts of the disk are highly influenced by the stellar radiation potentially ionizing the upper layers of the disk. As a first step to implement photoheating processes caused by stellar EUV, FUV, and X-ray photons, we calculate the radial column density locally to predict at what scale height the EUV and X-ray–FUV heated layers are most dominant. The radial column density is calculated by integrating the hydrogen nuclei column density along the radial line of sight (6)

where R* is the star radius and nρ is the number density of hydrogen nuclei. For this first paper in the series, since we do not consider detailed chemistry or the frequency of UV photons, we determine the boundaries in terms of Σr by post-processing. We focus on the definition of the radial column density to locate the boundary of the ionized region (EUV heated region) and the neutral layer (X-ray–FUV heated layer) which is empirically determined from the simulations of Nakatani et al. (2018a,b). We select the boundary between the ionized or EUV heated region, namely the H–H+ boundary, to be located at Σr =1019 cm−2, assuming the EUV stellar luminosity is ~ 1030 erg s−1. The EUV heated region typically has temperatures of orders of magnitude between 103 and 104 K and the hydrogen is fully ionized here. The boundary of the X-ray–FUV heated layer is defined at Σr =1022 cm−2, where the gas is heated to 103 K and it dominates in the neutral layers. While soft X-rays are considered to be an important component for photoevaporation, the locations of FUV and X-ray heated layers are nearly identical Σr ≤1021−22 cm−2 (Gorti & Hollenbach 2008; Nakatani et al. 2018b). Therefore, setting the boundary of the neutral photo-heated layer at Σr ≤1021−22 cm−2 is valid even in cases where X-rays effects are included. It should be noted that we neglect scattered light in our ray-tracing. For more information about the effects of stellar photons and the launch of winds, see Sects 4.2 and 4.3.

3 Results

3.1 Dynamical evolution

In order to investigate the kinematics we start by looking at the time evolution of the kinetic energy. Figure 1 shows the average kinetic energy at different grid cell resolutions. Our prescription of the global kinetic energy follows (7)

where vZ = vr cos(θ) – vθ sin(θ) is the vertical component of the velocity and Ekin is normalized with respect to (8)

where vk is the Keplerian velocity. Due to our large global domain, it becomes necessary to transform our spherical velocity components to cylindrical when describing the vertical kinetic energy. For the resolution analysis we simulate a large parameter space, starting from 12 cells per scale height until 203 cells per scale height. At the lowest resolution we observed that the VSI is not resolved as the kinetic energy is not growing (see Fig. 1). We observe that the initial perturbations are damped and the kinetic energy is decreasing leading to a fully laminar disk. For resolutions higher than 12 cells per scale height, the growth seen at ≥25 orbits is characteristic of the initial increase phase of the kinetic energy as the disk become unstable to the VSI, as found by Nelson et al. (2013). The 25 cells per scale height case is often used in the literature as a standard resolution to perform global simulations. However, we note that this resolution does not fully solve the VSI in accretion disks. It reaches an energy saturation slightly lower than 10−5, and after that it starts to increase monotonically without reaching a steady state, at least within the first 200 orbits. Because it is characteristic for all resolutions, the first energy saturation ~10−5 is a consequence of the linear phase of the instability spreading throughout the system. The second phase of increasing kinetic energy might be related to the slower growth of the body modes that start to become visible after 25 orbits, when the quickly growing finger modes have already saturated. During this second phase, the body modes grow in the disk until they reach saturation after around 90 orbits. Our results show that at a resolution of about 100 cells per scale height the kinetic energy levels converge. For even higher resolutions, which are computationally very expensive, we predict the same convergence trend.

Our total kinetic energy follows the same trend as the perturbed kinetic energy from Nelson et al. (2013). The resolution for the VSI used by Cui & Bai (2020) lies between 96 and 108 cells per scale height, which is in good agreement with our resolution study as well. Figure 2 shows a snapshot of the three main physical quantities in our analysis: the number density, the vertical velocity, and the kinetic energy in physical units. Overplotted are the contour lines in which the EUV heated region and the X-ray–FUV heated boundary layer are defined, based on our definition of Σr (Sect. 2.2).

The density (Fig. 2a) has a smooth profile until it reaches the upper layers. Our density floor starts to become important at HR = 13, which lies in a region where we do not expect VSI activity. As expected, the FUV photons penetrate deeper than EUV photons within the first 5 AU (Fig. 2). The X-ray–FUV heated layer and EUV heated region reach HR = 6.2 and HR = 9.7, respectively.At this photoevaporation region we predict that HI regions become optically thick against EUV with a Σr = 1019 cm−2 to Σr = 1022 cm−2 (Fig. A.1).

The profile of the vertical velocity (Fig. 2b) shows clearly the upward motion shaded in blue and the downward motion shaded in red. In this region, the vertical velocity reaches a fraction of the sound speed.

The kinetic energy profile (Fig. 2c) shows the highest concentration present as corrugated wavy structures resembling the body modes in the disk. In the next section we investigate this in detail from a more local perspective.

thumbnail Fig. 1

Time evolution of the kinetic energy for different resolutions and space averaged at R0 = 1 AU. The kinetic energy is normalized with the kinetic energy from pure Keplerian velocity. The different resolutions are color-coded (see inset).
thumbnail Fig. 2

Snapshot of the three main global variables in the simulation. Panel a: number density in cm−3. Panel b: vertical velocity normalized with respect to the sound speed. For better visibility of the midplane region, we show only an extent of ± 53° with respect to the midplane. Panel c: vertical kinetic energy in g cm−1 s−2. Overplotted are contour solid lines where the column density of the EUV-heated region is defined by Σr =1019 cm−2 (in purple), and the column density of the X-rays/FUV-heated boundary layer is defined by Σr = 1022 cm−2 (in dark green). The number density of hydrogen nuclei at the H/H+ boundary is ~109 cm−3 and at the X-rays/FUV-heated layer is ~1012 cm−3.

3.2 Local kinetic energy evolution over height and time

We show the time-evolution of the kinetic energy over HR at R0 in Fig. 3. We transform the velocity vectors to cylindrical and, in order to see the local kinetic energy in a more continuous spatial structure, we then use a cubic method to interpolate the data on a new grid space with a Z direction between −1 and 1 AU at R0. Figure 3 shows that the kinetic energy within few first orbits represents the fast growth of the finger modes in the upper layersof the disk, which is expected. At the midplane, there is no VSI activity since the vertical shear vanishes right at the midplane. Once the body modes grow, the large-scale velocity perturbations threaten the upper and lower hemisphere. At this stage the kinetic energy at the midplane also grows.

Two important features are depicted from all resolution cases in Fig. 3. The first is the presence at the beginning of vertical lanes of kinetic energy. We refer to these vertical kinetic energy lanes as finger modes, and they are consistent with the linear energy increase phases from Fig. 1. The second feature corresponds to the appearance of the body modes. As seen in Fig. 3, the growth rate of the body modes increases with resolution. For the 203 cells per scale height, the body modes saturate at around ~70 orbits. At HR = 0 (midplane), the kinetic energy is very low but even more detailed analysis in the midplane is performed in the last paragraph.

Figure 4 shows the time averaged vertical profile of the kinetic energy. The vertical profile follows a trend similar to that of the density profile. Only for the lowest resolution case do we observe a plateau feature, which is due to the insufficient resolution to solve for the VSI. The VSI seems to operate in the full vertical extent. We expect a physical boundary for the VSI to be the high-temperature region caused by the FUV radiation. We also expect this position to be where the X-ray–FUV heated boundary layer is located at HR = 6.2 (see Figs. 4 and 5). Figure 5 shows that the time averaged vertical velocity increases above 60% of the sound speed at the expected wind base at HR around 6.2. In the midplane layers the time averaged vertical motions remain small, at a level of around 1%.

Figure 6 shows the time-evolution of the local kinetic energy at the midplane at R0. The local kinetic energy at the midplane replicates the same two growth regimes as the averaged global kinetic energy seen in Fig. 1. The fast growth from the finger modes can even be seen at the vertical velocity at the midplane, and the slow growth of the body modes. During the growth phase of the generation of the body modes, the local kinetic energy in the midplane in all three resolutions follow a linear growth rate of ~0.13 orbits. This growth phase persists within the first ~70 orbits; it corresponds to the slow growth rate of the body modes until reaching a saturation level of ~10−4 and it is 2 times slower compared to the general growth rate found by Nelson et al. (2013). The 203 and 101 cells per scale height generate the body modes faster than the 50 cells per scale height, which is the same tendency as found in Fig. 3. We conclude that we need more than 50 cells per scale height to solve the VSI in accretion disks.

thumbnail Fig. 3

Time-evolution of the local kinetic energy at R0 = 1 over height for different resolutions (from top to bottom: 25, 50, 101, 203 cells per scale height). The local kinetic energy is normalized to the global kinetic energy (see Eq. (8)).
thumbnail Fig. 4

Averaged profile of the local kinetic energy normalized with respect to the kinetic energy in the midplane for the 25, 50, 101, and 203 cells per HR at R0 = 1 AU. Thevertical colored lines represent HR = 6.2, where the X-ray–FUV heated boundary layer is (in dark green), and the HR = 9.7, where the EUV heated region is located (in purple) (same as in Fig. 2). The models for 101 and 203 cells per HR are in good agreement. The density floor is located HR = 13, and therefore does not influence our result.
thumbnail Fig. 5

Averaged profile of the vertical component of the velocity normalized with respect to the sound speed for the 203 cells per HR. The vertical colored lines are the same as in Fig. 4. At HR = 6.2 the vertical velocity reaches 0.3cs. The hot ionized region is expected to start beyond HR > 6.2.
thumbnail Fig. 6

Time evolution of the local kinetic energy in the midplane for the 50, 101, and 203 cells per scale height resolutions at R0. The black dashed line represent the slope of the linear growth rate phase of the local energy in the midplane. The vertical gray lane represent the approximate time where the 101 and 203 cells per scale height resolutions reach saturation and the vertical dashed line represents the time where the 50 cells per scale height converges to stability.

3.3 Small-scale vortices

In Fig. 7 we show the result of the line integral convolution (LIC) method (e.g., Cabral & Leedom 1993) to visualize thevelocity magnitude vector field associated with velocity perturbations from the VSI, for the lowest and highest resolutions. The velocity magnitude includes the velocity field in every component and follows the ordinary prescription of , where vR = vr sin(θ) + vθ cos(θ). From the LIC method, we find small-scale vortices in both resolution cases; however, these small-scale vortices are not featured in the lowest resolution case.

The possibility of the presence of vortices formed by VSI was first found by Richard et al. (2016) in the r− ϕ plane. These vortices can be long-lived locally around 500 orbits (Manger & Klahr 2018) with an HR = 0.1 and with an azimuthal and radial extension of about 40 and 4 au (Flock et al. 2020), respectively. In our case the vortices appear in the r− Z plane, mostly inside the region between two large-scale upward and downward motions. They might be related to Kelvin-Helmholtz instabilities appearing between these large-scale motions, as was described by Latter & Papaloizou (2018). Moreover, these small-scale vortices happen in an HR < 5, and are therefore more important for dust evolution.

4 Implications for planet formation and winds in disk

The VSI is a large-scale global instability that potentially plays an important role in the evolution and dynamics of protoplanetary disks. In the midplane the dust settling overcomes the VSI, and as the dust concentrates and grows to millimeter size and greater, these particles settle against VSI and eventually undergo SI to favor planet formation (Lin 2019). However, the VSI in combination with the SI leads to radial dust concentration in long-lived accumulations which are significantly denser areas than those formed by the SI alone (Schäfer et al. 2020).

We observe symmetric features, which have not been observed before, after 100 orbits in the time evolution of the kinetic energy (Fig. 3). Furthermore, we refer to the LIC plots (Fig. 7) that show the presence of small vortices along the midplane, where grains can be concentrated. We encourage further analysis of this manner and its effect on planet formation in the future.

4.1 Importance of stellar photon layers and the launch of winds

The VSI can maintain small particles (<100 μm) floating (Lin 2019) until they eventually become ionized by FUV photons at the disk surface and contribute to heating in the atmosphere.Both the photoionization effects and the high temperature in this region make the surface of the disk and atmosphere an optimal scenario for photoevaporation processes, resulting in the launch of winds. In the disk atmosphere the MHD winds are likely favorable (Cui & Bai 2020) at less than 5 AU. Cui & Bai (2020) showed that VSI and MHD winds can coexist, while magnetized winds being the dominant process to transport angular momentum. However, this would depend on (1) the so-called wind base, which depends mainly on the position of the EUV and FUV layers, and (2) a realistic thermal profile. By considering the initial flux and the extinction along the radial line of sight we can determine the location and extension of the wind in the future work. The detailed role of the VSI in the mechanism of wind launching has to be investigated, especially because the magnetic coupling to the gas remains low in the FUV shielded region (HR < 6.2). Recent improvements by Wang et al. (2019) and Gressel et al. (2020), aimed at including more realistic scenarios, have included thermochemistry and irradiation in their global magnetohydrodynamical simulations.

4.2 Frequency-integrated cross sections

Perhaps the most robust way to define the influence of EUV and X-ray–FUV photon layers is by taking into consideration the frequency-dependent cross sections for each photon. The radial optical depth of each photon depends on the frequency range that is used in the integration of the ionization rates and the associated heating rates. This is very important to consider when characterizing the ionized region in the disk atmosphere, which is influenced by EUV photons. The atomic hydrogen abundance in this ionized region is negligible in the ionization front where we expect a radial optical depth of the EUV equal to 1 (Hollenbach et al. 1994; Tanaka et al. 2013; Nakatani et al. 2018a,b). This boundary between the ionized region and the neutral layer is instead located at an atomic hydrogen column density of 1018 cm−2. This atomic hydrogen column density can be translated into a hydrogen column density of 1019−20 cm−2 depending on the stellar EUV luminosity, and therefore is used here to define the EUV heated region.

Given that FUV photons can promote heating via photoelectric heating (Bakes & Tielens 1994) in the disk atmosphere, these photons are attenuated by dust once the hydrogen column density is between 1021 and 1022 cm−2 (Gorti et al. 2009; Nakatani et al. 2018a). When considering the radial optical depth of the FUV, the FUV opacity evaluated for graphite is about ~104 cm2 g−1, which is in the wavelength regime between 0.2 and 0.09 μm (6 eV < hν < 13.6 eV). Assuming a dust-to-gas mass ratio of 0.01, this gives a total FUV opacity of 100 cm2 g−1. While it is more likely that the dust-to-gas mass ratio is smaller than 0.01 in the upper layers, perhaps a more nominal range could be somewhat from 0.001 to 0.0001, shifting the radial optical depth of the FUV equal to 1 closer to the midplane. This FUV opacity prescription is consistent with the dust population used by Flock et al. (2016) to describe the physical properties of the inner rim. In realistic disk scenarios, which we consider for the second part of this project, the EUV and FUV flux relies on the wavelength-dependent initial flux we assume and the extinction along the radial line of sight at each location in the disk. Meanwhile, we expect the EUV to provide an important effect on the extension of the wind in the disk atmosphere.

thumbnail Fig. 7

Snapshot of vmag for the 203 cells per scale height (panel a) compared to the lowest resolution of 25 cells per scale height (panel b) at 200 orbits. The colored contour plot shows vmag and the streamlines show the velocity flow pattern when applying the LIC method. The vmag range (from 0 to 0.2 in code units) is narrowed down to highlight the low-velocity perturbations close to the midplane associated with the VSI.

4.3 Electron density

The temperature and electron number density can be inferred from thermally excited forbidden lines (e.g., [OI] 5577/6300 line ratio) in regions where high temperatures (>5000 K) lead to collisional excitations of electrons (Simon et al. 2016). Including or not the effects of X-rays, the electron relative abundance between the EUV heated layer and the X-ray–FUV heated layer is still ~ 10−4, which comes from atomic carbon ionization by FUV. The electron density in this region can be as high as ~ 104–105 cm−3 at ~1 AU (Nakatani et al. 2018a,b). We expect that at <1 AU, the X-ray, EUV, and FUV fluxes are higher; therefore, photons penetrate deeper in the disk and the electron density can be higher. On the other hand, farther from the star the electron abundance decreases, and as does the electron density.

5 Conclusions

We performed 2.5D (axisymmetric) hydrodynamical simulations for an isothermal accretion disk with different resolutions, up to 203 cells per scale height to fully resolve for the VSI. Resolution study for the VSI is important to assess its effect in different regions of the disk. We determine that the standard resolution of 50 cells per scale height is the lower limit to resolve the VSI. This determination is based on the total kinetic energy saturation level. Our key points can be summarized as follows:

  • 1.

    We found the VSI to operate throughout the entire disk until reaching our density floor 0.003 cm−3 at around HR = 13 at R0 = 1 for resolutions ≥50 cells per scale height. We expect turbulent scenarios where the X-ray–FUV heated boundary layer at HR = 6.2 is located, whereas the EUV heated region reaches HR = 9.7. Nevertheless, we encourage the use of high resolution (≥50 cells per scale height) to properly conduct the analysis of VSI.

  • 2.

    We found at a resolution of 25 cells per HR the appearance of a plateau around the midplane when looking at the local kinetic energy over height. This could serve as a diagnostic to check whether or not the VSI is fully resolved.

  • 3.

    We observe a clear dependence on the growth rate of the body modes with resolution. For the highest resolution the growth of the body modes saturates after 70 orbits.

  • 4.

    We diagnose the presence of small-scale vortices between the large-scale motion of the VSI when performing the LIC method. These small-scale vortices reach HR < 5.0; therefore they are more important for dust evolution. The vertical kinetic energy shows a more symmetric oscillation for the upper and lower hemisphere for the highest resolution. In the future we will fully investigate the dynamical influence of the heated layers on the VSI and the disk evolution. A resolution of 50 cells per scale HR seems to capture the main characteristics of the VSI, and so is suitable for future studies.

Acknowledgements

This work is supported by the European Research Council (ERC) project under the European Union’s Horizon 2020 research and innovation programme number 757957. The LIC images were generated by Thomas Müller from the Haus dier Astronmomie. We would like to thank the matplotlib team (Hunter 2007) for the good quality tool in order to better visualize thedata.

Appendix A Radial column density

We conduct a surface density analysis for the determination of the EUV heated region and X-ray–FUV heated boundary layer (see Sect. 2.2). Figure A.1 shows the contour plot of the radial column density; overplotted are the two regions of interest.

thumbnail Fig. A.1

Radial column density in cm−2. Overplotted are solid lines where Σr = 1019 cm−2 represents the EUV heated region (in purple), and Σr = 1022 cm−2 represents the X-ray–FUV heated boundary layer (in dark green).

References

  1. Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006a, MNRAS, 369, 216 [NASA ADS] [CrossRef] [Google Scholar]
  2. Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006b, MNRAS, 369, 229 [NASA ADS] [CrossRef] [Google Scholar]
  3. Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson, AZ: University of Arizona Press), 475 [Google Scholar]
  4. Arlt, R., & Urpin, V. 2004, A&A, 426, 755 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822 [NASA ADS] [CrossRef] [Google Scholar]
  6. Ballabio, G., Alexander, R. D., & Clarke, C. J. 2020, MNRAS, 496, 2932 [CrossRef] [Google Scholar]
  7. Banzatti, A., Pascucci, I., Edwards, S., et al. 2019, ApJ, 870, 76 [NASA ADS] [CrossRef] [Google Scholar]
  8. Barker, A. J., & Latter, H. N. 2015, MNRAS, 450, 21 [NASA ADS] [CrossRef] [Google Scholar]
  9. Cabral, B., & Leedom, L. C. 1993, in Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’93 (New York, NY, USA: Association for Computing Machinery), 263 [Google Scholar]
  10. Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485 [NASA ADS] [CrossRef] [Google Scholar]
  11. Cui, C., & Bai, X.-N. 2020, ApJ, 891, 30 [NASA ADS] [CrossRef] [Google Scholar]
  12. Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114 [NASA ADS] [CrossRef] [Google Scholar]
  13. Ercolano, B., & Pascucci, I. 2017, R. Soc. Open Sci., 4, 170114 [Google Scholar]
  14. Ercolano, B., Drake, J. J., Raymond, J. C., & Clarke, C. C. 2008, ApJ, 688, 398 [NASA ADS] [CrossRef] [Google Scholar]
  15. Fang, M., Pascucci, I., Edwards, S., et al. 2018, ApJ, 868, 28 [NASA ADS] [CrossRef] [Google Scholar]
  16. Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2016, ApJ, 827, 144 [NASA ADS] [CrossRef] [Google Scholar]
  17. Flock, M., Turner, N. J., Mulders, G. D., et al. 2019, A&A, 630, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155 [CrossRef] [Google Scholar]
  19. Fricke, K. 1968, Z. Astrophys., 68, 317 [NASA ADS] [Google Scholar]
  20. Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 [Google Scholar]
  21. Gorti, U., & Hollenbach, D. 2008, ApJ, 683, 287 [NASA ADS] [CrossRef] [Google Scholar]
  22. Gorti, U., Dullemond, C. P., & Hollenbach, D. 2009, ApJ, 705, 1237 [NASA ADS] [CrossRef] [Google Scholar]
  23. Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126 [CrossRef] [Google Scholar]
  24. Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654 [NASA ADS] [CrossRef] [Google Scholar]
  25. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  26. Klahr, H., Pfeil, T., & Schreiber, A. 2018, Handbook of Exoplanets (Cham: Springer), 138 [Google Scholar]
  27. Latter, H., & Papaloizou, J. 2018, MNRAS, 474, 3110 [CrossRef] [Google Scholar]
  28. Lin, M.-K. 2019, MNRAS, 485, 5221 [NASA ADS] [CrossRef] [Google Scholar]
  29. Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17 [NASA ADS] [CrossRef] [Google Scholar]
  30. Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001 [NASA ADS] [CrossRef] [Google Scholar]
  31. Manger, N., & Klahr, H. 2018, MNRAS, 480, 2125 [NASA ADS] [CrossRef] [Google Scholar]
  32. Mignone, A. 2007, ApJ, 170, 228 [Google Scholar]
  33. Mignone, A. 2014, Comput. Phys., 270, 784 [Google Scholar]
  34. Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018a, ApJ, 857, 57 [CrossRef] [Google Scholar]
  35. Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018b, ApJ, 865, 75 [NASA ADS] [CrossRef] [Google Scholar]
  36. Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, A&A, 435, 2610 [Google Scholar]
  37. Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512 [NASA ADS] [CrossRef] [Google Scholar]
  38. Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415 [NASA ADS] [CrossRef] [Google Scholar]
  39. Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13 [NASA ADS] [CrossRef] [Google Scholar]
  40. Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880 [NASA ADS] [CrossRef] [Google Scholar]
  41. Pascucci, I., Sterzik, M., Alexander, R. D., et al. 2011, ApJ, 736, 13 [NASA ADS] [CrossRef] [Google Scholar]
  42. Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150 [NASA ADS] [CrossRef] [Google Scholar]
  43. Pfeil, T., & Klahr, H. 2020, ApJ, submitted [arXiv:2008.11195] [Google Scholar]
  44. Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571 [NASA ADS] [CrossRef] [Google Scholar]
  45. Rüdiger, G., Arlt, R., & Shalybkov, D. 2002, A&A, 391, 781 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190 [CrossRef] [EDP Sciences] [Google Scholar]
  47. Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781 [Google Scholar]
  48. Simon, M. N., Pascucci, I., Edwards, S., et al. 2016, ApJ, 831, 169 [NASA ADS] [CrossRef] [Google Scholar]
  49. Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Tanaka, K. E. I., Nakamoto, T., & Omukai, K. 2013, ApJ, 773, 155 [NASA ADS] [CrossRef] [Google Scholar]
  51. Turner, N., Fromang, S., Gammie, C., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. Klessen, C. Dullemond, & T. Henning, (Tucson, AZ: University of Arizona Press), 411 [Google Scholar]
  52. Wang, L., Bai, X.-N., & Goodman, J. 2019, ApJ, 874, 90 [NASA ADS] [CrossRef] [Google Scholar]
1

PLUTO 4.3 is an open source code available for download: http://plutocode.ph.unito.it/download.html

All Tables

Table 1

Grid setup and initial parameters for the different runs.

In the text

All Figures

thumbnail Fig. 1

Time evolution of the kinetic energy for different resolutions and space averaged at R0 = 1 AU. The kinetic energy is normalized with the kinetic energy from pure Keplerian velocity. The different resolutions are color-coded (see inset).
In the text
thumbnail Fig. 2

Snapshot of the three main global variables in the simulation. Panel a: number density in cm−3. Panel b: vertical velocity normalized with respect to the sound speed. For better visibility of the midplane region, we show only an extent of ± 53° with respect to the midplane. Panel c: vertical kinetic energy in g cm−1 s−2. Overplotted are contour solid lines where the column density of the EUV-heated region is defined by Σr =1019 cm−2 (in purple), and the column density of the X-rays/FUV-heated boundary layer is defined by Σr = 1022 cm−2 (in dark green). The number density of hydrogen nuclei at the H/H+ boundary is ~109 cm−3 and at the X-rays/FUV-heated layer is ~1012 cm−3.
In the text
thumbnail Fig. 3

Time-evolution of the local kinetic energy at R0 = 1 over height for different resolutions (from top to bottom: 25, 50, 101, 203 cells per scale height). The local kinetic energy is normalized to the global kinetic energy (see Eq. (8)).
In the text
thumbnail Fig. 4

Averaged profile of the local kinetic energy normalized with respect to the kinetic energy in the midplane for the 25, 50, 101, and 203 cells per HR at R0 = 1 AU. Thevertical colored lines represent HR = 6.2, where the X-ray–FUV heated boundary layer is (in dark green), and the HR = 9.7, where the EUV heated region is located (in purple) (same as in Fig. 2). The models for 101 and 203 cells per HR are in good agreement. The density floor is located HR = 13, and therefore does not influence our result.
In the text
thumbnail Fig. 5

Averaged profile of the vertical component of the velocity normalized with respect to the sound speed for the 203 cells per HR. The vertical colored lines are the same as in Fig. 4. At HR = 6.2 the vertical velocity reaches 0.3cs. The hot ionized region is expected to start beyond HR > 6.2.
In the text
thumbnail Fig. 6

Time evolution of the local kinetic energy in the midplane for the 50, 101, and 203 cells per scale height resolutions at R0. The black dashed line represent the slope of the linear growth rate phase of the local energy in the midplane. The vertical gray lane represent the approximate time where the 101 and 203 cells per scale height resolutions reach saturation and the vertical dashed line represents the time where the 50 cells per scale height converges to stability.
In the text
thumbnail Fig. 7

Snapshot of vmag for the 203 cells per scale height (panel a) compared to the lowest resolution of 25 cells per scale height (panel b) at 200 orbits. The colored contour plot shows vmag and the streamlines show the velocity flow pattern when applying the LIC method. The vmag range (from 0 to 0.2 in code units) is narrowed down to highlight the low-velocity perturbations close to the midplane associated with the VSI.
In the text
thumbnail Fig. A.1

Radial column density in cm−2. Overplotted are solid lines where Σr = 1019 cm−2 represents the EUV heated region (in purple), and Σr = 1022 cm−2 represents the X-ray–FUV heated boundary layer (in dark green).
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.