Radiation Hydrodynamics Simulations of Photoevaporation of Protoplanetary Disks by Ultra Violet Radiation: Metallicity Dependence

TL;DR

This study uses radiation hydrodynamics simulations to explore how metallicity influences the photoevaporation and lifetime of protoplanetary disks, revealing a complex dependence with lower metallicity leading to faster disk dispersal.

Contribution

It presents the first comprehensive simulations incorporating self-consistent radiative transfer, chemistry, and dust physics across a wide metallicity range to analyze disk photoevaporation.

Findings

Photoevaporation rate decreases with increasing metallicity in certain ranges due to dust shielding.

At very low metallicities, photoevaporation is driven mainly by HI photoionization heating.

Disk lifetime is shorter at sub-solar metallicities, consistent with observations.

Abstract

Protoplanetary disks are thought to have lifetimes of several 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 $\sim 10000$ years, and the metallicity dependence of mass-loss rate. Our simulations follow hydrodynamics, extreme and far ultra-violet 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 $10^{-4}~ Z_\odot \leq Z \leq 10 ~Z_\odot$. The photoevaporation rate is lower with higher metallicity in the range of $10^{-1} \,Z_\odot \lesssim Z \lesssim 10 \,Z_\odot$, 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 $10^{-2} \,Z_\odot \lesssim Z \lesssim 10^{-1}\,Z_\odot$, because FUV photoelectric heating becomes less effective than dust-gas collisional cooling. The temperature in the neutral region decreases, and photoevaporative flows are excited only in an outer region of the disk. At $10^{-4}\,Z_\odot \leq Z \lesssim 10^{-2}\,Z_\odot$, HI 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 $Z=0.3~Z_\odot$ than at $Z = Z_\odot$, being consistent with recent observations of the extreme outer galaxy.

Tables

Table 1. Table 1 Properties of the model

Stellar parameters
Stellar mass	$0.5M_{\odot }$
Stellar radius	$2R_{\odot }$
FUV luminosity	$3\times {10}^{32}\mathrm{erg}{s}^{-1}$
EUV luminosity	$6\times {10}^{41}{s}^{-1}$
Gas/dust properties
Species	H I, H II, \ceH2, CO, O I, C II, \cee-
Carbon abundance	$0.927\times {10}^{-4}\times Z/Z_{\odot }$
Oxygen abundance	$3.568\times {10}^{-4}\times Z/Z_{\odot }$
Dust to gas mass ratio	$0.01\times Z/Z_{\odot }$

Species	$j\to i$	${\nu }_{ij}$[Hz]	$A_{ij}$ [/s]	${\gamma }_{ij}^e$ [${\mathrm{cm}}^3/s$]	${\gamma }_{ij}^{\mathrm{HI}}$ [${\mathrm{cm}}^3/s$]	reference
C II	2$\to$ 1	$1.9\times {10}^{12}$	$2.4\times {10}^{-6}$	$2.8\times {10}^{-7}{\left(T/100K\right)}^{-0.5}$	$8.0\times {10}^{-10}{\left(T/100K\right)}^{0.07}$	1,2
O I	2$\to$ 1	$4.7\times {10}^{12}$	$8.9\times {10}^{-5}$	$1.4\times {10}^{-8}$	$9.2\times {10}^{-11}{\left(T/100K\right)}^{0.67}$	1,3
O I	3$\to$ 1	–	$1.0\times {10}^{-10}$	$1.4\times {10}^{-8}$	$4.3\times {10}^{-11}{\left(T/100K\right)}^{0.80}$	1,3
O I	4$\to$ 1	–	$6.3\times {10}^{-3}$	$1.0\times {10}^{-10}$	$1.0\times {10}^{-12}$	1,3
O I	5$\to$ 1	–	$2.9\times {10}^{-4}$	$1.0\times {10}^{-10}$	$1.0\times {10}^{-12}$	1,3
O I	3$\to$ 2	$2.1\times {10}^{12}$	$1.7\times {10}^{-5}$	$5.0\times {10}^{-9}$	$1.1\times {10}^{-10}{\left(T/100K\right)}^{0.44}$	1,3
O I	4$\to$ 2	–	$2.1\times {10}^{-3}$	$1.0\times {10}^{-10}$	$1.0\times {10}^{-12}$	1,3
O I	5$\to$ 2	–	$7.3\times {10}^{-2}$	$1.0\times {10}^{-10}$	$1.0\times {10}^{-12}$	1,3
O I	4$\to$ 3	$4.7\times {10}^{14}$	$7.3\times {10}^{-7}$	$1.0\times {10}^{-10}$	$1.0\times {10}^{-12}$	1,3
O I	5$\to$ 3	–	$0$	$1.0\times {10}^{-10}$	$1.0\times {10}^{-12}$	1,3
O I	5$\to$ 4	$5.4\times {10}^{14}$	$1.2$	$0$	$0$	1,3

Label	Reaction	Rate Coefficient  a	Reference b
k1	\ceH + e -¿ H+ + 2e	$\exp [-32.71396786$	1
		$+13.536556\ln T_{\mathrm{eV}}$
		$-5.73932875{\left(\ln T_{\mathrm{eV}}\right)}^2$
		$+1.56315498{\left(\ln T_{\mathrm{eV}}\right)}^3$
		$-0.2877056{\left(\ln T_{\mathrm{eV}}\right)}^4$
		$+3.48255977\times {10}^{-2}{\left(\ln T_{\mathrm{eV}}\right)}^5$
		$-2.63197617\times {10}^{-3}{\left(\ln T_{\mathrm{eV}}\right)}^6$
		$+1.11954395\times {10}^{-4}{\left(\ln T_{\mathrm{eV}}\right)}^7$
		$-2.03914985\times {10}^{-6}{\left(\ln T_{\mathrm{eV}}\right)}^8]$
k2	\ceH+ + e -¿ H + $\gamma$	$\exp [-28.6130338$	1
		$-0.72411256\ln T_{\mathrm{eV}}$
		$-2.02604473\times {10}^{-2}{\left(\ln T_{\mathrm{eV}}\right)}^2$
		$-2.38086188\times {10}^{-3}{\left(\ln T_{\mathrm{eV}}\right)}^3$
		$-3.21260521\times {10}^{-4}{\left(\ln T_{\mathrm{eV}}\right)}^4$
		$-1.42150291\times {10}^{-5}{\left(\ln T_{\mathrm{eV}}\right)}^5$
		$+4.98910892\times {10}^{-6}{\left(\ln T_{\mathrm{eV}}\right)}^6$
		$+5.75561414\times {10}^{-7}{\left(\ln T_{\mathrm{eV}}\right)}^7$
		$-1.85676704\times {10}^{-8}{\left(\ln T_{\mathrm{eV}}\right)}^8$
		$-3.07113524\times {10}^{-9}{\left(\ln T_{\mathrm{eV}}\right)}^9]$
k12	\ceH2 + e -¿ 2H + e	$4.4\times {10}^{-10}{T}^{0.35}\exp \left(-1.02\times {10}^5/T\right)$	1
k13	\ceH2 + H -¿ 3H	$k_H{\left(k_L/k_H\right)}^a$,	1
		$k_L\equiv 1.12\times {10}^{-10}\exp \left(-7.035\times {10}^4/T\right)$,
		$k_H\equiv 6.5\times {10}^{-7}/\sqrt{T}\exp \left(-5.2\times {10}^4/T\right)\left(1-\exp \left(-6000/T\right)\right)$
		$\log n_{\mathrm{cr}}\equiv 4-0.416\log \left(T/1.0\times {10}^4\right)-0.327{\left(\log \left(T/1.0\times {10}^4\right)\right)}^2$,
		$a\equiv {\left(1+n_{\text{H}}/n_{\mathrm{cr}}\right)}^{-1}$
k19	\ce3H -¿ H2 + H	$5.5\times {10}^{-29}/T$	1
k20	\ce2H + H2 -¿ 2H2	$R_{\mathrm{k19}}/8$	1
k21	\ce2H2 -¿ 2H + H2	$k_H{\left(k_L/k_H\right)}^a$,	1
		$k_L\equiv 1.18\times {10}^{-10}\exp \left(-6.95\times {10}^4/T\right)$,
		$k_H\equiv 8.125\times {10}^{-8}/\sqrt{T}\exp \left(-5.2\times {10}^4/T\right)\left(1-\exp \left(-6000/T\right)\right)$,
		$\log n_{\mathrm{cr}}\equiv 4.845-1.3\log \left(T/1.0\times {10}^4\right)+1.62{\left(\log \left(T/1.0\times {10}^4\right)\right)}^2$,
		$a\equiv {\left(1+n_{\text{H}}/n_{\mathrm{cr}}\right)}^{-1}$
k22	\ce2H -¿ H+ + e + H	$1.7\times {10}^{-4}{R}_{\mathrm{k1}}$	1
k23	\ce2H -¿C[dust] H2	$6.0\times {10}^{-17}\sqrt{T/300}{f}_a\left(Z/Z_{\odot }\right)$	1
		$\times {\left[1.0+4.0\times {10}^{-2}\sqrt{T+T_{\mathrm{dust}}}+2.0\times {10}^{-3}T+8.0\times {10}^{-6}{T}^2\right]}^{-1}$,
		$f_a\equiv {\left[1.0+\exp \left(7.5\times {10}^2\left(1/75-T_{\mathrm{dust}}^{-1}\right)\right)\right]}^{-1}$
p1	\ceH + $\gamma$ \ce-¿ H+ + e	$R_{\mathrm{Ionize}}$  (cf. Eq. (A5))	-
p2	\ceH2 + $\gamma$ \ce-¿ 2 H	$R_{\text{ce}H2,\mathrm{diss}}$  (cf. Appendix B.1)	2
p3	\ceCO + $\gamma$ \ce-¿ C+ + O	$R_{\text{ce}CO,\mathrm{diss}}$  (cf. Appendix B.2)	3
k24	\ceC+ + O -¿ CO	$R_{\mathrm{CO},\mathrm{form}}$  (cf. Appendix B.3)	4

Equations

$$\mathcal{DG}=0.01\times Z/Z_{\odot }.$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho \mathbf{v}=0,$$

$$\frac{\partial \rho v_r}{\partial t}+\nabla \cdot \left(\rho v_r\mathbf{v}\right)=-\frac{\partial P}{\partial r}-\rho \frac{G{M}_{\ast }}{r^2}+\rho \frac{v_{\theta }^2+v_{\varphi }^2}{r},$$

$$\frac{\partial \rho v_{\theta }}{\partial t}+\nabla \cdot \left(\rho v_{\theta }\mathbf{v}\right)=-\frac{1}{r}\frac{\partial P}{\partial \theta }-\rho \frac{v_{\theta }{v}_r}{r}+\frac{\rho v_{\varphi }^2}{r}\cot \theta ,$$

$$\frac{\partial \rho v_{\varphi }}{\partial t}+{\nabla }^l\cdot \left(\rho v_{\varphi }\mathbf{v}\right)=0,$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left(H\mathbf{v}\right)=-\rho v_r\frac{G{M}_{\ast }}{r^2}+\rho \left(\Gamma -\Lambda \right),$$

$$\frac{\partial n_{\text{H}}{y}_i}{\partial t}+\nabla \cdot \left(n_{\text{H}}{y}_i\mathbf{v}\right)=n_{\text{H}}{R}_i.$$

$$\nabla \cdot \mathbf{F}=$$

$${\nabla }^l\cdot \mathbf{F}=$$

$$e$$

$$P$$

$$\gamma =1+\frac{y_{\text{HI}}+y_{\text{HII}}+y_{\text{ceH2}}+y_{\text{e}}}{\frac{3}{2}{y}_{\text{HI}}+\frac{3}{2}{y}_{\text{HII}}+\frac{5}{2}{y}_{\text{ceH2}}+\frac{3}{2}{y}_{\text{e}}},$$

$$E$$

$$H$$

$$n_{\text{H}}=n_0{\left(\frac{R}{1\mathrm{AU}}\right)}^{-9/4}\exp \left[-\frac{z^2}{2h^2}\right],$$

$${\dot{M}}_{\mathrm{ph}}={\int }_{S}d\mathbf{S}\cdot \rho \mathbf{v}=r_S^2{\int }_{S}d\theta d\varphi \sin \theta \rho v_r,$$

$$\eta =\frac{1}{2}{\mathbf{v}}^2+\frac{\gamma }{\gamma -1}{c}_s^2-\frac{G{M}_{\ast }}{r},$$

$$z=f(R,Z)=a(Z)R^2+b(Z)R,$$

$$A_{\text{V}}={\Sigma }_{\mathrm{d}}{N}_{\text{H}}={\Sigma }_{\mathrm{d}}{\int }_{R_{\ast }}^{r}d{r}^{\prime }{n}_{\text{H}},$$

$$n_{\text{H}}\sim \frac{1}{2{\Sigma }_{\mathrm{d}}(Z/Z_{\odot })r}.$$

$$T_{\mathrm{fit}}=$$

$$T_0(Z)=$$

$$\alpha (Z)=$$

$$\eta$$

$$r>r_{\text{min}}\equiv r_0{\left[\frac{\gamma -1}{(1+2{\mathcal{M}}^{-2})\gamma -1}\frac{\mu m_{u}G{M}_{\ast }}{{\mathcal{M}}^2{r}_{0}k{T}_0}\right]}^{1/(1-\alpha )}.$$

$$a=$$

$$b=$$

$${\dot{M}}_{\text{FUV}}^{\text{ana}}={\int }_{\eta >0}ds2\pi R\rho v_{\mathrm{p}}\sin \beta ,$$

$${\dot{M}}_{\text{FUV}}^{\text{ana}}=$$

$$=$$

$${\dot{M}}_{\text{FUV}}^{\text{ana}}\simeq$$

$$\simeq$$

$${\dot{M}}_{\text{FUV}}^{\text{sim}}(r)=r^2{\int }_{S_{\mathrm{I}}(r)}d\theta d\varphi \sin \theta \rho v_r,$$

$${\dot{M}}_{\text{FUV}}^{\text{ana,fit}}\simeq$$

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{F}_{\nu }\right)=-n_{\text{HI}}{\sigma }_{\nu }{F}_{\nu },$$

$${\sigma }_{\nu }=6.3\times 10^{-18}{\left(\frac{h\nu }{h{\nu }_1}\right)}^{-3}\mathrm{c}{\mathrm{m}}^2,$$

$$F_{\nu }(r,\theta ,t)=\frac{{\Phi }_{\nu }(R_{\ast })}{4\pi r^2}\exp \left[-{\sigma }_{\nu }{N}_{\text{HI}}\right],$$

Full text

1. INTRODUCTION

Protoplanetary disks are geometrically thin Keplerian disks surrounding pre-main-sequence stars e.g., (Shu et al., 1994). They are considered to be the birth places of planets, and thus 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 $H-K$ 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 $10\%$ for the cluster age of $\gtrsim 6\,{\rm Myr}$ (Haisch et al., 2001; Hernández et al., 2007; Meyer et al., 2007; Mamajek, 2009; Fedele et al., 2010; Ribas et al., 2014). Hence the typical disk lifetime is estimated to be $\sim 3-6\,{\rm Myr}$ 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 $\lesssim 10\%$ within the cluster age of $\lesssim 1\,{\rm Myr}$. 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 ${\rm AU}$, 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 $\Sigma\propto t^{-p}$. 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 (e.g., 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; $6{\rm\,eV}<h\nu<13.6{\rm\,eV}$), extreme ultraviolet (EUV; $13.6{\rm\,eV}<h\nu<0.1\,{\rm keV}$), and X-rays ($h\nu>0.1\,{\rm keV}$). 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 $\dot{M}_{\text{ph}}\sim 10^{-10}\,M_{\odot}\,{\rm yr}^{-1}$. In contrast, recent 2D radiative transfer calculations of Tanaka et al. (2013) suggest that the direct component of EUV is dominant. They derive $\dot{M}_{\text{ph}}\sim 10^{-9}\,M_{\odot}\,{\rm yr}^{-1}$. FUV can effectively heat denser regions of a disk than EUV, because FUV is attenuated at a higher column density ($\sim 10^{21}\,{\rm cm^{-2}}$) 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 $\sim 10^{-8}\,M_{\odot}\,{\rm yr}^{-1}$ 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 $\dot{M}_{\text{ph}}\sim 10^{-8}\,M_{\odot}\,{\rm yr}^{-1}$ (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 $\dot{M}_{\text{ph}}$. 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 necessary 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, a 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 $10^{-4}~{}Z_{\odot}\leq Z\leq 10~{}Z_{\odot}$. We assume a central $M_{*}=0.5~{}M_{\odot}$ star with constant EUV photon number luminosity $\Phi_{\rm EUV}=6\times 10^{41}~{}{\rm s}^{-1}$ and FUV luminosity $L_{\rm FUV}=3\times 10^{32}~{}{\rm erg}~{}{\rm s}^{-1}$. 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. X-rays can also drive photoevaporation (Alexander et al., 2004; Ercolano et al., 2008, 2009; Gorti & Hollenbach, 2008, 2009; Owen et al., 2010, 2012), but we do not include X-rays in the present study. Here, we focus on the metallicity dependence of UV-driven photoevaporation.

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 $Z/Z_{\odot}$. 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 $Z=Z_{\odot}$. Hence, we give the dust to gas mass ratio $\mathscr{DG}$ by

[TABLE]

The gas-phase elemental abundances of carbon and oxygen are $y_{\text{\rm C}}=0.927\times 10^{-4}~{}Z/Z_{\odot}$ and $y_{\text{\rm O}}=3.568\times 10^{-4}~{}Z/Z_{\odot}$, respectively (Pollack et al., 1994; Omukai, 2000). 111 The abundance of species $i$ is defined as the ratio of its number density to hydrogen nuclei number density: $y_{i}\equiv n_{i}/n_{\text{\rm H}}~{}.$ We adopt chemical symbol notation for elemental abundances and Romanian notation for chemical abundances. For example, $y_{\text{\rm C}}$ and $y_{\text{\rm CI}}$ denote the elemental abundance of carbon and the chemical abundance of neutral carbon atoms, respectively. This is also the case with density and column density.

The parameters used in our model are listed in Table 1.

2.2. Basic Equations

We use two dimensional spherical polar coordinates $(r,~{}\theta)$, taking into account the time evolution of gas density, all the three components of velocity $\bm{v}=(v_{r},~{}v_{\theta},~{}v_{\phi})$, gas energy including relevant heating/cooling sources, and chemical abundances including advection and chemical reactions. The basic equations are

[TABLE]

In the above equations, $\rho,~{}\bm{v}$, and $P$ are gas density, velocity, and pressure, respectively, and $G$ 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, $M_{\rm disk}/M_{*}\sim 0.01$. We denote the total energy and enthalpy per unit volume of gas as $E$ and $H$, respectively, and $\Gamma$ is a heating rate per unit mass (specific heating rate), and $\Lambda$ is a cooling rate per unit mass (specific cooling rate). We denote the fractional abundance of each of the seven chemical species as $y_{\text{\rm HI}},~{}y_{\text{\rm\text{H{\cal II}}}},~{}y_{\text{\rm\ce{H2}}},~{}y_{\text{\rm CO}},~{}y_{\text{\rm OI}},~{}y_{\text{\rm\text{C{\cal II}}}},~{}y_{\text{\rm e}}$. Chemical reaction rates $R_{i}$ 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

[TABLE]

where $\bm{F}$ 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:

[TABLE]

where $e$ is specific energy of gas, $\gamma$ is adiabatic index, $k$ is the Boltzmann constant, $T$ is gas temperature, $\mu$ is mean molecular weight, and $m_{u}$ is the atomic mass unit. The ratio of specific heat $\gamma$ is defined as

[TABLE]

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

[TABLE]

The computational domain is set to be on $r=[1,\,400]\,{\rm AU}$ and $\theta=[0,~{}\pi/2]{\rm~{}rad}$. We need to use a sufficiently large radial outer boundary so that the computational domain contains transonic points of photoevaporative flows at the metallicities of interest (see the discussion in Section 4.5). The sink region $(\leq 1\,{\rm AU})$ is out of the computational domain, but disk materials exist there in reality and shield stellar photons. We take into account this effect approximately by assuming the sink density distributions are radially uniform and the densities are given by those of the innermost cells in the computational domain. Thus, the sink column densities are calculated as $N_{i}^{\text{sink}}=sn_{\text{i, m}}$, where $i$ is a label of the chemical species, $s$ is the sink size ($1\,{\rm AU}$), and $n_{\text{i, m}}$ is the density of the chemical species in the innermost cell. The stellar photon fluxes are reduced by the sink column densities in our simulations. We assume axisymmetry around the rotational axis $(\theta=0)$ and mid-plane symmetry $(\theta=\pi/2)$ of a disk. We use $128$ grid cells logarithmically spaced in the radial direction. In the meridional direction, we use different resolutions in two domains divided by $\theta=1$. In each domain, we use 80 uniform grid cells. The high resolution in $1\leq\theta\leq\pi/2$ allows to resolve the scale height of a disk and the launch points of photoevaporation flows, which are called photoevaporation bases.

The effective gravitational radius for an ionized gas $(T=10^{4}{\rm\,K})$ is $\simeq 1.4(M_{*}/M_{\odot})\,{\rm AU}$ (Liffman, 2003). Our inner extent of the computational domain is larger than the effective gravitational radius for a $0.5~{}M_{\odot}$ star. Therefore, our calculations might miss the contribution of mass-loss from the region near the effective radius. However, the resulting base density profile of an ionized gas is expected to show, and actually has, a scaling of $\propto R^{-1.5}$, where $R$ is the cylindrical radius (Tanaka et al., 2013). In this case, the mass-loss is dominated by the contribution from outer regions of a disk, and the contribution from the region near the effective gravitational radius is sufficiently small. We have run simulations with small inner boundaries of $r_{\text{inner}}=0.1\,{\rm AU},~{}0.35\,{\rm AU},~{}0.5\,{\rm AU}$, to confirm that the resulting photoevaporation rate is almost the same as that of a simulation with $r_{\text{inner}}=1\,{\rm AU}$. The contributions from $R\leq 10\,{\rm AU}$ is only about a few percent of the total. Hence, we use $r_{\text{inner}}=1\,{\rm AU}$ for the inner boundary of our computational domain.

We note that the absorption of direct EUV photons by the inner $(<1\,{\rm AU})$ disk could be important. If all the direct stellar photons are absorbed by the inner disk and its atmosphere, only diffuse photons emitted through recombination can reach the outer region. In the simulations with small inner boundaries of $r_{\text{inner}}=0.1\,{\rm AU},~{}0.35\,{\rm AU},~{}0.5\,{\rm AU}$, we find that the density of the ionized atmosphere is sufficiently small not to shield the EUV photons, and that the direct photons actually reach $r>1$ AU. In the outer region, the heating rate and ionization rate are almost the same as in the simulation with a boundary of $r_{\text{inner}}=1\,{\rm AU}$. We thus obtain essentially the same photoevaporation rate from the two simulations. Therefore, we justify using the computational domain of $r=[1,~{}100]\,{\rm AU}$.

2.3. Cooling/Heating

We implement photoionization heating caused by EUV and photoelectric heating caused by FUV. We use the analytic formula presented by Bakes & Tielens (1994) to calculate photoelectric heating. Bakes & Tielens (1994) assume the MRN distribution (Mathis et al., 1977) for the dust model to derive the formula. The same size distribution is assumed for small carbon grains, polycyclic aromatic hydrocarbons (PAHs). Note that the observed PAH abundances around T Tauri stars are typically several tens times smaller than the ISM value (Gorti & Hollenbach, 2008, 2009). We examine the effect of the PAH abundance on disk photoevaporation rates in Section 4.3.

We also implement radiative recombination cooling of H II (Spitzer, 1978), dust-gas collisional cooling (Yorke & Welz, 1996), Ly${\rm\alpha}$ 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). Other collisional excited lines (CELs) can be important cooling in H II regions but they are neglected in this study for simplicity. We discuss the validity of this simplification in Section 4.7.

We do not include O I photoionization explicitly in our calculations. To treat O I cooling in the H II region approximately, while saving computational time, we set the O I abundance as $y_{\text{\rm OI}}\left(1-y_{\text{\rm\text{H{\cal II}}}}\right)$. This approximation is based on the fact that O I ionization energy is close to H I ionization energy. Although a more detail treatment of O I photoionization would be necessary to model the fine structure of the O I and O II regions, we simplify the O I chemistry because O I cooling remains subdominant in the H II region, compared with adiabatic cooling (see also our discussion in the above). 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, whose ionization energy is close to the dissociation energy of CO, is quickly converted to C II  following CO photodissociation. In practice, we assume that the CO dissociation front is located at the same position of the C II ionization front. 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 a chemical species. These approximations above greatly save computational cost. We have checked the validity of this approximation by performing post-process calculations with solving C I photoionization consistently. The results show that the C I region is geometrically thin, with at most a $\sim 10\%$ thickness of the C II and CO regions, and otherwise the structures of C II/CO regions are hardly affected after the post-processing. The details of the chemical reactions are described in Appendix B.

2.5. Radiative Transfer

We solve radiative transfer to calculate photo-chemical reaction rates, photo-heating rates, and dust temperatures consistently. Gas and dust 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 plays 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. The dust absorption of EUV photons is not dominant in our computational domain with the assumed EUV luminosity. 222 The EUV luminosity yields the maximum density of the ionization front to be $n_{\text{\rm H}}\sim 10^{6}\,{\rm cm^{-3}}$ in the innermost region of the computational domain. With this density, H I becomes optically thick against EUV within the length of $\sim 0.1\,{\rm AU}$ near the ionization front. The corresponding H I column density is $N_{\text{\rm HI}}\sim 10^{18}\,{\rm cm^{-2}}$. EUV absorption by dust grains is effective at much higher column densities ($N_{\text{\rm H}}\sim 10^{21}\,{\rm cm^{-2}}$). Hence, the assumption of effectively optically thin dust is valid for $N_{\text{\rm H}}\sim 10^{18}\,{\rm cm^{-2}}$.

FUV radiative transfer is also solved by ray-tracing in order 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 the grain temperatures by solving radiation 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 (FLD) 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 an optically thick disk (Kuiper & Klessen, 2013). Although the FLD approximation does not strictly hold in disk wind regions ($A_{\text{V}}<1$), the region is directly illuminated by the stellar irradiation, and hence, the local radiation field is dominated by the stellar irradiation component rather than the diffuse radiation component. In the region near the photoevaporation base ($A_{\text{V}}\sim 1$), the direct field is attenuated to some extent, but the region is optically thin for the diffusion component. The dust temperatures in this region are largely determined by the direct irradiation as the wind region. We have explicitly checked that the dust temperatures derived with and without including the diffusion component agree with each other well. The difference is $\sim 2\%$ on average, and at most $6\%$. Thus, we conclude that our calculations provide accurate dust temperatures in wind regions. We note that our radiation transfer model has been applied in many studies of massive star formation and feedback effects (Kuiper et al., 2010a, 2012; Kuiper & Yorke, 2013a; Kuiper et al., 2015, 2016), massive accretion disks (Kuiper et al., 2011; Meyer et al., 2017b, a), 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 Conditions

The disk is assumed to consist of an initially neutral gas. There, all of hydrogen nuclei are assumed to be in molecular (\ceH2) form and all of carbon nuclei are in CO at $t=0$.

The initial grain and gas temperatures are set to be $T=T_{\rm dust}=100{\rm~{}K}\left({R}/{1\,{\rm AU}}\right)^{-1/2}$ (e.g., Kenyon & Hartmann, 1987) except the case of $Z=10^{-4}~{}Z_{\odot}$. In the case of $Z=10^{-4}~{}Z_{\odot}$, we first calculate the thermo-chemical structure without updating the density structure for $\sim 1\,{\rm Myr}$ 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,

[TABLE]

where $R$ and $z$ are positions in two dimensional cylindrical polar coordinates $(R,~{}z)=(r\sin\theta,~{}r\cos\theta)$, $h$ is the scale height of a disk, which is defined as $h\equiv c_{s}/\Omega_{\rm K}$, where $c_{s}$ is isothermal sound speed and $\Omega_{\rm K}$ is the Keplerian angular velocity. We denote $n_{0}$ as the mid-plane density of a disk at $1\,{\rm AU}$, and we set $n_{0}=10^{14}\,{\rm cm^{-3}}$. In Eq. (15), the surface density $\Sigma(\simeq\sqrt{2\pi}h\rho_{\rm m};~{}\rho_{\rm m}\text{ is the mid-plane density structure})$ is assumed to have the profile of $\Sigma\propto R^{-1}$. 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) is important in neutral (H I , \ceH2) regions. These two processes drive disk photoevaporation in our simulations. 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 a semi-analytic model to interpret our numerical results (Section 3.4).

3.1. Structure of a Solar Metallicity Disk

3.1.1 Density, Velocity, and Temperature Structures

Figure 1 shows photoevaporative flows from both H II regions and neutral regions. Gas flows from neutral regions are excited by FUV heating. We perform a test simulation in which the FUV heating is initially included but is switched off at the time $t=t_{\rm c}\equiv 100\,{\rm AU}/1{\rm\,km\,s^{-1}}\simeq 4.74\times 10^{2}\,{\rm yr}$. The neutral flows disappear soon after the FUV heating is switched off. We have thus confirmed that FUV is the main driver of the neutral photoevaporative flows in our simulations (Gorti & Hollenbach, 2009; Owen et al., 2012). Note that X-rays, which are not included here, can also drive neutral flows (Alexander et al., 2004; Ercolano et al., 2008, 2009; Gorti & Hollenbach, 2008, 2009; Owen et al., 2010, 2012).

FUV radiation is attenuated by dust once the hydrogen column density $N_{\text{\rm H}}\gtrsim 10^{21}\,{\rm cm^{-2}}$ in the case of $Z=Z_{\odot}$, while EUV radiation is strongly attenuated once H I column density becomes $N_{\text{\rm HI}}\gtrsim 10^{17}\,{\rm cm^{-2}}$. Therefore, FUV photons typically reach and heat the denser regions of a disk than EUV photons. The typical density of the neutral flows, $n_{\text{\rm H}}\sim 10^{5}-10^{7}\,{\rm cm^{-3}}$, is much larger than the typical density of the H II region flow, $n_{\text{\rm H}}\sim 10^{3}-10^{4}\,{\rm cm^{-3}}$, 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 due to gas expansion rather than radiative recombination cooling. The recombination timescale, $t_{\rm rec}\sim 10^{2}\,{\rm yr}~{}(n_{\text{\rm H}}/10^{4}\,{\rm cm^{-3}})^{-1}$, is longer than the sound-crossing time in the ionized gas, $t_{\text{II}}\simeq(100\,{\rm AU}/30{\rm\,km\,s^{-1}})\sim 16\,{\rm yr}$. Hence, the gas flows out of the disk system before recombining.

The heating/cooling processes bring the gas temperature to $\sim 10^{4}{\rm\,K}$ in this region. The corresponding sound speed is $c_{s}\sim 10{\rm\,km\,s^{-1}}$. The gas is accelerated outward by the local pressure gradient. The poloidal velocity ${v}_{\rm p}=\sqrt{v_{r}^{2}+v_{\theta}^{2}}$ reaches a few times of the sound speed ($\sim 30{\rm\,km\,s^{-1}}$) 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 in the region between the H II ionization front and the \ceH2 photodissociation front, while the dominant process is \ceH2 line cooling in the \ceH2 region. Dust-gas collisional cooling becomes dominant among the three coolants in regions with much larger densities. Similar features are observed in previous studies (e.g., Nomura & Millar, 2005; Nomura et al., 2007), but \ceH2 cooling is not included in these studies. Our simulations show 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 $\sim 10^{2}-10^{3}{\rm\,K}$ ($c_{s}\sim 1-3{\rm\,km\,s^{-1}}$). The gas is accelerated by the pressure gradient and achieves $\sim 1-5{\rm\,km\,s^{-1}}$ in the neutral region while it expands.

3.1.2 Distribution of Hydrogen-Bearing Species

\ce

H2 photoevaporative flows are excited through the following processes (Figure 1). \ceH2 advection associated with photoevaporation replenishes \ceH2 molecules into the neutral gas. It makes the height of the H/\ceH2 boundary large. Heinzeller et al. (2011) argue that the H/\ceH2 boundary above a protoplanetary disk can move upward owing to the advection with winds, but hydrodynamics are not directly incorporated in their study. Our hydrodynamical simulations confirm that the H/\ceH2 boundary is actually raised 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 is lower than the replenishing rate of \ceH2. The self-shielding becomes effective when \ceH2 column density is $N_{\text{\rm\ce{H2}}}\gtrsim 10^{14}\,{\rm cm^{-2}}$. We give the self-shielding function of \ceH2 as $f_{\rm shield}={\rm min}[1,~{}(N_{\text{\rm\ce{H2}}}/10^{14}\,{\rm cm^{-2}})^{-0.75}]$ (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 photodissociation by self-shielding rather than dust shielding.

It has been proposed, in the study 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 $z\simeq 70\,{\rm AU}$ at $R=50\,{\rm AU}$ (see Figure 1), while Woitke et al. (2009) shows that it is $z\sim 15-20\,{\rm AU}$ at $R\simeq 50\,{\rm AU}$. Thus, hydrodynamics significantly affects the chemical structure of protoplanetary disks, and the actual chemical structure is different from the results of a hydrostatic calculation.

In the upper regions above the H/\ceH2 boundary in Figure 1, FUV photons are unshielded, and the \ceH2 abundance is determined by the balance between the strong (unshielded) photodissociation and the \ceH2 formation on dust grains. In the lower regions below the H/\ceH2 boundary, gas advection effectively replenishes \ceH2 molecules in addition to the \ceH2 formation on grains. In the H I region, weak (shielded) photodissociation and \ceH2 formation on dust grains determine the \ceH2 abundance. The typical \ceH2 abundance is $y_{\text{\rm\ce{H2}}}\lesssim 10^{-5}$ and remains roughly constant.

3.1.3 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 where the dust shielding factor $\Theta_{3}({A}_{\rm V})$ is large. This indicates that dust is the most important 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 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 (Figure 2).

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 $\theta$ 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 $z\simeq 20\,{\rm AU}$ at $R=50\,{\rm AU}$, while Woitke et al. (2009) shows $z\sim 15\,{\rm AU}$ at $R\simeq 50\,{\rm AU}$. 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

3.2.1 Structure of Photoevaporationing Flow

Figure 5 presents the structure of the photoevaporative flow with different metallicities, $Z=10^{0.5}Z_{\odot}$, $Z=10^{-0.5}Z_{\odot}$, and $Z=10^{-4}Z_{\odot}$ from the top to bottom panel. Although the photoevaporative flow is excited for all these cases, the dense neutral flow only appears with $Z=10^{0.5}Z_{\odot}$ and $Z=10^{-0.5}Z_{\odot}$. Remarkably, we also see that the neutral flow for $Z=10^{-0.5}Z_{\odot}$ is more denser than that for $Z=10^{0.5}Z_{\odot}$. In fact, the typical density of the neutral flow is $n_{\text{\rm H}}\sim 10^{5}-10^{6}\,{\rm cm^{-3}}$ for $Z=10^{0.5}Z_{\odot}$ and $n_{\text{\rm H}}\sim 10^{5}-10^{7}\,{\rm cm^{-3}}$ for $Z=10^{-0.5}Z_{\odot}$. We have confirmed that the density at the base of the neutral flow is almost proportional to $Z^{-1}$ in our simulations. The figure suggests that the metallicity of $Z\gtrsim 10^{-0.5}Z_{\odot}$ 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 $A_{\text{V}}\propto N_{\text{\rm H}}~{}Z$, 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 $Z=10^{-0.5}Z_{\odot}$ disk is much higher than that in the $Z=10^{0.5}Z_{\odot}$ disk (Figure 5). We conclude that, for $Z\gtrsim 10^{-0.5}Z_{\odot}$, 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 $Z\lesssim 10^{-0.5}Z_{\odot}$ and almost ceases at $Z=10^{-4}Z_{\odot}$. 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 a 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).

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 $Z=10^{0.5}~{}Z_{\odot}$, $Z=10^{-0.5}~{}Z_{\odot}$, and $Z=10^{-4}~{}Z_{\odot}$ 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 $Z\geq 10^{-0.5}~{}Z_{\odot}$. These cooling rates decrease with metallicity, and become so small that the adiabatic cooling dominates in both the H II and H I regions at the lowest metallicity $Z=10^{-4}~{}Z_{\odot}$.

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 $\sim Z^{-1}$ and $Z$, 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 $Z^{-1}$ 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, or 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. In the lowest metallicity range of $Z\lesssim 10^{-2}~{}Z_{\odot}$, FUV heating does not even excite neutral photoevaporation.

3.2.2 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.

3.2.3 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 $\Theta_{3}({A}_{\rm V})$ with $Z=10^{0.5}~{}Z_{\odot}$ and $Z=10^{-0.5}~{}Z_{\odot}$, while it is \ceH2 shielding factor $\Theta_{2}(N_{\text{\rm\ce{H2}}})$ with $Z=10^{-4}~{}Z_{\odot}$. 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 $\dot{M}_{\text{ph}}$ by integrating the mass flux component normal to a spherical surface $S$:

[TABLE]

where $d\bm{S}$ is an infinitesimal surface element vector orthogonal to the spherical surface and $r_{S}$ is the radius of $S$. If the specific enthalpy

[TABLE]

is negative at the boundary $S$, we regard that the gas remains bound in the disk. Therefore, we sum up only the gas with $\eta>0$ in Eq. (16). Without this condition, the bound disk which has a very high density can give a large contribution to Eq. (16) even with very small velocity.

Figure 9 shows the resulting $\dot{M}_{\text{ph}}$ of the different metallicity disks estimated by Eq. (16) with $r_{S}=100\,{\rm AU},~{}150\,{\rm AU},~{}200\,{\rm AU},~{}250\,{\rm AU}$. We give each of the dots as the time-averaged value of $\dot{M}_{\rm ph}$ from $t=0$ to $16\,t_{\text{c}}\simeq 7.58\times 10^{3}\,{\rm yr}$.

The figure shows the photoevaporation rates increase with the measuring radius $r_{S}$. This trend indicates $\eta>0$ is satisfied in the outer region than the gravitational radius, where $\eta=0$ (Liffman, 2003).

The gravitational radius is inversely proportional to the gas temperature (Hollenbach et al., 1994; Liffman, 2003). In other words, the gas temperature $T_{\text{esc}}$ necessary for escape is inversely proportional to the radius. If the base temperature $T_{\text{base}}$ decreases more smoothly than $T_{\text{esc}}(\propto r^{-1})$, there is a radius where $T_{\text{base}}=T_{\text{esc}}$. This is regarded as photoevaporative flows are excited from anywhere in the further region. In this case, $\dot{M}_{\text{ph}}$ calculated by Eq. (16) increases with $r_{S}$ and do not converge. In our simulation, the base temperatures decrease as $T_{\text{base}}\propto r^{-\alpha}~{}(\alpha<0.5)$, while $T_{\text{esc}}$ decreases more rapidly following $r^{-1}$. Thus, $\dot{M}_{\text{ph}}$ generally increases with $r_{S}$ as Figure 9 shows at least for $r_{s}\lesssim 250\,{\rm AU}$. We do not consider the further large $r_{S}$ because other effects such as external photoevaportion (e.g., Adams et al., 2004; Facchini et al., 2016) dominate the evolution in such an outer part in typical environments. Including such effects is beyond the scope of the current work.

In the metallicity range of $10^{-4}~{}Z_{\odot}\leq Z\leq 10^{-2}~{}Z_{\odot}$, $\dot{M}_{\text{ph}}$ is roughly independent of metallicity, as shown in the top panel of Figure 9. As discussed in Section 3.2, FUV heating is less efficient than dust cooling in this metallicity range, and neutral flows are not driven. Ionized gas flows are still driven by EUV heating, of which the rate is independent of metallicity, and hence give a roughly constant and dominant contribution to the photoevaporation.

The resulting EUV photoevaporation rate is $\dot{M}_{\rm ph,EUV}\simeq 1.0\times 10^{-9}\,M_{\odot}\,{\rm yr}^{-1}$ 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 $(2.7-7.3)\times 10^{-10}\,M_{\odot}\,{\rm yr}^{-1}$ for $\Phi_{\rm EUV}=6\times 10^{41}~{}{\rm s}^{-1}$ and $M=0.5M_{\odot}$, which is smaller than our $\dot{M}_{\rm ph,EUV}$. 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 those estimated by the 1+1D radiative transfer. The EUV photoevaporation rates seem 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 photoevaporation rates of 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 $10^{-1}~{}Z_{\odot}\leq Z\leq 10~{}Z_{\odot}$, both neutral flows and ionized flows are constantly excited as shown in Figure 5. As discussed in Section 3.2, FUV photons can reach and heat dense parts of the disk when the dust opacity is small. The density of the excited neutral photoevaporative flow is higher for lower metallicity, and the resulting $\dot{M}_{\text{ph}}$ is larger. In the context of massive star formation, a similar dust attenuation effect to regulate the EUV photoevaporation rates is reported in Tanaka et al. (2017).

As presented in EC10, the X-ray photoevaporation rate also increases as metallicity decreases. The metallicity dependence of the photoevaporation rates is approximated as $\dot{M}_{\text{ph}}\propto Z^{\delta}$. The slope $\delta$ is $-0.77$ in the range of $10^{-2}\,Z_{\odot}\leq Z\leq 2\,Z_{\odot}$ in EC10 while that of our study is $-0.85\pm 0.07$ in the metallicity range of $10^{-0.5}\,Z_{\odot}\leq Z\leq 10\,Z_{\odot}$. Despite the model differences between EC10 and our study, these slopes are both negative and take similar values. Clearly, the opacity is an important factor that determines the photoevaporation rate.

In the range of $10^{-2}~{}Z_{\odot}\leq Z\leq 10^{-1}~{}Z_{\odot}$, FUV heating becomes inefficient and cannot balance dust-gas collisional cooling, as discussed in Section 3.2. Then the neutral gas temperature decreases, and the low-temperature gas can evaporate out of the disk only in the outer region, where the central star’s gravity is weak. This means that the minimum radius $r_{\rm min,n}$ where neutral photoevaporative flows are excited gets larger as metallicity decreases. The mass loss rate beyond $r_{S}$ is not counted if $r_{\text{min,n}}>r_{S}$. Therefore, $\dot{M}_{\text{ph}}$ decreases almost suddenly at $Z=10^{-1.2}\,Z_{\odot}$ for $r_{S}=100\,{\rm AU}$ and at $Z=10^{-1.5}\,Z_{\odot}$ for $r_{S}=150\,{\rm AU},~{}200\,{\rm AU},~{}250\,{\rm AU}$.

3.4. Semi-Analytic Model

In this section, we develop a semi-analytic model to interpret our numerical results. As discussed in Sections 3.2 and 3.3, $\dot{M}_{\text{ph}}$ is largely determined by FUV-driven neutral flows with a strong metallicity dependence. We focus on modeling the FUV-driven photoevaporation rate $\dot{M}_{\text{FUV}}^{\text{ana}}$ with different metallicities of $10^{-2}~{}Z_{\odot}\lesssim Z\lesssim 10~{}Z_{\odot}$. Regarding the photoevaporation via the EUV irradiation, we simply assume a constant rate $\dot{M}_{\text{EUV}}^{\text{ana}}=1.0\times 10^{-9}\,M_{\odot}\,{\rm yr}^{-1}$. The EUV photoevaporation rate is taken from our calculation in Section 3.3.

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

The disk system is in a steady state. 2. 2.

Evaporative flows are launched from the regions where ${A}_{\rm V}\sim 1/2$. 3. 3.

All hydrogen are in the molecular form, but CO molecules are completely photodissociated at the base. 4. 4.

Evaporative flows are launched at the speed $\mathscr{M}c_{s}$, where $\mathscr{M}(R,Z)$ is the Mach number and $R$ is distance in cylindrical coordinates. 333 We use the lower-case letter of $z$ as height in cylindrical coordinates to distinguish it from metallicity, which is denoted by the upper-case letter of $Z$. 5. 5.

The azimuthal velocity is given by $v_{\phi}\sim\sqrt{GM_{*}/r}$ at the base. 6. 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. 7.

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

The profile of the base is approximated by a quadratic function

[TABLE]

where the coefficients $a(Z)$ and $b(Z)$ are provided later.

We use $\sim 2A_{\text{V}}$ as the exponent of the dust shielding factor for FUV heating. Therefore, the points where $A_{\text{V}}\sim 1/2$ approximately correspond to the boundary where FUV can reach in the disk. The visual extinction is defined as

[TABLE]

where $N_{\text{\rm H}}$ is hydrogen nuclei column density and $\Sigma_{\rm d}=5.34\times 10^{-22}~{}{\rm mag}~{}{\rm cm}^{2}$ 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. (19) can be approximately rewritten to $N_{\text{\rm H}}\sim n_{\text{\rm H}}r$. Hence, in our semi-analytical model, we approximate the visual extinction as $A_{\text{V}}\sim\Sigma_{\rm d}(Z/Z_{\odot})n_{\text{\rm H}}r$, and the base number density is given by

[TABLE]

Hydrogen dominates the gas mass in our chemistry model, and thus the base density is approximately given by $\rho_{\rm b}={m}_{\rm H}/(2\Sigma_{\rm d}Z/Z_{\odot}r)$.

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. (20). 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 $T_{\rm fit,d}=120{\rm\,K}(r/10\,{\rm AU})^{-0.35}$ 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

[TABLE]

where $r_{0}=100\,{\rm AU}$. Note that $\alpha$ takes a value in the range of $0.28<\alpha<0.40$ with $10^{-2}\,Z_{\odot}\leq Z\leq 10\,Z_{\odot}$.

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

[TABLE]

Thus, the condition, $\eta>0$ (the seventh assumption), corresponds to

[TABLE]

Therefore, the seventh assumption is equivalent to the assumption that photoevaporation is excited from the region where $r>r_{\text{min}}$.

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

[TABLE]

By using all the elements above, we can finally derive $\dot{M}_{\text{FUV}}^{\text{ana}}$. 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

[TABLE]

where $ds$ is a line element of the base and given by $ds=dR\sqrt{1+{f^{\prime}}^{2}}$, and $\beta=\beta(R,Z)$ is the angle of the poloidal velocity $\vec{v}_{\rm p}$ relative to the line element $d\vec{s}$. In our model, Eq. (29) is rewritten to

[TABLE]

where $R_{\text{max}}$ is the upper limit of the integration, and it is set to be the real root of $R_{\text{max}}^{2}+f(R_{\text{max}})^{2}=r_{\text{max}}^{2}=r_{S}^{2}$ in order to compare the model with the simulation results in Section 3.3. The analytical rate $\dot{M}_{\text{FUV}}^{\text{ana}}$ is set to zero if $R_{\text{min}}>R_{\text{max}}$, where $R_{\text{min}}(>0)$ is defined by the real root of $R_{\text{min}}^{2}+f(R_{\text{min}})^{2}=r_{\text{min}}^{2}$. We set $\bar{\mathscr{M}}=0.6$ and $\bar{\beta}=\pi/6~{}{\rm rad}$. These values are determined from the simulation results in the regions where $\eta>0$ 444Though $\mathscr{M}$ and $\beta$ depend on both metallicity and radius in general, we simply take their averages in metallicity and radius.. We approximate the gradient of the base as $\bar{f^{\prime}}=[f(R_{\text{max}})-f(R_{\text{min}})]/(R_{\text{max}}-R_{\text{min}})$. Then, Eq. (30) is rewritten as

[TABLE]

where $T_{\text{X}}\equiv T_{\rm fit}(r_{\text{max}})$.

The model photoevaporation rate $\dot{M}_{\text{model}}=\dot{M}_{\text{FUV}}^{\text{ana}}+\dot{M}_{\text{EUV}}^{\text{ana}}$ is shown by the red line in the top panel of Figure 11.

It is clear that our model explains well the metallicity dependence of the photoevaporation rate derived from our simulations. The discrepancy between the photoevaporation rates of the model and simulations are relatively large in $Z\lesssim 10^{-1}\,Z_{\odot}$. In this metallicity range, adiabatic cooling is comparable to or dominates over the other cooling/heating processes and thus primarily determines temperature in the neutral region. The base temperature is calculated to be higher in our model than in the simulations owing to the absence of adiabatic cooling. This suggests that hydrodynamical simulations are necessary to derive photoevaporation rates when the characteristic dynamical time is comparable to or shorter than the characteristic cooling time.

Eq. (31) also gives the $r$-dependence of $\dot{M}_{\text{FUV}}^{\text{ana}}$ by replacing $r_{\text{max}}$ with $r~{}(r\geq r_{\text{min}})$.

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

[TABLE]

In the equation, $S_{\rm I}(r)$ is the regions where $y_{\text{\rm HII}}<0.5$ in the spherical surface at $r$. We use the condition $y_{\text{\rm HII}}<0.5$ to calculate the contribution of the neutral photoevaporative flow to the photoevaporation rates. Figure 12 compares the $r$-dependence of the analytic photoevaporation rate with that of the simulation results. The model photoevaporation rate of Eq. (31) can explain not only the metallicity dependence of the photoevaporation rates but also the $r$-dependence of the photoevaporation rates.

In Eq. (31), we can give the approximate forms of $r_{\text{min}}$ as functions of metallicity $r_{\text{min}}\simeq 12.4~{}(Z/Z_{\odot})^{-0.55}\,{\rm AU}$. Also, we can approximate the last factor of Eq. (31) to $\sim 0.5$ with an error of less than four percent. With these quantities, Eq. (31) can be further approximated to the form which explicitly depends on $r$ and $Z$:

[TABLE]

where $\tilde{r}_{\text{2}}\equiv r/(10^{2}\,{\rm AU})$ and $\tilde{Z}\equiv(Z/Z_{\odot})$. In Figure 12, $\dot{M}_{\text{FUV}}^{\text{ana,fit}}$ (Eq. (33)) is shown by the green solid line and compared with $\dot{M}_{\text{FUV}}^{\text{ana}}$ (Eq. (31)) denoted by the red solid line.

4. DISCUSSION

4.1. Comparison with X-ray Photoevaporation

In the present study, $\dot{M}_{\text{ph}}$ increases with decreasing metallicity in the range of $10^{-1}~{}Z_{\odot}\leq Z\leq 10~{}Z_{\odot}$, which is similar to that of EC10, while $\dot{M}_{\text{ph}}$ decreases with metallicity in the range of $10^{-2}~{}Z_{\odot}\leq Z\leq 10^{-1}~{}Z_{\odot}$, which is different from that of EC10, who derive EUV/X-ray photoevaporation rates. Evidently, it is worth investigating the metallicity dependence of FUV/EUV/X-ray photoevaporation. In future work, we plan to incorporate X-ray radiative transfer in our photoevaporation model, and derive the metallicity dependence of photoevaporation excited by EUV/FUV/X-ray.

4.2. Disk Lifetime

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 $\Sigma\propto R^{-p}$ to $p=1$, the formula is given by $T_{\rm life}\propto\dot{M}_{\rm ph}^{-2/3}$ . 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.

Figure 13 compares the estimated lifetimes from the photoevaporation rates of our study with the metallicity dependence of the observational lifetimes. We choose $3\,{\rm Myr}$ at $Z=Z_{\odot}$ (Haisch et al., 2001) and $1\,{\rm Myr}$ at $Z=0.2~{}Z_{\odot}$ (Yasui et al., 2009, 2010) as the observational lifetimes. The formula merely gives the simple relation between the lifetime and photoevaporation rate. Therefore, it needs to be normalized to give a specific lifetime from a photoevaporation rate. In order to enable the comparison between the observational lifetimes and the estimated lifetimes, we normalize the formula so that it gives the same lifetime as the observational lifetime at $Z=Z_{\odot}$. The estimated lifetime is $\sim 2\,{\rm Myr}$ with $Z=0.2~{}Z_{\odot}$ in Figure 13, while the observational lifetime is $1\,{\rm Myr}$. Thus, the metallicity dependence of the estimated lifetimes is less steep than that of the observational lifetimes between $0.2~{}Z_{\odot}\lesssim Z\lesssim Z_{\odot}$.

We have only two data points from observations in Figure 13. A reasonable and meaningful comparison with our model requires more observational data. In addition, disk lifetimes are influenced by the accretion process, which has not been fully understood yet. However, at least, we can report here that the slope of the estimated lifetimes in the range of $10^{-0.5}~{}Z_{\odot}\lesssim Z\lesssim 10~{}Z_{\odot}$, $(-0.85)\times(-2/3)=0.57$, is quite consistent with that of the observational lifetimes, $0.68$. 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.

It has been observationally shown that the gas-giant occurrence decreases with the host star’s metallicity at $Z\gtrsim 10^{-0.5}~{}Z_{\odot}$, which is called “planet-metallicity correlation” (e.g., Gonzalez, 1997; Johnson et al., 2010; Mortier et al., 2013). The apparent correlation is thought to reflect the fact that planet formation is inefficient in a low-metallicity disk. Interestingly, EC10 shows that in the context of core accretion scenario the higher planet occurrence is attributed mainly to the faster core growth due to a larger amount of solids in a higher-metallicity disk rather than the reduced lifetimes due to X-ray photoevaporation. On the other hand, Wang & Fischer (2015) conclude the terrestrial planet occurrence is not as strongly dependent on metallicity as the gas-giant occurrence. This observational result would suggest that core growth of planets is not so strongly dependent on metallicity. In that case, metallicity dependence of FUV and/or X-ray photoevaporation could have effects on metallicity dependence of gas giant occurrence.

4.3. Grain Effects on FUV Photoevaporation

Photoelectric heating generally depends on both the local dust-to-gas mass ratio and the local size distribution of dust/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 (Takeuchi et al., 2005; Owen et al., 2011; Hutchison et al., 2016a, b). 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.

The PAH abundance significantly affects the gas-grain photoelectric heating rate and hence the resulting FUV photoevaporation rate (Gorti & Hollenbach, 2008, 2009). In our fiducial model, we adopt the ISM value, which may be larger than the PAH abundances around T Tauri stars (Geers et al., 2007; Oliveira et al., 2010; Vicente et al., 2013). Although there remains large uncertainties in the observationally determined PAH abundances, it is worth examining the overall impact of the assumed PAH abundance on our results.

According to Bakes & Tielens (1994), about a half of the total photoelectric heating rate is contributed by the grain species with sizes smaller than $\lesssim 15\,\text{\AA}~{}(N_{\text{C}}\lesssim 1500)$. By reducing the FUV heating rate to a half of that given by Eq. (A7), we can approximate an effective photoelectric heating rate without PAH contribution.

We perform additional simulations with using the reduced FUV heating rate. The resulting photoevaporation rates are shown in Figure 14.

Clearly, the FUV-driven flows contribute to $\dot{M}_{\text{ph}}$ in the range of $Z\gtrsim 0.1\,Z_{\odot}$ even in the case there is no PAH contribution.

The halved FUV heating makes base temperatures lower. The lower base temperatures yield the result that photoevaporative flows are excited in the outer region of the disks. In $Z\gtrsim Z_{\odot}$, the base temperatures are still high to excite photoevaporative flows in the large part of the disk even if PAHs do not contribute to the FUV heating. Consequently, the abundance and the size of PAH do not significantly change the photoevaporation rates. In $0.1\,Z_{\odot}\lesssim Z\lesssim Z_{\odot}$, the dust-gas collisional cooling is effective enough to suppress the excitation of photoevaporative flows even at $Z\sim 10^{-0.3}\,Z_{\odot}$. As a result, $\dot{M}_{\text{ph}}$ drops at higher metallicity in the case PAHs are absent than in the case PAHs exist.

Although the small PAH abundances result in smaller FUV photoevaporation rates as we demonstrate, other grain effects could also affect FUV photoevaporation rates (Gorti et al., 2015). For example, if dust growth and settling are incorporated, disk opacity for UV photons would be reduced. In this case, photoevaporation rates are increased because UV photons reach the higher density interior of the disk. Actually, a low visual extinction $A_{\text{V}}\sim 0.1-0.2$ is observed for high column density regions with $N_{\text{\rm H}}\sim 10^{22}\,{\rm cm^{-2}}$ (Vicente et al., 2013). This suggests that the effects of grain growth/settling might deplete dust grains with the size of $\sim 0.1\,{\rm\mu m}$ in the neutral region, and the effects would reduce the disk opacities for UV photons. Hence, for a comprehensive modeling of FUV photoevaporation, we need to take account of not only the reduced PAH abundance but also other effects such as grain growth, destruction/fragmentation, and settling.

4.4. MHD Wind

Magneto-hydrodynamics (MHD) driven disk wind has been proposed as another important mechanism for disk evolution (Suzuki & Inutsuka, 2009; Suzuki et al., 2010; Armitage et al., 2013; Bai & Stone, 2013a, b; Fromang et al., 2013; Lesur et al., 2013; Simon et al., 2013a, b; Gressel et al., 2015; Bai, 2016; Suzuki et al., 2016; Bai & Stone, 2017). Magnetorotational instability (MRI) excites turbulence that can drive a wind from disk surfaces (Suzuki & Inutsuka, 2009; Suzuki et al., 2010). Recent non-ideal MHD studies show that MRI is mostly suppressed because of low ionization degree in the interior of a disk, but magneto-centrifugal winds can be launched from disk surfaces (Bai & Stone, 2013b; Gressel et al., 2015). The winds extract the disk angular momentum and can promote the accretion onto the central star.

The MHD effects on disk evolution and photoevaporation have been studied independently, but the interplay between them is an important question for realistic modeling of the dispersal. Bai (2016) and Bai et al. (2016) examine evolution of protoplanetary disks with incorporating MHD and external thermal heating (irradiation). They show that the wind mass loss rate is actually affected by both the strength of magnetic field and thermal heating, and that it is characterized by the ratio of sound speed and Alfven speed at the base.

The base position and density, and the ionization degree there depend critically on metallicity. Therefore, for a complete picture of metallicity dependence of disk dispersal, it would be necessary to study the global evolution of protoplanetary disks with both MHD and photoevaporation with a detailed treatment of the relevant thermal processes.

4.5. 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 photoevaporative flow and the outer boundary of computational domain, where the flow is not yet accelerated up to $\mathscr{M}>1$. 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 at the sound speed. It can make photoevaporation profiles and rates have some features which vary with the timescale of the order of the crossing time. The effects would be significant when the outer boundary is so small that the computational domain does not contain transonic points of photoevaporative flows. In order to examine the effect of the bogus reflection, we carry out simulations with a smaller outer boundary $L=100\,{\rm AU}$, keeping the numbers of the cells the same as the settings described in Section 2.1.

Figure 15 shows the time-evolution of $\dot{M}_{\text{ph}}$ derived by Eq. (16) with $r_{S}=80\,{\rm AU}$. We find $\dot{M}_{\rm ph}$ varies periodically in runs with $Z\geq 10^{-1.5}\,Z_{\odot}$. Evaporative flows are driven by FUV in this metallicity range, but they are not yet accelerated to $\mathscr{M}>1$ at the outer boundary. The flows are spuriously reflected, causing oscillational trends in the photoevaporation rates. With $10^{-4}~{}Z_{\odot}\leq Z\leq 10^{-1.5}~{}Z_{\odot}$, the time evolution of $\dot{M}_{\text{ph}}$ is almost independent of metallicity and always similar to that of $Z=10^{-2}~{}Z_{\odot}$ indicated by the magenta line in Figure 15. Only the ionized flows contribute to $\dot{M}_{\text{ph}}$ in this metallicity range. The flow velocity usually exceeds the sound speed soon after they are launched at the base, and thus the oscillation does not appear in $\dot{M}_{\text{ph}}$.

We have also performed simulations with $L=200\,{\rm AU}$ and confirmed that the oscillation is actually caused by the bogus reflection. The top panel of Figure 16 shows the time-dependent $\dot{M}_{\text{ph}}$ of the $Z=Z_{\odot}$ and $Z=0.1~{}Z_{\odot}$ disks, where $r_{S}=80\,{\rm AU}$ again.

The oscillation of $\dot{M}_{\text{ph}}$ for $Z=Z_{\odot}$ is damped with time in the case of $L=200\,{\rm AU}$. 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 $r\sim 100\,{\rm AU}$, where the neutral gas has a subsonic velocity and contributes significantly to $\dot{M}_{\text{ph}}$. In the case of $Z=0.1\,Z_{\odot}$, the oscillation is still found with $L=200\,{\rm AU}$. We run simulations with the outer boundaries of $L=300\,{\rm AU}$ and $L=400\,{\rm AU}$. In these cases, the numerical oscillations disappear and the photoevaporation rates converge after $\sim 1000$ years calculation (bottom panel of Figure 16). With the enlarged outer boundaries, the transonic point of the photoevaporative flow lies inside the computational domain. The outgoing gas has a supersonic velocity at the boundary, and so the artificial reflection does not occur. Clearly, in order to eliminate or mitigate the outer boundary effect and to obtain converged photoevaporation rates, the transonic points of streamlines in photoevaporative flow should be included in the computational domain.

We perform simulations with $L=100,200,300,400\,{\rm AU}$ with all the metallicities we consider here. Figure 17 shows the resulting time-averaged photoevaporation rates.

The bogus reflection affects the resulting photoevaporation rates especially in the sub-solar metallicity range, where inefficient photoelectric heating reduces gas temperature and yields slower flow velocity. The photoevaporation rates converges if sufficiently large outer boundaries are set, as the red and green lines show in Figure 17. Note that the photoevaporation rates are in better agreement with the analytical one. This result just reflects the fact that excluding the bogus reflection by using large outer boundaries allows the simulations to reach a steady state, as assumed in the analytical model.

4.6. Photoevaporation Estimate

As discussed in Section 3.3, $\dot{M}_{\text{ph}}$ increases with $r_{S}$. In order for $\dot{M}_{\text{ph}}$ to converge, the mass flux should have a dependence $\rho_{\text{base}}v_{\text{base}}\propto R^{p}$ with $p<-2$ because ${\rm d}\dot{M}_{\text{ph}}\propto\rho_{\text{base}}v_{\text{base}}R^{2}~{}{\rm d}(\log R)$. The radius dependence of the base velocity $v_{\text{base}}$ is generally not so strong as $\rho_{\text{base}}$, so $\rho_{\text{base}}$ needs to be $\rho_{\text{base}}\propto R^{p^{\prime}}$ with $p^{\prime}\lesssim-2$. In our simulations, the base density has $p^{\prime}\geq-(1.2-1.5)$. Therefore, $\dot{M}_{\text{ph}}$ does not converge until $r_{S}$ reaches the disk edge. Similar results are reported by Tanaka et al. (2013) for EUV photoevaporation, and also can be inferred from the figures of Gorti & Hollenbach (2009) and the figure 4 of Owen et al. (2010), where the cumulative X-ray photoevaporation rate does not converge with radius up to $70\,{\rm AU}$. In summary, in order to obtain the total photoevaporation rate, the computational domain should contain the whole disk.

4.7. CELs in an H II regions

CELs such as O II (3730 Å, 3727 Å), N II (6585 Å, 6550 Å), O III (88.36 $\mu$m, 51.81 $\mu$m, 5008 Å, 4960 Å), Ne II (12.81 $\mu$m), S II (6733 Å, 6718 Å), and S III (33.48 $\mu$m, 18.71 $\mu$m, 9071 Å, 9533 Å) can be important cooling sources in an H II region especially when the gas metallicity is higher than the solar metallicity (Draine, 2011). The total CEL cooling rate is estimated to be $\sim 10^{3}~{}(n_{\text{\rm H}}/10^{3}\,{\rm cm^{3}})(Z/Z_{\odot})~{}{\rm erg}~{}{\rm g}^{-1}~{}{\rm s}^{-1}$ in a typical H II region (Draine, 2011). The rate of adiabatic cooling, which is shown to be an important cooling in photoevaporative winds that we study here, is given as $P(d/dt)(1/\rho)=(P/\rho)\nabla\cdot\bm{v}\sim c_{s}^{3}/r\sim 10^{6}~{}(c_{s}/30{\rm\,km\,s^{-1}})^{3}(r/1\,{\rm AU})^{-1}~{}{\rm erg}~{}{\rm g}^{-1}~{}{\rm s}^{-1}$. Thus, the CEL cooling can dominate over the adiabatic cooling in the inner, high density part of an H II region. If we choose $n_{\text{\rm H}}\sim 10^{6}\,{\rm cm^{-3}}~{}(r/1\,{\rm AU})^{-1.5}$ as a typical base density (Hollenbach et al., 1994; Tanaka et al., 2013), the CEL cooling dominates over the adiabatic cooling in $r\lesssim 1\,{\rm AU}$, where the base temperature can be lowered by a few to several tens percent. Note that this is outside our computational domain. The ionized gas will be then more strongly bound there and the contribution from $r\lesssim 1\,{\rm AU}$ to the photoevaporation rate is sufficiently small, as discussed already. In the wind region (atmosphere), the density is much smaller than the base density, and thus the CELs are less important. Overall, the impact of the CELs to the net photoevaporation rate is unimportant in our simulations.

5. SUMMARAY

We have performed radiation-hydrodynamical simulations of photoevaporation of protoplanetary disks with self-consistent modeling of multi-species chemistry. In particular, we have considered a broad range of metallicities from $10^{-4}~{}Z_{\odot}$ to $10~{}Z_{\odot}$ to examine the metallicity dependence, if any, of the disk lifetime. Our findings are summarized as follows:

• •

As metallicity decreases, dust shielding effect is reduced and FUV photons reach and heat denser regions of the disk. Thus $\dot{M}_{\text{ph}}$ increases with decreasing metallicity in the range of $10^{-1}\,Z_{\odot}\lesssim Z\lesssim 10\,Z_{\odot}$.

• •

As metallicity decreases, FUV photoelectric heating becomes less efficient than cooling in neutral regions. The temperature decreases, so that a large portion of the disk gas is gravitationally bound. This reduces the contribution of the FUV-driven neutral flow to photoevaporation at $10^{-2}\,Z_{\odot}\lesssim Z\lesssim 10^{-1}\,Z_{\odot}$.

• •

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

• •

In the metallicity range of $10^{-4}\,Z_{\odot}\lesssim Z\lesssim 10^{-2}\,Z_{\odot}$, EUV photons primarily drive photoevaporative flows if X-ray is not considered. Hence, $\dot{M}_{\text{ph}}$ is nearly independent of metallicity in this extremely low-metallicity environment.

• •

We develop a semi-analytical model of disk photoevaporation that describe accurately both the metallicity dependence of $\dot{M}_{\text{ph}}$ (See Figure 9) and the outflow profile $\dot{M}_{\text{FUV}}^{\text{sim}}(r)$ (See Figure 12).

• •

Generally, $\dot{M}_{\text{ph}}$ cumulatively increases with the radius where they are measured. Hence, $\dot{M}_{\text{ph}}$ depends on the disk radius. Global simulations are necessary to derive total photoevaporation rates.

• •

In numerical simulations, the bogus reflection at the outer boundary affects photoevaporative flow profiles. It can even lead to a wrong conclusion regarding photoevaporation rates. A sufficiently large outer boundary should be used so that it can contain transonic points of photoevaporative flows, or one would need to use a non-reflecting boundary condition.

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. Their photoevaporation rates are derived by hydrostatic calculations, and thus are subjected 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/EUV/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, Kei Tanaka, and Xuening Bai for helpful discussions and insightful comments on the paper. We also thank the anonymous referee for giving practical comments to improve the manuscript. 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

[TABLE]

where $\nu$ is a frequency of EUV photons, $F_{\nu}$ is the specific number flux of the direct EUV field, $n_{\text{\rm\text{H{\cal I}}}}$ is number density of H I, and $\sigma_{\nu}$ is the absorption cross section of H I. We use the approximate absorption cross section

[TABLE]

(e.g., Osterbrock & Ferland, 2006). In Eq. (A2), $\nu_{1}$ is the frequency at the Lyman limit ($h\nu_{1}\simeq 13.6{\rm\,eV};\lambda_{1}\equiv c/\nu_{1}\simeq 91.2{\rm nm}$, $c$ is the light speed). Eq. (A1) can be solved analytically

[TABLE]

where $\Phi_{\nu}(R_{*})$ is the specific photon number luminosity of the EUV emitted from the stellar surface, and $N_{\text{H{\cal I}}}=N_{\text{H{\cal I}}}(r,\theta,t)$ is the column density of hydrogen atoms between the stellar surface and a certain point in the computational domain:

[TABLE]

With Eq. (A3), the photoionization rate and the specific photoionization heating rate are given as

[TABLE]

respectively. We assume that the spectral energy distribution of the EUV photons is given by a black body spectrum with an effective temperature of $T_{\rm eff}=10^{4}~{}{\rm K}$. The corresponding total stellar EUV luminosity is then $\Phi_{\rm EUV}\simeq 1.5\times 10^{41}(R_{*}/R_{\odot})^{2}~{}{\rm s}^{-1}$. We use 81 frequency bins for the numerical integrations.

For the photoelectric heating rate, we use the analytic formula presented in BT94. BT94 derives the photoelectric heating rate theoretically by using the dust size distribution of the MRN dust model (Mathis et al., 1977). The analytic formula is

[TABLE]

where $\epsilon_{\rm pe}$ is the photoelectric effect efficiency of the grains, which corresponds to the ratio of the gas heating rate to FUV absorption rate of the grains, $\gamma_{\rm pe}$ is the ratio of the dust/PAH photoionization rate to the dust/PAH recombination rate with electrons, which is given by $\gamma_{\rm pe}\equiv G_{\rm FUV}\sqrt{T}/n_{\text{\rm e}}$. $G_{\rm FUV}$ is the FUV flux ($6{\rm\,eV}<h\nu<13.6{\rm\,eV}$) at the local point in the unit of the averaged interstellar flux $F_{\text{ISRF}}=1.6\times 10^{-3}~{}{\rm erg~{}cm^{-2}~{}s^{-1}}$, and given by $G_{\rm FUV}=L_{\text{FUV}}e^{-1.8A_{\text{V}}}/(4\pi r^{2}F_{\text{ISRF}})$. The last factor $Z/Z_{\odot}$ of Eq. (A7) is multiplied to take into account the effect of the grain amount on the heating rate.

A.2. Dust-Gas Collisional Cooling

Dust grains act as a cooling/heating agent for a gas via collisional heat transfer between gas and dust. We adopt the dust-gas collisional cooling function of Yorke & Welz (1996), which is given by

[TABLE]

where $a_{\rm dust},~{}\rho_{\rm dust},~{}m_{\rm dust}$, and $T_{\rm dust}$ are the mean dust size, dust mass density, mean dust mass, and dust temperature, respectively. We use the dust parameters of Yorke & Welz (1996); $a_{\rm dust}=5\times 10^{-6}{\rm~{}cm}$ and $m_{\rm dust}=1.3\times 10^{-15}{\rm~{}g}$.

A.3. Atomic/Molecular Line Cooling

We implement radiative recombination cooling of H II, Ly${\rm\alpha}$ cooling of H I, C II line cooling, O I line cooling, \ceH2 line cooling, and CO line cooling as line cooling sources of gas.

When a hydrogen ion recombines with a free electron in the H II region, approximately two-thirds of the electron energy $\sim kT$ is lost by radiative recombination (e.g., Spitzer, 1978). We adopt the radiative recombination cooling rate

[TABLE]

where $R_{\rm k2}$ is the reaction coefficient of H II recombination (the reaction labeled “k2” in Table B.3).

A neutral hydrogen atom is excited by collision and then de-excited by emitting a Ly${\rm\alpha}$ photon. We simply refer to this cooling process as Ly${\rm\alpha}$ cooling. We use the Ly${\rm\alpha}$ cooling function presented in Anninos et al. (1997):

[TABLE]

where $T_{\rm K}$ is gas temperature in Kelvin.

C II and O I  have fine-structure transitions, and they work as line cooling sources by spontaneous emissions. The total line cooling rate of each atom is given by the equation

[TABLE]

The label X indicates C II or O I. The indexes, $i,j~{}(=1,~{}2,~{}3,~{}...)$, are the label of an energy level, $x_{j}$ is population of level $j$, $A_{ji}$ is the Einstein A coefficient of the transition $j\rightarrow i$, and $\Delta E_{ji}$ is its corresponding energy difference, respectively. Each of the level populations is derived by solving the equations of statistical equilibrium simultaneously:

[TABLE]

where

[TABLE]

The collisional excitation (de-excitation) rate of the transition $i\rightarrow j$ with a collisional counterpart ${\rm\lambda}$ is represented as $\gamma_{ij}^{\rm\lambda}$. In this study, we treat these line emissions as optically thin for simplicity. Therefore, an escaping probability, absorption of external radiation, or induced radiation is ignored in Eq. (A15).

\ce

H2 and CO molecules have rovibrational transitions, and they also work as line cooling sources. These cooling rates could also be directly calculated by Eq. (A13), but we use the analytic formula presented in Galli & Palla (1998) for \ceH2 line cooling and the tabulated values presented in Omukai et al. (2010) for CO line cooling.

Appendix B Chemical Reactions

We take into account all the chemical reactions listed in Table B.3. We include photo-chemical reactions such as H I photoionization, \ceH2 photodissociation, and CO photodissociation as well as collisional reactions. We summarize these chemical reactions in this section.

B.1. Photodissociation of \ceH2

\ce

H2 is photodissociated by FUV photons in the energy range of $11.2{\rm\,eV}\lesssim h\nu\lesssim 13.6{\rm\,eV}$ as follows: \ceH2 is pumped up to an upper electronic bound state by absorbing an FUV photon, it goes back to an excited vibrational state of the ground electronic state with fluorescence emission, or it goes to a continuum vibrational state of the ground electronic state and then is photodissociated with fluorescence emission. Photodissociation occurs to $\sim 10$ - $15\%$ of pumped up \ceH2 molecules.

FUV photons are shielded by dust and \ceH2 molecules themselves. We adopt the photodissociation rate function and the self-shielding function presented in Draine & Bertoldi (1996). The photodissociation rate per unit volume is given by

[TABLE]

where $\tau_{\rm d,1000}$ is optical depth of dust at the wavelength of $1000$Å, and $I_{\rm diss}\simeq 4\times 10^{-11}~{}G_{\rm FUV2}~{}{\rm s^{-1}}$ is the unshielded photodissociation rate. The definition of $G_{\rm FUV2}$ is similar to Appendix A and given by $G_{\rm FUV2}=L_{\text{FUV}}/(4\pi r^{2}F_{\text{ISRF}})$. The self-shielding function $f_{\rm shield}$ is

[TABLE]

where $N_{0}\equiv 10^{14}\,{\rm cm^{-2}}$.

B.2. Photodissociation of \ceCO

CO is also photodissociated by processes similar to \ceH2 photodissociation. CO photodissociation is also a line process like \ceH2 photodissociation. CO shields FUV photons once the column density of CO becomes large. In addition, CO is shielded by \ceH2 molecules owing to line overlap. Note that \ceH2 is a more abundant species than CO.

We adopt the CO photodissociation function and shielding functions presented in Lee et al. (1996):

[TABLE]

where $p_{\rm diss}=1.03\times 10^{-10}{\rm s^{-1}}$ is the unshielded photodissociation rate of CO. The factors $\Theta_{1},\Theta_{2}$, and $\Theta_{3}$ are the self-shielding factor, the \ceH2 shielding factor, and the dust shielding factor, respectively. These quantities are tabulated in the table of Lee et al. (1996).

B.3. Carbon Chemistry

We assume that C I is converted to C II as soon as it is produced by the photodissociation of CO, as presented in Richling & Yorke (2000). In other words, we assume that the CO dissociation front coincides with the C II ionization front.

As the reverse process of CO photodissociation, we use the simplified chemistry model for CO formation described in Nelson & Langer (1997). In the model, CO formation is initiated by the reaction, \ceC+ + H2 -¿ CH2+ + $\gamma$ (reaction rate ; $k_{0}=5\times 10^{-16}{\rm~{}cm^{3}~{}s^{-1}}$). \ceCH2+ ions rapidly convert to \ceCH and \ceCH2 by dissociative recombination with electrons and ion-molecule reactions with \ceH2 molecules. These hydrocarbon radicals react with oxygen atoms to form CO (reaction rate; $k_{1}=5\times 10^{-10}{\rm~{}cm^{3}~{}s^{-1}}$) or are photodissociated to form ions (dissociation rate; $\Gamma_{\ce{CH}_{\rm x}}$). As a result, the effective CO formation rate per unit volume is given,

[TABLE]

(for the detailed derivation of this formula, see, e.g., Nomura & Millar (2005)). We set the dissociation rate of the hydrocarbons $\Gamma_{\ce{CH}_{\rm x}}$ to five times of the CO dissociation rate for simplicity.