Hydrodynamic Photoevaporation of Protoplanetary Disks with Consistent Thermochemistry

TL;DR

This study uses advanced hydrodynamic simulations with thermochemistry to analyze how ultraviolet and X-ray radiation drive photoevaporative winds in protoplanetary disks, revealing detailed wind structures and chemical reformation processes.

Contribution

It introduces a comprehensive simulation approach combining hydrodynamics, radiative transfer, and thermochemistry to study disk photoevaporation with consistent chemistry modeling.

Findings

Most models show a three-layer disk structure with a hot wind at ~30 km/s.

Reformation of molecules like CO, OH, and H2O occurs in the intermediate layer.

Mass-loss rates depend on radiation intensity and modeling assumptions.

Abstract

Photoevaporation is an important dispersal mechanism for protoplanetary disks. We conduct hydrodynamic simulations coupled with ray-tracing radiative transfer and consistent thermochemistry to study photoevaporative winds driven by ultraviolet and X-ray radiation from the host star. Most models have a three-layer structure: a cold midplane, warm intermediate layer, and hot wind, the last having typical speeds $\sim 30~\mathrm{km\ s}^{-1}$ and mass-loss rates $\sim 10^{-9}~M_\odot~\mathrm{yr}^{-1}$ when driven primarily by ionizing UV radiation. Observable molecules including CO, OH and H2O re-form in the intermediate layer and survive at relatively high wind temperatures due to reactions being out of equilibrium. Mass-loss rates are sensitive to the intensity of radiation in energy bands that interact directly with hydrogen. Comparison with previous works shows that mass loss rates are also sensitive to the treatment of both the hydrodynamics and the thermochemistry. Divergent results concerning the efficiency of X-ray photoevaporation are traced in part to differing assumptions about dust and other coolants.

Tables

Table 1. Table 1 Properties of the fiducial model

Item	Value
Radial domain	$2\text{au}\le r\le 100\text{au}$
Latitudinal domain	$0\le \theta \le \pi /2$
Resolution	$N_{\log r}=256$, $N_{\theta }=128$
Stellar mass	$1.0M_{\odot }$
$M_{\mathrm{disk}}$	$0.03M_{\odot }$
Mid-plane density	${10}^{10}{\left(R/100\text{au}\right)}^{-2.25}{\mathrm{cm}}^{-3}$
Mid-plane temperature	$20{\left(R/100\text{au}\right)}^{-0.5}K$
Luminosities [photon $s^{-1}$]
$7\mathrm{eV}$ (“soft” FUV)	$4.5\times {10}^{42}$
$12\mathrm{eV}$ (LW)	$1.6\times {10}^{40}$
$25\mathrm{eV}$ (EUV)	$5.0\times {10}^{40}$
$1\mathrm{keV}$ (X-ray)	$1.6\times {10}^{39}$
Initial abundances [$n_X/n_H$]
$H_2$	0.5
He	0.1
$H_{2}O$	$1.8\times {10}^{-4}$
CO	$1.4\times {10}^{-4}$
S	$2.8\times {10}^{-5}$
Si	$1.7\times {10}^{-6}$
Fe	$1.7\times {10}^{-7}$
Gr	$1.0\times {10}^{-7}$
Dust/PAH properties
$r_{\mathrm{dust}}$	$5\AA$
${\rho }_{\mathrm{dust}}$	$2.25g{\mathrm{cm}}^{-3}$
$m_{\mathrm{dust}}/m_{\mathrm{gas}}$	$7\times {10}^{-5}$
${\sigma }_{\mathrm{dust}}/H$	$8\times {10}^{22}{\mathrm{cm}}^2$

Table 2. Table 2 Models exploring parameter space

No.	Description	${\dot{M}}_{\mathrm{wind}}$	${\dot{M}}_{\mathrm{GH}}$	Total heating	Efficiency	$⟨v_r⟩$
		$\left({10}^{-9}{M}_{\odot }{\mathrm{yr}}^{-1}\right)$	$\left({10}^{-9}{M}_{\odot }{\mathrm{yr}}^{-1}\right)$	$\left({10}^{30}\mathrm{erg}{s}^{-1}\right)$		$\left(\mathrm{km}{s}^{-1}\right)$
(1)	(2)	(3)	(4)	(5)	(6)	(7)
0	Fiducial	2.5 $\pm$ 0.2	11.6	4.4	0.67	39
1	No LW photons	2.5 $\pm$ 0.3	9.3	4.1	0.67	38
2	100$\times$ LW photons	17.6 $\pm$ 2.1	61.3	9.1	0.60	18
3	No ”soft” FUV	1.1 $\pm$ 0.1	2.7	2.3	0.53	58
4	”Soft” FUV only	0.0	0.2	1.0	-	-
5	No EUV	0.0	3.7	1.6	-	-
6	10$\times$ EUV photons	9.4 $\pm$ 0.7	107.8	26.7	0.74	33
7	No X-ray	2.1 $\pm$ 0.2	6.9	2.6	0.80	38
8	10$\times$ X-ray photons	9.1 $\pm$ 0.4	55.6	14.1	0.42	24
9	$r_{\mathrm{dust}}=10\AA$	2.8 $\pm$ 0.6	10.5	3.5	0.68	30
10	$\mathrm{OECA10}{\mathrm{analog}}^{†}$	11.2 $\pm$ 4.2	105.0	0.8	0.54	5
11	Convergence test	2.7 $\pm$ 0.6	16.1	3.0	0.58	32

Table 3. Table 3 Approximate line luminosity

Model	O I $63\mu m$	O I $6300\AA$	S I $25\mu m$	Si II $35\mu m$	$H_2$ ro-vib	$\mathrm{OH}$/$H_{2}O$ ro-vib	$\mathrm{CO}$ ro-vib
0	$-4.31$	$-7.36$	$-4.00$	$-5.65$	$-4.67$	$-4.29$	$-3.89$
1	$-4.34$	$-8.17$	$-3.96$	$-5.73$	$-4.60$	$-4.36$	$-3.90$
2	$-3.75$	$-4.69$	$-4.17$	$-4.73$	$-4.11$	$-3.57$	$-3.91$
3	$-5.05$	$-7.72$	$-5.05$	$-6.48$	$-5.26$	$-3.61$	$-5.68$
4	$-5.05$	-	$-4.68$	-	$-5.61$	$-4.30$	$-4.26$
5	$-4.22$	$-7.31$	$-3.97$	$-5.58$	$-5.03$	$-4.63$	$-3.88$
6	$-4.43$	$-7.13$	$-3.88$	$-5.36$	$-3.65$	$-3.56$	$-3.78$
7	$-4.74$	$-8.02$	$-4.61$	$-6.30$	$-4.72$	$-4.26$	$-4.23$
8	$-3.57$	$-5.99$	$-3.11$	$-4.69$	$-4.15$	$-3.38$	$-3.37$
9	$-4.47$	$-4.94$	$-4.39$	$-5.94$	$-4.39$	$-4.07$	$-4.59$
10	$-4.38$	$-6.88$	-	$-4.55$	-	-	-
11	$-4.49$	$-5.24$	$-4.32$	$-5.98$	$-4.79$	$-4.94$	$-4.83$

Equations

$$F_{\mathrm{eff}}(\nu )=F_0(\nu )\left\{\frac{1-\exp [-\delta l/\lambda (\nu )]}{\delta l/\lambda (\nu )}\right\},$$

$$\begin{array}{rl}\frac{\mathrm{d}{n}^i}{\mathrm{d}t}& ={\mathcal{A}}_{jk}^i{n}^j{n}^k+{\mathcal{B}}_j^i{n}^j;\\ \frac{\mathrm{d}ϵ}{\mathrm{d}t}& =\Gamma -\Lambda ;\end{array}$$

$$\mathcal{B}\equiv \frac{v^2}{2}+\frac{\gamma p}{(\gamma -1)\rho }+\Phi ,$$

$${\mathrm{H}}_2+h\nu (\mathrm{LW})\to \{\begin{array}{ll}2\mathrm{H},& \sim 10\text{per cent};\\ {{\mathrm{H}}_2}^{\ast },& \sim 90\text{per cent}.\end{array}$$

$$\begin{array}{rl}{k}_{\mathrm{de}}(\mathrm{H})& \simeq 1.8\times 10^{-13}{\mathrm{cm}}^3{\mathrm{s}}^{-1}\times {\left(\frac{T}{\mathrm{K}}\right)}^{1/2}\exp \left(-\frac{1000\mathrm{K}}{T}\right);\\ k_{\mathrm{de}}({\mathrm{H}}_2)& \simeq 2.3\times 10^{-13}{\mathrm{cm}}^3{\mathrm{s}}^{-1}\times {\left(\frac{T}{\mathrm{K}}\right)}^{1/2}\exp \left(-\frac{18000\mathrm{K}}{T+1200\mathrm{K}}\right).\end{array}$$

$$\eta \sim \{\begin{array}{ll}1,& T_{\mathrm{dust}}<25\mathrm{K};\\ 0.6,& 25\mathrm{K}<T_{\mathrm{dust}}<80\mathrm{K};\\ 0.33,& 80\mathrm{K}<T_{\mathrm{dust}}<900\mathrm{K};\\ 0,& T_{\mathrm{dust}}>900\mathrm{K};\end{array}$$

$${\mathrm{Gr}}^{+}+e^{-}\to \mathrm{Gr};\mathrm{Gr}+e^{-}\to {\mathrm{Gr}}^{-}.$$

$${\mathrm{Gr}}^{-}+{\mathrm{X}}^{+}\to X+\mathrm{G};\mathrm{Gr}+X^{+}\to X+{\mathrm{Gr}}^{+}.$$

$$\mathrm{Gr}+h\nu \to {\mathrm{Gr}}^{+}+e^{-};{\mathrm{Gr}}^{-}+h\nu \to \mathrm{Gr}+e^{-}.$$

$${\Lambda }_{\mathrm{dust}}=\sum _{\mathrm{sp}}{n}_{\mathrm{sp}}{\left(\frac{8k_{\mathrm{B}}T}{\pi m_{\mathrm{sp}}}\right)}^{1/2}{\sigma }_{\mathrm{dust}}\alpha \times 2k_{\mathrm{B}}(T-T_{\mathrm{dust}}),$$

$$0=\frac{\mathrm{d}{ϵ}_{\mathrm{dust}}}{\mathrm{d}t}=\max \left[4{\sigma }_{\mathrm{B}}{T}_{\mathrm{ah}}^4(R){\sigma }_{\mathrm{dust}}q(T_{\mathrm{dust}}),\sum _{h\nu }{F}_{\mathrm{eff}}(h\nu )\sigma (h\nu )\right]-4{\sigma }_{\mathrm{B}}{T}_{\mathrm{dust}}^4{\sigma }_{\mathrm{dust}}q(T_{\mathrm{dust}}),$$

$$\beta \simeq \{\begin{array}{ll}{\left[{\tau }_0{\pi }^{1/2}\left(1.2+\frac{\sqrt{\ln {\tau }_0}}{1+{\tau }_0/10^5}\right)\right]}^{-1},& {\tau }_0\gtrsim 1;\\ \frac{1-\exp (-2{\tau }_0)}{2{\tau }_0},& {\tau }_0\lesssim 1.\end{array}$$

Full text

11footnotetext: Princeton University Observatory, Princeton, NJ 08544

1. Introduction

Protostellar/protoplanetary disks (hereafter PPDs) surrounding low-mass T Tauri stars are the birthplaces of planets and have typical lifetimes $\sim 10^{6}-10^{7}\ {\rm yr}$ lifespan (e.g. 1995Natur.373..494Z; 2001ApJ...553L.153H). Along with accretion onto the star, sequestration of mass in planets, and perhaps magnetized disk winds, photoevaporation by hard photons likely contributes to the dispersal of PPDs (1994ApJ...428..654H).

Hard photons in different energy bands experience different microscopic physics and have differing effects on PPDs. Following 2009ApJ...690.1539G, we use the term “far-UV (FUV)” for photon energies $6~{}{\rm eV}<h\nu<13.6~{}{\rm eV}$, “extreme-UV (EUV)” for $13.6~{}{\rm eV}<h\nu<0.1~{}{\rm keV}$, and “X-ray” for $h\nu>0.1~{}{\rm keV}$. While EUV may be blocked by the wind from the disk surface (e.g. 2005MNRAS.358..283A), FUV and X-ray radiation are more penetrating. All of these heat, dissociate, or ionize the gas via a plethora of mechanisms. In order to model photoevaporation of PPDs, therefore, one is required to take the richness of the microphysics into account, as well as its interaction with the hydrodynamics.

Evolving a hydrodynamic system coupled with thermochemistry to (quasi-) steady state could be prohibitively expensive if a large chemical reaction network were included. Past work on PPD photoevaporation has compromised (at least) one of the two aspects: hydrodynamics or thermochemistry. 2006MNRAS.369..216A; 2006MNRAS.369..229A modeled EUV photoevaporation in hydrodynamic simulations with minimum thermochemistry. On the other hand, calculations with detailed thermochemistry usually adopt semi-analytic prescriptions for the wind mass-loss rate rather than simulate multidimensional hydrodynamics e.g. 2008ApJ...683..287G; 2009ApJ...690.1539G (hereafter GH08, GH09). Some recent works conduct hydrodynamic simulations with interpolation tables for gas temperature drawn from hydrostatic scenarios (e.g. 2010MNRAS.401.1415O). More recently, 2012MNRAS.420..562H, 2016MNRAS.463.3616H, and 2017MNRAS.468L.108H have coupled hydrodynamics and thermochemistry in simulations of externally irradiated disks and pre-stellar cores; their code is three-dimensional, but their applications have been confined mostly to simplified geometries (spherical or cylindrical) for easier comparison to semi-analytic work.

This work focuses on a consistent combination of hydrodynamic simulation with a moderate-scale chemical network (24 species, $\sim 10^{2}$ reactions). We include the species and reactions that are relevant to photoevaporation, especially heating and cooling mechanisms. Full hydrodynamic simulations are carried out in 2.5-dimensions (axisymmetry), coupled with radiation, thermodynamics, and chemistry, by solving time-dependent differential equations in every zone throughout the simulation domain. Compared to simulations with interpolation tables for thermochemistry, this approach is able to deal with non-equilibrium processes, as when some chemical and hydrodynamic timescales are comparable. The long-term goal of our exploration is to predict observables, especially emission and absorption-line profiles and strengths of important atomic and molecular species, thereby constraining our wind models and the parameters that go into them (e.g. abundances, dust properties, EUV luminosities). We aim eventually to incorporate MHD processes, and expect that the combination of photoevaporative and magnetic effects will lead to higher mass-loss rates than each process acting alone. The hydrodynamic simulations presented here are first steps toward these goals.

This paper is structured as follows. In §2, we briefly summarize our numerical methods and physical approximations. Additional details concerning our treatment of thermochemical processes are given in the Appendices. §3 introduces the parameter choices underlying our fiducial model. §4 presents the main results of our calculations for this model, and for several other models that differ from the fiducial one in one or more parameters, with the goal of exploring the effects of these parameters on gross properties of the flow, especially the mass-loss rate. In §LABEL:sec:discussions, we discuss the role that different bands of radiation play, and also compare and contrast our results with those of 2009ApJ...690.1539G and 2010MNRAS.401.1415O. §LABEL:sec:conclusion concludes and summarizes the paper.

2. Methods

This section summarizes our methods. The computational scheme for hydrodynamics is first described, followed by our methods for radiative transfer and thermochemistry.

2.1. Hydrodynamics

Our modeling of PPD photoevaporation systems involves full hydrodynamic calculations. We use the grid-based, general-purpose, astrophysical code Athena++ (2016ApJS..225...22W; Stone et al., in preparation) in spherical coordinates $(r,\theta,\phi)$ but neglect all dependence on $\phi$: our simulations are axisymmetric. Magnetic fields are neglected in the present work, although Athena++ is fully capable of MHD (indeed optimized for it). We use the HLLC Riemann solver and van Leer reconstruction with improved order of accuracy using the revised slope limiter (see 2014JCoPh.270..784M). Consistent Multi-fluid Advection (CMA) is used to ensure strict conservation of chemical elements and species (e.g. 2010MNRAS.404....2G).

2.2. Radiative transfer

Absorption processes dominate scattering for most of the radiation that we consider: FUV, EUV, and $1~{}\mathrm{keV}$ X-rays (Draine_book; Verner+etal1996). An exception would be $\text{Ly}\alpha$ photons, which may dominate the FUV luminosity, and whose scattering into nonradial directions helps them to penetrate more deeply into the disk (e.g. 2011ApJ...739...78B). We find, however, that unscattered soft FUV photons penetrate the intermediate layer anyway, and more deeply than $\text{Ly}\alpha$. Like $\text{Ly}\alpha$, these photons dissociate $\mathrm{H_{2}O}$ and $\mathrm{OH}$, which can be important coolants, but not $\mathrm{H_{2}}$ or $\mathrm{CO}$(e.g. 1978ApJ...224..841S). The scattering of harder X-rays can be important for ionization and hence magnetic coupling of the upper layers of the disk (Igea+Glassgold1999; Bai+Goodman2009, e.g.), but we are neglecting magnetic fields here.

Therefore, in this paper, scattering is neglected, and radiative transfer consists only of radial ray tracing, the sources of all hard photons being assumed to lie at the origin ($r=0$). This is facilitated by our choice of spherical coordinates, although our algorithm can trace rays in nonradial directions also (Wang 2017, in preparation).

One ray is assigned to each radial column. Its luminosity is adjusted as it propagates through each cell according to the photoreactions within that cell. Some cells can be individually optically thick. Hence for photochemistry, we adopt as the effective flux at photon frequency $\nu$,

[TABLE]

where $F_{0}$ is the flux impinging on the inner face of the current cell, $\lambda(\nu)$ is the local absorption mean free path of photons at frequency $\nu$, and $\delta l$ is the chord length of the ray across the cell. (For radial ray tracing, $\delta l$ is simply the radial width of the cell.) Eq. (1) yields $F_{\mathrm{eff}}\rightarrow F_{0}$ as $(\delta l/\lambda)\rightarrow 0$.

2.3. Chemistry and Thermodynamics

In each cell, a coupled set of ordinary differential equations (ODEs) is solved to update the abundances of all $\mathcal{N}$ chemical species $\{n^{i}\}$ and internal energy density $\epsilon$. These equations read, nominally,

[TABLE]

in which the terms involving $\{\mathcal{A}^{i}_{\;jk}\}$ describe two-body reactions, while those in $\{\mathcal{B}^{i}_{\;j}\}$ represent photoionization and photodissociation. $\Gamma$ and $\Lambda$ are the heating and cooling rates per unit volume, respectively. $\{\mathcal{A}^{i}_{\;jk}\}$, $\Gamma$ and $\Lambda$ are usually functions of temperature $T$. The thermal energy density $\epsilon=c_{V}(\{n^{i}\})T$, where $c_{V}$ is the heat capacity of the gas at constant volume. (Thermochemistry and hydrodynamics are solved in separate substeps, whence we use $c_{V}$ instead of $c_{P}$ here.) The ODEs (2) are solved in conjunction with the hydrodynamics by operator splitting. That is, they are advanced one time step after each hydrodynamic step, which has included advection of the chemical species, while holding the masses of all elements fixed within each cell. Photoreactions are included using the radiative fluxes computed as described in §2.2. The updated internal energy $\epsilon$ and number densities $\{n^{i}\}$ of all species are then used to initialize the next hydrodynamic step.

The ODEs (2) are usually stiff and hence numerically difficult. We use a standard ODE solver with adaptive implicit modules, CVODE (see hindmarsh2005sundials). The solution of these equations dominates our total computation time, typically by a factor $\gtrsim 10$ compared to the hydrodynamics. Nonetheless, this brute-force approach rewards us by being able to deal with non-equilibrium conditions, as will be discussed later in this paper.

Guided by GH08, GH09, and our own numerical experiments, we adopt 24 species that are most relevant to heating and cooling processes involved in PPD photoevaporation: $\mathrm{e}^{-}$ (free electrons), $\mathrm{H}^{+}$, H, $\mathrm{H_{2}}$, $\mathrm{\mathrm{H_{2}}}^{*}$ (using the $v=6$ vibrational state as a proxy for $\mathrm{H_{2}}$ in all excited states, see Appendix A.1 and TH85), He, $\mathrm{He}^{+}$, O, $\mathrm{O}^{+}$, $\mathrm{O}^{*}$ (the ${}^{1}D$ state of atomic oxygen as a proxy for all neutral excited states, see Appendix A.2), OH, $\mathrm{H_{2}O}$, C, $\mathrm{C}^{+}$, CO, S, $\mathrm{S}^{+}$, Si, $\mathrm{Si}^{+}$, Fe, $\mathrm{Fe}^{+}$, Gr, $\mathrm{Gr}^{+}$, $\mathrm{Gr}^{-}$. Here Gr and Gr± denote neutral and singly-charged dust grains, respectively.

We extract the reactions involving these species from the UMIST astrochemistry database (UMIST2013). However, the interstellar radiation fields and matter densities to which the standard UMIST database is usually applied are rather different from those of PPDs. We therefore exclude all reactions involving photons and dust grains in the UMIST library; instead, we evaluate those reaction rates separately.

Photoionization and photodissociation are critical mechanisms that affect photoevaporation. At each photon energy, the ionization cross section of each atomic species is evaluated using the data in Verner+Yakovlev1995; Verner+etal1996. For molecular species that can react with FUV photons, namely $\mathrm{H_{2}}$, CO, OH, and $\mathrm{H_{2}O}$ here, we adopt the FUV-induced photochemcial reaction rate based on 1985ApJ...291..722T for $\mathrm{H_{2}}$, 2009A&A...503..323V for CO (note that this photodissociation cross section is $\sim 10\times$ the value in TH85), and 2014ApJ...786..135A for $\mathrm{H_{2}O}$ and OH. The photochemical processes related to $\mathrm{H_{2}}$, C and CO may be subject to considerable self-shielding and cross-shielding. Using the radial column density data that are obtained by integrating along radial rays, we evaluate the impact of the self-/cross-shielding by adopting the analytic formulae in 2009A&A...503..323V (for CO) and TH85 (for C), and 1996ApJ...468..269D (for $\mathrm{H_{2}}$). It is worth noting that the FUV-induced processes in parallel with photodissociation of $\mathrm{H_{2}}$, $\mathrm{H_{2}O}$ and OH, e.g. FUV pumping of $\mathrm{H_{2}}$ onto its excited states and its subsequent effects, can have considerable thermodynamic effects. We refer the reader to Appendices A.1 and A.2 for detailed discussion.

Heating and cooling processes are directly associated with chemical reactions. While the amount of energy deposited into and removed from the gas by photoionization and recombination can be estimated straightforwardly (see also Draine_book, eqs. 27.3, 27.23), the thermodynamic effects of other chemical reactions need elaboration, which is provided in Appendices A.1 through A.3. There are other radiative mechanisms that remove energy from the gas, especially collisionally pumped ro-vibrational transitions of molecules, and fine-structure transitions of atoms. We briefly summarize those mechanisms and our method for evaluating them in Appendix A.4.

Dust grains are usually crucial in PPD photoevaporation. Following the arguments in GH08, as well as 2006A&A...459..545G; 2007A&A...476..279G; 2011ApJ...727....2P, we suggest that polycyclic aromatic hydrocarbons (PAHs) overwhelm dust grains of other sizes in terms of the following effects, thanks to their dominant contribution to total dust surface area: photoelectric heating of gas, dust-gas collisional energy transfer, recombination with free electrons, dust-assisted molecular hydrogen formation, and neutralization of positive ions. We include the processes listed above as outlined in Appendix A.3.

3. Choice of Fiducial model

This section presents the setup of our fiducial model, whose main properties are listed in Table 1. Other models, each differing from the fiducial in one parameter, are described in §4.2.

The simulation domain is axisymmetric, extending from $2~{}\,\mbox{\sc au}$ to $100~{}\,\mbox{\sc au}$ in radius ($r$) and [math] to $\pi/2$ in colatitude ($\theta$). All models are presumed to be symmetric about the equatorial plane, so that, for example, quoted mass-loss rates include outflows at $\theta>\pi/2$. All dependence on the azimuthal coordinate ($\phi$) is ignored. Outflow boundary conditions with a radial flow limiter (which inhibits radial inflow) are imposed at $r=2~{}\,\mbox{\sc au}$ and $r=100~{}\,\mbox{\sc au}$, and reflecting boundary conditions at $\theta=0$ and $\theta=\pi/2$. Our standard resolution is $256$ radial by $128$ latitudinal zones, the radial zones being logarithmically spaced, and the latitudinal zones equally spaced.

The gravitational field is that of a $1~{}M_{\odot}$ star located at the origin. The disk, whose self-gravity is neglected, is initialized in hydrostatic and centrifugal balance, except for slight imbalances due to numerical discretization. The disk density and temperature profile follow the steady state solution in 2013MNRAS.435.2610N, in which we set the midplane density as $n=10^{10}~{}{\rm cm}^{-3}$ and temperature $T=20~{}{\rm K}$ at $r=10^{2}~{}\,\mbox{\sc au}$, with radial power index being $(-2.25)$ for density and $(-0.5)$ for temperature–this profile yields a disk mass $\approx 0.03M_{\odot}$ within $100~{}\,\mbox{\sc au}$. The density and temperature profiles roughly agree with GH09, but the latter are not quite hydrostatic.

All radiation emanates from the origin of spherical polar coordinates. Our simulation domain does not cover the origin, and the rays are not attenuated before they reach the inner boundary. The source is isotropic, but those rays that reach the midplane region at the inner boundary are discarded (we also test not discarding those rays, finding negligible differences in the dissociation layer and in the wind). Each ray has four discrete energy bins, representing four important bands of photon energy: $h\nu=7~{}{\rm eV}$ for FUV photons that do not interact appreciably with hydrogen molecules (“soft FUV” hereafter), $12~{}{\rm eV}$ for Lyman-Werner (“LW” for short) band photons, $25~{}{\rm eV}$ for EUV photons, and $1~{}{\rm keV}$ for X-ray photons. $\text{Ly}\alpha$ photons are neglected, as discussed above. The number of photons radiated in each energy bin per unit time follows the luminosity model described in GH08 and GH09: (1) a $9000~{}{\rm K}$ black body spectral profile for FUV ($6~{}{\rm eV}<h\nu<13.6~{}{\rm eV}$) with total luminosity $L_{\mathrm{FUV}}=10^{31.7}~{}{\rm erg}~{}{\rm s}^{-1}$; (2) an additional EUV-photon emission rate 111GH08 and GH09 assumed different EUV luminosities for their fiducial models. Here we adopt the value specified in GH09. $\Phi_{\mathrm{EUV}}=10^{40.7}~{}{\rm s}^{-1}$; (3) and an X-ray luminosity $L_{X}=10^{30.4}~{}{\rm erg}~{}{\rm s}^{-1}$.

The initial elemental abundances are determined by the values in Table 1 (wherein $n_{\mathrm{H}}$ is the number density of hydrogen nuclei). These choices generally follow a subset of those in GH08, with the additional assumption that elements appear in chemical compounds if possible. These initial abundances are uniform throughout the simulation domain.

Our assumptions about the dust turn out to be important for our results. GH08 and GH09 treated two populations of grains: (i) an MRN-type power-law distribution with a minimum grain radius of $50~{}\mathrm{\AA}$, maximum of $20~{}{\rm\mu m}$, and a dust-to-gas ratio by mass of $10^{-2}$; and (ii) PAH grains with abundance $8.4\times 10^{-8}$ per hydrogen nucleus. The first population has a total geometrical cross section of $2\times 10^{-22}~{}{\rm cm}^{2}$ per hydrogen nucleus ($\sigma_{\rm dust}/\mathrm{H}$). The authors do not state the radius of their PAH grains explicitly, but they refer to 2001ApJ...554..778L, and we interpret this to mean that their PAHs can be approximated by spheres of radius $6~{}\mathrm{\AA}$. It would follow that the contribution of their PAHs to the cross section is $\sigma_{\rm dust}/\mathrm{H}\approx 9.5\times 10^{-22}$, i.e. several times larger than that of their MRN population, although the contribution to the dust-to-gas mass ratio is only $\sim 10^{-4}$. As noted above, A.3, the principal effects of dust, especially heating and absorption of radiation, are expected to be dominated by the smallest grains—PAHs.

For simplicity, we prefer to work with a single-sized grain population. We therefore neglect MRN grains and take the approximate relative abundance for our PAH-like grain species (Gr) as $10^{-7}$ per hydrogen atom, slightly greater than that of GH08, and a PAH radius of $5~{}\mathrm{\AA}$, i.e. approximately $60$ carbon atoms per PAH: see 2001ApJS..134..263W). The dust-to-gas mass ratio is then $0.7\times 10^{-4}$, and $\sigma_{\rm dust}/\mathrm{H}=8\times 10^{22}~{}{\rm cm}^{2}$.

Although variable dust abundance is fully allowed by our code, for the sake of simplicity we set the relative abundance of Gr to be uniform and assume that the dust comoves with the gas.

We run the simulation for $1.2\times 10^{4}~{}{\rm yr}$ with microphysics enabled but the central radiation sources turned off until the disk structure is fully numerically relaxed, and the temperature profile converges to that set by the artificial heating profile (§A.3.3, which is sufficiently close to the initial profile. The chemical abundances do not change during this relaxation process except by passive advection. We confirm after this process that the disk is indeed in hydrodynamic equilibrium and has no outflow. Then, at $t=3.6\times 10^{3}~{}{\rm yr}$, irradiation is turned on and remains on for the rest of simulation; this lasts $\gtrsim 500~{}{\rm yr}$, sufficiently long compared to the radial flow timescale $\tau\sim(100~{}\,\mbox{\sc au})/(30~{}{\rm km\ s^{-1}})\approx 16~{}{\rm yr}$ so as to reach an approximate quasi-steady state. On Princeton University’s local computer cluster perseus, $500~{}{\rm yr}$ of simulated time takes $\sim 100~{}\mathrm{\ hrs}$ of wall-clock time on 128 CPUs. About 95 per cent of the time is consumed by the thermochemical calculations for the fiducial model, the hydrodynamic and ray-tracing steps being relatively quick.

We also calculate several models that differ from the fiducial in one or more parameters, as described in §4.2.

4. Results

In this section, we first present and elaborate the fiducial simulation (see §3), then compare the the results of the variant models shown in 2 (see §4.2).

4.1. Fiducial Model

Fig. 1 displays meridional plots of the structure of our fiducial model averaged over the final $100~{}{\rm yr}$ of the simulation. The white curves shown in the top row of panels are streamlines, the integral curves of the vector field $\rho\boldsymbol{v}_{p}$ ($\boldsymbol{v}_{p}$ is the poloidal velocity), spaced by constant mass-loss rate $10^{-10}~{}M_{\odot}~{}{\rm yr}^{-1}$: that is to say, this is the mass flux between neighboring streamlines when integrated over azimuth and multiplied by two to include the reflection of the computational region below the equatorial plane. Streamlines that meet the outer boundary with a negative value of the Bernoulli parameter

[TABLE]

are not plotted, and the outflow along such streamlines is omitted from the computation of the total mass-loss rate. Here $v$ is the magnitude of fluid velocity vector, $p$ the gas pressure, $\gamma\approx 5/3$ the adiabatic index, and $\Phi$ the gravitational potential. With this mask we get rid of (very slow) radial flows near the mid-plane: since the density there is six orders of magnitude higher than the wind, a tiny radial velocity fluctuation could otherwise give a spurious contribution to the mass-loss rate. As displayed in Fig. 1, the streamlines terminate on the disk at the surface where $\mathcal{B}$ becomes negative. We consider this surface to be the base of the wind. (As discussed in §LABEL:sec:GH09_comparison, this definition of the wind base differs from that of GH09.)

Fig. 2 shows several flow variables along two representative streamlines originating from cylindrical radii $R\equiv r\sin\theta=5~{}\,\mbox{\sc au}$ and $R=15~{}\,\mbox{\sc au}$.

The density and temperature profiles shown by Fig. 1 can be divided into three relatively distinct regions:

• •

Midplane layer: $0<(z/R)\lesssim 0.3$ ($R=r\sin\theta$ being cylindrical radius), $T<10^{2}{\rm K}$. The structure here is basically unchanged from the initial conditions.

• •

Intermediate layer: $0.3\lesssim(z/R)\lesssim 0.6$, $10^{-19}~{}{\rm g}~{}{\rm cm}^{-3}\lesssim\rho\lesssim 10^{-16}~{}{\rm g}~{}{\rm cm}^{-3}$, $10^{2}\lesssim T\lesssim 10^{3}~{}{\rm K}$. The total mass in this layer is $\sim 10^{-6}~{}M_{\odot}$. EUV photons scarcely penetrate this region, whose properties are controlled by FUV and X-ray processes: photodissociation and photoelectric heating, as well as radiative cooling by collisionally excited molecular and/or atomic transitions. Most $\mathrm{H_{2}}$ molecules and a lot of $\mathrm{CO}$ molecules survive in this region because of significant self- and cross-shielding of Lyman-Werner photons. Soft FUV photons that do not interact much with molecular hydrogen are relatively unshielded and pervade the intermediate layer, photodissociating $\mathrm{OH}$ and $\mathrm{H_{2}O}$), penetrating to the bottom of the layer, or escaping through the outer radial boundary.

• •

Wind layer: $(z/R)\gtrsim 0.6$, $\rho\lesssim 10^{-19}~{}{\rm g}~{}{\rm cm}^{-3}$, $T\gtrsim 10^{4}~{}{\rm K}$. This region is filled with mostly ionized gas, flowing outwards at radial velocity $v_{r}\sim 30~{}{\rm km\ s^{-1}}$. Photoionization heating and adiabatic expansion dominate the thermodynamics of this region.

If we integrate the $\mathcal{B}$-masked radial mass flux at the the $r=100~{}\,\mbox{\sc au}$ boundary (and its reflection at $z<0$) and average over the last $100~{}{\rm yr}$ of our fiducial run, we obtain a total mass-loss rate $\dot{M}_{\mathrm{wind}}\simeq 3.4\times 10^{-9}~{}M_{\odot}~{}{\rm yr}^{-1}$, corresponding to a disk dispersal timescale $\sim 10^{7}~{}{\rm yr}$. The mass-loss rate is lower than that of GH09 (see §LABEL:sec:GH09_comparison for further discussion). However, our mass-loss rate undergoes significant fluctuations, and is uncertain to at least $\sim 10$ per cent. Fig. 3 plots the mass-loss rate for the last $200~{}{\rm yr}$ of the (lower resolution) fiducial run (Model 0). They correlate with what appears to be a thermal instability of the outer disk, whereby it swells vertically, intercepts more radiation, and then swells further but also migrates at a few ${\rm km\ s^{-1}}$ through the outer boundary, temporarily increasing $\dot{M}$. This behavior is smoothed over by the time averages used to make Fig. 1. These swellings, being slower and denser than the general wind, partly shield themselves against photodissociation of some molecules, especially $\mathrm{H_{2}}$ and $\mathrm{CO}$, so that those molecules survive farther into the outflow than they would otherwise.

Even outside these swellings, there are also molecules surviving in regions with rather high temperature ($\sim 10^{3}~{}{\rm K}$, or even up to $\sim 10^{4}~{}{\rm K}$). $\mathrm{H_{2}O}$ and $\mathrm{OH}$ molecules exist at the surface of the intermediate layer, detached from the midplane (last panel of Fig. 1). The reformation rates of $\mathrm{H_{2}O}$ and $\mathrm{OH}$ are comparable to photodissociation at that surface. At the cooler temperatures below it, inside the intermediate layer, reformation is less efficient but photodissociating FUV is still present. The wind region, on the other hand, does not have sufficient $\mathrm{H_{2}}$ (reactions that are most efficient in forming $\mathrm{H_{2}O}$ and $\mathrm{OH}$ need $\mathrm{H_{2}}$ as reactants, while the reactions that convert atomic H to OH and $\mathrm{H_{2}O}$ are very slow).

In Fig. 4, we plot the distribution of $\mathrm{CO}$, $\mathrm{OH}$ and $\mathrm{H_{2}O}$ in the wind region and intermediate layer, in the plane of by $\log_{10}T$ and $v_{r}$. (The temperature $T$ here represents the kinetic temperature of the local (mostly $\mathrm{H}$ & $\mathrm{He}$) gas, not the vibrational or even rotational excitation temperature of the molecules.) For those molecules, a tail on the high temperature ($T\sim 10^{3}\mbox{--}10^{4}~{}{\rm K}$) and intermediate radial velocity ($v_{r}\sim 5\mbox{--}10~{}{\rm km\ s^{-1}}$) end of the 2-D distribution indicate their survival at the bottom in the wind region. Such hot molecular gas would be less prominent had we assumed local thermochemical equilibrium. For the luminosity in the LW and EUV bands of our fiducial model, it can be estimated that the timescale of CO photodissociation is $\sim 0.1~{}{\rm yr}$ at $r\sim 10~{}\,\mbox{\sc au}$. Given the speed of photoevaporative outflow, this timescale is sufficient for some CO to survive $\sim 0.1-1~{}\,\mbox{\sc au}$ into the hot wind. These timescales are sensitive to radial distance (from the radiation sources), to the way photodissociation is modeled (see §2.3), and to the LW and EUV band luminosity. Observational constraints on such molecules could be an important check on these models, and might diagnose the role of UV in driving PPD winds.

4.2. Exploring the Parameter Space

To explore the effects of our input parameters, we have run a number of additional simulations, most differing from the fiducial run in one parameter. These models and some synoptic results are listed in Table 2. We discuss some of these models here, and others in §§LABEL:sec:GH09_comparison-LABEL:sec:compare-oeca10 in relation to the works by GH09 and OECA10.

In the fiducial model the luminosity in the Lyman-Werner band is tiny compared to that in soft ($h\nu<11.3~{}{\rm eV}$) FUV photons: around 0.35 per cent, using the $9000~{}{\rm K}$ black body SED. However, as observed by e.g. 2000ApJ...544..927G, the SED for FUV radiation is rather variable from object to object and often more luminous in the LW band than the black-body model adopted by GH08 and GH09. Hence we include a series of models, with 0 and 100 times the fiducial luminosity in the Lyman-Werner band, to cover this uncertainty. We also test 0 and 10 times EUV or X-ray luminosity to diagnose the impact of those photons that can ionize atomic hydrogen.

For very small grains such as our PAHs, the grain absorption cross section for FUV and EUV photons depends on total grain mass rather than grain area. We have a much smaller grain mass than GH09. Model 9 in Table 2 has double the dust radius ($r_{\mathrm{dust}}=10~{}\mathrm{\AA}$) and therefore eight times the dust mass at the same relative number density ($n_{\mathrm{Gr}}/n_{\mathrm{H}}=10^{-7}$).

To test our truncation errors, we repeat the fiducial run at resolution $128\times 64$, i.e. coarser by $\times 2$ in both latitude and radius. This convergence test is run for much longer time period ($\sim 2000~{}{\rm yr}$) to better characterize fluctuations around the mean state.

Fig. 5 illustrates the hydrodynamic structure of a few representative models. These plots are based on time averages over $100~{}{\rm yr}$, so that the flow field is in approximate steady state. In the runs with $100\times$ LW photons (Model 2) or $10\times$ X-ray photons (Model 8), a thick neutral atomic layer exists at the top of the intermediate layer. In this layer, the temperature and sound speed reach a local maximum with respect to height or latitude, and significant outflows may occur. This causes the jagged shape of the sonic curves in the third and fourth panels.

Figs. LABEL:fig:heating and LABEL:fig:cooling show the vertical distributions of heating and cooling mechanisms at $R=15~{}\,\mbox{\sc au}$, a typical location where the outflow streamlines originate. The three layer structure (§4.1) is obvious in most of the models. Details of those structures vary with model parameters, with implications for the mechanisms responsible.

The panels of Fig. LABEL:fig:heating convey some general impressions about the heating mechanisms. The vertical heating profile usually has two peaks: one at the bottom of the intermediate layer, the other at the top of it. Photoionization heating by the harder (EUV and X-ray) photons dominates, unless these photons are absent or are overwhelmed by photons in other bands (e.g. Model 2, $100\times$ LW photons; see discussions below). On the cooling side (Fig. LABEL:fig:cooling), the $\mathrm{OH}$ and/or $\mathrm{H_{2}O}$ ro-vibrational transitions and S I $25~{}{\rm\mu m}$ transition dominate at the bottom of intermediate layer, the Si II $35~{}{\rm\mu m}$ and O I $63~{}{\rm\mu m}$ transitions in the middle of that layer, and $\mathrm{H_{2}}$ ro-vibrational cooling near the top. In the “wind” region, cooling and heating are dominated by recombination and the photoelectric effect. Using the integrated cooling rate, we have estimated some of the important line luminosities (Table LABEL:table:line-luminosity).