Anatomy of Photoevaporation Base: Linking the Property of the Launched Wind to Irradiation Flux

Ultraviolet and X-rays from radiation sources disperse nearby gas clumps by driving winds due to heating associated with the photochemical processes. This dispersal process, photoevaporation, constrains the lifetimes of the parental bodies of stars and planets. To understand this process in a universal picture, we develop an analytical model that describes the fundamental physics in the vicinity of the wind-launching region. The model explicitly links the density and velocity of photoevaporative winds at the launch points to the radiation flux reaching the wind-launching base, using a jump condition. The model gives a natural boundary condition for the wind-emanating points. We compare the analytical model with the results of radiation hydrodynamic simulations, where a protoplanetary disk is irradiated by the stellar extreme-ultraviolet, and confirm good agreement of the base density and velocity, and radial profiles of mass-loss rates. We expect that our analytical model is applicable to other objects subject to photoevaporation not only by extreme-ultraviolet but by far-ultraviolet/X-rays with suitable modifications. Future self-consistent numerical studies can test the applicability.


Introduction
In a variety of astrophysical environments, we find a situation where a clump of gas is exposed to ultraviolet (UV) and X-ray radiation fields. For instance, interstellar medium (ISM) and molecular clouds are irradiated by the ultraviolet (UV) radiation from, if present, nearby massive stars (e.g., Kahn 1954;Oort et al. 1955;Goldsworthy 1958Goldsworthy , 1961Axford 1961;Elmergreen 1976;Reipurth 1983;Bertoldi 1989;Bertoldi & McKee 1990;Lefloch & Lazareff 1994;O'dell & Handron 1996). Protoplanetary disks (PPDs) and the atmospheres of planets are also subject to the UV and X-rays from the host star (e.g., Watson et al. 1981;Hollenbach et al. 1994;Shu et al. 1994;Lammer et al. 2003;Yelle 2004;Murray-Clay et al. 2009). At high redshifts, strong ionizing sources such as Pop III stars and quasars illuminate intergalactic medium and (mini)halos (e.g., Shapiro 1986;Donahue & Shull 1987;Shapiro & Giroux 1987;Shapiro et al. 2004;Iliev et al. 2005;Park et al. 2016;Nakatani et al. 2020). One of the important outcomes from the irradiation is heating associated with the photochemical processes, which can result in driving winds. This wind-driving process, termed photoevaporation, has affected a wide variety of objects through the history of the universe. Irradiated objects lose mass via photoevaporative winds. It limits the available time for star/planet formation and determines the fate of the parental bodies.
Since the underlying physics is overall similar regardless of irradiated objects or radiation sources, it is expected possible to develop a general theory of photoevaporation by understanding the basics. In most of relevant studies, however, the main focus has been placed more on resulting quantities due to photoevaporation (e.g., mass-loss rates, wind velocity, and lifetimes of a specific object) than on the physical process itself. Thus, there is room for further investigation into the fundamentals. It still remains unclear, for example, how the photon energy relates to the velocity of launched winds, what quantities determines the mass flux of photoevaporative winds at the wind-emanating points, what conditions set the escapability of photoheated gas against gravity, etc.
To address these issues, as the first step, we study the windlaunching region (base) to understand the relation between irradiation and photoevaporative winds. We primarily focus on the physical link between the irradiation flux and the base property (namely, launch density and velocity) in the present study. To this end, we develop an analytical model describing the basic physics operating in the vicinity of the windlaunching region in a first-principles approach.
The paper is organized as follows. In Section 2, we first introduce previous analytical models, and then we present our newly constructed model. We compare the analytical model with self-consistent radiation hydrodynamic simulations of photoevaporating PPDs in Section 3 to test the applicability of the model. We present a discussion in Section 4, and summarize and conclude in Section 5.

Analytical Model
Our analytical model for photoevaporative winds builds on the classical theories for the ionization front propagation in the ISM. We therefore first introduce the representatives of such models in Section 2.1. Then, we present our analytical model in Section 2.2.

Preliminaries
Suppose the situation where ionizing radiation is incident on neutral medium (H I region) from one side (Figure 1). It forms a sharp ionization front that advances into the neutral gas while leaving the ionized medium behind it (H II region). The neutral gas ahead of the ionization front has density ρ 1 , pressure p 1 , and velocity v 1 in the rest frame of the advancing front. The neutral gas is immediately altered to an ionized gas with density ρ 2 and pressure p 2 . The ionized gas leaves the front at a relative velocity of v 2 .
By approximating the ionization front to a discontinuous boundary, Kahn (1954) introduced a jump condition using the conservation laws of mass, momentum, and energy, where h 1 and h 2 are the specific enthalpy for the neutral and ionized gas, respectively, and F is the total energy flux impinging on the ionization front. The total energy flux F is expressed as ò n = n n ¥ F Fd th with the specific energy flux F ν and the ionization potential hν th . The total energy flux F is not the intrinsic flux but only counts the excessive photon energy beyond the ionization potential (see also below). Here, the unknown variables are the three quantities in the ionized region (indicated by the subscript "2"), and the others are model parameters.
As one ionizing photon is responsible for one ionization, 3 the total number flux of ionizing photons J is equal to the flux number of atomic hydrogen flowing into the front, where n i is number density of neutral or ionized gas. The mass density ρ i and number density n i are related via gas mass per particle, m, as ρ i = mn i . The specific energy flux relates to the specific photon number flux J ν as F ν = h(ν − ν th )J ν . Equations (1)-(3) reduce to a quadratic equation with respect to ρ 2 /ρ 1 , and the roots are , which correspond to either The ionization front is classified according to v 1 (or 1  ) and the magnitude of ρ 1 /ρ 2 . 4 We do not go into the further details about the classification because it is not very relevant to the present study.
The above model lets the ionized gas have an arbitrarily high temperature if a high F is given. In practice, however, the temperature increase is limited by efficient cooling processes like H II radiative recombination cooling and other metal line cooling. An alternate jump condition that can avoid this issue (Axford 1961;Spitzer 1978). The effects of coolants are effectively taken into account by dropping the energy equation but, as in the jump condition for isothermal shock, giving the isothermal sound of speed ahead of and behind the front (c 1 and c 2 , respectively) as model parameters, The isothermal sound speeds are determined by the equilibrium temperatures in this model (usually, c 1 /c 2 = 1 is assumed). Physical quantities are defined in the rest frame of the advancing ionization front. Neutral gas in the upstream region flows into the front at velocity v 1 and with density ρ 1 and pressure p 1 . The gas is ionized while passing through the front and gets v 2 , ρ 2 , and p 2 . Ionizing photons reach the boundary with the photon number flux J.
Note that the unknown variables are ρ 2 and v 2 here. Solving the above equations with respect to the density ratio ρ 2 /ρ 1 gives roots similar to that of Equation (4), and v R correspond to D  and R  in Kahn's model. The ionization front is classified according to v 1 and the signs in Equation (9). We refer the readers to the original papers for the classifications.

Analytic Model for Photoevaporative Winds
The previous models are developed for the problems of H II regions where strong radiation sources form a sharp ionization front advancing into H I clouds. This picture is qualitatively different from that of steadily photoevaporating objects, in which winds emerge from the surface of a static/stationary bulk region. The surface does not change its location significantly while the object retains a sufficient mass. Hence, to apply the previous models to the problems of photoevaporation, suitable modifications are needed to the laws of conservation. We make such modifications to develop a similar model available for photoevaporating objects in this section.

Model Assumptions
Photoevaporation is a process that net deposited energy is distributed to thermal and kinetic energy of the gas. It indicates the importance of modeling the energy partition. We therefore use the conservation of energy as well as those of mass, momentum, as in Kahn's model.
The relation J = n 1 v 1 used in the previous models originates from the law of conservation for the number of the ionized hydrogen across the ionization front. This relation is obtained by assuming that the downstream region is fully ionized while the upstream region is fully neutral, and that the integrated recombination rate within the ionization front is ignored (see also Appendix for reference). The integrated recombination rate can be ignored when it is dominated by the advection term n 2 v 2 . It is always justified in the previous models for the assumption of an infinitesimally thin ionization front. However, this is not generally the case for photoevaporating objects, since the integrated recombination term can be essential to form a stationary ionization front. We thus do not use the relation J = n 1 v 1 in our model. Alternatively, we treat v 1 as another unknown variable in addition to ρ 2 and v 2 . The other variables (ρ 1 , c 1 , and c 2 ) are given as model parameters.
Dropping the J = n 1 v 1 relation corresponds to excluding the conservation law of the chemical species that is the main absorber of the heating photons, from the determination of the unknown variables. This approach lets us construct a model without any other chemistry-related unknown variables or model parameters, namely, the degree of ionization in the upstream/ downstream regions and the integrated recombination rate within the ionization front. Such unknown variables or model parameters are unfavorable for the simplicity of the model. This approach is also advantageous to explicitly describe and understand the relation between the injected energy and the resulting wind property, since the fluid's velocity is derived via energy conservation rather than chemical conservation.
In our model, we introduce a geometrically thin, finite-width processing layer where the bulk of the photoenergy is absorbed (the cyan layer in Figure 2). Having this layer is more realistic than simply assuming a discontinuous front and is essential for a standing photoevaporation front. The thin-layer assumption lets us ignore any geometrical effects, such as gas expansion transverse to the layer's plane, the differences between the upstream and downstream regions in the gravitational potential, and the angular momentum of the gas. Hereafter, we also refer to this layer as the photoevaporation front, photoevaporation base, or, simply, base.
Assuming a finite-width processing layer results in introducing a term of the integrated cooling rate in the energy conservation Equation (3). This term brings another parameter, which is unknown a priori, to our model and needs to be taken into account for the general applicability. Nevertheless, in the A dense gas clump is illuminated by the radiation field that excites photoevaporative winds from the surface of the cold gas. The cold gas can be a gravitationally bound object. The photoheated layer is optically thin to the heating photons. (Bottom panel) Zoomed-in view near the boundary between the photoheated and cold regions. We consider a processing layer where the cold gas is converted into a hot, photoheated gas. The bulk of the excess energy is absorbed in this processing layer and goes to gas' thermal and kinetic energy. The processing layer's width is assumed sufficiently geometrically thin and is of the order of s á ñn 1 1 ( ) . In the present study, we also refer to this processing layer as the photoevaporation base, and we call v 2 the launch velocity. present study, we ignore this term to avoid complications. This means that we consider the cases where the bulk of photoenergy goes into the gas kinetic and thermal energy while only a negligible fraction is lost by radiative cooling. Constructing a model that covers the cooling-dominated cases should be addressed in future works.
Neglecting the integrated cooling term implies that our model is applicable only to a partial region in the parameter space spanned by the photoevaporating object's property and ionizing radiation strength. Using a generalized analytical model for steadily evaporating clouds, Bertoldi (1989) quantitatively showed that the assumptions of a very thin front and ionization/thermal equilibrium in the H II region are justified when the ionizing radiation is sufficiently strong, i.e., when the photoheating and ionization timescales are shorter than the dynamical timescale. Cooling is not negligible in this case. We consider that this conclusion holds for other photoevaporating objects. In fact, in the context of planetary atmospheric escape, it has been known that a high ionizing radiation forms a sharp ionization front, and balance between photoionization and recombination gives the estimation for the density at the ionization front (so-called recombination-limited escape; Murray-Clay et al. 2009). Thus, our assumption appears to break for strongly irradiated objects. An alternative model is necessary for this case.
At a lower irradiating flux, there is a regime where the deposited energy is mainly processed via adiabatic expansion within the ionization front rather than being radiatively lost. This regime is similar to so-called energy-limited escape of planetary atmosphere (Watson et al. 1981;Murray-Clay et al. 2009) but not exactly the same; the ionization front is supposed to be geometrically thick in the energy-limited escape, while our model assumes a thin front. Our model assumptions are likely met when the dynamical timescale is longer than the ionization timescale while being shorter than the photoheating timescale. Thus, the application range depends on the irradiation spectrum in addition to the object property and radiation strength. We expect our model to work the most for low deposited energy per ionization.

Conservation Laws
Since our interest is in photoevaporation from a standing photoevaporation front, we also set ansatz that the gas is subsonic ahead of the front (v 1 = c 1 ). To summarize, the adopted model assumptions are: where g is the external conservative body force vector, ds is a line element vector normal to the front, Δh ≡ h 2 − h 1 is the difference in enthalpy between the photoheated and cold layers, and F is the photoheating flux reaching the processing layer. In our model, ρ 2 , v 2 , and v 1 are unknown variables, while ρ 1 , c 1 , and c 2 are given parameters. The photoheating flux may be slightly attenuated up to the photoheated layer and not exactly equal to the intrinsic photoheating flux of the radiation source. Considering a finite-width front naturally leads us to incorporate work done by the external body force in the momentum conservation equation. The integration on the righthand side is taken within the processing layer with a width of the order of s á ñn 1 1 ( ) . 5 Thus, where sˆis a unit vector for s, and Δτ (≈1) is the optical thickness of the processing layer. Here, we have assumed that g s ·ˆcan be approximated to constant within the processing layer ( This relation holds regardless of the density profile within the layer, which is usually a monotonically increasing function from the photoheated layer to the cold layer. Since we are interested in decelerating forces here, we assume < 0  in this paper (note, however, that our model holds despite the sign of ). Since we have not set restrictions to the photoheating source other than that the flux is attenuated as dF = − n〈σ〉Fds, this model would also be applicable to FUV-and X-ray-driven winds as far as the model assumptions are met. The base is not necessarily an ionization front. It can be, for instance, an H 2 dissociation front, the photosphere of a small-sized dust disk, etc.
The momentum equation requires as a necessary condition for v 2 to be real. The necessary condition also sets a lower limit to the cold layer density ρ 1 . For a stratified cold gas where the pressure gradient balances with the external body force, this expression is approximately rewritten to The left-hand-side is the processing layer's width, and the right-hand side is the pressure scale length of the stratified gas just ahead of the base. Equation (17) can also be interpreted as the maximum mean pressure gradient across the base s á ñp n 1 1 1 ( ) needs to be larger than the pressure gradient of the cold layer in the hydrostatic equilibrium. If a system breaks Equation (17), it has a processing layer not sufficiently thin to ignore geometrical effects, in conflict with the model assumption (a). Thus, Equation (17) sets a limit to the applicability of the model. Equations The actual width is thicker than s á ñn 1 1 ( ) by several times since the density is lower than n 1 in the processing layer.
where M º v c 2 2 2 and a dimensionless parameter This parameter is positive owing to the requirement of Equation (16). Hereafter, we refer to f as pressure photoheating parameter. 6 In Equation (18), we have used an approximation , where c p is the specific heat capacity at constant pressure in the photoheated layer (c p = 5/2 for a monatomically ideal gas, and c p = 7/2 for a diatomically ideal gas). The discriminant of Equation (18) (18) has therefore only one real root in M 0 2  for any f or c p . This is also evident from Equation (18) by looking at f as a function of The discriminant D < 0 guarantees a positive value within the square root. The value within the absolute-value sign is positive for 3c p /2 > f 2 ,otherwise, negative or zero. Note that the solution is independent of M 1 . The cold gas velocity contributes to Equation (20) as higher-order terms (the lowest order is supposed to be M 1 2 ( )  ). We plot M 2 as a function of f in Figure 3.
Once M 2 is determined, the other unknowns are immediately computed Equation (21) sets the upper limit of c 1 /2c 2 to M 1 , validating the assumption of M < 1 1 in Equations (13) and (14) and selfconsistency of our model. The upper limit of the density ratio ρ 2 /ρ 1 is c c 1 2 2 2 . The density of cooler, accelerating gas is larger than that of the evaporating gas by a factor of c c 2 2 1 2 at least. The launch velocity in the reference frame is thus almost equal to that in the rest frame of the base, i.e., (v 2 − v 1 ) ≈ v 2 . In the limit of M  1 2 , M 2 and M 1 are approximately proportional to F. In this case, the density ratio is essentially independent of F, meaning an approximate thermal pressure balance. Therefore, the effective adiabatic index is nearly zero within the processing layer.
We can also derive an approximate root of Equation (18)  . This relation does not include information on the cold gas. This implies that the cold gas works as a nearly static wall (v 1 ≈ 0) for the expanding gas. The processing layer gas expands as if expansion is allowed only on the photoheated-layer side. Equation (13) relates the momentum flux injected into the photoheated layer (rv 2 2 ) to the pressure gradient and the external force. As for the pressure photoheating parameter f, it characterizes to which of Δh or v 2 2 2 the photoheating energy mainly goes (see Equation (14)). The launch velocity is enthalpy-flux limited for f < (c p + 1/2)/2. In this regime, the left-hand side of the cubic equation, Equation (18), is relatively dominant over the right-hand side in the determination of v 2 . The injected energy flux F largely goes into enthalpy flux ρ 2 v 2 Δh, and therefore v 2 is set by the expansion velocity. The resulting v 2 is accordingly subsonic ). In contrast, for f > (c p + 1/2)/2, the launch velocity is dynamics limited (or kinetic-energy-flux limited). The contribution of the momentum conservation is relatively strong to determine v 2 (see Equation (18)). A large portion of the injected energy flux F is converted to the kinetic energy flux r v 2 2 2 3 in terms of the momentum conservation. The launch velocity v 2 is supersonic and originates from acceleration by the pressure gradient force rather than thermal expansion due to photoheating. The density ratio is correspondingly stronger in this regime (Equation (22)).
Equation (20) appears as if an arbitrarily high v 2 would be achievable by injecting a very high F into the base. However, the launch velocity must be determined in a manner that is consistent with the large-scale structure of photoevaporative winds. For instance, in a spherically symmetric isothermal system, the wind structure will be that of the Parker wind. The velocity profile is uniquely set as a function of distance from the center, according to the boundary condition (or the value of the conservative quantity along the streamlines). The launch velocity and the base location need to be consistent with this velocity profile. A similar discussion can be made for other configurations, e.g., disk winds. Self-similar disk wind models explain that the power-law index of the base density profile limits an allowable range of v 2 in order for the streamline to reach infinity (Clarke & Alexander 2016;Sellek et al. 2021). Self-similar thermal winds are known to generally take the maximum of the allowed M 2 at the base, otherwise, the whole R-z plane cannot be filled with the self-similar winds (Sellek et al. 2021). In this way, the large-scale profiles of the disk winds also matter in determining the base property.
Our model originally has six parameters (ρ 1 , ρ 2 , v 1 , v 2 , c 1 , c 2 ) and imposes three conditions (Equations (12)-(14)). Thus, the solutions have three degrees of freedom. The large-scale profiles and other external information, which we have not incorporated in our model (e.g., the conservation law of chemical species, large-scale profiles), may be responsible for further reducing the degree of freedom to determine the specific values of the parameters uniquely.

Model Limitations and Caveats
The assumptions (a)-(c) will be the basic premises the systems need to meet in using our model. The processing layer can be geometrically thin (Assumption (a)) when heating is sufficiently rapid compared to the dynamical timescale of the gas. Radiative cooling can be neglected (Assumption (b)) if mechanical cooling dominates loss of the internal energy. It requires a shorter gas expansion timescale than radiative cooling timescale. The ratio of these timescales depends on the processing layerʼs width and the main radiative cooling source (s). We expect that since a higher base density is achieved for strong irradiation, radiative cooling would dominate over mechanical cooling as the irradiating flux increases. Hence, qualitatively, the flux needs to be intermediate not to break Assumptions (a) or (b). Quantifying this criterion is an issue dependent on objectʼs property. We will address it in future works.
We have assumed a subsonic upstream region (Assumption (c)) to ignore M 1 2 terms in the momentum and energy conservation (Equations (13) and (14)). It is valid if the upstream region is static or stationary, but there are cases when this assumption breaks down. For example, in externally UVirradiated clouds or disks, the ionization front develops within FUV-excited photoevaporative flows (e.g., Johnstone et al. 1998;Richling & Yorke 2000). Ignoring the M 1 2 terms is not justified there.
Another caveat is that we have modeled the region nearby the base without taking into account the large-scale flow structure. As discussed above, streamlines causally connect from the base to infinity (for escaping winds). The launch velocity v 2 and other quantities at the base must be consistent with the large-scale profiles, as mentioned in the previous section. The significance of our model is that it reduces degrees of freedom for the problem. It is rephrased as the model sets boundary conditions the base needs to meet.
Equations (12)- (14) build on the premise that photoenergy can freely go into the kinetic energy of the gas. In general, hydrostatic solution (v 1 = v 2 = 0) is possible even under the effect of irradiation. It would be realized, e.g., when the gravity is sufficiently strong to bind the photoheated gas. The deposited photoenergy must be drained by radiative cooling in such cases. This is also understood from Equations (12)-(14) since F = 0 is necessary with v 1 = v 2 = 0. It implies that radiative cooling cannot be ignored in this wind-prohibited limit.
In order for our model to predict the mass-loss rates, the downstream gas is required to have kinetic plus thermal energy exceeding net gravitational binding energy. The predictability is lost in a gravitationally bound region. There, when irradiation is strong, the photoheated gas would form a nearly hydrostatic atmosphere, and the deposited photoenergy is mostly lost through radiative cooling in a similar manner to the wind-prohibited limit above. When radiation is weak and cooling is negligible, the mass-loss rates would be determined by the photoheating rate and the binding energy per unit mass, as in the energy-limited escape of planetary atmosphere (Murray-Clay et al. 2009). A continuous energy injection is needed for the flowing-out gas to escape from gravitational binding in this limit.
We have ignored cooling effects within the processing layer in Equation (14), assuming that they are insignificant compared to photoheating (or adiabatic cooling). The net effect of coolants is reducing F by the total cooling rate integrated over the processing layer. Incorporating this effect brings another unknown parameter to the model, as mentioned in Section 2.2.1. If the correcting term is suitably included, our model is available by replacing F with an effective photoheating flux F eff . This parameter must be F eff > ρ 1 v 1 Δh to yield v 2 > 0; that is, the photoheating term must be larger than the cooling term by at least ρ 1 v 1 Δh, given the presence of excited winds.
The gravitational term has not appeared in Equation (14) owing to the assumption of a thin processing layer. For the systems where this assumption breaks, one needs to consider the effects of gravity on F eff as well as cooling effects.
Regardless of the limitations, we have obtained these findings through our model: (i) M 2 is characterized by f and c p and is a monotonically increasing function with respect to f (Figure 3), and (ii) the launch velocity is categorized into the enthalpy-flux-and dynamics-limited regimes.
We stress that the purpose of this study is to understand the physical link between irradiation and excited photoevaporative winds. Our analytical model successfully describes the fluidʼs response to the injected radiation energy into an almost static layer within the range of application. Since we do not set any restrictions on the irradiated gas geometry or photoheating process in our model, it is expected applicable to any objects and with any photoheating processes.

Comparison with Numerical Simulations
We compare physical quantities at the photoevaporation base derived by our analytical model with those directly computed by self-consistent radiation hydrodynamic simulations. We examine the validity of the model assumptions and the modelʼs predictability of ρ 2 and v 2 . Our analytical model can apply to any objects subject to photoevaporation, such as molecular clouds, PPDs, planetary atmosphere, and galactic (mini)halos. Nevertheless, in the present study we focus on photoevaporation of PPDs, which have a relatively complex geometry compared to spherical clouds/atmospheres, and therefore would be optimal to demonstrate the range of application.
Prior numerical works have shown that for PPDs, photoevaporation can be driven by FUV (6 eV  hν  13.6 eV; e.g., Richling & Yorke 1997;Johnstone et al. 1998;Störzer & Hollenbach 1999;Gorti et al. , 2015Wang & Goodman 2017;Nakatani et al. 2018aNakatani et al. , 2018bNakatani et al. , 2021Komaki et al. 2021), EUV (13.6 eV  hν  100 eV; e.g., Hollenbach et al. 1994;Shu et al. 1994;Yorke & Welz 1996;Richling & Yorke 1997;Font et al. 2004;Alexander et al. 2006), and X-rays (100 eV  hν  10 keV; e.g., Alexander et al. 2004;Ercolano et al. 2008Ercolano et al. , 2009Ercolano et al. , 2021Owen et al. 2010;Picogna et al. 2019;Wölfer et al. 2019;Monsch et al. 2021). FUV can heat the gas via photoelectric effects of grains and/or H 2 pumping, and EUV/X-ray can via photoionization of various chemical species whose ionization potential is exceeded by the photon energy. We adopt the methods of Nakatani et al. (2018a), and they are briefly summarized in Section 3.2. We refer the reader to Nakatani et al. (2018a) for more details regarding the numerical implementations. While our analytical model is expected applicable to various objects regardless of photonʼs energy band, we run simulations and make comparisons only for EUV-driven photoevaporation of PPDs in this study. This simplification can make the discussions simpler and save the computational cost significantly. This also suffices for the purpose of the investigation here. Thorough comparisons for FUV-and X-ray-driven photoevaporation with a wide variety of objects are ideal but may be unpractical in a single paper. It will be addressed in future numerical studies.

Application of the Necessary Condition for PPD Photoevaporation
Before we move on to the comparison between the analytical model and the simulation results, we estimate , the work per unit volume done by the external body force g. For steady PPDs, g includes the stellar gravity and the centrifugal force. We can calculate  analytically by assuming that the centrifugal force cancels out the radial component of gravity, leaving only the vertical component where G is the gravitational constant, M * is the stellar mass, θ is the polar angle at the point a ray hits the base, and β is the angle between the pole and the tangential line to the base (Figure 4). We have assumed the base is geometrically thin in the derivation of Equation (25). Then, using Equations (16) and (25), we obtain a rough threshold radius above which our model assumption is met,  (26) is an implicit relation with respect to r since the parameters other than M * are r dependent.

Simulation Methods
We suppose an axisymmetric disk irradiated by the direct EUV field from the host star with mass of M * . The disk is composed of gas and dust with the dust-to-gas mass ratio of Figure 4. Schematic of the model with which we derive M 2 and ρ 2 analytically to compare with the simulation results in Section 3.3.2. The stellar UV ray hits the disk surface (base) with an polar angle θ. The tangential line to the base has a polar angle β. The ray is incident on the surface with θ − β with respect to the base. Photoevaporative winds are assumed to launch nearly vertically to the base. 0.01. We use a reduced thermochemistry model where three chemical species (H, H + , and e) are included, and three chemical reactions are taken into account: EUV photoionization H + γ → H + + e, radiative recombination H + + e → H + γ, and collisional ionization H + e → H + + 2e. The excess energy of the thermalizing electrons associated with EUV photoionization goes into the heat for the gas and drive photoevaporative winds.
The radiation hydrodynamic simulations are performed in 2D spherical polar coordinates (r, θ). We use a publicly available hydrodynamic simulation code PLUTO (Mignone et al. 2007) suitably modified by implementing radiative transfer and thermochemistry (Kuiper et al. 2010;Kuiper & Klessen 2013;Nakatani et al. 2018aNakatani et al. , 2018bKuiper et al. 2020). The code solves the time evolution of gas density ρ, the three components of velocity v = (v r , v θ , v f ), total gas energy density, and chemical abundances {y i }, where p is the gas pressure; Γ is the total specific heating rate, and Λ is the total specific cooling rate; R i is the total reaction rate for the corresponding chemical species i (=H, H + , e); E and H are the total gas energy density and the total enthalpy density, respectively, with the specific heat ratio γ. We adopt the equation of state for an ideal gas, p = ρk B T/μm H , where k B is Boltzmann constant, and μ is mean molecular weight. The divergence for the azimuthal component of the Euler equation, Equation (30), is written in the angular momentum conserving form, i.e., the specific angular momentum of a fluid element conserves. The divergence ∇ l is defined as For the photoheating process, we incorporate EUV-induced photoionization heating. The UV flux is computed in a frequencydependent manner with ray-tracing radiative transfer. We do not incorporate the diffuse EUV but instead use the case B recombination rate. The EUV emission rate of the star is set to Φ EUV = 0.5 × 10 41 s −1 (e.g., ) in the fiducial run, and the EUV spectrum is assumed to be the blackbody with the effective temperature T eff = 10 4 K. The cross section of atomic hydrogen for EUV photoionization is given by s n n =ń --6.3 10 cm 18 1 2.8 2 ( ) (e.g., Osterbrock & Ferland 2006), where ν is he photon's frequency and ν 1 is that at the Lyman limit (hν 1 ≈ 13.6 eV ≈ 912 Å). For coolants, we incorporate radiative recombination cooling of H II (Spitzer 1978), Lyα cooling of H I (Anninos et al. 1997), and dust-gas collisional heat transfer (Yorke & Welz 1996). Dust temperatures are determined by solving 2D radiative transfer for both of the direct stellar irradiation and diffuse thermal emission of dust with a hybrid scheme (Kuiper et al. 2010(Kuiper et al. , 2020Kuiper & Klessen 2013). We assume ISMlike grains and use the dust opacity table in Draine & Lee (1984).
The computational domain is set to r = [0.01, 2] × r g ≈ [0.089, 18] au and θ = [0, π/2] rad, where r g is the gravitational radius (Hollenbach et al. 1994 The domain resolves both regions above and below the r exp (Equation (26)). The disk is axisymmetry around the rotational axis (θ = 0) and is symmetric with respect to the midplane (θ = π/2). The computational domain is logarithmically spaced with 384 cells in the radial direction. The meridional direction is divided into two domains at θ = 1 rad. The upper region (0 θ 1 rad) is spaced with 180 cells, and the lower region (1 rad θ π/2 rad ) is spaced with 300 cells. We use much higher spatial resolution for the lower region where the bulk of the disk mass is contained and the density/pressure gradients are steep. Our simulations simultaneously solve hydrodynamics, radiative transfer, and thermochemistry, and are thus selfconsistent. Self-consistency is indispensable to accurately trace the energy partition from the irradiation to the gas' kinetic plus thermal energy and to make the comparison with our analytical model meaningful.

Results
The disk system reaches a quasi-steady state within about 30 yr after the computation starts. The flow structure is almost steady on the astronomical unit scale but is time variable on the scale of the base's width, s D~á ñ » ( ) ( ) . There, the EUVheated gas intermittently evaporates from the opaque disk in an explosive manner. The variable timescale is roughly given by the sound crossing time over Δ. To make our analysis and comparison visually simpler, we time average the snapshots of the simulation over ∼10 yr and mainly use the data within R  7 au, which is fairly away from the computational outer boundary. We note that our model yields an error scaling with the thickness of the base if applied to a nonsteady system.

Base Property
We empirically define the photoevaporation base at which the column density of atomic hydrogen is N H I = 10 18 cm −2 and study the profiles of various physical quantities along lines normal to the base. Figure 5 shows an example of such analysis, where we see the profiles along the normal line located at R = 2.13 au. The color maps on the left show the density structure; the top panel is a global view, and the bottom panel is a zoomed view around the normal line (the light-green solid line).
The profiles of physical quantities along the normal line are shown by the three panels on the right side in Figure 5. As our analytical model explains in Section 2.2, the gas is nearly isobaric across the base; density and temperature change by orders of magnitude, while pressure is continuous. The abrupt change in density is also evident by the zoomed map; it is observed as a nearly discontinuous change in the colors at the base (the yellow dashed line). The bottom panel on the right-hand side shows that the velocity component normal to the base v n dominates the transverse component v t in the very vicinity of the base. This is also evident from the streamlines curving in the zoomed map. It is noteworthy that all these qualitative characters justify the assumptions and results of our analytical model.
To compare with the simulation results quantitatively, we extract the model parameters (c 1 , c 2 , and ρ 1 = mn 1 ) from the obtained profile. The taken values are n 1 = 1.7 × 10 7 cm −3 , c 1 = 1.3 km s −1 , c 2 = 8.5 km s −1 , and are represented by the colored dots in the top-right and middle-right panels. We also need to give  and F. We calculate  via Equation (25). Regarding F, given a base property, the upper limit of f is set by the non-attenuated heating flux F 0 as . The density contours for n H = 10 6 , 10 7 cm −3 and the base almost overlap. The arrows in the top represent the velocity field. The lengths are scaled by the magnitude. The velocity field is drawn by streamlines in the bottom panel. The light-green solid line is a normal line to the base at R = 2.13 au, along which we plot the quantities in the right panels. (Right panels) Profiles of the quantities along a normal line to the base (the light-green line in the density maps). The top, middle, and bottom panels show number density/gas temperature, isothermal sound speed/pressure, and velocities/Mach number, respectively. We show both the normal (v n ) and parallel (v t ) components of the velocity. The profiles are plotted with solid lines. The secondary x-axis of the top-right panel is the coordinate along the light-green lines in the color maps. We also show n 2 and M 2 derived from the analytical model (Equations (20) and (22)), using n 1 , c 1 , and c 2 indicated by the dots in the top and middle panels. (See the text for more details about the derivation of n 2 and M 2 .) The predicted values agree well with those calculated by the radiation hydrodynamic simulations, validating the predictability of our model. We note that the ionization front, where density and temperature sharply change, is properly resolved with about six computational cells.
where ΔE is the mean excess energy deposited to the gas per photon. The actual F reaching the base is slightly smaller than the non-attenuated value because of absorption along the ray's line of sight from the host star to the base. The reduction is found to be of the order of ∼10% in our simulation, corresponding to N H I ∼ 10 17 cm −2 . While Equation (33) provides a good estimate for the f parameter, we directly take the value of F from the simulation to eliminate the error due to flux attenuation. The obtained value is f = 0.49 in the specific case of Figure 5. (For reference, » f 1.2 max is obtained by substituting the taken values into Equation (33).) As we will see below, this f yields ρ 2 and v 2 consistent with those in the simulations. We also measure the angles θ and β from the simulation. The obtained value is (θ − β) = 0.10 rad in the specific case shown in Figure 5. Since the base has a finite width, the angle (θ − β) has uncertainty of the order of Using ρ 1 , c 1 , c 2 , and f, we derive M = 0.2 2 and ρ 2 /ρ 1 = 0.02 with the analytical model (Equations (20) and (22)). The derived M 2 and ρ 2 /ρ 1 are represented by the horizontal dashed lines in the top and bottom panels on the right-hand side of Figure 5. They are in good agreement with those obtained by the simulations. We also compare the vertical profiles between the model and simulations at different radii (R 7 au) using the same procedure as above. For R  1 au, M 2 is orders of magnitude smaller than unity, which means escaping photoevaporative winds are strongly gravitationally inhibited. The EUV heating balances with the sum of the line cooling in the photoheated and processing layers at these radii. Our model does not apply here as discussed in Section 2.2.3. (Indeed, the model predicts M 1 2 ( )much higher than the simulated value.) For R  1 au, we observe vertical acceleration across the base. The derived M 2 and ρ 2 /ρ 1 agree with those from the simulations within a factor of 2 and within several tens of percents, respectively. This demonstrates the validity and predictability of the analytical model. For 1 au  R  1.5 au, v n is immediately decelerated to nearly zero over a few times the processing layer's width.

Mass-loss Rates
We compare the mass-loss rates  M between the analytical model and the simulations. We measure  M of the simulations by integrating the mass flux over a spherical surface S with a radius of r, The integration is performed only for N H I 10 18 cm −2 to exclude the disk region. In a quasi-steady state,  M sim accounts for the integrated mass flux of the winds launched from a partial disk contained within S. Regarding the mass-loss rate of the analytical model, we integrate the mass flux normal to the base (Figure 4), where R in is the lower limit of the integration range, and v n is the velocity component normal to the base. We set R in to 1 au, since winds are visually identified from R  1 au (see Figure 5 and Section 3.3.1). The mass flux ρv n = ρ 2 v 2 is computed by the same procedure described in Section 3.3.1; ρ 1 , c 1 , c 2 , and other geometrical parameters are taken from the simulations, and f is computed with those taken quantities to calculate ρ 2 and v 2 . A fully analytical derivation of  M is ideal, but it requires modeling the R and Φ EUV dependence of these model parameters. This is beyond the scope of this study. We therefore use this semi-analytic approach to calculate  M ana . Comparison between  M sim and the semi-analytically derived  M ana is still meaningful to validate our model's predictability for the mass flux of photoevaporative winds. Figure 6 compares  M ana and  M sim for various EUV emission rates: Φ EUV = 1 × 10 39 , 5 × 10 40 , and 1 × 10 42 s −1 . In all cases,  M ana agree well with  M sim regardless of the different measurement techniques (Equations (35) and (34)). It indicates that our analytical model can predict the mass flux of photoevaporative winds accurately, confirming the validity of the model again. We have also checked that Equation (16) is met at the base (even in R < 1 au), being consistent with our model.

Discussion
We have shown the validity of our analytical model by comparison with the results of self-consistent radiation hydrodynamic simulations in Section 3. The analytical model requires the cold layer density ρ 1 and the isothermal sound speeds in the cold layer c 1 and the photoheated layer c 2 as model parameters. In Section 4.1, we discuss how these quantities are determined, which explains why these quantities need to be treated as parameters rather than derivable quantities in an analytical manner.

Input Parameters
Our model derives the launch velocity v 2 and density ρ 2 of photoevaporative winds as a response to an injected photoheating energy into the processing layer. We treat the isothermal sound speeds, c 1 and c 2 , and the cold gas density ρ 1 as input parameters. The sound speeds are set by thermochemistry, including adiabatic expansion/compression in the photoheated and cold layers. Similarly, the cold gas density ρ 1 (and the photoheating flux F) depends on radiative transfer, chemistry, and advection over a much larger scale than the width of the processing layer ( s á ñn ; 1 1 ( ) see also Appendix). Determination of the parameters is therefore another issue to be discussed through different physics from the analytical model that focuses on the processing layer scale, but it this is beyond the scope of the present study. One may need numerical computations to get the accurate values of these parameters, as in Section 3.3. For reference, we demonstrate how ρ 1 and F are derived using a simple analytical model described in Appendix. Nevertheless, our model is also available with approximated values for the input parameters, as has often been done in relevant studies. For instance, c 2 ≈ 10 km s −1 is often chosen as a typical sound speed for ionized gas at the equilibrium temperature T ∼ 10 4 K (e.g., Spitzer 1978;Hollenbach et al. 1994). Equation (23)  1 , respectively; see Figure 6), and this rough equation explicitly shows the overall physical determinants for the mass-loss rates. The errors come from the facts that c 2 and (θ − β) are dependent on Φ EUV and R, and that radiative cooling becomes more effective for higher Φ EUV . With weak EUV radiation, the photon heats an upper layer of the disk, and the base has a more flared structure. This results in a relatively large (θ − β). As for the photoheated gas temperature, the heated gas flows out before it reaches an equilibrium temperature (∼ 10 4 K) and thus has c 2 < 10 km s −1 . In our simulations, typical values at R  10 au are found to be c 2 ≈ 10 km s −1 and (θ − β) ∼ 0.01-0.1 with Φ EUV = 10 42 s −1 , and c 2 ≈ 5-10 km s −1 and (θ − β) ∼ 0.1-0.3 with Φ EUV = 10 39 s −1 . The estimated  ¢ M ana gets even closer to  M sim when taking into account these corrections; the corrected one corresponds to  M ana (the red lines in Figure 6).
Hence, hydrodynamics effects on the gas internal energy can be significant for a weak irradiation (long heating time). In such cases, using c 2 including the effects of adiabatic expansion/ compression is preferred to accurately calculate the physical quantities of the launched gas with our model. Performing selfconsistent radiation hydrodynamic simulations is one way to derive an accurate c 2 . Alternatively, pure thermochemistry models can also give a reasonable c 2 if adiabatic cooling is properly included. One needs information on flow geometry to compute adiabatic cooling without performing hydrodynamic simulations. For disk photoevaporation, recent analytical models of thermally driven winds (Clarke & Alexander 2016;Sellek et al. 2021) may be available for that purpose. Analytic modeling of c 2 (or flow structure) is a next step in future work. It is also essential for a more comprehensive understanding of the basics regarding photoevaporation.
Generally, c 2 and F are not completely independent parameters but can be tightly correlated especially when adiabatic cooling is significant compared to other coolants in the photoheated layer. The isothermal sound speed c 2 is set by the photoheating flux in the photoheated layer, which is close to F near the base, and cooling rates of coolants and adiabatic expansion in the photoheated layer. This is an implicit relation with respect to the photoheated layer's temperature, and  (16)), is satisfied at all radii in the planes. therefore the relation between F and c 2 may not be analytically derivable even if flow geometry is given. One may use an empirical relation or just set c 2 as an independent parameter to derive the launch velocity and density self-consistently from conservation laws. Our model uses the latter approach. Treating c 2 as an independent model parameter is in essence equivalent to giving the flow geometry (or adiabatic cooling) and cooling rates in the photoheated layer. As discussed in Section 2.2.2, the values of the model variables (ρ 2 , v 2 , v 1 ) are set so that they are consistent with the large-scale structure of the system. Overall, to derive the specific values, information on the large-scale structure would be needed to reduce degree of freedom for the problem.

Summary and Conclusions
Photoevaporation is an essential process to determine the fates of various objects in a wide range of scales. It has been known that ultraviolet and X-rays are capable of driving the winds. To understand this process in a unified picture, we have analyzed the layers in the vicinity of the wind-launching region by developing an analytical model. We have modeled the region as a photoheated layer and a cold layer divided by a processing layer (Figure 2). We have defined the density, velocity, and sound speed of the photoheated layer (ρ 2 , v 2 , c 2 ) and the cold layer (ρ 1 , v 1 , c 1 ). Using the conservation laws of mass, momentum, and energy (Equations (12)-(14)), we have derived the analytical forms of ρ 2 , v 2 , and v 1 as a function of the other parameters involved in this problem (radiation flux, ρ 1 , v 1 , c 2 , etc). Our model builds on these assumptions: (a) the processing layer is geometrically thin, (b) radiative cooling is negligible within the processing layer, and (c) the gas is subsonic in the cold layer. These assumptions set the range of the application of our model (see also Section 2.2.3).
Our findings are summarized as follows: 1. The properties of the launched winds are characterized by a parameter f (Equation (19)) and the specific heat at constant pressure c p . 2. The launch velocity v 2 is a monotonically increasing function of f ( Figure 3). The launched gas is subsonic for f < (c p + 1/2)/2. The radiation energy mainly goes into the enthalpy flux, and v 2 is determined by the gas thermal expansion (enthalpy-flux limited). If f > (c p + 1/2)/2, the radiation energy largely goes into the kinetic energy in terms of momentum conservation (dynamics limited). The launch velocity is supersonic in this regime. 3. The density ratio of the photoheated layer to the cold layer ρ 2 /ρ 1 decreases for higher f (Equation (22)). It is essentially determined by pressure balance for f = 1.
Our model treats the upstream density (ρ 1 ) and upstream/ downstream sound speeds (c 1 and c 2 ) as input parameters. Their values depend on various physical processes not only at the base but also at a larger scale, e.g., radiative transfer from the radiation source to the base, flow geometry, and chemical reactions along the ray. Our model explains that the problem of deriving the launch density and velocity reduces to the problem of determining these parameters.
We have compared our analytical model with the results of self-consistent radiation hydrodynamic simulations for photoevaporating PPDs. We have confirmed that our model well explains the launched wind property. Since our model does not limit the geometry of irradiated objects or photoheating processes, we expect it to be applicable to various objects irradiated not only by EUV but also by FUV and X-rays. Comparisons with self-consistent radiation hydrodynamic simulations are necessary to explore the applicability in future works.

(
) . The derived quantities (ρ 2 , v 1 , v 2 ) are in principle valid for any combination of F and ρ 1 . In practice, on the other hand, these two quantities can relate to each other. A high upstream density can yield high-density flows that shield the photons. It may be thus unrealistic to give an arbitrarily high ρ 1 for a given F. Whether the photons can reach the layer with ρ 1 is a problem to be solved by radiative transfer within the region where the bulk of the photons are absorbed. It depends on the main photon absorbers, their production rates, multidimensional flow geometry, etc. and is determined by larger scale physics and chemistry than the scale of the processing layer. Numerical computations would be thus required to get ρ 1 accurately. Nevertheless, in this appendix we demonstrate how F and ρ 1 link to each other, using a simplified 1D model. The discussion provides us with a way to roughly estimate the upstream density ρ 1 and with a reference for general situations.
As in the case of EUV photoionization heating, we assume that a chemical species is responsible for the absorption of the photons and yields a product as a result of the photochemical reaction. The chemical abundance of the product x is approximately zero in the upstream region and increases to x 2 in the downstream region. In a steady state, advection and (photo)chemical reactions balance as where v is gas velocity, N is the photon number flux, and R x is the net chemical reaction rate per unit volume (R x ∝ n 2 if dominated by two-body chemical reactions). Note that ∇ · (nv) = 0 in a steady state. We further assume that the bulk of the photons is absorbed within a length of  s D~á ñ --L n x s n dln d 1 1 1 | ( ) | ( ) in the downstream region. This scale length is characterized by photochemistry and characteristic lengths of the system, e.g., curvature of the irradiated object and the scale height for a gravitationally stratified gas. We integrate Equation (A1) over ΔL from the surface just behind the processing layer to calculate the flux actually reaching the layer. Equation (A1) reduces to where N 0 is the unattenuated photon number flux, N is the photon number flux reaching the base, and Δx 2 is the abundance difference within the downstream region and normally takes positive values. The second term on the righthand side accounts for the absorption to increase the abundance by Δx 2 within ΔL. The third term is absorption by continuously produced atoms/molecules via chemical reactions. Equation (A2) explains that the flux reaching the processing layer is residual after absorbed by originally existing and reproduced atoms/molecules within ΔL. Equation (A2) also implies two regimes depending on which of the second and third terms on the right-hand side is dominant, as has similarly been discussed in the context of molecular clump photoevaporation (e.g., Bertoldi 1989; Henney 2001; Nakatani & Yoshida 2019). If the second term dominates the third term, the regime is called advection dominated; we call the other regime chemistry dominated. The photoheating flux F is described by the product of N and the mean excess energy per photon ΔE.
For the advection-dominated regime, Equation (A2) is approximated to F ≈ F 0 − n 2 v 2 Δx 2 ΔE. Substituting this into Equation ( is practically in the chemistrydominated regime for a high-density region, since the third term on the right-hand side is essentially proportional to n 2 if two-body chemical reactions are the most relevant (it could be ∝n 3 for further high-density regions). We define a mean reaction coefficient as The flux is necessary to be positive at the processing layer to drive photoevaporative winds. This requirements sets an upper limit for n 2 and thereby for n 1 through Equation ( The flow is launched at a subsonic velocity if the cold gas density is high (n 1 > n c ).