Hydrodynamical simulations of protoplanetary disks including irradiation of stellar photons. I. Resolution study for Vertical Shear Instability (VSI)

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 magnetic-driven winds. We focus on diagnosing the the vertical extent of the VSI at 203 cells per scale height and allude at what resolution per scale height we obtain convergence. Finally, we determine the regions where EUV, FUV and X-Rays are dominant in the disk. 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. We determine the 50 cells per scale height to be the lower limit to resolve the VSI. For higher resolutions, greater than 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 the radial column density to be 10$^{19}$ cm$^{-2}$ located at $H_\mathrm{R}\sim9.7$, and the FUV/X-Rays-heated boundary layer defined by 10$^{22}$ cm$^{-2}$ located at $H_\mathrm{R}\sim6.2$, making it necessary to introduce the need of a hot atmosphere. For the first time, we report the presence of small scale vortices in the r-Z plane, between the characteristic layers of large scale vertical velocity motions. Such vortices could lead to dust concentration, promoting grain growth. Our results highlight the importance to combine photoevaporation processes in the future high-resolution studies of the turbulence and accretion processes in disks.


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 for the gas kinematics and the overall gas evolution in protoplanetary disks. In theory, the presence of the VSI (Goldreich, P. and Schubert, G. 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, R. P. and Gressel, O. and Umurhan, O. M. (2013); Stoll & Kley (2014); Barker, A.J. and Latter, H.N. (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, R. P. and Gressel, O. and Umurhan, O. M. 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, understanding the effect of the VSI together with realistic thermal profiles in the disk are not well known. The evolution of the VSI has being tested in the context of planet formation and disk evolution. Recent simulations conducted by Lin (2019) have found 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 (e.g., Shu et al. (1994); Hollenbach et al. (1994); Clarke et al. (2001); Owen et al. (2012); Alexander et al. (2006aAlexander et al. ( ,b, 2014; Ercolano & Pascucci (2017)) that photoevaporation by energetic photons is a potential agent of disk dispersal. Investigating the effect of high energy radiation onto protoplanetary disks can be quite complex, but recent efforts (e.g., Ercolano et al. (2008); Gorti et al. (2009) ;Owen et al. (2010Owen et al. ( , 2011) have attempted to explain the effect of high energy radiation onto the disk through observations of winds. A prominent signature of photoevaporation processes is the detection of forbidden lines in the protoplanetary disk spectrum. The detection of NeII line from the surface of TW Hya at around 10 AU (Pascucci et al. 2011), shows the important influence of high energetic photons onto the surface layers of the disk, however, the vertical extension of the wind is quite inconclusive yet.
Article number, page 1 of 9 arXiv:2010.06711v1 [astro-ph.EP] 13 Oct 2020 Moreover, Ballabio et al. (2020) provided 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, possible magneto driven winds (Fang et al. 2018;Banzatti et al. 2019). Photoevaporative process in the inner disk (¡ 5 AU) has been tested by Wang et al. (2019) and Gressel et al. (2020) which 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 magneto-centrifugal mechanism, are super-Keplerian while Wang et al. (2019) found to be sub-Keplerian. The mass loss in Gressel et al. (2020) is ∼ 10 −7 M yr −1 whereas Wang et al. (2019) The coexistence between the VSI and photoevaporation processes by stellar photons is of particular interest, especially in regions where we expect steep temperature profile and fast thermal relaxation timescales which few have found (e.g., Klahr et al. (2018); Lyra & Umurhan (2019); Pfeil & Klahr (2019) that the VSI might operate down to ∼1 AU. The effect of viscous heating in the disk remains relatively small, making it plausible to use passive irradiated disks models to accurately represent the thermal structure (Flock et al. 2019). As a first step towards to include photoevaporation processes in our HD simulations, our aim in this first part 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 follow: in §2 we present the methodology in which we describe the disk model and boundary conditions to resolve the VSI using a global 2D isothermal accretion disk configuration. In §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-Rays photons can have in the inner parts of the disk. In §4 we discuss and present a perspective of future work about the implications to planet formation and winds in disk. Finally, §5 we present our conclusions.

Methods
We performed our HD simulations using PLUTO 4.3 1 by Mignone (2007). The dynamic of 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: 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 2 s ρ, where c s is the isothermal sound speed. The code solve 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(θ) (e.g. Ormel et al. (2015)) for the meridional domain in order to have a better resolution

Disk model and Boundary conditions
Protoplanetary disks, rotating with Keplerian angular velocity and vertical stratification are dynamically stable according to Solberg-Høiland criteria (e.g. Rüdiger, G. et al. (2002)). The accretion disk setup in equilibrium in cylindrical coordinates (R, Z) (Nelson, R. P. and Gressel, O. and Umurhan, O. M. 2013) are defined: =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 Table1), R 0 is the reference radius, and H 0 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 Ω k = GM/R 3 is the Keplerian frequency. The disk scale height in terms of the radius where the reference scale height ratio is H 0 /R 0 = 0.1 (see Table 1). c s is related to the temperature through the power-law 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 . Overplotted are contour solid lines where the column density of the EUV-heated region is defined by Σ r =10 19 cm −2 (in purple), and the column density of the X-rays/FUV-heated 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 X-rays/FUV-heated layer is ∼10 12 cm −3 .
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 nonuniform grid spacing based on the approach from Mignone (2014). We use Harten, Lax, Van Leer (HLL) Riemann solver. For the time integration we adopt second-order Runge-Kutta and we set the Courant-Friedrichs-Lewy (CFL) number to 0.25. Table.1 summarize the model parameter and the different runs for each resolution. For the 101 and 203 cells per H R , the θ inner limit is shifted slightly inwards due to high computational precision issue when increasing the resolution at the θ boundary. We investigated the PPM3 spatial reconstruction using 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 not yet stable. Therefore, we only show our analysis and discussion for the VSI using the linear reconstruction method for this paper.

EUV,FUV and X-Rays 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-rays photons, we calculate the radial column density locally to predict at what scale height the EUV-and Xrays/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 where R * is the star radius, n ρ is the number density of hydrogen nuclei. For this first part, since we do not consider detailed chemistry or frequency of UVs, 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-rays/FUVheated layer) which is empirically determined from the simula- tions of Nakatani et al. (2018a,b). We select the boundary between the ionized or EUV-heated region namely 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 to 10 4 K and the hydrogen is fully ionized here. The boundary of the Xrays/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. Though, soft X-Rays are considered to be an important component for photoevaporation, the locations of FUV-and X-Rays-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 photo-heated layer at Σ r ≤ 10 21−22 cm −2 is valid even in cases where X-Rays effects are included. Note that we neglect scattered light from our ray-tracing. For more information about the effects of stellar photons and the launch of winds, see §4.1.1 and §4.1.2.

Dynamical evolution
In order to investigate the kinematics we start by looking at the time evolution of the kinetic energy. Fig.1 shows the averaged where v Z = v r cos(θ)v θ sin(θ) is the vertical component of the velocity and E kin is normalized with respect to 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 the 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, similar as found by Nelson, R. P. and Gressel, O. and Umurhan, O. M. (2013). The 25 cells per scale height case is often used in the literature as a standard resolution to perform global simulations. However, we emphasize that this resolution does not fully solve the VSI in accretion disks. It reaches an energy saturation slightly less than 10 −5 and after that, it starts to increase monotonically without reaching to 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  Fig.4. At H R =6.2, the vertical velocity reach to 0.3c s . We expect the hot ionized region to start beyond H R > 6.2.
orbits, when the fast growing finger modes have already saturated. During this second phase, the body modes are growing 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 level converge. For even higher resolutions, which are computationally very expensive, we predict to have the same convergence trend.
Comparing our total kinetic energy, it follows the same trend as the perturbed kinetic energy from Nelson, R. P. and Gressel, O. and Umurhan, O. M. (2013). The resolution for the VSI used by Cui & Bai (2020) lies in between 96 and 108 cells per scale height which is in good agreement to our resolution study as well. Fig.2 shows a snapshot of our 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-rays/FUVheated boundary layer are defined based on our definition of Σ r ( §2.2).
The density (Fig.2 panel (a)) has a smooth profile until reaching 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 X-rays/FUV-heated layer and EUV-heated region reach a H R = 6.2 and H R = 9.7, respectively. At this photoevaporation region we predict that HI regions becomes 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.2 panel (b)) shows clearly the upward motion shaded in blue and the downward motion shaded in red. In this region, the vertical velocity reach a fraction of the sound speed.
The kinetic energy profile (Fig.2 panel (c)) shows the highest concentration present as corrugated-shape like structures resembling the body modes in the disk. In the next section, we investigate this in details in a more local perspective.

Local kinetic energy evolution over height and time
We show the time-evolution 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 betwen -1 and 1 AU at R 0 . Fig.3 shows that the kinetic energy is growing first the finger modes in the upper layers of 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 are threatening the upper and lower hemisphere. At this stage also the kinetic energy at the midplane grows.
Two important features are depicted from all resolution cases in Fig.3. The first one is the presence of vertical lanes of kinetic energy initially. We refer 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. Fig.4 shows the time averaged vertical profile of the kinetic energy. The vertical profile follow a similar trend as the density profile. Only for the lowest resolution case, 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 which is caused by the FUV radiation. We also expect this position to be where the X-rays/FUV-heated boundary layer is located at H R =6.2 (see Fig.4 & Fig.5). Fig.5 shows that the time averaged vertical velocity increases to levels of several tens of percent at the expect wind base at H R around 6.2. In the midplane layers, the time averaged vertical motions remain small, on a level around 1%. Fig.6 shows the time-evolution 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 that 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 and 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, R. P. and Gressel, O. and Umurhan, O. M. (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 higher than 50 cells per scale height to solve the VSI in accretion disks.

Small scale vortices
In Fig.7, we show the result of the Line Integral Convolution (LIC) method (e.g. Cabral & Leedom (1993)) to visualize the velocity magnitude vector field associated to velocity perturbations from the VSI, for the lowest and highest resolutions. The velocity magnitude includes the velocity field in every compo- where v R =v r sin(θ) + v θ cos(θ). From the LIC method, we notice small scale vortices in both resolutions 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, Samuel and Nelson Richard P., and Umurhan Orkan M. (2016) in the r − φ plane. These vortices can be long-lived 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 large scale upward and downward motions. Moreover, these small scale vortices happen in an H R < 5, therefore, being more important for dust evolution.

Discussion of implications to planet formation and winds in disk
The VSI is a large scale global instability that potentially plays an important role for the evolution and dynamics of protoplanetary disks. In the midplane, the dust settling overcomes the VSI, and as the dust concentrates and grows to mm size and greater, these particles settle against VSI and eventually undergoes 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 shows 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.

Importance of stellar photons layers and the launch of winds
The VSI can maintain small particles (< 100 µm) lofted (Lin 2019) until eventually becoming ionized by FUV photons at the disk surface and contributing to heating in the atmosphere. Both, the photoionization effects and the high temperature in this region, makes 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 in fact VSI and MHD winds can coexist with magnetized winds the main driver of angular momentum transport. However, this would rather depend on 1) the so-called wind base which depend 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 deter-mine the location and extension of the wind in the future work. The detailed role of the VSI for the mechanism of wind launching has to be investigated, especially due to the fact that the magnetic coupling to the gas remains low in the FUV shielded region H R < 6.2. Recent improvements, towards including a more realistic scenarios, by Wang et al. (2019) and Gressel et al. (2020) have included thermochemistry and irradiation in their global magnetohydrodynamical simulations.

Frequency-integrated cross-sections
Perhaps the most robust way to define the influence of EUV and X-rays/FUV photons layers is by taking into consideration the frequency-dependent cross-sections for each photon respectively. 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, that is influenced by EUV photons. The atomic hydrogen abundance in this ionized region is negligible up in the ionization front, where we expect the 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 rather 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, used here to define the EUV-heated region. Given the fact that FUV photons can promote heating via photoelectric heating (Bakes & Tielens 1994) in the disk atmosphere, such photons are attenuated by dust once the hydrogen column density is in between 10 21 cm −2 to 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 that is in the wavelength regime between 0.2 to 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 cm 2 g −1 .
Though, it is more likely that the dust-to-gas mass ratio is rather 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 disks 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 for the extension of the wind in the disk atmosphere.

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 and the X-rays/FUV-heated layer is still ∼ 10 −4 that is coming 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 X-Rays, EUV, and FUV fluxes are higher, therefore, photons penetrate deeper in the disk and the electron density can be higher. On the other hand, the farther from the star, the electron abundance decreases and so the electron density.

Conclusions
We have 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 to be the lower limit to resolve the VSI. This determination is based on the total kinetic energy saturation level. Our key points can be summarize as 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 X-rays/FUV-heated boundary layer at H R = 6.2 is located, whereas, the EUV-heated region reach a H R = 9.7. Nevertheless, we encourage the use for high resolution (≥ 50 cells per scale height) to properly conduct 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 small scale vortices between the large scale motion of the VSI when performing the LIC method. These small scale vortices reach a H R < 5.0, therefore, being 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 suitable for future studies.