The Effect of Ultraviolet Photon Pumping of H2 in Dust-deficient Protoplanetary Disks

We perform radiation hydrodynamics simulations to study the structure and evolution of a photoevaporating protoplanetary disk. Ultraviolet and X-ray radiation from the host star heats the disk surface, where H2 pumping also operates efficiently. We run a set of simulations in which we varied the number of dust grains or the dust-to-gas mass ratio, which determines the relative importance between photoelectric heating and H2 pumping. We show that H2 pumping and X-ray heating contribute more strongly to the mass loss of the disk when the dust-to-gas mass ratio is D≤10−3 . The disk mass-loss rate decreases with a lower dust amount, but remains around 10−10−11 M ⊙yr−1. In these dust-deficient disks, H2 pumping enhances photoevaporation from the inner disk region and shapes the disk mass-loss profile. We thus argue that the late-stage disk evolution is affected by the ultraviolet H2 pumping effect. The mass-loss rates derived from our simulations can be used in the study of long-term disk evolution.


INTRODUCTION
Circumstellar disks around newly born stars evolve to be protoplanetary disks (PPDs) where planets are formed.Near infrared observations show that the disk fraction in a star-forming region decreases as the age of the central star, suggesting that PPDs disappear in 3-6 Myr (e.g., Haisch et al. 2001;Meyer et al. 2007;Hernández et al. 2007;Mamajek 2009;Bayo et al. 2012;Ribas et al. 2014).
Photoevaporation is suggested to be an important disk dispersal mechanism.High-energy photons such as far-ultraviolet (FUV; 6 eV ≲ hν < 13.6 eV), extremeultraviolet (EUV; 13.6 eV ≲ hν ≲ 100 eV), and X-ray (100 eV ≲ hν ≲ 10 keV) emitted from the central star heat a large portion of a PPD through various radiative and photo-chemical processes.It is generally found that EUV photons are mainly absorbed by abundant neutral hydrogen close to the central star, while FUV and Xray photons penetrate deep into the disk to drive dense ayano.komaki@phys.s.u-tokyo.ac.jp photoevaporative flows (e.g., Richling & Yorke 2000;Ercolano et al. 2009;Gorti & Hollenbach 2009).FUV photons are the cause of the photoelectric effect on dust grains in the disk, whereas EUV and X-ray contribute to the heating through ionization of hydrogen, helium, and other heavy elements.These processes altogether raise the gas temperature of the disk surface and generate photoevaporative flows.
Previous studies showed that the strength of photoevaporation is determined by a number of physical properties such as metallicity, stellar mass, luminosity and chemical composition (Ercolano et al. 2009;Nakatani et al. 2018a,b;Wölfer et al. 2019;Komaki et al. 2021).For example, Nakatani et al. (2018a) performed radiation hydrodynamics simulations with varying the gas metallicity, and showed that photoevaporation is suppressed in low-metallicity disks because of the low abundance of dust grains.The mass-loss rate and the overall disk evolution are highly sensitive to the effective disk heating rate, and thus it is important to perform direct numerical simulations in order to study the structure of a photoevaporating disk and to derive realistic physical and chemical profiles under a variety of environments.
The disk dust properties are characterised by the dustto-gas mass ratio and the size distribution.Since these quantities may be inhomogeneous in time and in space within the disk, the relative importance of the heating processes can also significantly change according to the disk evolution.Dust grains undergo physical processing during the disk evolution.Numerical calculations show that small dust grains entrain photoevaporative flows (Hutchison et al. 2016).Dust settling may occur in PPDs, as inferred from infrared observations and SED modelling (D'Alessio et al. 1999(D'Alessio et al. , 2006;;Grant et al. 2018).Dust sedimentation can also significantly reduce the local dust-to-gas mass ratio of the disk surface.Millimeter observations suggest that the dust mass decreases as the system age increases (Mathews et al. 2012;Ansdell et al. 2016;Pascucci et al. 2016;Ansdell et al. 2017).The total dust-to-gas mass ratio of a disk likely varies through its evolutionary phases.Gorti et al. (2015) perform a set of one-dimensional disk evolution simulations by incorporating viscous accretion, photoevaporation, and dust evolution.They show that the dust evolution results in the increase of small dust grains, which affects the dust attenuation, photoelectric heating, dust-gas collision, and formation of molecules.Clearly, dust grain growth is an important process for disk evolution, and the amount of dust can critically affect the disk structure and dispersal process.Unfortunately, previous studies do not address explicitly how the reduced amount of dust grains, i.e., low dustto-gas ratios, affect photoevaporation and disk dispersal efficiency.
In this paper, we study the detailed process of photoevaporation and the dispersal of dust-deficient disks.Following the basic methodology of Komaki et al. (2021), we perform a set of radiation hydrodynamics simulations.Our simulations incorporate the so-called UV pumping of H 2 molecules (H 2 pumping hereafter), which can be an important heating process in the region where a substantial amount of H 2 molecules are continuously formed and destroyed.H 2 pumping is a two-step process; Lyman-Werner (LW) photons in the energy range 11.2 eV ≲ hν < 13.6 eV excite H 2 molecules to the electronically excited state, followed by de-excitation to vibrationally excited states of the ground electronic state (Tielens & Hollenbach 1985).Collisional de-excitation to the vibrational ground states results in heating of the gas by distributing the de-excitation energy to the surrounding atoms and molecules.Note that, while the FUV photoelectric heating rate roughly scales with the local amount of small dust grains, the H 2 pumping operates even in regions with little dust content.In order to examine how the disk's thermo-chemical structure changes with dust amount, we perform radiation hydrodynamics simulations, for the first time, with varing the dust-to-gas mass ratio of the disk, D, in the range of 10 −6 -10 −1 .
H 2 pumping is incorporated to closely examine the effect of H 2 pumping on the disk temperature.We explain our methods in Section 2. We show our results in Section 3 and discussions in Section 4. In Section 5 we summarize.

NUMERICAL SIMULATIONS
We perform disk photoevaporation simulations including hydrodynamics, radiative transfer and nonequilibrium chemistry in a coupled and self-consistent manner.We use the open source code PLUTO (Mignone et al. 2007) suitably modified for simulations of PPDs.The details of the implemented physics are found in Nakatani et al. (2018a,b) and Komaki et al. (2021).We assume that the disk is axisymmetric around the rotational axis and adopt the 2D spherical coordinates (r, θ).We consider three components of the gas velocity, We first run our fiducial case, where the stellar mass M * = 1 M ⊙ and D = 10 −2 , with a detailed treatment of H 2 pumping.In the case of D = 10 −2 , we calculate the dust temperature in the same way as in Komaki et al. (2021), where we calculate and tabulate the dust temperature as a function of the distance from the central star and the column density, based on the results of more costly radiative transfer calculations.In our runs with other dust-to-gas mass ratios, we solve the dust temperature in a self-consistent manner by incorporating both the direct and diffuse radiation.We perform ray-tracing for the direct component and adopt the fluxlimited-diffusion (FLD) approximation for the diffused component (Kuiper et al. 2010;Kuiper & Klessen 2013;Kuiper et al. 2020).
The governing equations are given by The first five equations express fluid equations, and the last one coupled with the rate equations expresses the abundance evolution of chemical species.We use ρ, v, P , and M * to represent gas density, velocity, pressure, and stellar mass, respectively.The fifth equation(Eq.( 5)) is the energy equation where E and H express the total energy and enthalpy, including the kinetic energy per unit volume.In the equation, Γ and Λ are the specific heating rate and the specific cooling rate, which means the heating and cooling rate per unit mass.The azimuthal component of the Euler equation is described in the angular momentum conservation form(Eq.( 6)).We define y i as the abundance of chemical species i, n H as the number density of elemental hydrogen, and R i as the relevant chemical reaction rate.
We follow the chemical reactions tabulated in Nakatani et al. (2018a,b).We use the parallel simulation code PLUTO (Mignone et al. 2007) to solve hydrodynamics.We incorporate EUV/X-ray photoionization heating, FUV photoelectric heating (Bakes & Tielens 1994), and heating by H 2 photodissociation and by H 2 pumping.We also incorporate dust-gas collisional cooling (Yorke & Welz 1996), fine-structure cooling of C II and O I (Hollenbach & McKee 1989;Osterbrock 1989;Santoro & Shull 2006), molecular line cooling of H 2 and CO (Galli & Palla 1998;Omukai et al. 2010), hydrogen Lyman α line cooling (Anninos et al. 1997), and radiative recombination cooling (Spitzer 1978) as cooling sources.We calculate the C II and O I cooling rates assuming the LTE level populations.
The chemical reaction network follows the abundances of eleven chemical species: H I, H II, H -, He I, H 2 , H 2 + , H 2 * , CO, O I, C II and electrons in the run with M * = 1 M ⊙ , D = 10 −2 .We denote the vibrational excited state of H 2 (v = 6) as H 2 * and treat as a distinct chemical species (Tielens & Hollenbach 1985).H 2 has hundreds of excited states, but the pseudo-level of v = 6 approximately represents the excited levels altogether, which then reduces the computational cost significantly.We also use the pseudo-level of H 2 * to study the effect of H 2 pumping on the disk thermal structure.We assume that C I is quickly ionized to become C II after the dissociation of CO (Nelson & Langer 1997;Richling & Yorke 2000).
We incorporate H 2 * to calculate both the excitation rate and the de-excitation rate.We calculate the H 2 pumping rate, R pump as R pump = 3.4 × 10 −10 βG 0 e −2.5Av s −1 , following Tielens & Hollenbach (1985).We define β, G 0 and A v as the shielding factor, the intensity of the incident FUV field and the visual extinction respec-tively.We determine the shielding factor by H 2 following Draine & Bertoldi (1996).The H 2 molecules in the vibrational excited state de-excite and fluoresce.We incorporate the collisional de-excitation by H I or H 2 as the de-excitation reactions.We also consider the direct dissociation of H 2 * by FUV radiation and by spontaneous radiation (Tielens & Hollenbach 1985).We define k de (H), k de (H 2 ), R de , A(H * 2 ) as the reaction rates by the collision with H I, the collision with H 2 , the direct dissociation by FUV radiation and the spontaneous radiation respectively.The reaction rates are written as where T gas represents the gas temperature (Wang & Goodman 2017).Each collisional de-excitation process deposits 2.6 eV and each de-excitation by direct FUV radiation deposits 0.4 eV (Tielens & Hollenbach 1985).
We explicitly treat H 2 * as a distinct chemical species.We have run an additional simulation with a simpler implementation for comparison and testing purpose as follows.In the fiducial run with M * = 1 M ⊙ and D = 10 −2 , we calculate the heating rate by H 2 pumping without explicitly including the excited molecular hydrogen, as in Röllig et al. (2006), Gressel et al. (2020) and Nakatani et al. (2021).Specifically, we calculate the heating rate as follows.
where χ is the FUV radiation intensity in units of the local interstellar radiation field, P tot and n are the formation rate of H 2 * and the particle number density, respectively.The coefficients ∆E eff , A eff , D eff , γ eff express the average excitation energy, the effective spontaneous emission rate, the effective dissociation rate, and the effective collisional de-excitation rate assuming that molecular hydrogen has two levels, the ground state and the pseudo-excited level.We set P tot ∆E eff = 9.4 × 10 −22 erg s −1 , A eff = 1.9 × 10 −6 s −1 , D eff = 4.7 × 10 −10 s −1 and γ eff = 5.4 × 10 −13 T gas s −1 cm −3 following Röllig et al. (2006) and Gressel et al. (2020).Röllig et al. (2006) derived Equation ( 11) by fitting the heating rate calculated in the thermodynamical simulations that consider 100 energy levels of H 2 .We have tested and checked that using the fitting formula does not change the resulting mass-loss profiles appreciably, while reducing the computational time by ∼ 1.5 times.
We compare the result further in detail in Section 3.1.
In the following, we shall refer to the heating rate for H 2 pumping (Eq.[11]) as Röllig method.
The size distribution of dust grains is assumed to follow dn(a) ∝ a −3.5 da with a denoting dust radii (Draine & Lee 1984).Bakes & Tielens (1994) show that small grains are crucial for efficient photoelectric heating.In practice, dust dynamics such as sedimentation, radial drift and/or turbulent mixing can change the dust distribution across the disk.However, for simplicity, we set the dust-to-gas mass ratio to constant throughout the computational domain.The underlying assumption is that the dust grains responsible for photoelectric heating are well mixed within the gas so that we treat the mixture of dust and gas as a single fluid.We perform a set of simulations with different constant values of the dust-to-gas mass ratio.Note that the dust amount near the disk surface, where photoevaporative flows are launched, critically determines the disk mass-loss rate.We perform disk evolution simulations with varying the dust-to-gas mass ratio of the disk D in the range of 10 −6 -10 −1 as a parameter to examine how the effective heating and cooling processes vary with D. When we adopt lower D, we assume that the dust grains of all the sizes are depleted equally.Since photoelectric heating is largely caused by dust grains smaller than ∼ 20µm, the value of D in the present study is meant to reflect the amount of small dust grains that contributes to the photoelectric heating.In our simulations, the dust opacity, photoelectric heating rate, dust-gas collisional cooling rate, the rate of H 2 formation catalyzed by grains are scaled with (D/0.01).Our fiducial dustto-gas mass ratio is D = 10 −2 , which corresponds to the local ISM value.The gas-phase elemental abundance of carbon is set to y C = 0.927 × 10 −4 and that of oxygen to y O = 3.568 × 10 −4 (Pollack et al. 1994).We assume that y C and y O are constant even when we set different D.
We take into account FUV, EUV and X-ray radiations as the high-energy radiation from the central star.As in our previous paper Komaki et al. (2021), we choose the FUV, EUV and X-ray luminosities, ϕ EUV , L FUV , L X , following Table 1 in Gorti & Hollenbach (2009), which shows typical 1 Myr pre-main sequence star luminosities.
We define the computational domain of the polar angle to [0, π/2], and the radial coordinate to [0.1r g , 20r g ].The gravitational radius r g for a fully ionized gas with where G is the gravitational constant.The gravitational radius sets the estimate of the effective boundary where the photoevaporative flows can escape out of the potential well of the host star.However hydrodynamics simulations show that the photoevaporative flows are also excited inside the gravitational radius from the critical radius ∼ 0.2r g (Liffman 2003;Font et al. 2004).Motivated by findings, we set the computational range to include the region inside the graviational radius.We use r g as a reference physical length scale in this paper.We start our simulations setting the initial disk condition to be a 1 million year system, when the mass loss by accretion becomes lower (Clarke et al. 2001;Alexander et al. 2006;Owen et al. 2010;Suzuki et al. 2016;Kunitomo et al. 2020).We assume that the disk mass, M disk , is 3% of the central stellar mass (Andrews & Williams 2005).We assume the initial surface density profile, Σ and the initial temperature distribution, T ini , as Komaki et al. (2021), to be We run the simulations for 8.4×10 3 (M * / M ⊙ ) yr.This is 10 times of the gas Kepler rotational period around the central star at r = r g and this is long enough for the time-averaged mass-loss to converge.

RESULTS
In this section, we present the results of our simulations.We first describe the effect of H 2 pumping on the temperature and the chemical structure in the run with M * = 1 M ⊙ and D = 10 −2 in Section 3.1.We then calculate the mass-loss rate of this run in Section 3.2.Finally in Section 3.3, we describe how the disk thermal structure changes when varying the dust-to-gas mass ratio.

Thermochemical Structure
Figure 1 shows the time-averaged snapshot of the simulation with M * = 1 M ⊙ and D = 10 −2 and with incorporating H 2 * and H 2 pumping.To make this plot, we average the outputs from 840 yr to 8400 yr to avoid the initial transient and to smooth out (temporarily) fluctuating features.Photoevaporative flows are launched from the surface of the disk, where the column density satisfies N H2 ≈ 10 20 cm −2 (white line in the top panel).We define this surface as the base of the photoevaporative flows.EUV radiation heats the gas near the central star and at high latitudes, shaping the polar H II region where y HII > 0.5.The outgoing gas cools by adiabatic cooling and has a temperature of T gas ∼ 3000 K.
The bottom portion of Figure 1 shows the H 2 pumping heating rate Γ H * 2 and also the molecular fraction y H2 .The H 2 * is most abundant near the dissociation front with a number fraction ∼ 0.04, where Γ H * 2 reaches ∼ 10 4 erg s −1 g −1 .We note that the H 2 dissociation front locates at a lower polar angle than the base.The temperature of the photoevaporative flows is determined by the balance between FUV photoelectric heating and H 2 cooling in the neutral region.At the base of the outflow, OI cooling is the main cooling process.
In the dense inner region at r ≲ 30 au, H 2 * molecules get de-excited mainly by collisions with neutral hydro-gen atoms.There, the gas is largely heated by this H 2 pumping process.The critical density for H 2 * spontaneous radiation is ρ c ≃ 5.0 × 10 −20 g cm −3 assuming T gas = 3000 K.In our simulation, the density at r ≳ 30 au is lower than ρ c , and thus the spontaneous radiation primarily de-excites the H 2 * molecules there, and the resulting gas heating rate is small.We note that, in our simulations, the timescale for H 2 * de-excitation is shorter than the local hydrodynamical timescale by a factor of ∼ 10 2 , so the H 2 * molecules almost immediately get de-excited at the location where they are excited by the incident radiation.
We have run an additional simulation with adopting Röllig method, which does not explicitly treat H 2 * molecules.By closely examining the outputs of the two simulations, we find that the disk structure is essentially the same in regions where the gas temperature is 300 K ≤ T gas ≤ 3000 K.Such regions are found typically in the vicinity of the disk surface as shown in Figure 1, where strong photoevaporative flows are driven.Noticeable differences are found in low temperature regions with T gas ≤ 300 K, where radiative heating is unimportant.We thus conclude that our main results on disk photoevaporation are not significantly affected by the choice of the two calculation methods of H 2 pumping.

The Mass-loss Rate
The disk gas is heated by the radiation from the central star, and photoevaporateive flows are launched from the disk.The disk mass-loss rate is determined by the heating process around the base.In our fiducial run, FUV photoelectric heating is the main heating source at the base, and its rate is greater than the heating rate of H 2 pumping by a factor of ∼ 10 3 throughout the base.The resulting temperature distribution is almost the same as in the run without H 2 pumping at the base, indicating that the photoevaporative flows are driven mainly by FUV photoelectric heating.
We use the simulation outputs to calculate the massloss rate Ṁ = where dS represents the spherical surface unit area at r = 100 au.We define η as the Bernoulli function given by We define v p = v 2 r + v 2 θ as the poloidal velocity, γ as the specific heat ratio, and c s as the sound velocity.We assume that only the gas satisfying η > 0 has sufficient mechanical energy to escape from the system eventually.
The mass-loss rate slightly fluctuates in time mainly because of spurious reflection of subsonic gas at the outer computational boundary.The averaged mass-loss rate (from 840 yr to 8400 yr) is The magnitude of the fluctuation is from −13% to +10% compared to the averaged value.The mass-loss rate with H 2 pumping is almost the same within the fluctuation as that without H 2 pumping ( Ṁ = 2.6 × 10 −9 M ⊙ yr −1 ).This is again because H 2 pumping heats the gas mainly in the vicinity of the H 2 dissociation front and does not affect the temperature structure in the photoelectricheated region.Hence, the effect of H 2 pumping on driving photoevaporative flows is limited if H 2 -rich flow is excited by photoelectric heating.
We also calculate the mass loss rate for the simulation with Röllig method.The obtained time-averaged Ṁ differs only by 4% from our fiducial case.This is consistent with the already discussed finding that the accurate treatment of H 2 pumping is important only in low-temperature, neutral regions in the disk.In the following sections, we use simulations that include H 2 * molecules explicitly.

Dependence on the Dust-to-Gas Mass Ratio
We perform a series of photoevaporation simulations with varying the dust-to-gas mass ratio in the range D = 10 −6 -10 −1 .Figure 2 shows the specific heating and cooling rates at the base in each simulation.In the case of D = 10 −1 , since the dust optical depth is high, FUV photons cannot deeply penetrate into the disk compared to the fiducial case (D = 10 −2 ).Photoevaporation is caused by FUV photoelectric heating at N H2 = 10 19 cm −2 , which is one order of magnitude lower than the D = 10 −2 case.(We define this surface as the base for this run.)The mass-loss rate decreases because of the lower base density (Figure 3).
The FUV photoelectric heating rate approximately follows Γ FUV ∝ D at the base.For D = 10 −3 , it gets lower to be comparable to heating rates for X-ray heating and H 2 pumping.At the outflow base, the main cooling process is OI cooling for low dust content with D ≤ 10 −3 .In our simulations, we do not vary the gas metallicity when we decrease the dust-to-gas mass ratio, and thus the relative importance of gas-phase metals increases for lower D ≤ 10 −3 .The temperature of the H 2 region can not get sufficiently high to have η > 0. It is in contrast to the fiducial run where molecular flow is observed.The mass-loss rate is lower than that of the fiducial run accordingly.
For D = 10 −4 -10 −6 , the FUV photoelectric heating is ineffective, and the thermochemical structure differs from the higher-D runs.EUV photons heat the gas at high latitudes, whereas X-ray heating is dominant in the H I region instead of FUV photoelectric heating.The gas temperature is lowered to 300 K by adiabatic cooling and O I line emission there.H 2 pumping is generally dominant only near the dissociation front, and the other H 2 region is heated by X-ray up to N H ∼ 10 22 cm −2 .Figure 2 shows the specific heating and cooling rates along the base.The relative contributions of the heating processes vary with the distance from the host star.H 2 pumping and X-ray are dominant at the radii inner and outer than 5-7r g , respectively.The inner regions are favorable for H 2 pumping since the density exceeds the critical value.The gas temperature is 200-300 K in the molecular region at r = 10r g .It is lower by a factor of two than would have been achieved by photoelectric heating.The same is true at other radii, and thus the η = 0 boundary is present at two orders of magnitude smaller column density, N H2 = 10 18 cm −2 , than the D = 10 −2 case.The flow accordingly has a low density.It results in the significant decrease in the mass-loss rate compared to the D ≥ 10 −3 runs (see Figure 3).
In order to extract the contribution of EUV-driven flows from the total mass-loss rates in the low-D cases, we count only the mass flux of the ionized outflows (y HII ≥ 0.5) using Eq.12.The mass-loss rates of the ionized gas are found to be roughly an order of magnitude smaller than the total, meaning that the contribution of EUV-driven flows is limited.This small contribution is mainly due to the fact that H 2 pumping increases the local scale height of the neutral gas at the innermost radii, and the vertically extended gas shields EUV that would otherwise have reached the outer radii.As a result, the H II region is confined within the conical volume at low polar angle.The total mass-loss rates are measured in the same manner as in Section 3.2.
The mass-loss rate is approximated as Ṁ ∼ 6.4 × 10 −0.23(log 10 D) 2 −0.73(log 10 D)−10 M ⊙ yr −1 (13) in the cases with D ≥ 10 −3 .In the case with D = 10 −4 , 10 −5 , the gas temperature is lower than the dustabundant cases by the factor of ∼ 2, the mass-loss rate becomes low accordingly.In the case with D = 10 −6 , the gas is not heated by dust-gas collisional heat transfer at the midplane because of the poor dust abundance.The scale height becomes low and the gas gathers in the lowlatitude region consequently.The base density and the mass-loss rate becomes larger.Since the abundance of H 2 does not strongly depend on the amount of small dust grains, we expect that the mass-loss rate does not vary from the case with D = 10 −6 towards even lower values of D.
We note that the actual spectrum includes emission line features in addition to the continuum FUV field.Observations towards T Tauri stars suggest that hydrogen Lyman-α and C IV line photons cause H 2 pumping (e.g., Herczeg et al. 2002Herczeg et al. , 2004Herczeg et al. , 2006;;Yang et al. 2011).These photons have energies lower than 11.2 eV but can effectively pump H 2 * molecules in the vibrationally excited states.Hence strong emission lines can increase the direct photodissociation rate R de given by Equation (9).Since the direct dissociation rate and the associated heating rate are proportional to the FUV flux (G 0 ), we can examine the pumping effect by increasing R de or equivalently G 0 .To this end, we have performed a test calculation for the case of D = 10 −3 with a ten times larger R de .The result shows that the enhanced direct dissociation of H 2 * does not contribute significantly to the heating rate because de-excitation of H 2 * by collisions with hydrogen atoms is dominant in the inner disk region.The outer part is primarily heated by X-ray photons (Figure 2), and the boosted pumping is not important there.Lyman-α and C IV photons may affect the abundance of H 2 * in the vicinity of the central star, but do not affect the gas temperature nor the total mass-loss rate.In our future work, we will incorporate detailed radiation transfer with realistic spectra with emission lines for direct comparison with spectral observations.Clarke et al. (2001) suggest that the disk dispersal is driven by a combination of accretion and photoevaporation.In their simulation, the disk gas is dispersed quickly after a gap opens in the disk, and photoevaporation remains effective at the outer region.Kunitomo et al. (2020) use a set of disk evolution simulations which incorporate accretion, magnetohydrodynamic (MHD) winds and photoevaporation, to study the respective contribution to disk dispersal.It is necessary to perform disk evolution simulations with all the relevant processes implemented accurately, in order to compare the numerical results with observations.Our simulations allow us to use the detailed structure of evaporating disks for the study of long-term disk evolution, and thus enable comparison of the disk dispersal time with observationally inferred disk lifetimes (Komaki et al. 2023).Nakatani et al. (2018b) investigated the dependence of the mass-loss rate on the gas metallicity Z.They varied Z with setting the dust-to-gas mass ratio as D = 0.01 × (Z/Z ⊙ ).The gas-phase elemental abundances of metals are also scaled in proportion to (Z/Z ⊙ ).On the other hand, we fix the gas-phase elemental abundances of carbon and oxygen constant regardless of D in the present study.Our run with D = 10 −2 corresponds to the Z/Z ⊙ = 1 case of Nakatani et al. (2018b), which assumes ∼ 10 times higher luminosities for FUV and EUV, and ∼ 2.5 times higher X-ray luminosity.The overall stronger radiation yields a higher mass-loss rates in all the metallicity range.The difference between the low-Z runs of Nakatani et al. (2018b) and our low-D runs appears in the major cooling process.For example, with small D, O I line cooling is more effective than dustgas collisional cooling, which has not been observed in Nakatani et al. (2018b).
Wang & Goodman (2017) performed hydrodynamics simulations with the same disk mass, the same central stellar mass and the same stellar luminosities as in the present study.They assume that the disk size is 100 au and the abundance of small dust grains is fixed to be 7 × 10 −5 .The corresponding total dust-to-gas mass ratio is D = 2.2 × 10 −3 by assuming the MRN dust size distribution.In order to compare their results directly with our simulations, we estimate the mass-loss rate by using Eq. 13, to obtain Ṁ = 1.7 × 10 −9 M ⊙ yr −1 , which is consistent with the result of Wang & Goodman (2017).In their simulation, the disk surface is heated to T gas > 10 4 K by EUV photons.They assume that each EUV photon has energy of 25 eV, while we consider a spectral shape of high-energy radiation.Photons with higher energy tend to penetrate deeper into the gas and give higher energy to the gas.This leads to the distinct thermal structure with the broader H II region.

DISCUSSION
While the results presented in the previous section have a variety of interesting implications, there are still uncertainties originating from the model assumptions.In this section, we further discuss possible variations of our results and show caveats.

Dust Evolution
It is expected that dust sedimentation causes the local dust-to-gas mass ratio to vary in the vertical direction.Dust grains can also drift radially toward the inner region.In the present study, we set the dust-to-gas mass ratio D as a convenient parameter to quantify and examine the effect of dust abundance on disk photoevaporation.We have shown that the mass-loss rate and the mass-loss profile change sensitively depending on D.
Infrared and millimeter observations suggest that PPDs contain dust grains with various sizes (Acke et al. 2004;Pascual et al. 2016;Davies et al. 2018).The dust size distribution and its variation have considerable effects on disk heating and the resulting mass-loss rate.Small grains have large surface areas per mass, and thus contribute predominantly to the net photoelectric heating.If a PPD achieves a state deficient in small dust grains, it is less susceptible to FUV photoelectric heating.Hutchison et al. (2016) show that the dust size differs between wind regions and the disk because small dust grains samaller than 0.3 µm diameter are entrained by photoevaporating flows and the mass-loss rate is two magnitudes smaller than the gas mass-loss rate.We estimate the effect of dust entrainment in the case where the amount of dust grains decreases for ∼ 3 Myr.If we assume a mass loss rate of Ṁ ∼ 10 −9 M ⊙ yr −1 (Figure 3), the disk gas mass decreases by 10% by photoevaporation.If only dust grains smaller than 0.3 µm are entrained, the sum of dust surface will also decrease by 10%.Then the photoelectric heating rate also decreases by 10% because it is proportional to the dust surface area.Since the disk lifetime is around a few million years, it is unlikely that the dust transportation by photoevaporative flows affect the disk dispersal significantly more than the above estimate.However, mixing of gas and dust near the base may be a complicated process, and thus ideally detailed hydrodynamics simulations with explicit treatment of dust size evolution and dynamical gas-dust (de)coupling would be needed to fully address this issue.

X-ray Luminosity Dependence
Recent observations suggest that EUV, FUV, X-ray luminosities differ considerably even among the same spectral-type stars (Gullbring et al. 1998;Flaccomio et al. 2003;France et al. 2012;Vidotto et al. 2014;France et al. 2018).Kunitomo et al. (2021) consider the variation of the emissivity of the high-energy radiation together with stellar evolution.For example, for a central star with mass M * ≥ 1.5 M ⊙ , the X-ray luminosity decreases by a factor of ∼ 10 4 at the age of 1-10 Myr.To examine the impact of X-ray radiation, we perform an additional set of photoevaporation simulations with varying X-ray luminosities to L X = 2.51 × 10 29 erg s −1 and L X = 2.51×10 31 erg s −1 , which are 0.1 and 10 times our fiducial value, L X,f , respectively.We keep the stellar mass to M * = 1 M ⊙ and the dust-to-gas mass ratio to D = 10 −2 , 10 −6 .
The resulting mass-loss rates are listed in Table 1.In the low X-ray cases, the gas is heated by both H 2 pumping and X-ray radiation.As the X-ray luminosity decreases, the gas temperature in the neutral region becomes lower by the factor of ∼ 3.This sensitivity indicates that the mass-loss rate is dependent on X-ray luminosity.

Other effects
In our simulations, we set the FUV luminosity without considering the detailed spectral shape of the central star.Realistic stellar spectrum should include emission lines such as Lyα and C IV. Schindhelm et al. (2012) have shown that ∼ 81% of the stellar radiation in the FUV range is contributed by hydrogen Lyα photons.In order to examine the effect of the "boost" of FUV radiation by the line emission, we perform a simulation with 10 times higher L FUV .The result shows that Ṁ is enhanced by a factor of 3 than the fiducial case, which clearly suggests that the resulting Ṁ can be higher if a realistic spectrum with Lyα radiation is considered.
Finally, we discuss the effect of H 2 pumping in an extremely metal/dust-poor disk.An interesting question is how the H 2 pumping affects the dispersal of disks that consists of only hydrogen and helium.Our simulation with D = 10 −6 can be regarded as almost a metal-free case in the early universe (Figure 2).The result shows that the disk around a low-mass star is heated by H 2 pumping and X-ray radiation and loses its gas mass by photoevaporation at a rate of ∼ 2.2×10 −10 M ⊙ yr −1 .In our future study (Komaki et al., in preparation), we examine the structure, stability and the dispersal of zerometallicity disks in the context of the formation of the first stars.

SUMMARY
Recent millimeter observations suggest that the dust amount in a PPD decreases with disk evolution.The gas mass variation still remains uncertain, but the dust abundance and its long-term variation may be important in the process of disk photoevaporation, and thereby planet formation through disk dispersal.In the present paper, we have performed a set of radiationhydrodynamics simulations with varying the dust-to-gas mass ratio to quantify the effect of H 2 pumping on the photoevaporation of dust-deficient disks.In the standard case with M * = 1 M ⊙ and D = 10 −2 , H 2 pumping heats the evaporating gas driven by FUV photoelectric heating.H 2 pumping does not directly contribute to the total mass loss but affects the temperature structure of the outflow.However, we have also shown that H 2 pumping heating is effective at the inner region within several gravitational radii where the density is higher than the critical density for collisional de-excitation of H 2 molecules.
We have found that the major heating process varies with the dust-to-gas mass ratio of the disk.In relatively dust-rich cases with D ≥ 10 −3 , the disks are heated mainly by the FUV photoelectric process, and the outer molecular region evaporates efficiently.In the other, dust-deficient disks with D ≤ 10 −3 , the dominant heating process is caused by H 2 pumping and X-ray radiation.The derived mass-loss rate as a function of D is given by Eq. 13 for the cases of D ≥ 10 −3 .The massloss rate is relatively small at D < 10 −3 where the FUV photoelectric heating is inefficient.
Since the dust amount can vary on a long timescale as the disk evolves, the temperature structure, the massloss profile, and the chemical composition would also change on similar timescales.In future work, we shall study the long-term dispersal process of a PPD by incorporating the formation, growth, and destruction of dust grains.

Figure 1 .
Figure 1.(top) The time-averaged disk structure of the run incorporating H2 pumping.The color map on the left portion shows the gas temperature, Tgas, and the right portion shows the density distribution, ρ.The white line shows the location of the base, which satisfies NH 2 = 10 20 cm −2 .The navy contour lines on the left side represent where the gas temperature is 100 K, 300 K, 1000 K, and 3000 K.The streamline represents the poloidal velocity of the gas, which is defined as vp = v 2 r + v 2 θ .For clarity, we omit to plot the streamlines where the velocity is less than 0.1km s −1 , which are mostly seen in the optically thick disk region.(bottom) The time-averaged disk chemical structure of the run incorporating H2 pumping.The color map on the left portion shows the specific heating rate for H2 pumping.The right portion shows the abundance of H2.The right portion also shows the distribution of H-bearing species.The yellow line on the right side represents the ionization front where the abundance of H II is 0.5, and the brown line indicates the dissociation front where the abundance of H2 is 0.25.

Figure 2 .
Figure2.The time-averaged specific heating and cooling rates measured at the base in our runs with different D. The solid lines show FUV photoelectric heating (orange), H2 photodissociation heating (yellow), H2 pumping heating (brown), and X-ray heating (magenta).The dotted lines are for dust-gas collisional cooling (purple), H2 cooling (cyan), and O I cooling (green).

Figure 3 .
Figure 3.The time-averaged mass-loss rate for each dustto-gas mass ratio.The blue dots represent the simulation results.The bars indicate the range of typical variation of the instantaneous rate without time-averaging.The dotted line is a functional fit we derive for D = 10 −3 − 10 −1 where the main disk heating process is photo-electric heating.

Table 1 .
The mass-loss rates in the runs with different D and LX.