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/00046361/202039294  
Published online  30 November 2020 
Hydrodynamical simulations of protoplanetary disks including irradiation of stellar photons
I. Resolution study for vertical shear instability
^{1}
MaxPlanck Institute for Astronomy,
Königstuhl 17,
69117
Heidelberg,
Germany
email: flores@mpia.de; flock@mpia.de
^{2}
RIKEN Cluster for Pioneering Research,
21 Hirosawa,
Wakoshi,
Saitama
3510198,
Japan
email: ryohei.nakatani@riken.jp
Received:
30
August
2020
Accepted:
11
October
2020
Context. In recent years hydrodynamical (HD) models have become important to describe the gas kinematics in protoplanetary disks, especially in combination with models of photoevaporation and/or magnetically driven winds. Our aim is to investigate how vertical shear instability (VSI) could influence the thermally driven winds on the surface of protoplanetary disks.
Aims. In this first part of the project, we focus on diagnosing the conditions of the VSI at the highest numerical resolution ever recorded, and suggest at what resolution per scale height we obtain convergence. At the same time, we want to investigate the vertical extent of VSI activity. Finally, we determine the regions where extreme UV (EUV), farUV (FUV), and Xray photons are dominant in the disk.
Methods. We perform global HD simulations using the PLUTO code. We adopt a global isothermal accretion disk setup, 2.5D (2 dimensions, 3 components) which covers a radial domain from 0.5 to 5.0 and an approximately full meridional extension. Our simulation runs cover a resolution from 12 to 203 cells per scale height.
Results. We determine 50 cells per scale height to be the lower limit to resolve the VSI. For higher resolutions, ≥50 cells per scale height, we observe the convergence for the saturation level of the kinetic energy. We are also able to identify the growth of the “body” modes, with higher growth rate for higher resolution. Full energy saturation and a turbulent steady state is reached after 70 local orbits. We determine the location of the EUV heated region defined by Σ_{r} = 10^{19} cm^{−2} to be at H_{R} ~ 9.7 and the FUV–Xray heated boundary layer defined by Σ_{r} = 10^{22} cm^{−2} to be at H_{R} ~ 6.2, making it necessary to introduce a hot atmosphere. For the first time we report the presence of smallscale vortices in the r − Z plane between the characteristic layers of largescale vertical velocity motions. Such vortices could lead to dust concentration, promoting grain growth. Our results highlight the importance of combining photoevaporation processes in the future highresolution studies of turbulence and accretion processes in disks.
Key words: hydrodynamics / instabilities / methods: numerical / protoplanetary disks
© L. FloresRivera et al. 2020
Open 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 highenergy radiation on protoplanetary disks can be quite complex, but recent efforts have attempted to explain the effect of highenergy 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 highenergy 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 raytracing 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 magnetocentrifugal mechanism are superKeplerian, while Wang et al. (2019) found it to be subKeplerian. 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), farultraviolet (FUV), and Xray 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.3^{1} 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 $$\frac{\partial \rho}{\partial t}+\nabla \cdot (\rho v)=0,$$(1) $$\frac{\partial (\rho v)}{\partial t}+\nabla \cdot (\rho vv)+\nabla P=\rho \nabla \Phi ,$$(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 $P={c}_{\text{s}}^{2}\rho $, where c_{s} 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°.
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 SolbergHø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 $$\rho (R,Z)={\rho}_{\text{0}}{\left(\frac{R}{{R}_{0}}\right)}^{p}\mathrm{exp}\left(\frac{GM}{{c}_{\text{s}}^{2}}\left[\frac{1}{\sqrt{{R}^{2}+{Z}^{2}}}\frac{1}{R}\right]\right),$$(3) $$\Omega (R,Z)={\Omega}_{\text{k}}{\left[(p+q){\left(\frac{{H}_{\text{R}}}{R}\right)}^{2}+(1+q)\frac{qR}{\sqrt{{R}^{2}+{Z}^{2}}}\right]}^{1/2},$$(4)
where ${\rho}_{\text{0}}=\frac{{\Sigma}_{\text{0}}}{\sqrt{2\pi}{H}_{\text{0}}{R}_{\text{0}}}$ = 4.5 × 10^{−10}g cm^{−3} = 2.7 × 10^{14} cm^{−3} is the initial density at the midplane, and Σ_{0} is the initial surface density (see Table 1), R_{0} is the reference radius, H_{0} is the reference scale height, G is the gravitational constant, M is the mass of the star, p is the powerlaw fitting exponent of the density profile, and ${\Omega}_{\text{k}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\sqrt{GM/{R}^{3}}$ is the Keplerian frequency. The disk scale height in terms of the radius is $${H}_{\text{R}}={H}_{\text{0}}{\left(\frac{R}{{R}_{\text{0}}}\right)}^{(q+3)/2},$$(5)
where the reference scale height ratio is H_{0} /R_{0} = 0.1 (see Table 1) and c_{s} is related to the temperature through the powerlaw fitting exponent q. The radial profile of c_{s} is given based on Eq. (5) so that c_{s} = H_{R}/Ω_{k}. The scale length of the vertical mode is proportional to H_{0}.
In order to compute the fluxes in each cell, we employ secondorder piecewise linear spatial reconstruction. This reconstruction method is appropriate for both uniform and nonuniform 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 secondorder RungeKutta algorithm and we set the CourantFriedrichsLewy (CFL) number to 0.25. Table 1 summarizes the model parameter and the different runs for each resolution. For 101 and 203 cells per H_{R}, the θ inner limit is shifted slightly inwards due to the high computational precision issue when increasing the resolution at the θ boundary. We investigated the thirdorder piecewise parabolic method (PPM3) spatial reconstruction using thirdorder RungeKutta (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 Xray 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 Xray photons, we calculate the radial column density locally to predict at what scale height the EUV and Xray–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 $${\Sigma}_{r}(r)={\displaystyle {\int}_{{R}_{\text{*}}}^{r}{n}_{\rho}}\text{d}r,$$(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 postprocessing. 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 (Xray–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} =10^{19} cm^{−2}, assuming the EUV stellar luminosity is ~ 10^{30} erg s^{−1}. The EUV heated region typically has temperatures of orders of magnitude between 10^{3} and 10^{4} K and the hydrogen is fully ionized here. The boundary of the Xray–FUV heated layer is defined at Σ_{r} =10^{22} cm^{−2}, where the gas is heated to 10^{3} K and it dominates in the neutral layers. While soft Xrays are considered to be an important component for photoevaporation, the locations of FUV and Xray heated layers are nearly identical Σ_{r} ≤10^{21−22} cm^{−2} (Gorti & Hollenbach 2008; Nakatani et al. 2018b). Therefore, setting the boundary of the neutral photoheated layer at Σ_{r} ≤10^{21−22} cm^{−2} is valid even in cases where Xrays effects are included. It should be noted that we neglect scattered light in our raytracing. 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 $${E}_{\text{kin}}=\frac{1}{2}{\displaystyle {\int}_{v}\rho}{v}_{Z}^{2}\text{d}V,$$(7)
where v_{Z} = v_{r} cos(θ) – v_{θ} sin(θ) is the vertical component of the velocity and E_{kin} is normalized with respect to $${E}_{\text{norm}}=\frac{1}{2}{\displaystyle {\int}_{v}\rho}{v}_{k}^{2}\text{d}V,$$(8)
where v_{k} 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 Xray–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 H_{R} = 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 Xray–FUV heated layer and EUV heated region reach H_{R} = 6.2 and H_{R} = 9.7, respectively.At this photoevaporation region we predict that HI regions become optically thick against EUV with a Σ_{r} = 10^{19} cm^{−2} to Σ_{r} = 10^{22} 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.
Fig. 1 Time evolution of the kinetic energy for different resolutions and space averaged at R_{0} = 1 AU. The kinetic energy is normalized with the kinetic energy from pure Keplerian velocity. The different resolutions are colorcoded (see inset). 
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 EUVheated region is defined by Σ_{r} =10^{19} cm^{−2} (in purple), and the column density of the Xrays/FUVheated boundary layer is defined by Σ_{r} = 10^{22} cm^{−2} (in dark green). The number density of hydrogen nuclei at the H/H^{+} boundary is ~10^{9} cm^{−3} and at the Xrays/FUVheated layer is ~10^{12} cm^{−3}. 
3.2 Local kinetic energy evolution over height and time
We show the timeevolution of the kinetic energy over H_{R} at R_{0} 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 R_{0}. 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 largescale 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 H_{R} = 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 hightemperature region caused by the FUV radiation. We also expect this position to be where the Xray–FUV heated boundary layer is located at H_{R} = 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 H_{R} around 6.2. In the midplane layers the time averaged vertical motions remain small, at a level of around 1%.
Figure 6 shows the timeevolution of the local kinetic energy at the midplane at R_{0}. 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.
Fig. 3 Timeevolution of the local kinetic energy at R_{0} = 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)). 
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 H_{R} at R_{0} = 1 AU. Thevertical colored lines represent H_{R} = 6.2, where the Xray–FUV heated boundary layer is (in dark green), and the H_{R} = 9.7, where the EUV heated region is located (in purple) (same as in Fig. 2). The models for 101 and 203 cells per H_{R} are in good agreement. The density floor is located H_{R} = 13, and therefore does not influence our result. 
Fig. 5 Averaged profile of the vertical component of the velocity normalized with respect to the sound speed for the 203 cells per H_{R}. The vertical colored lines are the same as in Fig. 4. At H_{R} = 6.2 the vertical velocity reaches 0.3c_{s}. The hot ionized region is expected to start beyond H_{R} > 6.2. 
Fig. 6 Time evolution of the local kinetic energy in the midplane for the 50, 101, and 203 cells per scale height resolutions at R_{0}. 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 Smallscale 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 ${v}_{\text{mag}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\sqrt{{v}_{\text{Z}}^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{v}_{\text{R}}^{2}}$, where v_{R} = v_{r} sin(θ) + v_{θ} cos(θ). From the LIC method, we find smallscale vortices in both resolution cases; however, these smallscale 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 longlived locally around 500 orbits (Manger & Klahr 2018) with an H_{R} = 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 largescale upward and downward motions. They might be related to KelvinHelmholtz instabilities appearing between these largescale motions, as was described by Latter & Papaloizou (2018). Moreover, these smallscale vortices happen in an H_{R} < 5, and are therefore more important for dust evolution.
4 Implications for planet formation and winds in disk
The VSI is a largescale 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 longlived 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 socalled 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 (H_{R} < 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 Frequencyintegrated cross sections
Perhaps the most robust way to define the influence of EUV and Xray–FUV photon layers is by taking into consideration the frequencydependent 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 10^{18} cm^{−2}. This atomic hydrogen column density can be translated into a hydrogen column density of 10^{19−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 10^{21} and 10^{22} 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 ~10^{4} cm^{2} g^{−1}, which is in the wavelength regime between 0.2 and 0.09 μm (6 eV < hν < 13.6 eV). Assuming a dusttogas mass ratio of 0.01, this gives a total FUV opacity of 100 cm^{2} g^{−1}. While it is more likely that the dusttogas 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 wavelengthdependent 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.
Fig. 7 Snapshot of v_{mag} 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 v_{mag} and the streamlines show the velocity flow pattern when applying the LIC method. The v_{mag} range (from 0 to 0.2 in code units) is narrowed down to highlight the lowvelocity 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 Xrays, the electron relative abundance between the EUV heated layer and the Xray–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 ~ 10^{4}–10^{5} cm^{−3} at ~1 AU (Nakatani et al. 2018a,b). We expect that at <1 AU, the Xray, 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 H_{R} = 13 at R_{0} = 1 for resolutions ≥50 cells per scale height. We expect turbulent scenarios where the Xray–FUV heated boundary layer at H_{R} = 6.2 is located, whereas the EUV heated region reaches H_{R} = 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 H_{R} 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 smallscale vortices between the largescale motion of the VSI when performing the LIC method. These smallscale vortices reach H_{R} < 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 H_{R} 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 Xray–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.
Fig. A.1 Radial column density in cm^{−2}. Overplotted are solid lines where Σ_{r} = 10^{19} cm^{−2} represents the EUV heated region (in purple), and Σ_{r} = 10^{22} cm^{−2} represents the Xray–FUV heated boundary layer (in dark green). 
References
 Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006a, MNRAS, 369, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006b, MNRAS, 369, 229 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Arlt, R., & Urpin, V. 2004, A&A, 426, 755 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Ballabio, G., Alexander, R. D., & Clarke, C. J. 2020, MNRAS, 496, 2932 [CrossRef] [Google Scholar]
 Banzatti, A., Pascucci, I., Edwards, S., et al. 2019, ApJ, 870, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Barker, A. J., & Latter, H. N. 2015, MNRAS, 450, 21 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Cui, C., & Bai, X.N. 2020, ApJ, 891, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Ercolano, B., & Pascucci, I. 2017, R. Soc. Open Sci., 4, 170114 [Google Scholar]
 Ercolano, B., Drake, J. J., Raymond, J. C., & Clarke, C. C. 2008, ApJ, 688, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, M., Pascucci, I., Edwards, S., et al. 2018, ApJ, 868, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2016, ApJ, 827, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Turner, N. J., Mulders, G. D., et al. 2019, A&A, 630, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155 [CrossRef] [Google Scholar]
 Fricke, K. 1968, Z. Astrophys., 68, 317 [NASA ADS] [Google Scholar]
 Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 [Google Scholar]
 Gorti, U., & Hollenbach, D. 2008, ApJ, 683, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Gorti, U., Dullemond, C. P., & Hollenbach, D. 2009, ApJ, 705, 1237 [NASA ADS] [CrossRef] [Google Scholar]
 Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126 [CrossRef] [Google Scholar]
 Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Klahr, H., Pfeil, T., & Schreiber, A. 2018, Handbook of Exoplanets (Cham: Springer), 138 [Google Scholar]
 Latter, H., & Papaloizou, J. 2018, MNRAS, 474, 3110 [CrossRef] [Google Scholar]
 Lin, M.K. 2019, MNRAS, 485, 5221 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K., & Youdin, A. N. 2015, ApJ, 811, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001 [NASA ADS] [CrossRef] [Google Scholar]
 Manger, N., & Klahr, H. 2018, MNRAS, 480, 2125 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A. 2007, ApJ, 170, 228 [Google Scholar]
 Mignone, A. 2014, Comput. Phys., 270, 784 [Google Scholar]
 Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018a, ApJ, 857, 57 [CrossRef] [Google Scholar]
 Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018b, ApJ, 865, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, A&A, 435, 2610 [Google Scholar]
 Ormel, C. W., Shi, J.M., & Kuiper, R. 2015, MNRAS, 447, 3512 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880 [NASA ADS] [CrossRef] [Google Scholar]
 Pascucci, I., Sterzik, M., Alexander, R. D., et al. 2011, ApJ, 736, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Pfeil, T., & Klahr, H. 2020, ApJ, submitted [arXiv:2008.11195] [Google Scholar]
 Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G., Arlt, R., & Shalybkov, D. 2002, A&A, 391, 781 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190 [CrossRef] [EDP Sciences] [Google Scholar]
 Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781 [Google Scholar]
 Simon, M. N., Pascucci, I., Edwards, S., et al. 2016, ApJ, 831, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tanaka, K. E. I., Nakamoto, T., & Omukai, K. 2013, ApJ, 773, 155 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Wang, L., Bai, X.N., & Goodman, J. 2019, ApJ, 874, 90 [NASA ADS] [CrossRef] [Google Scholar]
PLUTO 4.3 is an open source code available for download: http://plutocode.ph.unito.it/download.html
All Tables
All Figures
Fig. 1 Time evolution of the kinetic energy for different resolutions and space averaged at R_{0} = 1 AU. The kinetic energy is normalized with the kinetic energy from pure Keplerian velocity. The different resolutions are colorcoded (see inset). 

In the text 
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 EUVheated region is defined by Σ_{r} =10^{19} cm^{−2} (in purple), and the column density of the Xrays/FUVheated boundary layer is defined by Σ_{r} = 10^{22} cm^{−2} (in dark green). The number density of hydrogen nuclei at the H/H^{+} boundary is ~10^{9} cm^{−3} and at the Xrays/FUVheated layer is ~10^{12} cm^{−3}. 

In the text 
Fig. 3 Timeevolution of the local kinetic energy at R_{0} = 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 
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 H_{R} at R_{0} = 1 AU. Thevertical colored lines represent H_{R} = 6.2, where the Xray–FUV heated boundary layer is (in dark green), and the H_{R} = 9.7, where the EUV heated region is located (in purple) (same as in Fig. 2). The models for 101 and 203 cells per H_{R} are in good agreement. The density floor is located H_{R} = 13, and therefore does not influence our result. 

In the text 
Fig. 5 Averaged profile of the vertical component of the velocity normalized with respect to the sound speed for the 203 cells per H_{R}. The vertical colored lines are the same as in Fig. 4. At H_{R} = 6.2 the vertical velocity reaches 0.3c_{s}. The hot ionized region is expected to start beyond H_{R} > 6.2. 

In the text 
Fig. 6 Time evolution of the local kinetic energy in the midplane for the 50, 101, and 203 cells per scale height resolutions at R_{0}. 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 
Fig. 7 Snapshot of v_{mag} 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 v_{mag} and the streamlines show the velocity flow pattern when applying the LIC method. The v_{mag} range (from 0 to 0.2 in code units) is narrowed down to highlight the lowvelocity perturbations close to the midplane associated with the VSI. 

In the text 
Fig. A.1 Radial column density in cm^{−2}. Overplotted are solid lines where Σ_{r} = 10^{19} cm^{−2} represents the EUV heated region (in purple), and Σ_{r} = 10^{22} cm^{−2} represents the Xray–FUV heated boundary layer (in dark green). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.