Photoevaporating Protoplanetary Disks

# Radiation hydrodynamics simulations of photoevaporation of protoplanetary disks: Metallicity dependence

## Abstract

Protoplanetary disks are thought to have lifetimes of million years in the solar neighborhood, but recent observations suggest that the disk lifetimes are shorter in a low metallicity environment. We perform a suite of radiation hydrodynamics simulations of photoevaporation of protoplanetary disks to study the disk structure and its long-term evolution of years, and the metallicity dependence of mass-loss rate. Our simulations follow hydrodynamics, extreme and far ultraviolet radiative transfer, and non-equilibrium chemistry in a self-consistent manner. Dust grain temperatures are also calculated consistently by solving the radiative transfer of the stellar irradiation and grain (re-)emission. We vary the disk gas metallicity over a wide range of . For our fiducial model with a 0.5 central star with solar metallicity, the time-averaged photoevaporation rate is . The photoevaporation rate is lower with higher metallicity in the range of , because dust shielding effectively prevents far-ultra violet (FUV) photons from penetrating into and heating the dense regions of the disk. The photoevaporation rate sharply declines at even lower metallicities in , because FUV photoelectric heating is not efficient any more to raise the gas temperature and to drive outflows. At , H i photoionization heating acts as a dominant gas heating process and drives photoevaporative flows with roughly a constant rate. The typical disk lifetime is shorter at than at , being consistent with recent observations of the extreme outer galaxy. Finally, we develop a semi-analytic model that accurately describes the structure of photoevaporative flows and the metallicity dependence of mass-loss rate.

protoplanetary disks – stars: formation – infrared: planetary systems – stars: pre-main-sequence – ultraviolet: stars

## 1. Introduction

Protoplanetary disks are geometrically thin Keplerian disks surrounding pre-main-sequence stars e.g., (Shu et al., 1994). The central star grows through accretion of the protoplanetary disk and thus the disk evolution largely determines the mass of the star. Also, protoplanetary disks are the birth places of planets. Studying the structure and evolution of a protoplanetary disk is crucial in understanding planet formation.

Observationally, a star surrounded by a circumstellar disk shows a larger excess than a star without a circumstellar disk because of the dust infrared (IR) emission (Lada & Adams, 1992). This IR-excess is a robust indicator of the presence of a protoplanetary disk. Applying this diagnostics to the members of a cluster, one can estimate the disk fraction, which is the ratio of the member stars with disks to the total number of the members. It is observationally known that the disk fraction of the nearby clusters exponentially decreases with increasing cluster age, and it typically falls below for the cluster age of (e.g., Haisch et al., 2001). Hence the typical disk lifetime is estimated to be for the nearby clusters (Alexander et al., 2014; Gorti et al., 2016; Ercolano & Pascucci, 2017).

Interestingly, recent observations of the extreme outer Galaxy, where the metallicity is significantly lower than in the solar neighborhood, suggest that the typical disk lifetime is short (Yasui et al., 2009, 2010, 2016a, 2016b). The disk fraction there declines steeply with increasing cluster age and becomes within the cluster age of . It appears that a protoplanetary disk in low metallicity environments disperses earlier and/or faster than that of solar metallicity.

Protoplanetary disks lose their mass mainly via stellar accretion associated with angular momentum transfer, especially at the early stage of disk evolution (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974). The evolutional timescale is estimated to be of the order of a million years at several tens of , and can be even longer at further outside regions (Hollenbach et al., 2000; Armitage, 2011). Thus, viscous evolution alone cannot explain the observationally inferred disk lifetimes. Furthermore, viscous evolution predicts that the surface density should decrease with time as . It would contradict the existence of observed transitional disks (Andrews & Williams, 2005) and a much shorter transitional timescale than a lifetime (e.g., Skrutskie et al., 1990; Kenyon & Hartmann, 1995; Alexander et al., 2014). Clearly, in order to explain the timescale of disk dispersal and transition, there must be some other important physical mechanism(s).

Several dynamical processes such as photoevaporation (e.g., Hollenbach et al., 1994), MHD wind (Suzuki & Inutsuka, 2009), stellar wind (e.g., Elmegreen, 1979), and giant planet formation (Rice et al., 2003) have been proposed so far. In particular, photoevaporation is proposed as a main driver of disk dispersal. Photoevaporation appears to produce transitional disks when the effect is included in simulations of viscous disk evolution (Clarke et al., 2001; Alexander et al., 2006; Owen et al., 2010).

Photoevaporation from a disk is thought to occur in the following manner. The circumstellar disk is irradiated by the central star and/or by a nearby star. In optically thin regions, the gas temperature increases through thermalization of the electrons which are ejected from atoms and dust grains by absorption of high energy photons such as far ultraviolet (FUV; ), extreme ultraviolet (EUV; ), and X-rays (). The “hot” gas escapes from the star-disk system, and flows out of the disk. This causes considerable disk mass loss.

According to Hollenbach et al. (1994), who performed 1+1D radiative transfer calculations, the diffuse EUV component is dominant for exciting photoevaporation and drive mass loss at a rate of . In contrast, recent 2D radiative transfer calculations of Tanaka et al. (2013) suggest that the direct component of EUV is dominant. They derive a photoevaporation rate of . FUV can effectively heat denser regions of a disk than EUV, because FUV is attenuated at a higher column density () than EUV in general. FUV photoevaporation rates are thus generally higher than EUV photoevaporation rates (Gorti & Hollenbach, 2009; Owen et al., 2012). Gorti & Hollenbach (2009) conclude that FUV photoevaporation rates are of the order of for typical young low-mass stars. Photoevaporation excited by X-ray irradiation from a young low-mass star has also been studied (Alexander et al., 2004; Ercolano et al., 2008, 2009; Gorti & Hollenbach, 2008, 2009; Owen et al., 2010, 2012). X-rays are also attenuated at a larger column density comparable to FUV, and thus X-ray photoevaporation also gives a typical photoevaporation rate of (Ercolano et al., 2008, 2009; Owen et al., 2010, 2012).

Ercolano & Clarke (2010) (hereafter, EC10) derives the metallicity dependence of X-ray photoevaporation rates by hydrostatic calculations and an analytic formula to estimate a disk lifetime for a given photoevaporation rate. By applying the formula to the hydrostatic disks, EC10 derives the metallicity dependence of protoplanetary disk lifetimes. The obtained lifetimes monotonically decreases with metallicity, which appears to be consistent with the observed trend that disk lifetimes decrease with metallicity.

Unfortunately, none of these previous studies has derived the metallicity dependence of photoevaporation rates by using hydrodynamical simulations. Simultaneous modeling of photoevaporation and dynamical disk evolution is needed to study how photoevaporative flows change the density structure of a disk, and from where and with which speed photoevaporative flows are launched. Moreover, the previous studies which calculate the hydrodynamics of photoevaporating protoplanetary disks irradiated by a central low-mass star do not solve radiative transfer self-consistently. It is also important to incorporate non-equilibrium chemistry. It allows to study chemical evolution coupled with hydrodynamics and to estimate accurately the relevant heating/cooling rates. Yorke & Welz (1996) and Richling & Yorke (2000) self-consistently solve hydrodynamics and radiative transfer, including non-equilibrium chemistry, to follow the evolution of a disk irradiated by a central B star and an external radiation source, respectively. In the present paper, we solve the hydrodynamics of photoevaporating protoplanetary disks with self-consistent EUV/FUV radiative transfer and non-equilibrium chemistry. Our chemistry solver includes molecular species such as \ceH2 and \ceCO as well as atomic species. Dust temperatures are also calculated self-consistently. We run a set of simulations of a disk with different metallicities. We calculate the photoevaporation rates in order to derive, if any, the metallicity dependence of EUV/FUV photoevaporation rates and estimate disk lifetimes.

The paper is organized as follows. In Section 2, we present methods and the problem settings of our simulation. In Section 3, we discuss the simulation results and present an analytic model of the photoevaporation rates. A final discussion and a summary are given in Section 4 and Section 5, respectively.

## 2. Numerical Simulations

In order to calculate fluid dynamics of photoevaporating protoplanetary disks, we make use of a modified version of the publicly available code PLUTO (version 4.1; Mignone et al., 2007). We also summarize the following physical processes we implement in the code: radiative transfer, non-equilibrium chemistry network, and relevant heating/cooling processes.

### 2.1. Method

We consider the photoevaporation of protoplanetary disks caused by the UV irradiation from a central star, covering a broad range of different metallicities . We assume a central star with constant EUV photon number luminosity and FUV luminosity . Although the stellar UV emissivities will vary with different metallicities, we ignore such potential variation to concentrate on the roles of heavy elements contained within the disk. We keep the above parameters fixed throughout our simulations.

Our multi-species chemistry model is based on Omukai (2000) and Omukai et al. (2005, 2010). We assume that the medium consists of gas and dust grains. The gas contains seven chemical species: H i, H ii, \ceH2, CO, O i, C ii, and electron (hereafter, we refer to H i, H ii, and \ceH2 as H-bearing species and CO, O i, and C ii as metal species). We assume that the amount of the gas-phase metal elements and grains are proportional to relative metallicity . The dust-to-gas mass ratio and the gas-phase elemental abundances of carbon and oxygen are set to the values of local interstellar clouds in the case of . Hence, we give the dust to gas mass ratio by

 DG=0.01×Z/Z⊙. (1)

The gas-phase elemental abundances of carbon and oxygen are and , respectively (Pollack et al., 1994; Omukai, 2000). 7

### 2.2. Basic Equations

We use two dimensional spherical polar coordinates , taking into account the time evolution of gas density, all the three components of velocity , gas energy including relevant heating/cooling sources, and chemical abundances including advection and chemical reactions. The basic equations are

 ∂ρ∂t+∇⋅ρv=0 , (2) ∂ρvr∂t+∇⋅(ρvrv)=−∂P∂r−ρGM∗r2+ρv2θ+v2ϕr , (3) ∂ρvθ∂t+∇⋅(ρvθv)=−1r∂P∂θ−ρvθvrr+ρv2ϕrcotθ , (4) ∂ρvϕ∂t+∇l⋅(ρvϕv)=0 , (5) ∂E∂t+∇⋅(Hv)=−ρvrGM∗r2+ρ(Γ−Λ) , (6) ∂n\rm Hyi∂t+∇⋅(n\rm Hyiv)=n\rm HRi . (7)

In the above equations, , and are gas density, velocity, and pressure, respectively, and is the gravitational constant. We do not include the gas self-gravity which is currently negligible with the typical mass ratio between the star and disk, . We denote the total energy and enthalpy per unit volume of gas as and , respectively, and is heating rate per unit mass (specific heating rate), and is cooling rate per unit mass (specific cooling rate). We denote the fractional abundance of each of the seven chemical species as . Chemical reaction rates include all the relevant reactions (cf. Table B.3).

PLUTO discretizes the azimuthal component of Euler equations in an angular momentum conserving form. The divergence operator of Eq. (5) is represented by a different form compared to those of the other equations, and these divergence operators are defined as

 ∇⋅F= 1r2∂∂rr2Fr+1rsinθ∂∂θsinθFθ , (8) ∇l⋅F= 1r3∂∂rr3Fr+1rsin2θ∂∂θsin2θFθ , (9)

where is an arbitrary vector. Also, we do not consider angular momentum transfer due to viscous friction. We solve time evolution within the dynamical timescale of a disk which is much smaller than the viscous timescale.

We use the equation of state for an ideal gas:

 e =kTμmu(γ−1) , (10) P =ρkTμmu , (11)

where is specific energy of gas, is adiabatic index, is the Boltzmann constant, is gas temperature, is mean molecular weight, and is the atomic mass unit. The ratio of specific heat is defined as

 γ=1+y\rm HI+y\rm\text{H{\cal II}}+y\rm% \ce{H2}+y\rm e32y\rm HI+32y% \rm\text{H{\cal II}}+52y\rm\ce{H2}+32y% \rm e , (12)

where the contributions of the small abundances of the metal species are neglected. With the equation of state, the total energy and enthalpy per unit volume are explicitly written as

 E =12ρv2+ρe=12ρv2+Pγ−1 , (13) H =E+P=12ρv2+γPγ−1 . (14)

The computational domain is set to be on and . We assume axisymmetry around the rotational axis and mid-plane symmetry of a disk. We use grid cells logarithmically spaced in the radial direction. In the meridional direction, we use different resolutions in two domains divided by . In each domain, we use 80 uniform grid cells. The high resolution in allows to resolve the scale height of a disk and the launch points of photoevaporation flows.

### 2.3. Cooling/Heating

We implement photoionization heating caused by EUV and photoelectric heating caused by FUV (Bakes & Tielens, 1994). We also implement radiative recombination cooling of H ii (Spitzer, 1978), dust-gas collisional cooling (Yorke & Welz, 1996), Ly cooling of H i (Anninos et al., 1997), fine-structure line cooling of O i and C ii (Hollenbach & McKee, 1989; Osterbrock, 1989; Santoro & Shull, 2006), and molecular line cooling of \ceH2 and CO (Galli & Palla, 1998; Omukai et al., 2010).

We note that we do not include O i photoionization in our non-equilibrium chemistry. In an H ii region, O i is photoionized by EUV, and thus O i cooling is reduced. In order to include this effect approximately, we set O i abundance by when we calculate O i cooling rates. The heating/cooling rates are described in detail in Appendix A.

### 2.4. Chemical Reactions

We incorporate the relevant chemical reactions of the seven chemical species tabulated in Table B.3. As well as collisional chemical reactions, we implement the photo-chemical reactions: photoionization of H i, photodissociation of \ceH2 (Draine & Bertoldi, 1996), and photodissociation of CO (Lee et al., 1996).

We follow Richling & Yorke (2000) and assume that C i is converted to C ii soon after CO photodissociation. Namely, we assume that the CO dissociation front is located at the same position of the C ii ionization front. This procedure greatly saves computational cost. As the reverse reaction of CO photodissociation, we adopt the simplified chemistry model of Nelson & Langer (1997). This model can treat the formation of CO molecules from C ii via the reactions of hydrocarbon radicals without explicitly including C i as chemical species. The details of the chemical reactions are described in Appendix B.

We solve radiative transfer to calculate photo-chemical reaction rates, photo-heating rates, and dust temperatures consistently. Column densities are updated at each time-step. EUV radiative transfer is solved by ray-tracing. The diffusion component is neglected in our simulation, and we use case B recombination. Compared with the diffusion component, the direct component would have a dominant role in EUV photoevaporation (Tanaka et al., 2013), as discussed in Section 3.3. Although EUV photons are absorbed by H i and dust in general, we ignore the absorption by dust. Its effect on the absorption of EUV photons can be much less dominant than H i in the case of the EUV luminosity and computational domain we set above. 8

FUV radiative transfer is also solved by ray-tracing to calculate photoelectric heating rates, \ceH2 photodissociation rates, and CO photodissociation rates. We include the absorption of FUV photons by \ceH2 and CO molecules. The details of EUV/FUV radiative transfer are described in Appendix A and Appendix B.

We calculate grain temperatures by solving the transfer of both direct and diffusion components. We use a hybrid scheme; the direct component (stellar irradiation) is solved by ray-tracing, while the diffusion component due to thermal (re-)emission is solved by flux-limited-diffusion approximation. For these processes, we use the radiation transport module presented in Kuiper et al. (2010b). The hybrid scheme allows us to accurately model shadows caused by a optically thick disk (Kuiper & Klessen, 2013). This module has been applied in many astrophysical studies of massive star formation and feedback (Kuiper et al., 2010a, 2012; Kuiper & Yorke, 2013a; Kuiper et al., 2015, 2016), massive accretion disks (Kuiper et al., 2011; Meyer et al., 2017), stellar evolution (Kuiper & Yorke, 2013b), the formation of primordial stars (Hosokawa et al., 2016) as well as planet formation (Marleau et al., 2017). In our simulation, we use the opacity table taken from Draine & Lee (1984).

### 2.6. Initial Condition

The disk is assumed to consist of an initially neutral gas. All of hydrogen are assumed to be \ceH2 molecules and all of carbon are assumed to be CO at .

The initial grain and gas temperatures are set to be (e.g., Kenyon & Hartmann, 1987) except the case of . In the case of , we first calculate the thermo-chemical structure without updating the density structure for and start the simulation after that, otherwise the gas temperature does not couple with the dust temperature within the timescale of interest in the region near the mid-plane, which might be unrealistic.

The initial density structure is set to be hydrostatic equilibrium,

 n\rm H=n0(R1 AU)−9/4exp[−z22h2] , (15)

where and are positions in two dimensional cylindrical polar coordinates , is the scale height of a disk, which is defined as , where is isothermal sound speed and is the Keplerian angular velocity. We denote as the mid-plane density of a disk at , and we set . In Eq. (15), the surface density is assumed to have the profile of . The initial density distribution is shown in the top panel of Figure 1.

## 3. Results

Photoionization heating (hereafter, EUV heating) plays a dominant role in H ii  regions, while photoelectric heating (hereafter, FUV heating) does in neutral (H i , \ceH2) regions. As a result of the heating, photoevaporation is driven in those regions. In this section, we first discuss physical quantities such as density, velocity field, temperature, and chemical structure of a photoevaporating disk with solar metallicity and then we show their metallicity dependence (Section 3.1 and Section 3.2). Next, we study how the resulting photoevaporation rates vary with different metallicities (Section 3.3). Finally, we develop an semi-analytic model to interpret our numerical results (Section 3.4).

### 3.1. Physical Structure of a Solar Metallicity Disk

#### Density, Velocity, and Temperature Structures

As Figure 1 shows, photoevaporation occurs in both H ii region and neutral region. The neutral flow is formed by FUV heating. In order to check this, we perform a test simulation in which the FUV heating is initially included but artificially switched off at the time . After switching off the FUV heating, the neutral flow disappears soon. Thus, with our self-consistent radiation hydrodynamical simulation, we have confirmed that FUV is the main driver of the neutral photoevaporative flow (Gorti & Hollenbach, 2009; Owen et al., 2012).

FUV is strongly attenuated by dust once column density of hydrogen nuclei becomes in the case of , while EUV is strongly attenuated once H i column density becomes . Therefore, FUV photons typically reach and heat the denser regions of a disk than EUV photons. The typical density of the neutral flow () is much larger than the typical density of the H ii region flow () as visualized in Figure 1.

As shown in Figure 2, in the H ii region,

the main heating source is EUV heating, and the main cooling source is adiabatic cooling associated with gas expansion rather than radiative recombination cooling. The recombination timescale, , is typically longer than the crossing time of the ionized flow, . Hence, it typically flows out of the computational domain before recombining.

The heating/cooling sources lead to a gas temperature of in this region. The corresponding sound speed is . The gas is accelerated and expands according to the corresponding pressure gradient. The poloidal velocity achieves a few times of the sound speed () in the H ii region of Figure 1, as is also presented by the previous hydrodynamical simulation of EUV photoevaporation (Font et al., 2004).

In the neutral region, FUV heating balances O i cooling, \ceH2 cooling, and dust-gas collisional cooling. The most effective cooling source is O i line cooling between the H ii ionization front and the \ceH2 photodissociation front, while it is \ceH2 line cooling in the higher temperature regions of the \ceH2 region. Dust-gas collisional cooling becomes dominant among the three coolants in the much dense region. The similar features are observed in previous studies (e.g., Nomura & Millar, 2005; Nomura et al., 2007), except that there is no region where \ceH2 line cooling is the dominant cooling source. This discrepancy is attributed to the simple fact that \ceH2 cooling is not included as a cooling source in their models. Hence, we conclude that \ceH2 line cooling can be an effective cooling source as well as O i cooling and dust-gas collisional cooling in the neutral region of disks.

Adiabatic heating/cooling is subdominant in the region where FUV heating is dominant in contrast to the H ii region. The resulting temperature is (). The gas is accelerated by the pressure gradient and achieves in the neutral region while it expands.

#### Distribution of Hydrogen-Bearing Species

\ce

H2 photoevaporative flow is excited as Figure 1 shows. It is \ceH2 advection associated with photoevaporation that replenishes \ceH2 molecules into the neutral flow. It makes the height of the H/\ceH2 boundary large, as in Heinzeller et al. (2011). The study concludes that the H/\ceH2 boundary above a protoplanetary disk may be vertically transported to the upper region owing to the advection caused by winds, though hydrodynamics are not directly incorporated. Our hydrodynamical simulation confirms that the H/\ceH2 boundary is actually raised to the upper region by (FUV) photoevaporative advection from the dense region, where \ceH2 molecules are abundant.

In order to excite \ceH2 flow in the atmosphere, FUV photons should be sufficiently attenuated by dust shielding and/or \ceH2 self-shielding so that the \ceH2 photodissociation rate will be lower than the replenishing rate of \ceH2. The self-shielding becomes effective when \ceH2 column density is . In this study, we give the self-shielding function of \ceH2 as (Draine & Bertoldi, 1996, cf. Appendix B). As shown in Figure 3,

the photodissociation front coincides with the boundary where the \ceH2 self-shielding factor (the blue line in the bottom panel of Figure 3) sharply declines, i.e. self-shielding becomes strongly effective. Thus, \ceH2 molecules replenished by photoevaporation protect themselves against photodissciation by self-shielding rather than dust shielding.

It has been proposed, in the context of the protoplanetary disk chemistry, that self-shielding protects \ceH2 molecules against photodissociation especially in outer region of the disk (e.g., Woitke et al., 2009; Walsh et al., 2012). The height of H/\ceH2 boundary is much larger than those of the previous studies. For example, the height of H/\ceH2 boundary in our study is at (see Figure 1), while Woitke et al. (2009) shows that it is at . Thus, the hydrodynamics of the disk gas can significantly change the chemical structure of protoplanetary disks from the chemical structure expected by a hydrostatic calculation. This implies that the distributions of coolants are changed by the effect of hydrodynamics.

In the upper-height regions above the H/\ceH2 boundary in Figure 1, FUV photons are unshielded, and the \ceH2 abundance is determined by the balance between the unshielded photodissociation and the \ceH2 formation on grains. In the lower-height regions below the H/\ceH2 boundary, the \ceH2 advection is present and also replenishes \ceH2 molecules in addition to the \ceH2 formation on grains. Thus, in the lower-height region, the \ceH2 abundance is determined by the balance between the shielded photodissociation and the \ceH2 advection instead of the \ceH2 formation on grains. In H i region, (unshielded) photodissociation balances \ceH2 formation on grains instead of \ceH2 advection. The typical \ceH2 abundance is with the density of the region. No other chemical reaction is more efficient to form \ceH2 molecules than the formation on grains. This yields the abundance of \ceH2 molecules to stay low as long as the advection of \ceH2 molecules is absent. Thus, if there is little \ceH2 advection, H i region is formed and continues to exist.

#### Distribution of Metal Species

CO molecules are protected from photodissociation by self-shielding, \ceH2 shielding, and dust shielding of FUV photons (see Appendix B.2 for details).

Figure 4 shows that CO molecules are photodissociated until the dust shielding factor gets sufficiently small. This indicates that dust is the most dominant shielding source for FUV photons among the three kinds of the shielding sources. Therefore, the position of the CO photodissociation front is determined by dust shielding. This is, as shown in Figure 2, why the CO photodissociation front is almost identical to the boundary where FUV heating is effective in contrast to \ceH2 photodissociation front which is determined by the self-shielding of \ceH2.

C ii ionization front is assumed to be identical to CO photodissociation front which is caused by FUV in our model. Therefore, in Figure 2, the position of the C ii ionization front is not identical to that of the H ii ionization front but is embedded in the higher density region (i.e. the larger region) than the H ii ionization front.

The CO photodissociation front is almost identical to the boundary above which FUV heating is effective. Therefore, the difference in the height of the \ceC+/CO boundary between our study and previous hydrostatic studies is smaller than that of the H/\ceH2 boundary. For example, our study shows the height of the \ceC+/CO boundary is at , while Woitke et al. (2009) shows at . Hence, the chemical structure of CO molecules is not significantly affected by photoevaporation in contrast to that of \ceH2 molecules in our model.

### 3.2. Variations with Different Metallicities

#### Structure of Photoevaporationing Flow

Figure 5 presents the structure of the photoevaporative flow with different metallicities, , , and from the top to bottom panel. Although the photoevaporative flow is excited for all these cases, the dense neutral flow only appears with and . Remarkably, we also see that the neutral flow for is more denser than that for . In fact, the typical density of the neutral flow is for and for . We have confirmed that the density at the base of the neutral flow is almost proportional to in our simulations. The figure suggests that the metallicity of is required to excite the FUV-driven neutral photoevaporative flow, but that its density is higher with the lower metallicity once launched.

Since the visual extinction is proportional to the column density of grains along a line of sight , FUV photons can reach the denser part of the disk with the lower metallicity. This explains why the density of the neutral flow in the disk is much higher than that in the disk (Figure 5). We conclude that, for , the neutral flow has the higher density with the lower metallicity because the FUV radiation can reach and heat the dense part of the disk.

Next, we consider why the neutral photoevaporative flow turns to become weak for and almost ceases at . With our assumed dust-to-gas mass ratio in proportional to the metallicity, the relative amount of grains to the gas decreases with metallicity. Therefore, the specific FUV heating rate becomes small as metallicity decreases. In addition, under our chemistry model, the electron abundance is set to be equal to the abundance of the ionized carbon generated by CO photodissociation in the neutral region. The recombination timescale of charged grains becomes long at a fixed gas density as metallicity decreases, and dust grains are easy to be charged positively. Because of the deep coulomb potential of the positively charged grains, electrons become hard to be ejected from dust grains by the photoelectric effect. This yields low efficiency of the photoelectric effect (cf., Eq. A7) and reduces the resulting heating rate in the low density part of the neutral region (the region close to the \ceH+/H boundary). This effect influences the profile of photoevaporation especially with in our simulation.

Likewise, the specific cooling rates also become generally small with decreasing the metallicity. This behavior is clearly shown by the second row of Figure 6, which summarizes the specific heating and cooling rates within the , , and disks from the left to right column. Whereas the main cooling source is adiabatic cooling in the H ii region, O i cooling, \ceH2 cooling, and dust-gas collisional cooling dominate in the neutral region for . These cooling rates decrease with decreasing the metallicity, and become so small that the adiabatic cooling dominates in both the H ii and H i regions at the lowest metallicity .

The specific FUV heating rate, O i cooling rate, and dust-gas collisional cooling rate all decrease with metallicity owing to the decreasing amount of grains and metal species. However, the temperature of the neutral region also falls with metallicity as Figure 6 shows. This implies that FUV heating becomes less effective than cooling in the region as metallicity decreases.

In the low density part of the neutral region where O i cooling and \ceH2 cooling are dominant, the FUV heating rate is reduced by the low photoelectric efficiency in addition to the small amount of grains, as metallicity decreases. The O i cooling rate is reduced only by the small amount of O i and the \ceH2 cooling rate does not explicitly depend on metallicity. Therefore, compared with these cooling sources, FUV heating becomes relatively ineffective as metallicity decreases. In the high density part of the neutral region, the temperature is determined by the balance between the FUV heating and dust-gas collisional cooling. The photoelectric efficiency which depends on the electron density does not strongly depends on metallicity in the region, because the hydrogen nuclei density and the electron abundance in this region are basically proportional to and , respectively. Therefore, the specific FUV heating rate is basically proportional to metallicity (the amount of grains). Whereas, the specific dust-gas collisional cooling rate depends on dust temperature and is proportional to metallicity and hydrogen nuclei density. Dust temperature is determined by the balance between absorption and (re-)emission whose opacities are proportional to metallicity, and thus dust temperature does not strongly depend on metallicity. The density is proportional to in the region, so the specific dust-gas collisional cooling does not have explicit metallicity dependence in this region. As a result, similar to the low density part of the neutral region, FUV heating becomes relatively ineffective than dust-gas collisional cooling in the high density part of the neutral region. Hence, as metallicity decreases, FUV heating is reduced more strongly than O i cooling, \ceH2 cooling, dust-gas collisional cooling in the neutral region, so that the temperature of the neutral region falls with metallicity.

As metallicity decreases, FUV heating becomes unable to give neutral gas a sufficient energy to escape from the gravitational binding of the central star. With , the inefficient FUV heating leads to the gravitational dragging dominating over the pressure gradient, and the radial velocity occasionally has even a negative value. In the lowest metallicity range of , FUV heating does not even excite neutral photoevaporation.

#### Distribution of Hydrogen-bering Species

As discussed in Section 3.1, the chemical structures of \ceH2 and H i are determined by the balance of photodissociation, the photoevaporative advection of \ceH2, and \ceH2 formation on grains. As the solar metallicity disk, \ceH2 molecules are protected against photodissociation by self-shielding with any metallicity as shown by Figure 7.

The H/\ceH2 boundary is determined by the balance of \ceH2 advection and unshielded photodissociation, and it depends on the radius where sufficient \ceH2 flow occurs. The resulting gas temperature of the neutral region becomes high with high metallicity due to the efficient FUV heating. This allows \ceH2 molecules to evaporate even from the inner regions of a disk where the central star’s gravitational binding is strong. The \ceH2 flow density is small with high metallicity due to a large attenuation of dust (See also Section 3.4 for more quantitative discussions). Thus, with high metallicity, low-density \ceH2 flow is excited even from the inner region of a disk, and the small density H/\ceH2 boundary is formed.

The density of the \ceH+/H boundary is determined by the balance of photoionization and recombination of ionized hydrogen, and therefore it is independent of metallicity. The density of H/\ceH2 boundary is small with high metallicity. Hence, the density of the H/\ceH2 boundary becomes close to that of the \ceH+/H boundary with high metallicity, and this leads to a geometrically thin H i region with high metallicity as Figure 5 shows.

#### Distribution of Metal Species

The amount of dust is small at small metallicity, and the dust shielding factor becomes subdominant among the three shielding factors of CO photodissociation as metallicity decreases. As Figure 8 shows,

the most dominant shielding factor is the dust extinction factor with and , while it is \ceH2 shielding factor with . Thus, as metallicity becomes low, the most dominant attenuation source turns from dust to \ceH2, whose abundance does not depend on metallicity. In contrast, it is similar to the solar metallicity disk from Section 3.1 that the CO photodissociation front is embedded in the dense regions of the disks at all metallicities, as shown in Figure 6.

### 3.3. Metallicity Dependence of Photoevaporation Rate

We calculate photoevaporation rates by integrating the mass flux component normal to a spherical surface :

 ˙Mph=∫SdS⋅ρv=r2∫Sdθdϕsinθρvr , (16)

where is an infinitesimal surface element vector which is orthogonal to the spherical surface. We choose as the radius of (80% of the outer radial boundary) in order to avoid including unphysical effects from the outer boundary. This is the similar method for measuring a photoevaporation rate to Owen et al. (2010).

The blue dots in the top panel of Figure 9 show the resulting photoevaporation rates of the different metallicity disks estimated by Eq. (16). The photoevaporation rate is periodically variable in time as shown by the bottom panel of Figure 9. We give each of the dots as the time-averaged value of from to , which is almost seven times longer than the typical period.

The FUV photoevaporation rate is previously calculated by Gorti & Hollenbach (2009) (hereafter, GH09). GH09 estimates the photoevaporation rate of a solar metallicity disk without hydrodynamical simulation by giving the surface mass loss rate analytically. The study conclude that the resulting FUV photoevaporation rate is for and , while our result shows the photoevaporation rate of the disk is for and . The difference between these photoevaporation rates may mainly be attributed to the difference of the computational domains. We set the outer boundary of the computational domain to be and estimate the mass loss rate within , while GH09 sets the outer boundary to be and estimate the mass loss rate within the whole region. In fact, the photoevaporation rate of GH09 within is , which is well consistent with our result. Moreover, we run the simulation with the outer boundary of and got the photoevaporation rate of , which is also consistent with GH09.

In the metallicity range of , the photoevaporation rates are almost constant, as shown in the top panel of Figure 9. As discussed in Section 3.2, FUV heating is ineffective because of a smaller amount of grains and a smaller photoelectric effect efficiency. In this metallicity range, neutral flow is not excited and does not contribute to the photoevaporation rate, while ionized flow excited by EUV heating and contributes to the photoevaporation rates. Thus, the photoevaporation rates are almost constant in this metallicity range.

The time evolution of the photoevaporation rate is also almost independent of metallicity because of the same reason with . In this metallicity range, the time evolution is similar to the time evolution of indicated by the magenta line in the bottom panel of Figure 9. EUV photoevaporation is excited at any metallicity with no metallicity dependence, and its rate is smaller than that of FUV photoevaporation. Therefore, the magenta line virtually shows the minimum of photoevaporation rates.

The resulting EUV photoevaporation rate is as shown by Figure 9. In previous studies such as Hollenbach et al. (1994) and Font et al. (2004), the typical EUV photoevaporation rate is given by for and , which is smaller than our EUV photoevaporation rate. The photoevaporation rates of the previous studies are derived on the basis of the idea that the diffusion component of EUV dominates the direct component of EUV in a disk system. However, Tanaka et al. (2013) recently shows that the direct component is more dominant than the diffuse component by solving 2D radiative transfer, which is clearly more adequate than the approximated 1+1D radiative transfer. The estimated photoevaporation rate is typically five times larger than the photoevaporation rate estimated by the 1+1D radiative transfer. The EUV photoevaporation rates seems to be underestimated in the previous studies such as Hollenbach et al. (1994) and Font et al. (2004), and this is the reason why the photoevaporation rates are smaller than that of our study, where the diffusion component of EUV is not incorporated. We note that the geometrical structure of a disk is also crucial to determine which of the EUV components is dominant and affects a resulting photoevaporation rate. These differences in the previous studies and our model might also reflect the differences in the geometrical structures of the disks. Hence, it is important to solve radiative transfer with including self-consistent flow structure and scale height in order to estimate a photoevaporation rate.

In the range of , both neutral flow and ionized flow are constantly excited as shown in Figure 5. This is also directly shown by the bottom panel of Figure 9, where the photoevaporation rates are larger than that of at any time. As discussed in Section 3.2, FUV can reach and heat the dense region of the disk with low metallicity because of a smaller attenuation of FUV photons. The resulting density of the excited neutral photoevaporation flow increases as metallicity decreases. Thus, the photoevaporation rates increase as metallicity decreases because of a small dust attenuation in the metallicity range.

As presented in EC10, X-ray photoevaporation rate is also increased as metallicity decreases because of the smaller attenuation, though EC10 does not take into account FUV heating. The slope of the photoevaporation rates is in the range of in EC10 of our simulation is in the metallicity range of in our study, while that of our study is in the metallicity range of in our study. These slopes take similar values in spite of the fact that the photoevaporation models are completely different between EC10 and our study. This implies that photoevaporation rates would generally show a negative correlation to metallicity in this metallicity range.

In the range of , the specific FUV heating becomes inefficient owing to a small amount of grains and low photoelectric effect efficiency, as discussed in Section 3.2. This leads to a low temperature of the neutral gas. With such low temperature, the photoevaporation of the neutral gas is excited occasionally, as the bottom panel of Figure 9 shows. This is in contrast to the high metallicity range (), where a neutral flow is constantly excited. Moreover, the neutral gas with such low temperature can evaporate only from the outer region, where the central star’s gravity is weak. This means that the minimum radius where the occasional neutral flow is excited tends to become large as metallicity decreases. We note that the photoevaporative mass loss rate is counted only from the region inner than by Eq. (16). The neutral mass loss is neglected if the minimum radius for photoevaporation is large than . Therefore, the photoevaporation rates decline almost discontinuously, partly owing to the measurement technique of the photoevaporation rates.

The decline of the photoevaporation rates with decreasing metallicity in the range of is contrary to that of X-ray photoevaporation, which shows the monotonic increase of the photoevaporation rates as metallicity decreases in the range of as presented in EC10. The main differences between EC10 and our study is that EC10 derives the X-ray photoevaporation rates by solving the hydrostatic and thermo-chemical equilibrium structures iteratively, while our study derives the FUV/EUV photoevaporation rates by hydrodynamical simulations. The differences of the resulting metallicity dependences of the photoevaporation rates are possibly attributed to the differences of the incorporated heating sources and the calculation methods. We note that X-ray photoevaporation rates might also have the maximum at a certain metallicity as our result. In the limiting case , the main absorbers of X-rays, or the main ejectors of primary electrons, become H i, \ceH2, and \ceHe rather than metal elements (included in gas- and grain-phase). The cross sections of H i, \ceH2, and \ceHe are typically smaller than metal elements especially in the high energy range (; Maloney et al., 1996; Wilms et al., 2000). The specific X-ray heating rate should become small as metallicity decreases. Hence, with , X-rays has a weak heating effect and FUV has no effect on gas heating, so that a photoevaporation rate would converge to a small EUV/X-ray photoevaporation rate which is independent of metallicity. In any case, it is important to model X-ray photoevaporation hydrodynamically and to derive the metallicity dependence.

With a metallicity around the peak in the top panel of Figure 9, FUV photons excite a dense flow as metallicity decreases, but the specific FUV heating rate itself starts to become inefficient to excite neutral photoevaporative flow constantly. The peak is generated as a result of these two competing effects.

### 3.4. Semi-Analytic Model

In this section, we develop a semi-analytic model to interpret our numerical results, which show the metallicity dependence of the photoevaporation rates. As discussed in Sections 3.2 and 3.3, it is the FUV heating that causes the metallicity dependence of the photoevaporation rates. We then only focus on modeling the FUV-driven photoevaporation rate with different metallicities of . Regarding the photoevaporation via the EUV irradiation, we simply assume a constant rate without modeling it.

We suppose a situation shown in the schematic picture Figure 10. We further adopt following assumptions to construct our model:

1. The disk system is in a steady state.

2. Evaporative flows are launched from the regions where (hereafter, we simply refer to the wind lauching region as “base”).

3. Hydrogen atoms are all in the molecular form, but CO molecules are completely photodissociated at the base.

4. Evaporative flows are launched at the speed , where is the Mach number and is distance in cylindrical coordinates. 9

5. The azimuthal velocity is given by at the base.

6. The gas temperature at the base is determined by the thermal balance between the dominant heating and cooling processes, i.e., photoelectric heating, \ceH2 cooling, dust-gas collisional cooling, and O i cooling.

7. The evaporation flow is launched only from the base where the gas has the positive specific enthalpy (Liffman, 2003).

8. The profile of the base is approximated by a quadratic function

 z=f(R,Z)=a(Z)R2+b(Z)R, (17)

where the coefficients and are provided later.

We use as the exponent of the dust shielding factor for FUV heating. Therefore, the points where approximately correspond to the boundary where FUV can reach in the disk. The visual extinction is defined as

 Extra close brace or missing open brace (18)

where is hydrogen nuclei column density and is visual extinction per hydrogen nucleon. We expect, from Figure 1 and Figure 5, that the density does not significantly vary along the line of sight from the central star, so the integral part of Eq. (18) can be approximately rewritten to . Hence, in our semi-analytical model, we approximate the visual extinction as , and the base number density is given by

 n\rm H∼12Σd(Z/Z⊙)r . (19)

Hydrogen dominates the gas mass in our chemistry model, and thus the base density is approximately given by .

The FUV flux is analytically given at each point of the base. The density of each chemical species is derived by the third assumption and Eq. (19). The dust temperature is determined by the balance between the absorption of stellar irradiation and (re-)emission at the base, so that it is basically independent of metallicity. We use as the base dust temperature in our analytical model, which we have derived from our simulation results. Under the assumptions above, we can calculate the base gas temperature by solving a single non-linear equation of thermal equilibrium (the sixth assumption) at any metallicity. The resulting temperature is well described by a fit

 Tfit= T0(Z)(rr0)−α(Z) , (20) T0(Z)= 3.15×10 (log(Z/Z⊙))2 (21) + 3.64×102 log(Z/Z⊙) + 5.28×102 , α(Z)= 1.42×10−2(log(Z/Z⊙))3 (22) + 6.56×10−2(log(Z/Z⊙))2 + 2.41×10−2log(Z/Z⊙) + 3.05×10−1 ,

where . Note that takes a value in the range of with .

Under the above assumptions, the specific enthalpy at the base can be written as

 η =12v2p+γγ−1c2s−(GM∗r−12v2ϕ) , (24) =12M2c2s+γγ−1c2s−GM∗2r.

Thus, the condition, (the seventh assumption), corresponds to

 Extra close brace or missing open brace (25)

Therefore, the seventh assumption is equivalent to the assumption that photoevaporation is excited from the region where .

The quadratic coefficients in Eq. (17) are obtained by fitting our simulations as

 a= [−0.303(log(Z/Z⊙)+7.92×10−2)2+0.534] (26) × (100 AU)−1 , b= [1.34×10−2(log(Z/Z⊙))3 (27) + 3.26×10−2(log(Z/Z⊙))2 + 4.46×10−3log(Z/Z⊙)+0.421] .

By using all the elements above, we can finally derive the FUV photoevaporation rate analytically. Note that our model is based on one-dimensional distributions of the relevant physical quantities along the base. The FUV photoevaporation rate is given by

 Extra close brace or missing open brace (28)

where is a line element of the base and given by , and is the angle of the poloidal velocity relative to the line element . In our model, Eq. (28) is rewritten to

 Extra close brace or missing open brace 2∫η>0dR√1+f′2 2πR ρMcssinβ = 2πΣd(Z/Z⊙)√mHkT0rα0μ (29) Extra close brace or missing open brace

where is the upper limit of the integration, and it is set to be the real root of in order to compare the model with the simulation results in Section 3.3. The analytical rate is set to zero if , where is defined by the real root of . We set and . These values are determined from the simulation results in the regions where 10. We approximate the gradient of the base as . Then, Eq. (29) is rewritten as

 Extra close brace or missing open brace Extra close brace or missing open brace Extra close brace or missing open brace ≃ Extra close brace or missing open brace (30) Extra close brace or missing open brace ×√1+¯f′2[2+¯f′(¯f′+b)](1−α/2) ,

where .

The model photoevaporation rate is shown by the red line in the top panel of Figure 9. It is clear that our model explains well the metallicity dependence of the photoevaporation rate derived from our simulations.

Eq. (30) also gives the -dependence of the model photoevaporation rate by replacing with .

In order to compare the -dependence of the model photoevaporation rates with that of the simulation results. We use Eq. (16) with small modification:

 Extra close brace or missing open brace (31)

In the equation, is the regions where in the spherical surface at . We use the condition to calculate the contribution of the neutral photoevaporative flow to the photoevaporation rates. Figure 11 compares the -dependence of the analytic photoevaporation rate with that of the simulation results. The model photoevaporation rate of Eq. (30) can explain not only the metallicity dependence of the photoevaporation rates but also the -dependence of the photoevaporation rates.

In Eq. (30), we can give the approximate forms of and as functions of metallicity and we give and . Also, we can approximate the last factor of Eq. (30) to with an error of less than four percent. With these quantities, Eq. (30) can be further approximated to the form which explicitly depends on and :

 Extra close brace or missing open brace 9.0×10−9 M⊙/yr ~Z−1 (32) ×[3.64 log~Z+5.28]1/2 Extra close brace or missing open brace

where and . In Figure 11, the fitted model photoevaporation rate of Eq. (32) is shown be the green solid line and compared with the model photoevaporation rate of Eq. (30) denoted by the red solid line.

## 4. Discussion

### 4.1. Comparison with X-ray Photoevaporation

The resulting metallicity dependence of photoevaporation rates increases with decreasing metallicity in the range of , which is similar to that of EC10, while it decreases with metallicity in the range of , which is different from that of EC10. As discussed in Section 3.3, this difference could be originated from the difference of the methods for estimating the metallicity dependence of a photoevaporation rate: EC10 derives X-ray/EUV photoevaporation rates hydrostatically, while our study derives FUV/EUV photoevaporation rates hydrodynamically.

The main absorbers of energetic X-rays (; Wilms et al. (2000)) are metal elements, and thus the specific heating rate decreases. X-ray heating rate would be less effective than cooling and the resulting temperature would be low as metallicity decreases, as shown in this study. Therefore, as the metallicity dependence of FUV/EUV photoevaporation in this study, a X-ray/EUV photoevaporation rate might decrease with metallicity in a certain metallicity range. In any case, in order to derive the metallicity dependence of X-ray photoevaporation and compare it with that of FUV photoevaporation, the hydrodynamical simulations of X-ray photoevaporation is required.

Furthermore, hydrodynamical simulations itself is important for deriving a photoevaporation rate. It directly gives the base position, the launching velocity, and the critical radius above which heated gas can escape from the gravitational binding of the central star. These quantities have to be assumed in hydrostatic models. In fact, these two methods are known to give different photoevaporation rates because of the uncertainties originated from the assumptions of the quantities above (e.g., Ercolano et al., 2009; Owen et al., 2010), though the order of the photoevaporation rate is almost the same. In future work, we plan to incorporate X-ray radiative transfer into our photoevaporation model, and derive the metallicity dependence of photoevaporation excited by EUV/FUV/X-ray.

The crossing time of photoevaporative flow is much shorter than the timescale of the lifetime. This implies that it is computationally expensive to simulate the photoevaporation of a protoplanetary disk until the disk disperses completely. Instead of calculating the global evolution, we can use the analytic formula presented by EC10 to estimate a lifetime by giving a photoevaporation rate. When we set the exponent of the initial surface density profile to , the formula is given by . By assuming that the initial mass and radius of protoplanetary disks are independent of metallicity, we can evaluate the metallicity dependence of lifetimes from the metallicity-dependent photoevaporation rates of our simulations.

We have only two data points from observations in Figure 12. A reasonable and meaningful comparison with our model requires more observational data. However, at least, we can report here that the slope of the estimated lifetimes in the range of , , is quite consistent with that of the observational lifetimes, . Hence, it is suggested that FUV photoevaporation also has the potential to explain the short lifetimes of the protoplanetary disks in low metallicity environments as X-ray photoevaporation.

The metallicity of a protoplanetary disk should directly affect the formation efficiency of planets/cores. Actually, it has been observationally revealed that planet occurrence decreases with the host star metallicity in , which is so called “planet-metallicity correlation” (e.g., Gonzalez, 1997). This correlation is considered to reflect the fact that planet formation is inefficient in low metallicity environments. The metallicity dependence of the amount of material for forming planets is thought to be one of the causes for the correlation. Likewise, the disk lifetime is also an important quantity for planet formation efficiency. A protoplanetary disk lifetime limits the formation timescales of planets. The derived lifetime decreases with metallicity in the range of in this study, which is consistent with the correlation. Therefore, as well as the metallicity dependence of the amount of material to form planets/cores, the metallicity dependence of the lifetimes attributed to photoevaporation can play an important role in the occurrence of planets.

In reality, a disk lifetime depends not only on the photoevaporation rate but also depends on other parameters such as the disk radius, the disk mass, and the radius of a gap. It would be oversimplified to derive lifetimes from photoevaporation rates by the formula.

### 4.3. Effects of Grain Distribution and MHD Wind

Photoelectric heating generally depends on both the local dust-to-gas mass ratio and the local size distribution of dust/polycyclic aromatic hydrocarbon (PAH) grains. Though we assume a constant dust-to-gas mass ratio and a constant size distribution in the whole computational domain, they are, in general, variable because of settling, grain growth, and entrainment into disk wind (e.g., Hutchison et al., 2016). In fact, a variable dust-to-gas mass ratio and a variable grain size distribution are observationally proposed in both radial and vertical directions (e.g., Pinte et al., 2016). Therefore, for the metallicity dependence derived in this study, the photoevaporation rate is further affected by spatial grain distribution, grain size distribution, and grain aerodynamics.

Furthermore, magnetically driven disk wind (MHD wind) has been proposed to be capable of dispersing a protoplanetary disk (Suzuki & Inutsuka, 2009; Suzuki et al., 2010, 2016). MHD winds causes significant mass loss especially at the early stage of protoplanetary disk evolution, where the surface density is large. MHD winds are triggered by magnetorotational instability which depends on the local charged particle density in the high-density region of a disk. Its mass loss rate would also depend on metallicity. Hence, for a complete picture of disk dispersal, it would be required to simulate the global evolution of the disks subjected to both photoevaporation and MHD wind with various metallicities.

### 4.4. Outer Boundary Effect

In general, the profiles of photoevaporation and the derived photoevaporation rate can be affected by the bogus reflection at the boundary of the computational domain to some extent, especially when out-going flow is subsonic. In this study, the reflection possibly happens in the region close to both the launching points of photoevaporation flow and the outer boundary of the calculation region, where the flow is not yet accelerated up to . The bogus reflection leads to smoothing the pressure gradient owing to the accumulation of the gas near the outer boundary. In this case, gravitational force dominates over pressure gradient, so that the gas is artificially decelerated in the region.

The bogus reflection propagates with the sound speed. It can make the profiles and rate of photoevaporation have, if any, some features which vary with the timescale of the crossing time. In order to check if there is any effect caused by the bogus reflection, we carry out the simulations with a large outer boundary . The numbers of the cells are the same as the settings described in Section 2.1. Figure 13 shows the time-dependent photoevaporation rates of the and disks, which is estimated by Eq. (16).

The oscillation of the photoevaporation rate of is damped with time in the case of . This can be interpreted as, by expanding the computational domain, the bogus reflection is disappeared in the region close to both the launching point and , where there is large contribution to the photoevaporation rate due to the high density. In the case of , the oscillation is not disappeared. This shows that it is a physical effect that the radial component of the velocity, , is decelerated by the gravity of the star and occasionally has even a negative value, when metallicity is low. Further, this negative should physically take place from a qualitative point of view because FUV heating becomes inefficient to excite photoevaporation and cooling timescale is long as metallicity decreases. Although the profiles and rates of photoevaporation are actually affected by the unphysical effects near the the boundary of the computational domain, the time-averaged photoevaporation rates are not largely changed even if the computational domain is set to be large. Hence, our results shows the physical metallicity dependence of an FUV/EUV photoevaporation rate. Also, we suggest that in order to estimate a photoevaporation rate within a region whose size is characterized by a certain radius, we should set the computational domain to be sufficiently large such as double the radius.

## 5. Summaray

We have performed radiation-hydrodynamical simulations of photoevaporation of protoplanetary disks with self-consistent modelling of multi-species chemistry. In particular, we have considered a broad range of metallicities from to to examine the metallicity dependence, if any, of the disk lifetime. Our findings are summarized as follows:

• FUV photons reach denser regions of a disk than EUV photons, and thus are the main driver of photoevaporative flows.

• In low metallicity disks, dust shielding effect is reduced and FUV photons reach and heat even denser regions of the disk. The photoevaporation rate increases with decreasing metallicity in the range of .

• At , the photoevaporation rate descreases with decreasing metallicity because of inefficient photoelectric heating owing to the smaller amount of grains. A large portion of the disk gas is gravitationally bound.

• The photoevaporation rate shows a peak at as a result of the combination of the above two effects.

• In the metallicity range of , EUV photons primarily drive photoevaporative flows. Hence the resulting photoevaporation rate is nearly independent of metallicity in this extremely low-metallicity environment.

• We develop a semi-analytical model of disk photoevaporation that desribes accurately both the metallicity dependence of the photoevaporation rate (See Figure 9) and the outflow profile (See Figure 11).

Ercolano & Clarke (2010) argue that X-ray photoevaporation also causes metallicity dependence of the photoevaporation rates, and that the result is roughly consistent with that of the observational disk lifetimes. Unfortunately, their photoevaporation rates are derived by hydrostatic calculations, and thus are subject to several critical assumptions on the dynamical process. Based on the findings in this paper, we argue that it is necessary to examine the metallicity dependence of X-ray photoevaporation by using hydrodynamical simulations. We will address this issue in our forthcoming paper (R. Nakatani et al., in prep.).

Our analytic model in Section 4.2, suggests that the FUV-driven photoevaporation can explain the short lifetimes of the disks in low metallicity environments. A complete model of the protoplanetary disk dispersal would need to incorporate FUV/FUV/X-ray radiative transfer and possibly the effect of magnetic fields. We aim to extend our work to simulate the long-term evolution of protoplanetary disks to derive their lifetimes and the metallicity dependence.

We thank David Hollenbach, Shu-ichiro Inutsuka, Takeru Suzuki, and Kei Tanaka for helpful discussions and insightful comments on the paper. RN has been supported by the Grant-in-aid for the Japan Society for the Promotion of Science (16J03534) and by Advanced Leading Graduate Course for Photon Science (ALPS) of the University of Tokyo. TH appreciates the financial supports by the Grants-in-Aid for Basic Research by the Ministry of Education, Science and Culture of Japan (16H05996). HN appreciates the financial supports by Grants-in-Aid for Scientific Research (25400229). RK acknowledges financial support via the Emmy Noether Research Group on Accretion Flows and Feedback in Realistic Models of Massive Star Formation funded by the German Research Foundation (DFG) under grant no. KU 2849/3-1. All the numerical computations were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

## Appendix A Cooling/Heating

In this section, we summarize the heating/cooling processes included in our simulations.

### a.1. Photo-heating

We implement the photo-heating processes by stellar EUV/FUV irradiation. We directly solve radiative transfer to calculate the photoionization heating (EUV heating) rate, while we simply use an analytic formula presented in Bakes & Tielens (1994) (hereafter, BT94) to obtain the photoelectric heating (FUV heating) rate.

We consider absorption of the direct EUV photons from the central star. We solve

 1r2∂∂r(r2Fν)=−n% H{\cal I}σνFν , (A1)

where is a frequency of EUV photons, is the specific number flux of the direct EUV field, is number density of H i, and is the absorption cross section of H i. We use the approximate absorption cross section

 σν=6.3×10−18(hνhν1)−3  cm2 , (A2)

(e.g., Osterbrock & Ferland, 2006). In Eq. (A2), is the frequency at the Lyman limit (