and in this case, the rapid ionization condition in the gravity-inhibited regime, $t_\mathrm{ioni} < \mathrm{min}\left[t_\mathrm{rec} (\cs  = \vg ), R/\mach_\mathrm{g}^\mathrm{II} \vg  \right]$, is
$t_\mathrm{ioni} < t_\mathrm{rec}$, i.e.,
$\varphi > \varepsilon  q x^{\beta + 1}$.
This is always satisfied in the Type~II gravity-inhibited regime.
On the other hand, the above ratio is less than unity in
\[
\left\{
(x, \varphi)~|~
x > \frac{\varepsilon ^{-1}}{2},
\varphi > \frac{\varepsilon  q x^{1+\beta}}{(2-x^{-1}\varepsilon ^{-1})^2},
\right\}
\]
and in this regime, the rapid ionization condition becomes
\[
\frac{t_\mathrm{ioni}}{t_\mathrm{rec}} < x^{-1}\varepsilon ^{-1} - \braket{1-\frac{t_\mathrm{ioni}}{t_\mathrm{rec}}},
\]
or equivalently $x < \varepsilon ^{-1}$.
% Note that this equation is used only in the gravity-inhibited regime to evaluate $\mach_\mathrm{g}^\mathrm{II}$.
After all, the rapid ionization condition in the gravity-inhibited regime, $t_\mathrm{ioni} < \mathrm{min}\left[t_\mathrm{rec} (\cs  = \vg ), R/\mach_\mathrm{g}^\mathrm{II} \vg  \right]$, reduces to $x < \varepsilon ^{-1}$.

% Now we analyze the other condition in Condition~II-2: $c_\mathrm{ch}^\mathrm{II} < \vg $.
% This is equivalent with $\varphi < \varphi_\mathrm{g}$ (\eqnref{eq:varphi_g}).
% The line $\varphi = \varphi_\mathrm{g}$ connects to the boundary line between Type~I and Type~II of this class ($\varphi = \varepsilon ^{3/2} x$) at $x = \varepsilon ^{-1}$.
% One can prove that $\varphi_\mathrm{g} > \varepsilon  x^{1/2}$ for $x < \varepsilon ^{-1}$, which guarantees \eqnref{eq:chII_available_condition:dimensionless} is satisfied at $\varphi = \varphi_\mathrm{g}$.
% For $x > \varepsilon ^{-1}$, $\varphi_\mathrm{g} < \varepsilon  x^{1/2}$ holds (with asymptotically approaching $\varphi = \varepsilon  x^{1/2}$ for $x \gg 1$), and it breaks the condition for \eqnref{eq:chII_original_def:dimensionless} to be available, \eqnref{eq:chII_available_condition:dimensionless}.
% This implies that the gravity-inhibited wind regime must be Type~I for $x > \varepsilon ^{-1}$, which is consistent with the discussions in \secref{sec:typeI:regime}.

The other condition in Condition~II-2, $c_\mathrm{ch}^\mathrm{II} < \vg $, is equivalent with $\varphi < \varphi_\mathrm{g}$ (see \eqnref{eq:varphi_g}).
To summarize Condition~II-2,
the gravity-inhibited regime for Type~II winds is defined by
\[
\left\{
(x, \varphi)~|~
x < \varepsilon ^{-1},
\varphi > \varepsilon  q x^{\beta + 1},
% \varphi_\mathrm{II} < x^{-(\beta + 2)}
\varphi < \varphi_\mathrm{g}
\right\}.
\]
% and an asymptotic form can be expressed by
% \[
% \left\{
%     (x, \varphi)~|~
%     x < \varepsilon ^{-1},
%     ~
%     \varphi > \varepsilon  q x^{\beta + 1},
%     ~
%     \varphi_\mathrm{II} < x^{-(\beta + 2)}
%     \right\}.
% \]
In this regime, winds have a sound speed of $\sim \vg $ and a Mach number of
\begin{equation}
\mach_\mathrm{g}^\mathrm{II}
\approx
\braket{\frac{\varphi q x^\beta}{\varepsilon  }}^{1/2}
\left[
x^{-1} \varepsilon ^{-1} - 1
+ \braket{\frac{\varepsilon  q  x ^{\beta + 1} }{\varphi}}^{1/2}
\right]^{-1}
\label{eq:mach_g:typeII}
\end{equation}
% The second term is guaranteed less than unity since $t_\mathrm{ioni} < t_\mathrm{rec}$ in the Type~II regime.
Mach number increases with $\varphi$ at individual distance owing to increasing heating with $\varphi$ and approaches unity at the $c_\mathrm{ch}^\mathrm{II} = \vg $ boundary (see \eqnref{eq:varphi_g_def}).

Finally, we evaluate Condition~II-3.
This condition clearly corresponds to the remaining parameter space that is not covered by Conditions~II-1 and II-2, that is,
\[
\left\{
(x, \varphi)~|~
\varphi < f_\mathrm{eq} x,
\varphi > \varphi_\mathrm{g},
\varphi > \varepsilon ^{3/2} x
\right\}.
\]
In this space, the heating timescale is shorter than the gravitational timescale for Type~II winds but is not short enough to heat the gas to the equilibrium temperature until traveling to the height of $\sim R$.
Therefore, $\cs  \sim c_\mathrm{ch}^\mathrm{II}$ and $\mach \sim 1$ in this regime.

% \mycomment{Asymptotic boundaries are actually not asymptotic lines of $\varphi_\mathrm{g}$ or $f_\mathrm{eq} x$. They seem completely useless, why? Probably it ignores the initially injected energy or additional energy. That underestimates both the characteristic sound speed significantly.}
% We can derive an asymptotic form of this parameter space by evaluating the condition, $c_\mathrm{ch}^\mathrm{II} > \ceq $, with the approximated $c_\mathrm{ch}^\mathrm{II}$ in the two limiting cases.
% When recombination is slow (\eqnref{eq:c2_in_slow_recombination_valid_paramspace}), $c_\mathrm{ch}^\mathrm{II} > \ceq $ is not met since $\varepsilon  < 1$ in this class.
% This means that additional energy deposition through reionizations is indispensable to form isothermal winds.
% In the case of fast recombination (\eqnref{eq:rapid_recombination_condition_paramspace}), $c_\mathrm{ch}^\mathrm{II}  > \ceq $ can be approximated to $c_2  > \ceq $, and this means
% \[
%     \varphi_\mathrm{II} >  x.
% \]
% Since $(\varepsilon  q)^{-1} > \mathrm{max}(1, q)\varepsilon $ and $(\varepsilon  q)^{-1} > \varepsilon ^{\beta +2} q^{-1}$ in this class,
% Condition~II-1 reduces to $\varphi_\mathrm{II} > x$.
% Isothermal, ionized winds develop in this parameter space.

% \mycomment{Asymptotic boundaries are actually not asymptotic lines of $\varphi_\mathrm{g}$ or $f_\mathrm{eq} x$. They seem completely useless, why? Probably it ignores the initially injected energy or additional energy. That underestimates both the characteristic sound speed significantly.}
% We can also derive an asymptotic form of the gravity-inhibited wind regime, $c_\mathrm{ch}^\mathrm{II} < \vg $, by considering the approximated $c_\mathrm{ch}^\mathrm{II}$ in the two limiting cases.
% In the slow recombination limit (\eqnref{eq:c2_in_slow_recombination_valid_paramspace}),
% $c_\mathrm{ch}^\mathrm{II} $ can be approximated to $ c_1 (= \varepsilon ^{1/2} \ceq )$, and the condition, $c_\mathrm{ch}^\mathrm{II} < \vg $, reduces to $x< \varepsilon ^{-1} $.
% Hence, the gravity-inhibited regime for Type~II winds with slow recombination corresponds to the parameter space where
% \[
%     \left\{
%     (x, \varphi)~|~
%     x < \varepsilon ^{-1},
%     ~
%     \varphi < \varepsilon ^{\beta + 2} q^{-1} x,
%     ~
%     \varphi > \varepsilon  q x^{\beta + 1}
%     \right\}.
% \]
% On the other hand, when recombination is fast (\eqnref{eq:rapid_recombination_condition_paramspace}),
% the condition, $c_\mathrm{ch}^\mathrm{II} < \vg $, reduces to
% \[
%     \varphi < (\varepsilon  q)^{-1} x^{-(\beta + 2)}.
% \]
% Thus, the gravity-inhibited regime for Type~II winds with rapid recombination corresponds to the parameter space where
% \[
%     \left\{
%     (x, \varphi)~|~
%     \varphi < (\varepsilon  q)^{-1} x^{-(\beta + 2)},
%     ~
%     \varphi > \varepsilon ^{\beta + 2} q^{-1} x
%     \right\}.
% \]
% This region is covered by $x < \varepsilon ^{-1}$,
% and hence the parameter space where $c_\mathrm{ch}^\mathrm{II} < \vg $ always satisfies the fast ionization condition.
% It is worth noting that in the slow recombination limit, where we only count the energy deposition until a launched hydrogen is photoinoized once (see \eqnref{eq:c1_def}) and neglect the additional energy deposition through reionizations,
% \eqnref{eq:typeII:mach_g:self_consistent_condition} cannot be met for $c_1 < \vg $.
% Then, the self-consistent condition to form steady escaping winds in this regime is
% \[
%     c_p \vg ^2 \approx
%     \frac{\avesigma  \aveenergy   \phieuv }{4\pi m r^2}
%     \frac{\mach_\mathrm{g}^\mathrm{II}\vg t_\mathrm{ioni}}{\mach_\mathrm{g}^\mathrm{II}\vg }.
% \]
% However, this condition cannot be satisfied with any $\mach_\mathrm{g}^\mathrm{II}$ for $\vg  > c_\mathrm{ch}^\mathrm{II} (= c_1) $.
% This indicates that neglecting the contributions from reionizations is invalid.
% i.e., $c_\mathrm{ch}^\mathrm{II} \approx \varepsilon ^{1/2}$ is too much approximated (zeroth order) to characterize flow speed.
% In other words, $c_1 < \vg $ and neglecting the second term on the RHS of \eqnref{eq:c2_reduced} are incompatible.
% This is because, in reality, energy deposition through reionizations contributes to heating, resulting in developing steady winds with $\cs  \sim \vg $, even when $c_1 < \vg $.
% Hence, for the cases of $c_1 < \vg $, the additional energy deposition through reionizations always contributes to forming steady Type~II winds in the gravity-inhibited regime.

The sound speed in this regime is obtained by solving the implicit equation Eq.\eqref{eq:chII_original_def:dimensionless} with respect to $c_\mathrm{ch}^\mathrm{II}/\ceq $.
The expression of the root would have a complex functional form with arguments of $\varepsilon $, $\varphi \varepsilon  q/x $, and $\varphi/\varepsilon  x$, which might be unpractical to be presented.
Instead, it would be convenient to show approximate forms of the root here.
In the Sa-class, the parameter space where $\cs  = c_\mathrm{ch}^\mathrm{II}$ always contains the slow recombination regime \eqnref{eq:c2_in_slow_recombination_valid_paramspace}.
% as $f_\mathrm{mid} > \varepsilon ^{3/2}$.
% Slow recombination is realized when $f_\mathrm{mid} \gg \varepsilon ^{3/2}$, which reduces to $q \ll 2^{\beta+3/2} \varepsilon ^{(2\beta+1)/2}$.
This is because in this class, the second term on the RHS of \eqnref{eq:chII_original_def:dimensionless} is exactly zero at the Type~I and Type~II boundary ($\varphi = \varepsilon ^{3/2}x$), and \eqnref{eq:c2_in_slow_recombination_valid_paramspace} always exists within the Type~II regime.
% This is realized when $q \ll \varepsilon ^{(2\beta + 1)/2}$, i.e., when $\varphi = \varepsilon ^{\beta+2}q^{-1}x$ (cf. \eqnref{eq:c2_in_slow_recombination_valid_paramspace}) is much larger than $\varphi = \varepsilon ^{3/2} x$ (cf. \eqnref{eq:typeII_regime_Saclass}).
% In this case, \eqnref{eq:c2_in_slow_recombination} gives an approximate value of $c_\mathrm{ch}^\mathrm{II}$ as
The approximated $c_\mathrm{ch}^\mathrm{II}$ is
\[
c_\mathrm{ch}^\mathrm{II} \approx c_1 = \varepsilon ^{1/2} \ceq .
\]
This approximation holds in the parameter space of
\[
\left\{
(x, \varphi)~|~
x > \varepsilon ^{-1},
\varphi \ll \varepsilon ^{\beta + 2} q^{-1} x,
% \varphi \ll f_\mathrm{mid} x,
\varphi > \varepsilon ^{3/2} x
\right\}.
\]
Here, the condition $x > \varepsilon ^{-1}$ originates from a requirement for winds to be escaping, $c_\mathrm{ch}^\mathrm{II} > \vg $.
The corresponding area is illustrated by the orange-shading bounded by the red and blue dashed lines in the top right panel in \fref{fig:classmaps}, as an example.
Having this area where the additionally deposited energy through reionizations is negligible is a unique characteristic of the Sa-class.
Note that in $x < \varepsilon ^{-1}$, gravity is so strong that the additional energy deposition through reionizations is always needed to form steadily escaping winds.

% \mycomment{edit here.}
% When $c_\mathrm{ch}^\mathrm{II} $ is smaller than $ v_g$, i.e., $ x  < \varepsilon  ^{-1}$, heating is very slow so that winds are inhibited by gravity.
% To form escape winds, the gas must be heated to $\cs  \sim \vg  $ at a height of $\sim R$.
% However, in this zeroth-order estimation where we neglect the energy deposition from reionizations, the total energy gain during the gas traveling to the height of $\sim R$ is always $c_p c_1^2 = c_p \varepsilon  \ceq ^2 < c_p \vg ^2$ regardless of the flow velocity.
% This is because, in this estimation, we only count the energy deposition until a launched hydrogen is photoinoized once (see \eqnref{eq:c1_def}).
% Therefore, $\cs $ cannot reach $\vg $ however long time the gas takes to be heated until reaching a height of $R$ (see also discussions below \eqnref{eq:c1_def}).
% This means that, if we neglected the additional energy deposition from reionizations, steady winds could not develop whenever the energy deposition of the first ionization is significantly small $c_1 < \vg $, or equivalently $x < \varepsilon ^{-1}$.

% This does not indicate, however, that steady winds do not develop when $c_1 < \vg $ but indicates that neglecting the contributions from reionizations is invalid.
% In other words, $c_1 < \vg $ and neglecting the second term on the RHS of \eqnref{eq:c2_reduced} are incompatible.
% This is because, in reality, energy deposition through reionizations contributes to heating, resulting in developing steady winds with $\cs  \sim \vg $, even when $c_1 < \vg $.
% Hence, for the cases of $c_1 < \vg $, we always need to take into account the energy deposition through reionizations.

In this class, the fast recombination regime (\eqnref{eq:rapid_recombination_condition_paramspace}) is not necessarily contained in the free wind regime where $\cs  = c_\mathrm{ch}^\mathrm{II}$.
This is physically understood as follows.
Suppose the case where $\varepsilon $ has a value slightly smaller than unity,
$c_\mathrm{ch}^\mathrm{II} $ needs only a small contribution from the additionally supplied energy through reionizations (the second term on the RHS of \eqref{eq:chII_original_def:dimensionless}) to reach $c_\mathrm{ch}^\mathrm{II} > \ceq $.
% This can be rephrased as $c_\mathrm{mid } = \sqrt{2}\varepsilon ^{1/2} > 1$ can be easily met when $\varepsilon $ has a value slightly smaller than unity.
Hence, in this class, it is possible that the fast recombination regime is not present in Regime~B-II (the free Type~II winds), and the slow recombination regime (\eqnref{eq:c2_in_slow_recombination_valid_paramspace}) directly connects to the boundary between Regime~A (isothermal winds) and Regime~B-II (free winds), at which $c_\mathrm{ch}^\mathrm{II} = \ceq $.

% Approximate forms of $\cs $ are given by Eqs.\eqref{eq:c2_in_slow_recombination} and \eqref{eq:cchII_def} in the slow and rapid recombination limits (Eqs.\eqref{eq:slow_recombination_condition} and \eqref{eq:rapid_recombination_condition}), respectively.
% Hence, in the zeroth-order, the sound speed of steadily heated, free Type~II winds can be given by
% \[
%     \cs
%     \sim
%     \left\{
%     \begin{array}{ll}
%     c_1
%     &
%     \text{for}~
%     \left\{
%     (x, \varphi)~|~
%     \varphi < \varepsilon ^{\beta + 2} q^{-1} x,
%     ~
%     \varphi > \varphi_\mathrm{g},
%     ~
%     \varphi > \varepsilon ^{3/2} x
%     \right\}.
%     \\
%     c_2
%     &
%     \text{for}~
%     \left\{
%     (x, \varphi)~|~
%     \varphi < f_\mathrm{eq} x,
%     ~
%     \varphi > \varphi_\mathrm{g},
%     ~
%     \varphi > \varepsilon ^{\beta + 2} q^{-1} x
%     \right\}.
%     \end{array}
%     \right.
% \]
The resulting parameter space for the Sa-class is summarized by the top right panel of \fref{fig:classmaps}.

\subsubsection{Type~II winds in Sr-class}

%In the Sr-class %($\{(\varepsilon  , q)| \varepsilon   < 1, \varepsilon  q < 1, q > 1 \}$),
The Type~II regime of the Sr-class is expressed by
\begin{equation}
\left\{
\begin{array}{ll}
\varphi >  \varepsilon   q  x ^{1+\beta}& \quad\mathrm{for}~ 1 \leq   x  < (\varepsilon  q)^{-2/(3+2\beta)} \\
\varphi >  \braket{\varepsilon  q}^{3/(3+2\beta)}  x  & \quad\mathrm{for}~  x  > (\varepsilon  q)^{-2/(3+2\beta)}
\end{array}
\right.
\label{eq:typeII_regime_Srclass}
\end{equation}
(cf. Eq.\eqref{eq:typeI_regime_Srclass}).
The fast ionization condition for Condition~II-1 is $\varphi > \mathrm{max}(1, q)\varepsilon  x$ in this class.
Since $f_\mathrm{eq} > 1$ for the soft-spectrum classes,
the Type~II regime meets the fast-heating condition, $\varphi > f_\mathrm{eq} x$.
Therefore, Condition~II-1 is equivalent with $\varphi > f_\mathrm{eq} x$ for the same reason as in the Sa-class.

Regarding Condition~II-2, we do a similar analysis to that we have done for Type~II winds in the Sa-class.
The condition $c_\mathrm{ch}^\mathrm{II} < \vg $ is again $\varphi < \varphi_\mathrm{g}$ (\eqnref{eq:varphi_g}),
and $\varphi_\mathrm{g}$ connects to the intersection between $\varphi = \varepsilon  q x^{1+\beta} $ and $\varphi  = (\varepsilon  q)^{3/(3+2\beta)} x$ at $x = (\varepsilon  q)^{-2/(3+2\beta)}$ (see \eqnref{eq:typeII_regime_Srclass}).
Therefore, the gravity-inhibited regime for Type~II winds is defined by
\[
\left\{
(x, \varphi)~|~
x < (\varepsilon  q)^{-2/(3+2\beta)},
\varphi < \varphi_\mathrm{g},
\varphi > \varepsilon   q x^{1+\beta}
\right\}
\]
In this regime, $\cs  \sim \vg $ and $\mach = \mach_\mathrm{g}$ (\eqnref{eq:mach_g:typeII}).

The other parameter space corresponds to the regime of free Type~II winds, and Condition~II-3 reduces to
\[
\left\{
(x, \varphi)~|~
\varphi > (\varepsilon  q)^{3/(3+2\beta)} x,
\varphi > \varphi_\mathrm{g},
\varphi < f_\mathrm{eq} x
\right\}
\]
The typical sound speed and Mach number are $\cs  \sim c_\mathrm{ch}^\mathrm{II}$ and $\mach \sim 1$, respectively.

We can derive an approximate form of $c_\mathrm{ch}^\mathrm{II}$ for this class instead of deriving the exact root from the implicit equation \eqref{eq:chII_original_def:dimensionless}.
In contrast to the Sa-class, the Sr-class does not cover the slow recombination regime, as \eqnref{eq:c2_in_slow_recombination_valid_paramspace} is outside of the Type~II regime (\eqnref{eq:typeII_regime_Srclass}).
This means additional deposited energy always contributes to the heating for Type~II winds in this class.
Hence, $c_\mathrm{ch}^\mathrm{II}$ can be approximated to  Eq.\eqref{eq:cchII_def},
\[
c_\mathrm{ch}^\mathrm{II} \sim c_2  = \ceq \braket{\frac{\varphi_\mathrm{II}}{x}}  ^{1/2(\beta+3)}
\]
in the parameter space,
\[
\left\{( x , \varphi)~|~
\varphi_\mathrm{II} <  x ,
\varphi > (\varepsilon  q)^{3/(3+2\beta)}  x ,
\varphi_\mathrm{II} >  x ^{-(\beta+2)}
\right\}
\]
This parameter space is represented by the area bounded by the gray solid and green dashed lines in \fref{fig:classmaps}.
One can observe that the gray line gives approximate boundaries with the gravity-inhibited and isothermal wind regimes.
For reference, we also plot this boundary in the Sa-class map of \fref{fig:classmaps}.
The gray line does not necessarily give approximate boundaries with the gravity-inhibited and isothermal wind regimes, since neglecting the contribution of the initially injected energy to heating is not always justified in the Sa-class.

% \input{subsections/Lcrit_removed_3}

\section{Comparison with Hydrodynamics Simulations}
\label{sec:comparison_HD_simulations}

To assess the predictability of our phenomenological model and the validity of the adopted assumptions, we briefly compare the analytic model with hydrodynamics simulations, where ray-tracing EUV radiative transfer and nonequilibrium thermochemistry are self-consistently treated.
We stress that self-consistent coupling between hydrodynamics and thermochemistry is required to make meaningful comparisons.
We provide the computational methods of the simulations in \secref{sec:sub:computational_methods} and compare the resulting disk structures with the analytical predictions of the model in \secref{sec:sub:comparison_with_simulations}.

\subsection{Computational Methods}
\label{sec:sub:computational_methods}
For the methodology, we largely follow Nakatani} {et~al.}, 2018, Nakatani} {et~al.}, 2018.
Here, we provide a brief overview of the methods employed, referring interested readers to the original papers for more comprehensive details on the computational techniques.

In our simulations, we start with an initially hydrostatic disk (in the poloidal directions) that is exposed to EUV radiation from a central source.
The disk is assumed to be both axisymmetric and midplane symmetric, with an initially isothermal temperature profile in the vertical direction.

The initial density distribution of the hydrostatic disk is
\[
\nh = n_{\rm m}(R)
\exp \left[ -\frac{z^2}{2h^2} \right] , \label{eq:inidenstr}
\]
where $n_{\rm m} (R)$ represents the number density at the midplane ($z= 0$),
and $h$ denotes the pressure scale height defined as the ratio of $\cs$ to the orbital frequency $\Omega \equiv \sqrt{GM_*/R^3}$.
% By integrating \eqnref{eq:inidenstr} w.r.t. $z$, we obtain the surface density profile,
% \[
%             \Sigma (R) =  \sqrt{2\pi} h \rho_{\rm m} (R) ,
% \]
% where $\rho_{\rm m}$ is mass density on the midplane, $\rho_{\rm m} = m  n_{\rm m}$.

The midplane density is determined by specifying the initial surface density profile $\Sigma$, given by
\begin{equation}
\Sigma = \Sigma_0 \braket{\frac{R}{R_0}}^{-1}
\label{eq:sigma}
\end{equation}
with $R_0$ and $\Sigma_0$ being a reference radius and surface density at $R_0$, respectively.
% the steep exponent is reported for CO of a hybrid disk. See Tab.4 in Di Folco +20.
The product of $\Sigma_0$ and $R_0$ is determined by the total disk mass
\begin{equation}
M_{\rm disk} = \int_{R_{\rm min} }^{R_{\rm max}} 2\pi R \Sigma \, \dd R
= 2\pi  (R_{\rm max} - R_{\rm min}) R_0\Sigma_0,
\label{eq:diskmass}
\end{equation}
where $R_{\rm min}$ and $R_{\rm max}$ are the inner and outer truncated radii, respectively, set to the same values as the inner and outer computational boundaries, which will be described below.

The initial temperature profile is given by
\[
T = 100 \Kelvin \braket{\frac{R}{1\au}}^{-1/2}.
\]
This profile is selected based on our previous simulations.
While the choice of this profile may alter the base's geometry, it does not notably influence the comparisons conducted in this section.
The initial value is immediately adjusted to account for the local thermochemical processes within the wind region and remains nearly constant in the stable disk region.
% The results are thus insensitive to the initial condition for $T$.
% As we will see, the resulting mass-loss rates are primarily set by the stellar EUV emission rate and hardly depend on the initial condition for $T$.

The simulations are performed in 2D spherical polar coordinates $(r, \theta)$.
We solve the time evolution of gas density $\rho$,
velocities $\vec{v} = (v_r,~ v_\theta, ~ v_\phi)$,
total gas energy density,
and chemical abundances, using a modified version of PLUTO [Mignone} {et~al.}, 2007]:  % $\{y_i | i = \ce{H, H^+, e-}\}$.
\begin{gather}
\frac{\partial \rho}{\partial t} + \nabla \cdot \rho \vec{v}                                             =  0 ,                  \\
\frac{\partial \rho v_r}{\partial t} + \nabla \cdot \left( \rho v_r \vec{v} \right)              =  -\frac{\partial p}{\partial r}
-\rho \frac{GM_*}{r^2} + \rho \frac{v_\theta^2 + v_\phi^2}{r}                   ,                       \\
\frac{\partial \rho v_\theta}{\partial t} + \nabla \cdot \left( \rho v_\theta \vec{v} \right)    = - \frac{1}{r}\frac{\partial p}{\partial \theta }
- \rho \frac{v_\theta v_r}{r} + \frac{\rho v_\phi^2}{r} \cot \theta                     ,                       \\
\frac{\partial \rho v_\phi}{\partial t} + \nabla ^l \cdot \left( \rho v_\phi \vec{v} \right)     = 0     ,       \label{eq:euler_phi}            \\
\frac{\partial E}{\partial t} + \nabla \cdot \left(H \vec{v} \right)                                     = - \rho v_r \frac{ GM_* }{r^2} +\rho \left( \Gamma_\mathrm{EUV} -\Lambda    \right),                        \\
\frac{\partial \nh y_\mathrm{HI} }{\partial t} + \nabla \cdot \left( \nh y_\mathrm{HI} \vec{v} \right)               = - \nh y_\mathrm{HI} k_\mathrm{ioni} + \alpha\nh^2 y_\mathrm{e} y_\mathrm{HII}       .
\\
\frac{\partial \nh y_\mathrm{HII} }{\partial t} + \nabla \cdot \left( \nh y_\mathrm{HII} \vec{v} \right)               =  \nh y_\mathrm{HI} k_\mathrm{ioni} - \alpha\nh^2 y_\mathrm{e} y_\mathrm{HII}
\\
\text{and} \quad
\frac{\partial \nh y_\mathrm{e} }{\partial t} + \nabla \cdot \left( \nh y_\mathrm{e} \vec{v} \right)               = \nh y_\mathrm{HI} k_\mathrm{ioni} - \alpha\nh^2 y_\mathrm{e} y_\mathrm{HII}          .
\end{gather}
\eqnref{eq:euler_phi} is written in the angular momentum conserving form.
Here $p$ denotes the gas pressure;
and $y_\mathrm{HII}$ and $y_\mathrm{e}$ are the chemical abundances of \ce{H+} and \ce{e-}, respectively;
and $\Lambda$ is the total cooling rate per unit mass (specific cooling rate).
We incorporate radiative recombination cooling and Ly$\alpha$ cooling in the present study.
The total gas energy density and enthalpy density are defined as
\[
\begin{gathered}
E =  \frac{1}{2} \rho v^2 + \frac{p}{\gamma-1}\\
H =  E + p,
\end{gathered}
\]
respectively, with $\gamma$ being the specific heat ratio ($\gamma = 5/3$) and the ideal equation of state.
%        The incorporated chemical reactions are described in \tref{tab:chem_reac}.
% We follow the evolution until the mass-loss rates become quasi-steady.
% It is reached within much less time than the viscous timescale ($\gtrsim 10\Myr$),
% and thus we do not incorporate viscous friction in the conservation equations for energy and angular momentum.
We use $M_\mathrm{disk} = 0.01 \Msun$ and $M_* = 1.0 \Msun$ throughout our simulations.

We run four models varying the EUV spectral hardness and $\phieuv$, covering high- and low-$\phieuv$ ranges for both soft and hard spectra.
We employ $\delta$-function spectra of the form $\delta(h\nu - h\nu_0)$, where $h\nu_0$ is the EUV photon energy.
In the soft-spectrum runs, $h\nu_0 = 14.6\eV$, while in the hard-spectrum runs, $h\nu_0 = 72.9\eV$.
For the soft-spectrum runs, $h\nu_0 = 14.6\eV$, which corresponds to the average energy of a blackbody spectrum with $T_\mathrm{eff} = 10^4\Kelvin$ above the Lyman limit.
For the hard-spectrum runs, we chose $h\nu_0 = 72.9\eV$ arbitrarily for experimental purposes, corresponding to the average energy of a blackbody spectrum with $T_\mathrm{eff} = 3\times 10^5\Kelvin$.
We explore two orders of magnitude in the dimension of the EUV luminosity: $L_\mathrm{EUV} = 10^{28}\erg \sec^{-1}$ and $10^{30}\erg \sec^{-1}$ for the low- and high-luminosity cases, respectively.
The corresponding $\phieuv$ is determined by $L_\mathrm{EUV} = h\nu_0 \phieuv$.
A summary of the input EUV parameters is provided in \tref{tab:numerical_models}.
The Model labels reflect both the magnitude of $L_\mathrm{EUV}$ (high: H, low:L) and spectral hardness (soft: S, hard: H).
\begin{table}[]
\centering
\begin{tabular}{l c c c}
Model &  Emission Rate $\Phi_{40}$ & Photon Energy (eV) &  Luminosity $L_{30}$ \\
\hline
HS   &   4.3 & 14.6 & 1.0 \\
HH   &   0.86 & 72.9 & 1.0 \\
LS     &  $4.3\e{-2}$ & 14.6 & 0.01\\
LH     &  $0.86\e{-2}$ & 72.9 & 0.01 \\
\end{tabular}
\caption{EUV parameters for each run. Here, $\Phi_{40} \equiv \phieuv/10^{40}\sec^{-1}$ and $L_{30} \equiv L_\mathrm{EUV}/10^{30}\erg\sec^{-1}$. In the model labels, the first letter indicates either high (H) or low (L) luminosity, while the second letter signifies either soft (S) or hard (H) spectrum. }
\label{tab:numerical_models}
\end{table}
The EUV emission rate $\phieuv$ of Model~HS is the closest to the supposedly typical value for pre-main-sequence low-mass stars \citep[$\phieuv \approx 10^{40}$--$10^{41}\sec$; e.g.,][]{2014_Alexander}, a value commonly adopted in the literature.
We have also performed runs with luminosities ten times higher than Models~HS and HH. However, these models produced qualitatively similar structures to those in Models~HS and HH, as we will discuss in \secref{sec:sub:comparison_with_simulations}.

Our computational domain spans $r = [0.1, 20]\times R_{\rm g}( = [R_\mathrm{min}, R_\mathrm{max}])$ and $\theta = [0, ~\pi/2] {\rm \, rad}$ with a resolution of $N_r \times N_\theta = 256 \times 160$ cells.
The grid is spaced logarithmically in the radial direction, whereas in the polar direction, two different resolutions are employed. The domain is divided into two regions at $\theta = 1\rad$, with each region having 80 cells uniformly spaced.
We have arrived at this grid resolution after conducting a convergence check for Models~HS and LH. Changing $N_r \times N_\theta$ to $128 \times 160$ and $512\times 320$, respectively, have yielded overall similar structures.
% A high resolution is needed for the disk bulk region to resolve the sharp pressure profile in the vertical direction.

\subsection{Comparisons with the Phenomenological Model}
\label{sec:sub:comparison_with_simulations}
\begin{figure*}[htbp]
%  python3.9 ~/Dropbox/Liouhei_Ryohei/PythonScripts/EUVX.py fig LEUV1e30Rhole0Rmin0.01Rmax2.0Teff1e4Tmid100K/data/data.2000.vtk
% commit 1a751dc9fbe6288ef391e4698c4c663f8e10bce7
% python3.9 ~/Dropbox/Liouhei_Ryohei/PythonScripts/EUVX.py fig LEUV1e30Rhole0Rmin0.01Rmax2.0Teff1e4Tmid100K_highreso/data/data.3000.vtk
% commit cb99b5dde72c9185bf4762ca80f8e0914486457a
\centering
% \includegraphics[clip, width = \linewidth]{figs/ComparisonSnapshots_LEUV1e30Rhole0Rmin0_01Rmax2_0Teff1e4Tmid100K-2000MakeAverageVTSfile_2000-2300.pdf}
\includegraphics[clip, width = \linewidth]{ComparisonSnapshots_LEUV1e30Rhole0Rmin0_01Rmax2_0Teff1e4Tmid100K_highreso-3000MakeAverageVTSfile_2200-2700.pdf}
\caption{
Snapshots of Model~HS, showcasing various physical quantities:
(a--d) number density, gas temperature, isothermal sound speed, and \ion{H}{1} abundance from left to right, respectively.
The green dashed contour represents the isothermal sonic surface,
% while the blue dashed contour indicates where the total enthalpy gets zero, $H-GM_*/r = 0$,
while the light-gray dashed contours indicate where EUV flux is absorption-attenuated by half and 90\% from the original value, i.e., $\chi_\mathrm{e} =\chi_\mathrm{i} = 0.5$ and 0.1, respectively.
These contours roughly demarcate the wind region from the steady disk region.
In panel~(a), cyan contours represents $\nh = 10^3, 10^4, 10^5, 10^6 \cm^{-3}$; in panel~(b), $T = 1000, 5000\Kelvin$; and in panel~(c), $\cs = 5\kms$.
Streamlines are depicted in panel~(d).
Panels~(e) and (i) offer the same snapshots for $\cs$ and $y_\mathrm{HI}$, respectively, but in the $x$--$z/R$ coordinates for comparison with the analytically estimated values through the phenomenological model, presented in panels~(f) and (j).
Panel~(g) illustrates the difference between $\cs$ in the simulation and that predicted by the model, with the magenta contours representing where the analytic estimate falls within a factor of two from $\cs$ in the simulation.
The whitish region indicates where the analytic model accurately predicts $\cs$ (and thereby $T$) of the simulation.
For comparison with the canonical picture, where $\cs = \ceq$ is assumed, panel~(h) similarly displays the difference between $\cs$ in the simulation and $\ceq$.
A comparison between panels~(g) and (h) reveals that our phenomenological model improves the analytical predictions of temperature structures.
The light-gray dashed contours in panels~(e)--(j) also show where $\chi_\mathrm{e} = 0.1$ and 0.5, consistent with panels~(a)--(d).
Panel~(k) presents the difference between the Mach number in the simulation to that predicted through the phenomenological model, with the green dashed contour again representing the isothermal sonic surface.
Lastly, the phase diagram for the spectrum of Model~HS is shown by the lower-right corner panel, with
a black horizontal bar indicates the location of this model on the phase diagram within the range displayed by panels~(a)--(k), i.e., $100\au$.
\label{fig:model_H-S}
\end{figure*}
\begin{figure*}[htbp]
% python3.9 ~/Dropbox/Liouhei_Ryohei/PythonScripts/EUVX.py fig LEUV1e30Rhole0Rmin0.01Rmax2.0Teff3e5Tmid100K/data/data.3500.vtk
% commit cb99b5dde72c9185bf4762ca80f8e0914486457a
\centering
\includegraphics[clip, width = \linewidth]{ComparisonSnapshots_LEUV1e30Rhole0Rmin0_01Rmax2_0Teff3e5Tmid100K-3500.pdf}
\caption{Same as \fref{fig:model_H-S} but for Model~HH. }
\label{fig:model_H-H}
\end{figure*}
\begin{figure*}[htbp]
% python3.9 ~/Dropbox/Liouhei_Ryohei/PythonScripts/EUVX.py fig LEUV1e29Rhole0Rmin0.01Rmax2.0Teff1e4Tmid100K/data/data.2000.vtk
% commit 1a751dc9fbe6288ef391e4698c4c663f8e10bce7
% python3.9 ~/Dropbox/Liouhei_Ryohei/PythonScripts/EUVX.py fig LEUV1e28Rhole0Rmin0.01Rmax2.0Teff1e4Tmid100K_highreso/data/data.3300.vtk
% python3.9 ~/Dropbox/Liouhei_Ryohei/PythonScripts/EUVX.py fig LEUV1e28Rhole0Rmin0.01Rmax2.0Teff1e4Tmid100K_highreso/data/data.10000.vtk
% commit cb99b5dde72c9185bf4762ca80f8e0914486457a
\centering
% \includegraphics[clip, width = \linewidth]{figs/ComparisonSnapshots_LEUV1e29Rhole0Rmin0_01Rmax2_0Teff1e4Tmid100K-2000.pdf}
% \includegraphics[clip, width = \linewidth]{figs/ComparisonSnapshots_LEUV1e28Rhole0Rmin0_01Rmax2_0Teff1e4Tmid100K_highreso-3300.pdf}
\includegraphics[clip, width = \linewidth]{ComparisonSnapshots_LEUV1e28Rhole0Rmin0_01Rmax2_0Teff1e4Tmid100K_highreso-10000.pdf}
\caption{Same as \fref{fig:model_L-S} but for Model~LS.}
\label{fig:model_L-S}
\end{figure*}
\begin{figure*}[htbp]
% python3.9 ~/Dropbox/Liouhei_Ryohei/PythonScripts/EUVX.py fig LEUV1e28Rhole0Rmin0.01Rmax2.0Teff3e5Tmid100K_highreso/data/data.2800.vtk
% commit cb99b5dde72c9185bf4762ca80f8e0914486457a
\centering
% \includegraphics[clip, width = \linewidth]{figs/ComparisonSnapshots_LEUV1e29Rhole0Rmin0_01Rmax2_0Teff3e5Tmid100K-2000.pdf}
\includegraphics[clip, width = \linewidth]{ComparisonSnapshots_LEUV1e28Rhole0Rmin0_01Rmax2_0Teff3e5Tmid100K_highreso-2800.pdf}
\caption{Same as \fref{fig:model_L-H} but for Model~LH.
\label{fig:model_L-H}
\end{figure*}
Figures~\ref{fig:model_H-S}, \ref{fig:model_H-H}, \ref{fig:model_L-S}, and \ref{fig:model_L-H} show snapshots of the simulations (panels~a--d) along with comparisons with analytical predictions by our phenomenological model (panels~e--k) for Models~HS, HH, LS, and LH, respectively.
In these figures, we use the data from the gas structure when it has reached a nearly quasi-steady state.
However, in Model~HS, we employ time-averaged data to smooth out small fluctuations observed in the vicinity of the ionization front.
In the phase diagram located at the lower-right corner of each figure, a black horizontal bar represents the spatial range displayed by panels~(a)--(k) for each model.

These figures directly illustrate the wide variety in the structures of $\nh$, $T$, and $y_\mathrm{HI}$, depending on $\phieuv$ and the spectral hardness, even when they share a common total input energy $L_\mathrm{EUV}$.
Temperatures and isothermal sound speeds (panels~b and c) reach the equilibrium values ($T \approx10^4\Kelvin$ and $\cs\approx 10\kms$) only at relatively small $R$, indicating where the canonical picture of EUV photoevaporation applies.
Temperatures overall decrease at larger distances due to relatively long heating timescales.

Comparing Models~HS and HH, one can observe that H-class spectra tend to result in a larger volume of the high-temperature region, even though $\varphi$ decreases with spectral hardness.
This reflects the fact that a larger fraction of the photon energy is consumed for photoionization with soft spectra. To heat the gas to high temperatures, soft spectra require multiple energy depositions through reionizations.
However, this process is inefficient compared to energy deposition to an atomic gas.

Panels~(e) and (f) demonstrate that our analytic model effectively captures the radial-dependent trend of $\cs$ (and consequently $T$) over a wide range of distances.
In fact, the whitish region in panel~(g), where $\cs^2$ in the simulation is within the factor of two difference from $\cs^2$ predicted by the analytic model, covers a much broader region than the whitish region of panel~(h), which shows the ratio of $\cs^2$ in the simulation to the equilibrium value $\ceq^2$.
This indicates that our model significantly improves the analytical prediction of $\cs$ compared to simply assuming $\cs = \ceq$ as in the canonical picture.
This improvement holds true for a broad range of $\phieuv$ and the spectral hardness, as observed when comparing panels~(e) and (f) of Figures~\ref{fig:model_H-S}, \ref{fig:model_H-H}, \ref{fig:model_L-S}, and \ref{fig:model_L-H}.
The enhancement is more noticeable as the wind temperature decreases with lower $L_\mathrm{EUV}$.

Our phenomenological model also adeptly captures the ionization structures.
In panels~(d), (i), and (j), the greenish region represents the \ion{H}{2} region for each run, with its size varying significantly depending on $\phieuv$ and the spectral hardness (as evident by comparing the four models), as well as the radial distance to the source (as seen in panels~(i) and (j) in \fref{fig:model_L-S}).
As indicated by the analytical model, a harder spectrum tends to result in atomic winds even when $\phieuv$ is at the same levels.
This is implied by comparing, for instance, Figures~\ref{fig:model_H-H} and \ref{fig:model_L-S};
the soft-spectrum run (Model~LS) shows a much more extended \ion{H}{2} region despite having $\sim 20$ times smaller $\phieuv$.
% A larger fraction of photo-energy is consumed by ionizing the gas with soft spectra.

Panels~(i) and (j) show that our analytic model predicts the radial extension of the \ion{H}{2} region and wind types across a broad range of $\phieuv$ and the spectral hardness with reasonable accuracy.
Contrary to the canonical model, where the EUV-heated gas is typically assumed to be completely photoionized ($y_\mathrm{HI}\ll 1$) regardless of the distance to the radiation source and $\phieuv$, our numerical and analytic models highlight the necessity of incorporating the finite timescale of photoionization to discuss the ionization structure, especially for low-mass stars.

In panel~(k), we also compare the Mach number in the simulation, $\mach^\mathrm{sim}$, with the prediction from our phenomenological model, $\mach^\mathrm{ana}$.
% (Recall that in our model, $\mach\sim 1$ for Regimes~A and B, and $\mach \sim \mach_\mathrm{g}^\mathrm{I}$ or $\mach_\mathrm{g}^\mathrm{II}$ for Regime~C; in the canonical picture, $\mach \sim 1$ everywhere.)
Despite some deviations, there is generally a close agreement between the simulation and the analytic estimate, even in the gravity-inhibited subsonic regions (Regime~C) of Model~LH (see \fref{fig:model_L-H}), capturing an overall radial-dependent trend in $\mach$.
This suggests that our model provides a better estimate for $\mach$ across a wider parameter space, compared to the classical assumption of $\mach\sim 1$ everywhere beyond the critical radius.

Nonetheless, there is still room for improvement in estimating $\mach$.
Our model's treatment of a constant $\mach$ fails to fully predict the accelerating nature of photoevaporative winds.
Additionally, discrepancies arise where the ionization front becomes parallel to the radial direction (see, e.g., the $x\lesssim 5$ region in panel~(i) of \fref{fig:model_H-S}).
\footnote{The radiation flux needs to have a perpendicular component to the base to launch winds [Nakatani} \& {Takasao}, 2022].}
% This region belongs to the gravity-inhibited Type~II regime (Regime~C-II), and thus, even if the energy is continuously injected to launch winds, they experience a large pressure drop.
% It leads to a small \ion{H}{1} abundance, and then reduced secondary energy deposition through reionizations. Then gas cannot obtain sufficient energy to escape, resulting in a static atmosphere, especially when the initially injected energy alone is not enough to escape, i.e., $c_1 < \vg $, or equivalently $\varepsilon  x < 1$.
% hard, high, H, heating-oriented
This parallel feature may be artificial, possibly influenced by the computational boundaries; however, it is also observed across convergence-check runs with higher- ($N_r \times N_\theta = 512\times 320$) and lower-resolutions ($N_r \times N_\theta = 128 \times 160$) for Model~HS.
Should the parallel structure persist consistently across various numerical grid configurations, it would suggest that the structure has physically developed by radiation-hydrodynamical effects, at least within the adopted problem setup.
If this is the case, the parallel feature implies that diffuse EUV might dominate over direct EUV in heating the \ion{H}{1} layer here, in contrast to the results of hydrostatic radiative transfer models \citep[e.g.,][]{2013_Tanaka}.
% as the ray-tracing radiative transfer method fails to describe the photon injection in such a configuration.
% However, the parallel shape of the ionization front is likely a numerical effect.
Further comparisons with carefully handled simulations are needed to examine numerical effects and explore the analytical predictability of $\mach$ in such regions, which we plan to address in future studies.
% \begelman{}.

% The development of the gravity-inhibited region may also be important in the context of understanding magneto-thermal disk winds \citep{2016_Bai_a}.
% As gravity-inhibited winds experience a notable pressure drop, the gas's vertical expansion may cause the magnetic field to align more closely with the rotational axis.
% Consequently, this alignment could lower the effectiveness of the Blandford-Payne mechanism.

We have also conducted simulations with luminosities ten times higher ($L_\mathrm{EUV} = 10^{31}\erg\sec^{-1}$) than Models~HS and HH (corresponding $\phieuv$ of $4.3\e{41}\sec^{-1}$ and $0.86\e{41}\sec^{-1}$, respectively).
The resulting structures are qualitatively the same as those of Models~HS and HH other than forming larger isothermal regions (Regime~A winds).
This outcome is expected based on the phase diagrams of Figures~\ref{fig:model_H-S} and \ref{fig:model_H-H} (see the lower-bottom panels).
The disks of Models~HS and HH extend across Regimes~A and B, and increasing the luminosities does not alter this characteristic.

The mass-loss rates (within $R < 100\au$) obtained from the simulations exhibit a wide variation: $2\e{-10}\Msun\yr^{-1}$ for Model~HS, $2\e{-9}\Msun\yr^{-1}$ for Model~HH, $1\e{-11}\Msun\yr$ for Model~LS, and $2\e{-10}\Msun\yr^{-1}$ for Model~LH.
These rates can vary significantly even when the same energy input rate $L_\mathrm{EUV}$ is applied.
In general, models with soft spectra yield mass-loss rates approximately an order of magnitude smaller than those with hard spectra.
This reflects the greater energy consumption by photoionization with soft spectra.
While our current model is not specifically tailored for predicting mass-loss rates, it can be extended to cover this aspect.
To that end, determining the appropriate base density for a wide range of $\phieuv$ and spectral hardness is essential, a task we aim to tackle in future studies.

In summary, our comparisons here highlight a significant enhancement in analytical predictions for the physical structures across a wide range of $\phieuv$ and spectral hardness.
This signifies a broadened understanding of EUV-driven photoevaporation for accretion disks through our phenomenological model.
However, for robustness, similar comparisons over an even broader range of $L_\mathrm{EUV}$, spectral hardness, and spectral shape are warranted, considering uncertainties in recombination timescales and assumptions regarding attenuation factors ($\chi_\mathrm{e}$ and $\chi_\mathrm{i}$) in our model (see discussions in \secref{sec:discussion:model_limitation}).
These analyses would facilitate further model improvements and enhance its practical applications to general EUV spectra.
While such more comprehensive comparisons are deferred to future studies, the brief comparisons provided here suffice to illustrate the remarkable generality of our model.

% \mycomment{2D view. The representative temperature should be interpreted as an average. Temperature will decrease as winds go upwards.}

\section{On the Radial Extension of the H~II region }
\label{sec:limiting_radius}
% \mycomment{Make this section removable later, in case this paper gets too long. This means the other sections should not depend on this section. However, it does not mean this section is unimportant.}

Our model has shown that the canonical picture of EUV photoevaporation, where EUV photoheating drives isothermal ionized winds with temperature and velocity of $\sim 10^4\Kelvin$ and $\sim 10\kms$, respectively, holds only in Regime~A-II.
The corresponding parameter space is formulated as $\varphi > \mathrm{max}(1,q)\varepsilon x$ for the hard-spectrum class and $\varphi > f_\mathrm{eq} x$ for the soft-spectrum classes.
For soft-spectrum classes, ionized escaping winds are also feasible with lower temperatures (Regime~B-II).
These findings suggest the presence of a maximum radial distance up to which launched winds can be photoionized before reaching heights of $\sim R$.
Beyond this distance, the photoionization timescale is longer than the recombination timescale or advection timescale.
Consequently, hydrogen atoms remain only partially photoionized until they ascend to heights of $\sim R$ (cf. panels (d), (i), and (j) in \fref{fig:model_L-S}).
On the scale beyond the maximum distance, the \ion{H}{2} region appears vertically extended rather than radially.
% In essence, the distance limit points where the extension of the \ion{H}{2} region becomes more vertical than radial.
In this section, we quantify this distance limit in physical units and briefly discuss under which conditions the canonical picture of isothermal ionized winds remains valid.

We first consider the H-class.
The boundary between Type~I and Type~II regimes is exclusively present within Regime~A, meaning that the canonical picture is always applicable to ionized winds in this class.
The distance limit to which the \ion{H}{2} region extends is calculated from the Type~II condition (\eqnref{eq:typeII_condition_H_class}) as
\[
\begin{split}
R_\mathrm{lim}    & =   R_\mathrm{c} \frac{\varphi}{\mathrm{max}(1,q)\varepsilon}
=   \frac{\phieuv}{\mathrm{max}(1,q)}   \frac{\avesigma }{4\pi \ceq }
\\
&   \approx 53\au
\braket{\frac{\phieuv}{10^{40}\sec^{-1}}}
\braket{\frac{\avesigma}{10^{-18}\cm^2}}
\\
&   \quad   \braket{\frac{\ceq}{10\kms}}^{-1} \times  \mathrm{min}(1,q^{-1}) .
\end{split}
% In [9]: 1e40 * 1.e-18 / (4 * np.pi*  10e5) / ms.const_au
% Out[9]: 53.194246577984984
\]
The nominal values for the parameters are chosen to be typical for T~Tauri stars.
This suggests that the canonical picture does not necessarily apply to disks around low-mass stars, but there are also atomic winds at the equilibrium temperature (Regime~A-I) or lower (Regime~B-I) for $R > R_\mathrm{lim}$.
For massive stars, whose $\phieuv$ is several orders of magnitude higher than $10^{40}\sec^{-1}$, assuming the canonical picture of an EUV-irradiated disk is justified.
% as discussed in Hollenbach} {et~al.}, 1994.
Since $R_\mathrm{lim}$ is essentially set by the ionization timescale, it is independent of the stellar gravity and is essentially set by $\phieuv$ and the spectral hardness.

The distance limit for the Sa-class is derived similarly as
\[
\begin{split}
R_\mathrm{lim}    & =   R_\mathrm{c} \frac{\varphi}{f_\mathrm{eq}}
=   \frac{\braket{\sqrt{(\varepsilon ^{-1}-1)^2 + 4q} - (\varepsilon ^{-1}-1)}^2}{4 q}   \frac{\phieuv  \avesigma }{4\pi \ceq }
\end{split}
% In [9]: 1e40 * 1.e-18 / (4 * np.pi*  10e5) / ms.const_au
% Out[9]: 53.194246577984984
\]
When $(\varepsilon^{-1} -1 )^2 \gg 4q$, which is often the case for this class, $R_\mathrm{lim}$ is approximated to
\[
R_\mathrm{lim} \approx \frac{q}{(\varepsilon^{-1}-1)^2}\frac{\phieuv  \avesigma }{4\pi \ceq }
= \frac{1}{(\varepsilon^{-1}-1)^2}\frac{3C^2 \phieuv  \alpha_\mathrm{eq} }{4\pi \ceq^2 }
\]
In contrast to the H-class, this class has a regime where freely escaping ionized winds have temperatures lower than the equilibrium value (Regime~B-II).
This motivates us to define the distance limit for the extension of lower-temperature ionized escaping winds (cf. Eqs.\eqref{eq:Sa-class:TypeIIcondition}, \eqref{eq:typeI_regime_Saclass}, or \eqref{eq:typeII_regime_Saclass}),
\[
R_\mathrm{lim}^\mathrm{B}     =   R_\mathrm{c} \frac{\varphi}{\varepsilon^{3/2}}
=   \frac{1}{\varepsilon^{1/2}}   \frac{\phieuv \avesigma }{4\pi \ceq }.
\]
Note that these distance limits are available when $\phieuv$ exceeds the secondary critical EUV emission rate, i.e., $\varphi > f_\mathrm{eq}$.
Hence, for Sa-class spectra, the canonical picture of EUV photoevaporation applies up to $R = R_\mathrm{lim}$, and beyond that, ionized nonisothermal winds extend up to $R\leq R_\mathrm{lim}^\mathrm{B}$ (cf. panels~(d), (i), and (j) in \fref{fig:model_L-S}).

% \mycomment{add the values of the distance limits to \tref{tab:sound_spees}.}

The radius limit of isothermal ionized winds for the Sr-class is the same as that for the Sa-class.
The secondary radius limit is derived as
\[
R_\mathrm{lim}^\mathrm{B}     =   R_\mathrm{c} \frac{\varphi}{(\varepsilon q)^{3/(3+2\beta)}}
=   \frac{\varepsilon}{(\varepsilon q)^{3/(3+2\beta)}}   \frac{\phieuv\avesigma }{4\pi \ceq }.
\]
However, Sr-class might not be achieved with realistic EUV spectra, as mentioned in \secref{sec:sub:Sr-class}.

The presence of the two limiting distances for the soft-spectrum classes suggests the potential to constrain $\avesigma$ and $\aveenergy$ by spatially resolving the radial extensions of the isothermal (Regime~A-II) and nonisothermal (Regime~B-II) ionized winds:
\[
\frac{\varepsilon^{1/2} q}{(\varepsilon^{-1}-1)^2} \approx \frac{R_\mathrm{lim, obs}}{R_\mathrm{lim, obs}^\mathrm{B}}.
\]
This comparison allows for the potential constraint of $\phieuv$ using the constrained spectrum parameters and observed radial extensions.
Analyzing ionized wind tracers with different excitation temperatures could facilitate such analysis.
We further discuss this point in \secref{sec:discussions:wind_observations}.
%, a point we intend to explore further in our future studies.

\section{Work Efficiency}
\label{sec:work_efficiency}
% \mycomment{Make this section removable later, in case this paper gets too long. This means the other sections should not depend on this section. However, it does not mean this section is unimportant.}

In photoevaporating systems, photoheating is responsible for supplying energy to increase the mechanical and thermal energies of the gas while partly being lost through cooling.
The efficiency of this energy conversion process largely depends on the irradiating energy flux.
By employing our model, we can make an order-of-magnitude estimation for this work efficiency and observe its variation with $\phieuv$ by comparing the magnitudes of the photoheating with the work required to liberate the gas from gravitational binding. %, i.e., the power of gravity for the escaping gas.
This efficiency provides a physical reference for understanding the basic energy flow operating within photoevaporative winds in different regimes.

We define work efficiency $\eta$ as the ratio between the integration of the gravitational power per gas mass,
\[
\mathscr{P}_\mathrm{g} \equiv \frac{GM_*}{2r^2} \mach \cs
\]
and the integrated specific heating rate $\Gamma_\mathrm{EUV}$ (Eq.\eqref{eq:gamma_EUV:2}):
% \[
%     \Gamma_\mathrm{EUV } = \frac{y_\mathrm{HI}}{m}\frac{\phieuv}{4\pi r^2 } \avesigma \aveenergy.
% \]
\begin{equation}
\begin{split}
\eta & \equiv  \braket{\int \frac{GM_*}{2r^2}\mach \cs \frac{\dd s}{\mach \cs} } \braket{\int \Gamma_\mathrm{EUV} \frac{\dd s}{\mach \cs}}^{-1} .
% \\
% & \approx \frac{\vg^2}{c_\mathrm{ch}^2} \frac{\mach\cs }{c_\mathrm{ch}}
\end{split}
\label{eq:eta:original}
\end{equation}
Again, the integration is taken along a streamline to the height of $\sim R$.
The factor of $1/2$ in $\mathscr{P}_\mathrm{g}$ effectively accounts for the contribution of centrifugal force to the net gravity.

By the definition of the characteristic sound speed for Type~I winds (\eqnref{eq:c_ch_def}), $\eta$ is approximated to
\begin{equation}
\eta
\approx \eta^\mathrm{I}
\equiv \braket{\frac{\vg}{c_\mathrm{ch}^\mathrm{I}}}^2 \frac{\mach\cs }{c_\mathrm{ch}^\mathrm{I}}
= \varphi^{-1} \mach \frac{\cs}{\ceq}
\label{eq:eta:TypeI}
\end{equation}
The specific value of $\eta$ differs by the wind regimes (see \tref{tab:sound_spees}) and is summarized in \tref{tab:eta}.
\begin{table}[htpb]
\centering
\begin{tabular}{l l R}
Regime  &  State of winds & \eta \\
\hline\hline
A-I  &   Isothermal, atomic &    \varphi^{-1}\\
B-I  &   Steadily-heated, free; atomic & \varphi^{-2/3} x^{-1/3}\\
C-I  &   Gravity-inhibited, atomic&  1\\
\hline
A-II &   Isothermal, ionized& \eta_\mathrm{A}^\mathrm{II}\\
B-II &   Free, ionized& x^{-1}(c_\mathrm{ch}^\mathrm{II}/\ceq)^{-2}\\
C-II &   Gravity-inhibited, ionized &  1 \\
\end{tabular}
\caption{Order-of-magnitude estimation for work efficiency in each wind regime. }
\label{tab:eta}
\end{table}
% For winds in Regime~C-I, $\cs \approx \vg$, and $\mach \sim (c_\mathrm{ch}^\mathrm{I}/\vg)^3$, and thus $\eta \approx 1$.
The efficiency generally decreases as the EUV emission rate increases.
% This is because photo energy is growingly consumed to increase the kinetic plus thermal energy of the gas and is also drained by cooling due to a high temperature.
% The efficiency approaches zero in the extreme case of $\varphi\rightarrow \infty$, indicating that almost all deposited energy is lost via cooling.
% Conversely, work efficiency is the highest in the opposite limiting case $\varphi \rightarrow 0 $, signifying that the gas is heated very slowly, keeping a nearly hydrostatic structure.
Note that for Type~I winds, $\eta$ is set by the local variables, namely $\mathscr{P}_\mathrm{g}$ and $\Gamma_\mathrm{EUV}$, as we assume all the relevant quantities remain constant along the streamline, thus giving $\eta^\mathrm{I} = \mathscr{P}_\mathrm{g}/\Gamma_\mathrm{EUV}$.

For Type~II winds, the work efficiency is approximated to
\[
\begin{split}
\eta
\approx\eta^\mathrm{II}
\equiv  x^{-1}
\left(
\varepsilon   \left[1 - \braket{\frac{x \varepsilon   q}{ \varphi }}^{1/2}
\braket{\frac{\cs}{\ceq }}^{-\beta}\right]
\right.
\\  &   \quad
\left.
+ \braket{\frac{\varphi \varepsilon   q}{ x }}^{1/2}
\braket{\frac{\cs}{\ceq }}^{-(\beta+1)}
\mach^{-1}
\right)^{-1}
\end{split}
\]
(cf. \eqnref{eq:chII_original_def:dimensionless:another_form}).
In Regime~C-II (gravity-inhibited ionized winds), where $\cs\approx\vg$ and $\mach\approx\mach_\mathrm{g}^\mathrm{II}$, the work efficiency is again calculated to be unity, as in Regime~C-I.
In Regime~B-II, using \eqnref{eq:chII_original_def:dimensionless}, $\eta^\mathrm{II}$ is expressed as
\[
\eta^\mathrm{II} = x^{-1}\braket{\frac{c_\mathrm{ch}^\mathrm{II}}{\ceq}}^{-2}.
\]
Given that $c_\mathrm{ch}^\mathrm{II}$ increases with the EUV emission rate, $\eta^\mathrm{II}$ monotonically decreases with increasing $\varphi$.
For Regime~A-II, substituting $\cs\approx\ceq$ and $\mach\approx 1$ into $\eta^\mathrm{II}$, the work efficiency reduces to
\[
\eta^\mathrm{II} =
\eta_\mathrm{A}^\mathrm{II}
\equiv
\left(
\varepsilon
\left[1 - \braket{\frac{x \varepsilon   q}{ \varphi }}^{1/2}
\right]
+ \braket{\frac{\varphi \varepsilon   q}{ x }}^{1/2}
\right)^{-1} .
\]
In the high flux limit of $\varphi/x \gg 1$, the work efficiency approximately scales as $\eta^\mathrm{II}\approx x^{-1/2}\varphi_\mathrm{II}^{-1/2}$.
In comparison  to Type~I winds, the decrease in $\eta$ with respect to $\varphi$ is slower.
This reflects that in this limit, energy is supplied in terms of reionizations, which is limited by recombination timescale proportional to $\varphi^{-1/2}$ (cf. \eqnref{eq:trec_base_density}).

In summary, the work efficiency generally declines as the EUV emission rate increases, indicating that more photo-energy is directed toward the gas's kinetic energy and is also drained through cooling as the temperature rises sufficiently to activate coolants.
In the high flux limit of $\varphi/x \gg 1$, where winds attain an isothermal and ionized state (Regime~A-II), the efficiency approximately scales as $\propto \varphi^{-1/2}$ approaching zero.
This indicates that a significant portion of the energy is lost by cooling, with only a small fraction contributing to the gas's mechanical and thermal energy.
In the opposite limit, where the gas is gravity-inhibited, the efficiency becomes $\eta \approx 1$, meaning efficient conversion of supplied energy into the gas's mechanical and thermal energy.
The highest efficiency is achieved when the gas is heated quasi-statically, keeping a practically hydrostatic structure.

We stress that the derived $\eta$ here is an order-of-magnitude estimation of energy conversion efficiency.
The actual work efficiency defined in the form of \eqnref{eq:eta:original} should not reach unity even in the gravity-inhibited regime, as the deposited energy must partly be consumed to increase the enthalpy of the gas at least.
Nevertheless, the order of $\eta$ and its $\varphi$-dependent trend would remain consistent in more strict estimations.
The overall qualitative behaviors of the work efficiency would be applicable to other astrophysical photoevaporating objects in various scales, such as planetary atmosphere [Murray-Clay} {et~al.}, 2009, Owen} \& {Alvarez}, 2016], molecular clouds, blackhole's accretion disk on the Bondi scale, and galactic minihalos.

\section{Discussions}       \label{sec:discussion}

We have constructed a model primarily tailored for EUV photoevaporation.
However, given that the underlying physics are essentially the same, with suitable modifications, the model can in principle be extended to account for other phenomena such as FUV- and X-ray-driven photoevaporation, as well as the photoevaporation of gas composed purely of metal species.
As an initial demonstration of this versatility, in \secref{sec:discussion:critical_luminosities}, we approximately estimate the critical FUV and X-ray luminosities,
and based on these estimates, we discuss the parameter spaces corresponding to Regimes~A, B, and C of FUV- and X-ray-driven winds.
In \secref{sec:Xray_divergent_conclusion}, We explore how this discussion sheds light on a long-standing issue in the field of protoplanetary disk photoevaporation: the contradictory conclusions on the effectiveness of X-ray photoevaporation.
We also examine the potential time-dependent variations of EUV photoevaporation classes arising from stellar and disk evolution in \secref{sec:discussion:time_dependent_evolution}.
Finally, we address the model limitations and caveats of our model in \secref{sec:discussion:model_limitation}.

\input{FUVX}

\subsection{Possible Time-Dependent Class Variation by Stellar and Disk Evolution}
\label{sec:discussion:time_dependent_evolution}

The classification and dynamics of disk photoevaporation are subject to variations contingent upon the evolutionary stages of both the star and the disk.
Various factors stemming from the evolution of both stellar and disk components wield influence over the classification and dynamics of disk photoevaporation.
This section provides a qualitative exploration of representative factors that contribute to such variations.

The accretion of disk material onto the star can generate strong UV and X-ray emissions by liberating the gravitational energy of the accreting gas.
As for FUV, the accretion-generated component can dominate the stellar component at early stages ($\sim 1\Myr$) where the accretion rate is relatively high [Gorti} {et~al.}, 2009].
The FUV luminosities likely exceed the critical value (\eqnref{eq:critical_FUV_luminosity}) at this stage if the disk still retains sufficiently abundant PAHs/very small grains.
This leads to driving a vigorous FUV photoevaporation.
However, as the system evolves and the accretion rate diminishes, stellar emissions become more dominant.
The luminosities are relatively small compared to the early stages, correspondingly resulting in weaker photoevaporation.
Additionally, if small grains are depleted in later stages, FUV photoevaporation can further attenuate.
% On one hand, photoevaporation may weaken as the disk evolves.

On the other hand, strong accretion implies strong inner MHD winds capable of shielding UV and X-ray photons that could otherwise reach the wind region $x > 1$ \citep[e.g.,][]{2010_Turner, 2018_Takasao, 2022_Takasao, 2018_Fang, 2023_Fang}.
The column density of the inner wind is higher when the system is younger and has a stronger wind.
Only hard photons manage to penetrate the inner region with an attenuated total flux, leading to a hard spectral class and weakened photoevaporation.
% However, as the system evolves,
However, as the system evolves, the surface density reduces relatively quickly at inner radii since MHD winds are more efficient in removing mass there [Bai}, 2016, Suzuki} {et~al.}, 2016].
Then, inner winds and disks become more transparent to UV and X-ray radiation, the spectrum transitions to softer.
This transition facilitates photons reaching the outer region without significant attenuation in the inner region, potentially exciting vigorous winds.
This interaction between radiation and inner winds is supported by observational evidence [Pascucci} {et~al.}, 2020].
Hence, as well as the time-variation of accretion-powered UV and X-ray, accounting for spectral class and wind regime variations is crucial for understanding disk dispersal.

Moreover, stellar evolution, particularly for stars exhibiting time-varying surface magnetic activity, significantly influences the classification and dynamics of disk photoevaporation [Kunitomo} {et~al.}, 2021, Nakatani} {et~al.}, 2023].
For instance, intermediate-mass stars may possess a surface convective zone at ages below $\sim 1\Myr$, resulting in strong EUV and X-ray emissions from magnetic activity.
During this phase, the spectral class is likely H, and strong photoevaporation may ensue due to potentially high UV and X-ray luminosities.
However, as these stars lose their surface convective zones around ages $\sim 1$--$10\Myr$, EUV and X-ray emissions can significantly weaken.
It results in softer EUV emissions from moderately hot photospheres, and photoevaporation is weakened due to a low EUV emission rate. % (\tref{tab:specific_cases}).

In conclusion, the long-term variability of stellar luminosity, stellar activity, disk accretion rate, and grain abundance significantly shape the classification and dynamics of photoevaporation at different stages of disk dispersal.
This indicates that the wind velocity, temperature, and ionization state significantly vary according to the system's evolutionary state,
underscoring the importance of considering temporal variations in photoevaporation for understanding overall disk evolution and interpreting observational data on winds.

% EUV photoevaporation is expected to become effective in the later evolutionary stage and sets the final dispersal.

\subsection{Model Limitation and Caveats}   \label{sec:discussion:model_limitation}
Our model is inherently subject to uncertainties originating from the assumptions we have adopted for simplicity and therefore has caveats to interpret the results.
Here, we will address these points in this section.
We emphasize that, regardless of these model limitations and caveats, our model offers valuable insights into the fundamental physics underlying EUV-driven photoevaporative winds and systematic variations of wind characteristics across a broad parameter space.

\paragraph{Incorporation of Helium photoionization}
Our model focuses solely on EUV heating associated with hydrogen photoionization, neglecting potential contributions from helium photoionization.
While the ionization potentials of \ion{He}{1} and \ion{He}{2} are relatively high (approximately $24.6\eV$ and $54.4\eV$, respectively), their overall impact on the photoheating rate and ionization state of hydrogen in photoevaporative winds is expected to be modest.
Despite the about an order of magnitude higher absorption cross-section per \ion{He}{1}/\ion{He}{2} than the cross-section per \ion{H}{1} in the EUV range,
the elemental helium abundance of $\approx 10\%$, limits the total photoheating rate increment to only a few-fold at most when compared to the heating rate used in our model ($\Gamma_\mathrm{EUV}$; \eqnref{eq:gamma_EUV}).
Similarly, the electron consumption resulting from the recombination of ionic helium species is projected to have a perturbing impact on the ionization state of hydrogen within photoevaporative winds.
Consequently, the order estimation in the present study is expected to remain largely unaltered even with the incorporation of helium contributions.

\paragraph{Self-attenuation effects in photoevaporative winds} We have approximated the absorption attenuation factors, $\chi_\mathrm{e}$ and $\chi_\mathrm{i}$ (see Eqs.\eqref{eq:gamma_EUV} and \eqref{eq:kioni}), as near unity, thereby assuming optically thin photoevaporative winds.
However, especially in cases involving atomic flows, the photoevaporative winds in the inner region can shield to prevent EUV from reaching the outer region.
This self-attenuation effect can reduce the attenuation factors significantly.
In this case, the actual flux reaching the disk is smaller than the value used in our model for a given $\phieuv$, and thus the actual heating timescale can correspondingly become longer if we incorporate corrections due to self-attenuation.
This would lead to an increase in the critical EUV emission rate $\Phi_\mathrm{c}$ by several times and corresponding adjustments in the boundaries of wind regimes on the phase diagram.
Hence, the derived critical EUV emission rate should be interpreted as a lower limit, suggesting that the derived radius limits presented in \secref{sec:limiting_radius} should be considered upper limits.

Since the self-attenuation is expected to be increasingly pronounced at outer radii, the actual boundary between Regimes~A and B may not adhere to a simple linear relation as depicted in our model but rather may be described by a convex downward function.
A similar is true for the boundary between Regimes~B and C; the gravity-inhibited regime can be somewhat radially extended than the predictions of our model.
% Additionally, the actual boundaries between wind regimes may undergo radial changes depending on the magnitude of wind attenuation.

Additionally, in the 2D cross-sectional view of a photoevaporating disk, the wind region would exhibit a vertically stratified structure with Regimes~A, B, and C present in distinct layers.
This stratification arises due to the irradiating spectrum becoming harder as it is attenuated.
In the most attenuated layer, both the photon number flux and the average cross-section are small, resulting in the layer being in a gravity-inhibited regime.
In a moderately attenuated layer, the photoheated layer can exhibit characteristics of Regimes~A and B if $\phieuv$ is sufficiently high.
We will explore this aspect in more detail in \secref{sec:comparison_HD_simulations}.

% In future studies, comparisons with hydrodynamics simulations are necessary to thoroughly investigate the impacts of wind attenuation on our findings.
% \mycomment{
% The approximation of $\chi_\mathrm{e} \sim 1$ may be valid only when the spectrum has a delta-function-like shape as $\Phi_\nu = \phieuv \delta (\nu - \nu_2)$; otherwise, wind regions can include $\chi_\mathrm{e} \ll 1$ when they are atomic. % If , $\aveenergy = h (\nu_2-\nu_1)$.
% When the spectrum has a certain width, $\chi_\mathrm{e}$ decays in a power-law manner rather than exponentially.
% }

\paragraph{Inclusion of neutral coolants}
Our model primarily relies on coolants effective in \ion{H}{2} regions, such as forbidden lines of ionic metals, Ly$\alpha$, and radiative recombination, which results in an equilibrium temperature of $\approx10^4\Kelvin$ for both atomic and ionized gas.
However, as also discussed in \secref{sec:discussion:critical_luminosities}, atomic and molecular cooling, such as infrared emission from \ion{O}{1} and \ce{H2}, can be significant in the neutral layer of a disk.
This can lead to a substantially lower equilibrium temperature on the order of $\sim 10^2$--$ 10^3\Kelvin$.
Thus, the equilibrium temperature would, in principle, depend on whether the gas is in a neutral or ionized state, which is determined by the luminosity of the ionizing radiation.

This realization suggests the existence of a new Type~I wind regime characterized by winds exhibiting equilibrium temperatures typical of neutral gas.
This regime would fall within the Type~I regime and bridge connections to the Type~II regime at the high-flux parameter space.
The prevalence of this regime is expected, particularly for H-class EUV spectra, which yield a low $\avesigma$, creating the base within the deep interior of the neutral layer.
Therefore, an important next step for the phenomenological model is to incorporate the effects of neutral coolants to formulate the new Type~I regime and determine the corresponding parameter space within the phase diagram.

\paragraph{Investigation of the inner region}
While our study primarily focuses on the outer radius $x > 1$, where vigorous winds are anticipated due to weak gravity, investigating the inner region $x < 1$ is also necessary.
This disk corona, i.e., nearly hydrostatic gas in $x <1$, has the potential to produce weak winds whose escaping rate exponentially decays inward, potentially serving as a source for time variation.
Investigating this region could give valuable insights into achievable velocities and temperatures therein, which could provide surface mass-loss rates available for population synthesis studies and be used to discuss the observability of infrared wind tracers.
% Investigating this unexplored parameter space is also beneficial for studying possible achievable velocity and temperature of the inner region, which could be observable using infrared wind tracers.

\paragraph{Uncertainty in recombination timescale}
There remains uncertainty in the recombination timescale used in this study (\eqnref{eq:trec_base_density}), particularly regarding the applicability of the base density formulation (\eqnref{eq:base_density}) across different EUV emission rates and spectra.
If the base density has another formulation within certain parameter spaces, it could impact the boundary of Types~I and II regimes.
Although our key findings --- e.g., the presence of two wind types with Type~II being achieved at a relatively high $\phieuv$ --- qualitatively remain unchanged by such modifications, it is necessary to examine the general applicability of \eqnref{eq:base_density} by comparing the base density with numerical calculations.

Additionally, our model has not accounted for the density reduction due to the expansion of evaporating gas, leading to an overestimation of additional energy supply through reionizations within the ionized winds.
Addressing this is another crucial aspect that needs to be incorporated into an updated model.

\paragraph{Integration with MHD winds}
Magnetohydrodynamics (MHD) winds are another important disk-dispersal process, especially during the early stage of protoplanetary disk evolution when the disk is relatively massive \citep[e.g.,][see also \secref{sec:discussion:time_dependent_evolution}]{2014_Turner, 2022_Pascucci}.
Our model could be applied to the determination of MHD winds' thermochemical states if the winds have sufficiently low densities, allowing photochemical processes to dominate other heating sources \citep[e.g., ambipolar diffusion ][]{2001_Garciab, 2001_Garcia}.
The velocity achieved by MHD winds would set whether the gas can reach the equilibrium temperature.
In high-velocity winds, the gas may not have enough time to be heated, resulting in relatively low temperatures.
Thus, corresponding modifications are necessary to the definition of $c_\mathrm{ch}$ (Eq.\eqref{eq:c_ch_def}) for the temperature estimation with accounting for this effect.

Similarly, since the base density of MHD winds is not related to the photoionization and recombination rates in general, the estimation of the recombination timescale also needs to be modified.
With the updated recombination timescale, the same criteria (Eqs.\eqref{eq:typeI_condition} and \eqref{eq:typeII_condition}) would still be applicable to infer the ionization state of the winds.
Overall, while our model is not directly applicable to MHD winds, it could be adapted by treating other heating processes, MHD winds' velocities, and base densities.
% - streamlines are determined by MHD winds but the same temperature estimation techniques would be available although the flow speed may not be determined thermally

\subsection{Prospects for Potential Applications to Photoevaporative Wind Observations}
\label{sec:discussions:wind_observations}

Since the main scope of this paper is to broaden our understanding of EUV-driven photoevaporation in a physically consistent manner, we will leave a detailed discussion on the observational applications and tests of our model for future studies.
Nevertheless, we describe here the potential utility of our model in investigating photoevaporative winds and stellar parameters and outline the steps that need to be addressed to that end.

We have demonstrated that EUV-driven winds have a temperature gradient in the radial direction, even within the \ion{H}{2} region.
When the spectrum is soft, such as those softer than a blackbody spectrum with an effective temperature of $T_\mathrm{eff} \lesssim 5\e{4}\Kelvin$, the \ion{H}{2} region extends across the isothermal and non-isothermal wind regimes (Regimes~A-II and B-II).
On the other hand, with hard spectra, the thermal balance dictates the temperature of the \ion{H}{2} region, resulting in isothermal winds.
In both cases, there is a maximum radius beyond which the \ion{H}{2} region cannot extend (\secref{sec:limiting_radius}), and the winds transition to being atomic with lower temperatures.

These findings suggest that the observed size of the ionized region size can vary depending on the ionized gas tracers used, which have different excitation temperatures.
By comparing the sizes of the ionized region using different tracers, we can estimate the average EUV energy of the source and then determine the source's EUV emission rate $\phieuv$, following a similar procedure as outlined in \secref{sec:limiting_radius}.

The expected structure of emissions from various tracers can be outlined as follows.
When the EUV spectrum is relatively hard, high-temperature ($\approx 10^4 \Kelvin$) ionized gas tracers (e.g., [\ion{Ne}{2}]\,12.81$\mum$, [\ion{O}{2}] and [\ion{O}{3}] optical lines) exhibit a compact emission (Regime~A-II).
This region is surrounded by a more extended emission of neutral gas tracers (e.g., [\ion{O}{1}]\,6300\AA, [\ion{O}{1}]\,63$\mum$, and [\ion{C}{2}]\,158$\mum$) with a decrease in temperature towards the radial direction (Regimes~A-I and B-I).
In contrast, when the spectrum is soft, the high-temperature ionized region is encircled by more extended emissions from intermediate-temperature ($\approx 5000$--$10^4\Kelvin$) ionized gas tracers (e.g., [\ion{S}{2}] and [\ion{N}{2}] optical lines; Regime~B-II),
followed by lower-temperature neutral gas tracers extending further outward (Regime~B-I).
Note that $\phieuv$ must exceed the critical EUV emission rates (cf. Eqs.\eqref{eq:critical_luminosity} and \eqref{eq:secondary_critical_luminosity}) to yield the high-temperature ionized region.

Recent JWST observations have revealed relatively compact [\ion{Ar}{2}] emissions compared to [\ion{Ne}{2}] and [\ion{Ne}{3}] emissions for T~Cha [Bajaj} {et~al.}, 2024, Sellek} {et~al.}, 2024].
These differences in emission sizes, however, are concluded to be likely due to distinct photoionization processes for each tracer, specifically EUV ionization versus X-ray ionization.
Considering that low-mass stars commonly exhibit X-rays, generalizing our model to the X-ray wavelengths appears to be a crucial step to making it available for inferring spectra and luminosities of low-mass stars based on the different emission extensions of wind tracers.

Nevertheless, our model can further be expanded to develop an analytic model for X-ray photoevaporation, as discussed in \secref{sec:discussion:critical_luminosities}.
By integrating our EUV model with such an X-ray photoevaporation model, an analytical approach can be formulated to estimate the line intensities of ionized gas tracers, similar to the approach by Hollenbach} \& {Gorti}, 2009.
Such a model provides analytical predictions generalized by incorporating hydrodynamical effects on temperature and ionization state determination,
thereby updating the currently existing models that assume the equilibrium temperatures determined by thermochemical balance.

Additionally, those models are of great use for investigating free-free emissions observed towards young sources and estimating $\phieuv$ based on the observed free-free flux density [Pascucci} {et~al.}, 2012, Ricci} {et~al.}, 2021].
The resulting spatial distribution of free-free emissions is valuable for interpreting the spatially resolved data of future ngVLA observations and exploring the origins of these emissions.

% ultraviolet spectroscopy? -> usually needs a high temperature than the typical HII region, like > 30000K
% See my onenote: Wind Observation -> Observations -> Ultraviolte Lines

\section{Summary and Prospects}   \label{sec:summary}

In the canonical picture of EUV photoevaporation, the photoevaporative winds consist of ionized gas at equilibrium temperatures $\approx 10^4\Kelvin$, corresponding to wind speeds $\approx 10\kms$ beyond the gravitational radius [Hollenbach} {et~al.}, 1994] or critical radius [Liffman}, 2003].
We have delved into the conditions under which this canonical picture holds true, constructing a phenomenological model based on the analytic approach of Begelman} {et~al.}, 1983.
The resulting generalized picture applies across the whole range of conditions where photoevaporation is possible, and yields estimates for the temperature, speed, and ionization state that EUV-driven disk winds exhibit as functions of four parameters: EUV luminosity, spectral hardness, distance from the radiation source, and the source's mass.
These are the main findings:
\begin{enumerate}
\item   As in the model of Begelman} {et~al.}, 1983, photoevaporative winds occur in three hydrodynamical regimes: isothermal (Regime~A); free-flowing (Regime~B), and gravity-inhibited (Regime~C).
Isothermal winds occur when the heating timescale is shorter than both the wind crossing and gravitational timescales, whereas gravity-inhibited winds result when the heating timescale is longer than the gravitational timescale.
In all other cases, the winds transition into free winds.

\item The canonical picture applies only in a region of the parameter space characterized by high EUV fluxes.
This parameter region is denoted Regime~A-II and is represented by red shading without vertical gray stripes in the phase diagrams of Figures~\ref{fig:H-class}, \ref{fig:Sa-class}, \ref{fig:Sr-class}, and \ref{fig:classmaps}, as summarized in \tref{tab:sound_spees}.

\item   To form isothermal winds, the EUV emission rate must at least exceed the critical value (\eqnref{eq:critical_luminosity}), which is defined by the stellar mass, equilibrium temperature, and spectral hardness.

\item For lower EUV fluxes, the winds are cooler than the equilibrium temperature due to a longer heating timescale.
This results in two distinct regimes:
\begin{itemize}
\item At larger disk radii where gravity is weak, the low-temperature winds escape freely without significant gravitational influence (Regime~B; green and orange regions in the figures).
\item Close-in, gravity is strong and inhibits the winds, leading to a substantial pressure drop (Regime~C; blue regions in the figures).
\end{itemize}

\item The phase diagrams divide into three classes depending on the hardness of the EUV spectrum, as parameterized by $\varepsilon$ and $q$ (Eqs.~\eqref{eq:tildes} and \eqref{eq:q}).
\begin{itemize}
\item   {\it H-class}: Characterized by hard spectra, depositing a relatively large energy per photoionization.
Ionized winds in this class are always characterized by isothermal flows at the equilibrium temperature (Regime~A-II), with the possibility of isothermal atomic winds in a limited parameter space (Regime~A-I).
Spectra from young low-mass stars with strong surface magnetic activity are expected to fall into this class.

\item   {\it Sa-class}:
One of the two soft-spectrum classes, where advection is more effective in maintaining atomic hydrogen in the winds than recombination.
With little energy deposited per photoionization, multiple reionizations are necessary before the gas heats to the equilibrium temperature.
This results in cooler ionized winds that either escape freely (Regime~B-II) or are gravity-inhibited (Regime~C-II).
Spectra in this class may be rarer than the H-class spectra, with main-sequence intermediate-mass stars potentially falling in the Sa class.
% Young stars dominated by accretion-produced soft EUV emission could also be in this class.

\item   {\it Sr-class}: Same as Sa-class, but with recombination more effective than advection in maintaining atomic hydrogen in the winds.
The key distinction is that Sa-class has a parameter region where energy deposition through reionizations is negligible in ionized winds (see \secref{sec:sub:Sa-class}).
Sr-class appears to occur rarely if ever in the context of protoplanetary disk dispersal.
\end{itemize}
\end{enumerate}

We have evaluated the model's predictions by comparing them against hydrodynamics simulations treating the radiative transfer and disequilibrium thermochemistry (\secref{sec:comparison_HD_simulations}).
There is broad concurrence between the analytic model and the numerical results across orders of magnitude in the EUV emission rate and a wide range in spectral hardness.  Refining the model and ensuring its applicability across various EUV environments will nevertheless require more exhaustive assessments with diverse EUV spectra and flow geometries.

To demonstrate the practical utility of the derived typical temperature and speed of the winds,
we have estimated the limiting radius of the \ion{H}{2} region to which ionized winds can extend (\secref{sec:limiting_radius}).
We have shown that this limiting radius is determined by the EUV emission rate, average absorption cross-section, isothermal sound speed at the equilibrium temperature, and spectral hardness.
Our analysis suggests that the EUV emission rate can potentially be constrained by comparing the radial sizes of the \ion{H}{2} regions observed in ionized gas tracers with different excitation temperatures --- an avenue we intend to explore in future investigations.
Furthermore, our analysis has revealed that while applying the canonical picture of EUV photoevaporation to disks around massive stars is well justified, the same does not hold for disks around young low-mass stars.
Therefore, when interpreting observations of ionized wind tracers or conducting numerical simulations of photoevaporation for low-mass stars, it is crucial to take into account the finite timescales of photoheating and photoionization.

We have also estimated the work efficiency, or ratio of the work necessary to move the gas against gravity to the energy deposited in that gas through photoheating (\secref{sec:work_efficiency}).
We find that the efficiency is of order unity for gravity-inhibited winds and is highest when the gas is quasi-statically heated.
Conversely, the efficiency drops as the EUV emission rate increases and more energy is directed towards gas kinetic and thermal energy or is lost to radiative cooling.
At the highest fluxes, the efficiency scales as $\phieuv^{-1/2}$, where $\phieuv$ represents the rate at which the source emits EUV photons.
The scaling in the efficiency aligns with that found by numerical simulations in the context of planetary atmospheric escape [Murray-Clay} {et~al.}, 2009, Owen} \& {Alvarez}, 2016],
suggesting that our results could be applied to understanding EUV photoevaporation across various classes of astrophysical objects (Mitani et al. in prep).

The analytic model has the potential to be extended to FUV and X-ray photoevaporation with suitable modifications, providing a tool to explore the effectiveness of these energy bands in disk dispersal over an even wider parameter space (\secref{sec:discussion:critical_luminosities}).
Additionally, the model has potential to explain divergent conclusions in the literature regarding X-ray photoevaporation (\secref{sec:Xray_divergent_conclusion}).
We acknowledge that uncertainties persist in the analysis, leaving room for further refinements, including those discussed in \secref{sec:discussion:model_limitation}.
Specifically, it is crucial to validate the predictions by comparing them against numerical simulations that treat the radiative transfer and time-dependent photochemistry.
Additionally, it will be essential to update the model to accommodate cases where the spectrum covers wider ranges in photon energy.
With these caveats in mind, we have discussed prospected applications of our model to observations \secref{sec:discussions:wind_observations}.

Since the model is analytic, it can be extended to enable exploring the detectability of photoevaporative winds and deriving mass-loss rates across a wide parameter space, which we hope will yield fresh insights into disk dispersal from both observational and theoretical perspectives.
Moreover, the model provides a solid foundation for understanding photoevaporation in astrophysical objects including planets, molecular clouds, and galactic mini-halos.

% The model gives important principles for evaporating objects.
% application range of Hollenbach} {et~al.}, 1994 models.
% Wind detectability
% Mass loss rates.
% \section{Mass-Loss Profiles}
% \input{subsections/mass_loss_profiles}

% \input{subsections/base_pressure}
% \input{subsections/base_analytic}
% \section{Comparison with Numerical Simulations}
% \input{subsections/simulation}
% \input{subsections/gravitational_radius}

\acknowledgments %
We thank the reviewer for the practical comments.
R.N. is supported by the Japan Society for the Promotion of Science (JSPS), Overseas Research Fellowship.
% The numerical computations were carried out on the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.
S.T. was supported by the JSPS KAKENHI grant Nos. JP21H04487, JP22K14074, and JP22KK0043.
This research was performed in part at the Jet Propulsion Laboratory, California Institute of Technology, under contract 80NM0018D0004 with the National Aeronautics and Space Administration and with the support of the NASA Exoplanets Research Program through grant 17-XRP17\_2-0081 to N.J.T.
Numerical computations were carried out on the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.
% This work was strongly benefited from the Core2disk-III residential program of Institut Pascal at Universit\'e Paris-Saclay, with the support of the program ``Investissements d'avenir'' ANR-11-IDEX-0003-01.
\copyright 2024. All rights reserved.

\vspace*{1cm}
\software{Numpy [Harris {et~al.}, 2020], Matplotlib [Hunter, 2007], Astropy [Astropy Collaboration} {et~al.}, 2013, Astropy Collaboration} {et~al.}, 2018, Astropy Collaboration} {et~al.}, 2022], SciPy [Virtanen {et~al.}, 2020]
\vspace*{5cm}

% \mycomment{Do not be critical to other works. Do not use the words, good, bad, etc. Just describe and list objective facts. }

%\software{ CASA \citep{McMullin2007}}

\bibliographystyle{aasjournal}
% \begin{thebibliography}{}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\providecommand{\url}[1]{\href{#1}{#1}}
\providecommand{\dodoi}[1]{doi:~\href{http://doi.org/#1}{\nolinkurl{#1}}}
\providecommand{\doeprint}[1]{\href{http://ascl.net/#1}{\nolinkurl{http://ascl.net/#1}}}
\providecommand{\doarXiv}[1]{\href{https://arxiv.org/abs/#1}{\nolinkurl{https://arxiv.org/abs/#1}}}

[Alexander} {et~al.}, 2014]
{Alexander}, R., {Pascucci}, I., {Andrews}, S., {Armitage}, P., \& {Cieza}, L.
2014, in Protostars and Planets VI, ed. H.~{Beuther}, R.~S. {Klessen}, C.~P.
{Dullemond}, \& T.~{Henning}, 475,
\dodoi{10.2458/azu_uapress_9780816531240-ch021}

[Alexander} {et~al.}, 2004]
{Alexander}, R.~D., {Clarke}, C.~J., \& {Pringle}, J.~E. 2004, \mnras, 354, 71,
\dodoi{10.1111/j.1365-2966.2004.08161.x}

[Alexander} {et~al.}, 2005]
---. 2005, \mnras, 358, 283, \dodoi{10.1111/j.1365-2966.2005.08786.x}

[Astropy Collaboration} {et~al.}, 2013]
{Astropy Collaboration}, {Robitaille}, T.~P., {Tollerud}, E.~J., {et~al.} 2013,
\aap, 558, A33, \dodoi{10.1051/0004-6361/201322068}

[Astropy Collaboration} {et~al.}, 2018]
{Astropy Collaboration}, {Price-Whelan}, A.~M., {Sip{\H{o}}cz}, B.~M., {et~al.}
2018, \aj, 156, 123, \dodoi{10.3847/1538-3881/aabc4f}

[Astropy Collaboration} {et~al.}, 2022]
{Astropy Collaboration}, {Price-Whelan}, A.~M., {Lim}, P.~L., {et~al.} 2022,
apj, 935, 167, \dodoi{10.3847/1538-4357/ac7c74}

[Bai}, 2016]
{Bai}, X.-N. 2016, \apj, 821, 80, \dodoi{10.3847/0004-637X/821/2/80}

[Bajaj} {et~al.}, 2024]
{Bajaj}, N.~S., {Pascucci}, I., {Gorti}, U., {et~al.} 2024, \aj, 167, 127,
\dodoi{10.3847/1538-3881/ad22e1}

[Bakes} \& {Tielens}, 1994]
{Bakes}, E.~L.~O., \& {Tielens}, A.~G.~G.~M. 1994, \apj, 427, 822,
\dodoi{10.1086/174188}

[Begelman} {et~al.}, 1983]
{Begelman}, M.~C., {McKee}, C.~F., \& {Shields}, G.~A. 1983, \apj, 271, 70,
\dodoi{10.1086/161178}

[Bertoldi}, 1989]
{Bertoldi}, F. 1989, \apj, 346, 735, \dodoi{10.1086/168055}

[Bertoldi} \& {McKee}, 1990]
{Bertoldi}, F., \& {McKee}, C.~F. 1990, \apj, 354, 529, \dodoi{10.1086/168713}

[Ercolano} {et~al.}, 2009]
{Ercolano}, B., {Clarke}, C.~J., \& {Drake}, J.~J. 2009, \apj, 699, 1639,
\dodoi{10.1088/0004-637X/699/2/1639}

[Ercolano} {et~al.}, 2008]
{Ercolano}, B., {Drake}, J.~J., {Raymond}, J.~C., \& {Clarke}, C.~C. 2008,
\apj, 688, 398, \dodoi{10.1086/590490}

[Fang} {et~al.}, 2023]
{Fang}, M., {Pascucci}, I., {Edwards}, S., {et~al.} 2023, \apj, 945, 112,
\dodoi{10.3847/1538-4357/acb2c9}

[Fang} {et~al.}, 2018]
---. 2018, \apj, 868, 28, \dodoi{10.3847/1538-4357/aae780}

[Font} {et~al.}, 2004]
{Font}, A.~S., {McCarthy}, I.~G., {Johnstone}, D., \& {Ballantyne}, D.~R. 2004,
\apj, 607, 890, \dodoi{10.1086/383518}

[Garcia} {et~al.}, 2001]
{Garcia}, P.~J.~V., {Cabrit}, S., {Ferreira}, J., \& {Binette}, L.
2001{\natexlab{a}}, \aap, 377, 609, \dodoi{10.1051/0004-6361:20011146}

[Garcia} {et~al.}, 2001]
{Garcia}, P.~J.~V., {Ferreira}, J., {Cabrit}, S., \& {Binette}, L.
2001{\natexlab{b}}, \aap, 377, 589, \dodoi{10.1051/0004-6361:20011145}

[Geers} {et~al.}, 2007]
{Geers}, V.~C., {van Dishoeck}, E.~F., {Visser}, R., {et~al.} 2007, \aap, 476,
279, \dodoi{10.1051/0004-6361:20078466}

[Gorti} {et~al.}, 2009]
{Gorti}, U., {Dullemond}, C.~P., \& {Hollenbach}, D. 2009, \apj, 705, 1237,
\dodoi{10.1088/0004-637X/705/2/1237}

[Gorti} \& {Hollenbach}, 2004]
{Gorti}, U., \& {Hollenbach}, D. 2004, \apj, 613, 424, \dodoi{10.1086/422406}

[G{\"u}del} {et~al.}, 2007]
{G{\"u}del}, M., {Briggs}, K.~R., {Arzner}, K., {et~al.} 2007, \aap, 468, 353,
\dodoi{10.1051/0004-6361:20065724}

[Habing}, 1968]
{Habing}, H.~J. 1968, \bain, 19, 421

[Harris {et~al.}, 2020]
Harris, C.~R., Millman, K.~J., van~der Walt, S.~J., {et~al.} 2020, Nature, 585,
357, \dodoi{10.1038/s41586-020-2649-2}

[Hollenbach} \& {Gorti}, 2009]
{Hollenbach}, D., \& {Gorti}, U. 2009, \apj, 703, 1203,
\dodoi{10.1088/0004-637X/703/2/1203}

[Hollenbach} {et~al.}, 1994]
{Hollenbach}, D., {Johnstone}, D., {Lizano}, S., \& {Shu}, F. 1994, \apj, 428,
654, \dodoi{10.1086/174276}

[Hollenbach} {et~al.}, 1993]
{Hollenbach}, D., {Johnstone}, D., \& {Shu}, F. 1993, in Astronomical Society
of the Pacific Conference Series, Vol.~35, Massive Stars: Their Lives in the
Interstellar Medium, ed. J.~P. {Cassinelli} \& E.~B. {Churchwell}, 26

[Hunter, 2007]
Hunter, J.~D. 2007, Computing in Science \& Engineering, 9, 90,
\dodoi{10.1109/MCSE.2007.55}

[Iliev} {et~al.}, 2005]
{Iliev}, I.~T., {Shapiro}, P.~R., \& {Raga}, A.~C. 2005, \mnras, 361, 405,
\dodoi{10.1111/j.1365-2966.2005.09155.x}

[Kahn}, 1954]
{Kahn}, F.~D. 1954, \bain, 12, 187

[Komaki} {et~al.}, 2021]
{Komaki}, A., {Nakatani}, R., \& {Yoshida}, N. 2021, \apj, 910, 51,
\dodoi{10.3847/1538-4357/abe2af}

[Kunitomo} {et~al.}, 2021]
{Kunitomo}, M., {Ida}, S., {Takeuchi}, T., {et~al.} 2021, \apj, 909, 109,
\dodoi{10.3847/1538-4357/abdb2a}

[Liffman}, 2003]
{Liffman}, K. 2003, PASA, 20, 337, \dodoi{10.1071/AS03019}

[Maeder} \& {Meynet}, 1987]
{Maeder}, A., \& {Meynet}, G. 1987, \aap, 182, 243

[Maloney} {et~al.}, 1996]
{Maloney}, P.~R., {Hollenbach}, D.~J., \& {Tielens}, A.~G.~G.~M. 1996, \apj,
466, 561, \dodoi{10.1086/177532}

[McElroy} {et~al.}, 2013]
{McElroy}, D., {Walsh}, C., {Markwick}, A.~J., {et~al.} 2013, \aap, 550, A36,
\dodoi{10.1051/0004-6361/201220465}

[Mignone} {et~al.}, 2007]
{Mignone}, A., {Bodo}, G., {Massaglia}, S., {et~al.} 2007, \apjs, 170, 228,
\dodoi{10.1086/513316}

[Murray-Clay} {et~al.}, 2009]
{Murray-Clay}, R.~A., {Chiang}, E.~I., \& {Murray}, N. 2009, \apj, 693, 23,
\dodoi{10.1088/0004-637X/693/1/23}

[Nakatani} {et~al.}, 2020]
{Nakatani}, R., {Fialkov}, A., \& {Yoshida}, N. 2020, \apj, 905, 151,
\dodoi{10.3847/1538-4357/abc5b4}

[Nakatani} {et~al.}, 2018]
{Nakatani}, R., {Hosokawa}, T., {Yoshida}, N., {Nomura}, H., \& {Kuiper}, R.
2018{\natexlab{a}}, \apj, 857, 57, \dodoi{10.3847/1538-4357/aab70b}

[Nakatani} {et~al.}, 2018]
---. 2018{\natexlab{b}}, \apj, 865, 75, \dodoi{10.3847/1538-4357/aad9fd}

[Nakatani} \& {Takasao}, 2022]
{Nakatani}, R., \& {Takasao}, S. 2022, \apj, 930, 124,
\dodoi{10.3847/1538-4357/ac63a0}

[Nakatani} {et~al.}, 2023]
{Nakatani}, R., {Turner}, N.~J., {Hasegawa}, Y., {et~al.} 2023, \apjl, 959,
L28, \dodoi{10.3847/2041-8213/ad0ed8}

[Nakatani} \& {Yoshida}, 2019]
{Nakatani}, R., \& {Yoshida}, N. 2019, \apj, 883, 127,
\dodoi{10.3847/1538-4357/ab380a}

[Nomura} {et~al.}, 2007]
{Nomura}, H., {Aikawa}, Y., {Tsujimoto}, M., {Nakagawa}, Y., \& {Millar}, T.~J.
2007, \apj, 661, 334, \dodoi{10.1086/513419}

[Owen} \& {Alvarez}, 2016]
{Owen}, J.~E., \& {Alvarez}, M.~A. 2016, \apj, 816, 34,
\dodoi{10.3847/0004-637X/816/1/34}

[Owen} {et~al.}, 2010]
{Owen}, J.~E., {Ercolano}, B., {Clarke}, C.~J., \& {Alexander}, R.~D. 2010,
\mnras, 401, 1415, \dodoi{10.1111/j.1365-2966.2009.15771.x}

[Pascucci} {et~al.}, 2023]
{Pascucci}, I., {Cabrit}, S., {Edwards}, S., {et~al.} 2023, in Astronomical
Society of the Pacific Conference Series, Vol. 534, Protostars and Planets
VII, ed. S.~{Inutsuka}, Y.~{Aikawa}, T.~{Muto}, K.~{Tomida}, \& M.~{Tamura},
567, \dodoi{10.48550/arXiv.2203.10068}

[Pascucci} {et~al.}, 2012]
{Pascucci}, I., {Gorti}, U., \& {Hollenbach}, D. 2012, \apjl, 751, L42,
\dodoi{10.1088/2041-8205/751/2/L42}

[Pascucci} {et~al.}, 2014]
{Pascucci}, I., {Ricci}, L., {Gorti}, U., {et~al.} 2014, \apj, 795, 1,
\dodoi{10.1088/0004-637X/795/1/1}

[Pascucci} {et~al.}, 2020]
{Pascucci}, I., {Banzatti}, A., {Gorti}, U., {et~al.} 2020, \apj, 903, 78,
\dodoi{10.3847/1538-4357/abba3c}

[Picogna} {et~al.}, 2019]
{Picogna}, G., {Ercolano}, B., {Owen}, J.~E., \& {Weber}, M.~L. 2019, \mnras,
487, 691, \dodoi{10.1093/mnras/stz1166}

[Ricci} {et~al.}, 2021]
{Ricci}, L., {Harter}, S.~K., {Ercolano}, B., \& {Weber}, M. 2021, \apj, 913,
122, \dodoi{10.3847/1538-4357/abf5d8}

[Sellek} {et~al.}, 2022]
{Sellek}, A.~D., {Clarke}, C.~J., \& {Ercolano}, B. 2022, \mnras, 514, 535,
\dodoi{10.1093/mnras/stac1148}

[Sellek} {et~al.}, 2024]
{Sellek}, A.~D., {Bajaj}, N.~S., {Pascucci}, I., {et~al.} 2024, \aj, 167, 223,
\dodoi{10.3847/1538-3881/ad34ae}

[Shapiro} {et~al.}, 2004]
{Shapiro}, P.~R., {Iliev}, I.~T., \& {Raga}, A.~C. 2004, \mnras, 348, 753,
\dodoi{10.1111/j.1365-2966.2004.07364.x}

[Shoda} \& {Takasao}, 2021]
{Shoda}, M., \& {Takasao}, S. 2021, \aap, 656, A111,
\dodoi{10.1051/0004-6361/202141563}

[Shu} {et~al.}, 1993]
{Shu}, F.~H., {Johnstone}, D., \& {Hollenbach}, D. 1993, \icarus, 106, 92,
\dodoi{10.1006/icar.1993.1160}

[Suzuki} {et~al.}, 2016]
{Suzuki}, T.~K., {Ogihara}, M., {Morbidelli}, A., {Crida}, A., \& {Guillot}, T.
2016, \aap, 596, A74, \dodoi{10.1051/0004-6361/201628955}

[Takasao} {et~al.}, 2018]
{Takasao}, S., {Tomida}, K., {Iwasaki}, K., \& {Suzuki}, T.~K. 2018, \apj, 857,
4, \dodoi{10.3847/1538-4357/aab5b3}

[Takasao} {et~al.}, 2022]
---. 2022, \apj, 941, 73, \dodoi{10.3847/1538-4357/ac9eb1}

[Tanaka} {et~al.}, 2013]
{Tanaka}, K.~E.~I., {Nakamoto}, T., \& {Omukai}, K. 2013, \apj, 773, 155,
\dodoi{10.1088/0004-637X/773/2/155}

[Turner} {et~al.}, 2010]
{Turner}, N.~J., {Carballido}, A., \& {Sano}, T. 2010, \apj, 708, 188,
\dodoi{10.1088/0004-637X/708/1/188}

[Turner} {et~al.}, 2014]
{Turner}, N.~J., {Fromang}, S., {Gammie}, C., {et~al.} 2014, in Protostars and
Planets VI, ed. H.~{Beuther}, R.~S. {Klessen}, C.~P. {Dullemond}, \&
T.~{Henning}, 411--432, \dodoi{10.2458/azu_uapress_9780816531240-ch018}

[Vicente} {et~al.}, 2013]
{Vicente}, S., {Bern{\'e}}, O., {Tielens}, A.~G.~G.~M., {et~al.} 2013, \apjl,
765, L38, \dodoi{10.1088/2041-8205/765/2/L38}

[Virtanen {et~al.}, 2020]
Virtanen, P., Gommers, R., Oliphant, T.~E., {et~al.} 2020, Nature Methods, 17,
261, \dodoi{10.1038/s41592-019-0686-2}

[Wang} \& {Goodman}, 2017]
{Wang}, L., \& {Goodman}, J. 2017, \apj, 847, 11,
\dodoi{10.3847/1538-4357/aa8726}

[Wilms} {et~al.}, 2000]
{Wilms}, J., {Allen}, A., \& {McCray}, R. 2000, \apj, 542, 914,
\dodoi{10.1086/317016}

[Woods} {et~al.}, 1996]
{Woods}, D.~T., {Klein}, R.~I., {Castor}, J.~I., {McKee}, C.~F., \& {Bell},
J.~B. 1996, \apj, 461, 767, \dodoi{10.1086/177101}

\end{thebibliography}

\begin{thebibliography}{}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\providecommand{\url}[1]{\href{#1}{#1}}
\providecommand{\dodoi}[1]{doi:~\href{http://doi.org/#1}{\nolinkurl{#1}}}
\providecommand{\doeprint}[1]{\href{http://ascl.net/#1}{\nolinkurl{http://ascl.net/#1}}}
\providecommand{\doarXiv}[1]{\href{https://arxiv.org/abs/#1}{\nolinkurl{https://arxiv.org/abs/#1}}}

[Alexander} {et~al.}, 2014]
{Alexander}, R., {Pascucci}, I., {Andrews}, S., {Armitage}, P., \& {Cieza}, L.
2014, in Protostars and Planets VI, ed. H.~{Beuther}, R.~S. {Klessen}, C.~P.
{Dullemond}, \& T.~{Henning}, 475,
\dodoi{10.2458/azu_uapress_9780816531240-ch021}

[Alexander} {et~al.}, 2004]
{Alexander}, R.~D., {Clarke}, C.~J., \& {Pringle}, J.~E. 2004, \mnras, 354, 71,
\dodoi{10.1111/j.1365-2966.2004.08161.x}

[Alexander} {et~al.}, 2005]
---. 2005, \mnras, 358, 283, \dodoi{10.1111/j.1365-2966.2005.08786.x}

[Astropy Collaboration} {et~al.}, 2013]
{Astropy Collaboration}, {Robitaille}, T.~P., {Tollerud}, E.~J., {et~al.} 2013,
\aap, 558, A33, \dodoi{10.1051/0004-6361/201322068}

[Astropy Collaboration} {et~al.}, 2018]
{Astropy Collaboration}, {Price-Whelan}, A.~M., {Sip{\H{o}}cz}, B.~M., {et~al.}
2018, \aj, 156, 123, \dodoi{10.3847/1538-3881/aabc4f}

[Astropy Collaboration} {et~al.}, 2022]
{Astropy Collaboration}, {Price-Whelan}, A.~M., {Lim}, P.~L., {et~al.} 2022,
apj, 935, 167, \dodoi{10.3847/1538-4357/ac7c74}

[Bai}, 2016]
{Bai}, X.-N. 2016, \apj, 821, 80, \dodoi{10.3847/0004-637X/821/2/80}

[Bajaj} {et~al.}, 2024]
{Bajaj}, N.~S., {Pascucci}, I., {Gorti}, U., {et~al.} 2024, \aj, 167, 127,
\dodoi{10.3847/1538-3881/ad22e1}

[Bakes} \& {Tielens}, 1994]
{Bakes}, E.~L.~O., \& {Tielens}, A.~G.~G.~M. 1994, \apj, 427, 822,
\dodoi{10.1086/174188}

[Begelman} {et~al.}, 1983]
{Begelman}, M.~C., {McKee}, C.~F., \& {Shields}, G.~A. 1983, \apj, 271, 70,
\dodoi{10.1086/161178}

[Bertoldi}, 1989]
{Bertoldi}, F. 1989, \apj, 346, 735, \dodoi{10.1086/168055}

[Bertoldi} \& {McKee}, 1990]
{Bertoldi}, F., \& {McKee}, C.~F. 1990, \apj, 354, 529, \dodoi{10.1086/168713}

[Ercolano} {et~al.}, 2009]
{Ercolano}, B., {Clarke}, C.~J., \& {Drake}, J.~J. 2009, \apj, 699, 1639,
\dodoi{10.1088/0004-637X/699/2/1639}

[Ercolano} {et~al.}, 2008]
{Ercolano}, B., {Drake}, J.~J., {Raymond}, J.~C., \& {Clarke}, C.~C. 2008,
\apj, 688, 398, \dodoi{10.1086/590490}

[Fang} {et~al.}, 2023]
{Fang}, M., {Pascucci}, I., {Edwards}, S., {et~al.} 2023, \apj, 945, 112,
\dodoi{10.3847/1538-4357/acb2c9}

[Fang} {et~al.}, 2018]
---. 2018, \apj, 868, 28, \dodoi{10.3847/1538-4357/aae780}

[Font} {et~al.}, 2004]
{Font}, A.~S., {McCarthy}, I.~G., {Johnstone}, D., \& {Ballantyne}, D.~R. 2004,
\apj, 607, 890, \dodoi{10.1086/383518}

[Garcia} {et~al.}, 2001]
{Garcia}, P.~J.~V., {Cabrit}, S., {Ferreira}, J., \& {Binette}, L.
2001{\natexlab{a}}, \aap, 377, 609, \dodoi{10.1051/0004-6361:20011146}

[Garcia} {et~al.}, 2001]
{Garcia}, P.~J.~V., {Ferreira}, J., {Cabrit}, S., \& {Binette}, L.
2001{\natexlab{b}}, \aap, 377, 589, \dodoi{10.1051/0004-6361:20011145}

[Geers} {et~al.}, 2007]
{Geers}, V.~C., {van Dishoeck}, E.~F., {Visser}, R., {et~al.} 2007, \aap, 476,
279, \dodoi{10.1051/0004-6361:20078466}

[Gorti} {et~al.}, 2009]
{Gorti}, U., {Dullemond}, C.~P., \& {Hollenbach}, D. 2009, \apj, 705, 1237,
\dodoi{10.1088/0004-637X/705/2/1237}

[Gorti} \& {Hollenbach}, 2004]
{Gorti}, U., \& {Hollenbach}, D. 2004, \apj, 613, 424, \dodoi{10.1086/422406}

[G{\"u}del} {et~al.}, 2007]
{G{\"u}del}, M., {Briggs}, K.~R., {Arzner}, K., {et~al.} 2007, \aap, 468, 353,
\dodoi{10.1051/0004-6361:20065724}

[Habing}, 1968]
{Habing}, H.~J. 1968, \bain, 19, 421

[Harris {et~al.}, 2020]
Harris, C.~R., Millman, K.~J., van~der Walt, S.~J., {et~al.} 2020, Nature, 585,
357, \dodoi{10.1038/s41586-020-2649-2}

[Hollenbach} \& {Gorti}, 2009]
{Hollenbach}, D., \& {Gorti}, U. 2009, \apj, 703, 1203,
\dodoi{10.1088/0004-637X/703/2/1203}

[Hollenbach} {et~al.}, 1994]
{Hollenbach}, D., {Johnstone}, D., {Lizano}, S., \& {Shu}, F. 1994, \apj, 428,
654, \dodoi{10.1086/174276}

[Hollenbach} {et~al.}, 1993]
{Hollenbach}, D., {Johnstone}, D., \& {Shu}, F. 1993, in Astronomical Society
of the Pacific Conference Series, Vol.~35, Massive Stars: Their Lives in the
Interstellar Medium, ed. J.~P. {Cassinelli} \& E.~B. {Churchwell}, 26

[Hunter, 2007]
Hunter, J.~D. 2007, Computing in Science \& Engineering, 9, 90,
\dodoi{10.1109/MCSE.2007.55}

[Iliev} {et~al.}, 2005]
{Iliev}, I.~T., {Shapiro}, P.~R., \& {Raga}, A.~C. 2005, \mnras, 361, 405,
\dodoi{10.1111/j.1365-2966.2005.09155.x}

[Kahn}, 1954]
{Kahn}, F.~D. 1954, \bain, 12, 187

[Komaki} {et~al.}, 2021]
{Komaki}, A., {Nakatani}, R., \& {Yoshida}, N. 2021, \apj, 910, 51,
\dodoi{10.3847/1538-4357/abe2af}

[Kunitomo} {et~al.}, 2021]
{Kunitomo}, M., {Ida}, S., {Takeuchi}, T., {et~al.} 2021, \apj, 909, 109,
\dodoi{10.3847/1538-4357/abdb2a}

[Liffman}, 2003]
{Liffman}, K. 2003, PASA, 20, 337, \dodoi{10.1071/AS03019}

[Maeder} \& {Meynet}, 1987]
{Maeder}, A., \& {Meynet}, G. 1987, \aap, 182, 243

[Maloney} {et~al.}, 1996]
{Maloney}, P.~R., {Hollenbach}, D.~J., \& {Tielens}, A.~G.~G.~M. 1996, \apj,
466, 561, \dodoi{10.1086/177532}

[McElroy} {et~al.}, 2013]
{McElroy}, D., {Walsh}, C., {Markwick}, A.~J., {et~al.} 2013, \aap, 550, A36,
\dodoi{10.1051/0004-6361/201220465}

[Mignone} {et~al.}, 2007]
{Mignone}, A., {Bodo}, G., {Massaglia}, S., {et~al.} 2007, \apjs, 170, 228,
\dodoi{10.1086/513316}

[Murray-Clay} {et~al.}, 2009]
{Murray-Clay}, R.~A., {Chiang}, E.~I., \& {Murray}, N. 2009, \apj, 693, 23,
\dodoi{10.1088/0004-637X/693/1/23}

[Nakatani} {et~al.}, 2020]
{Nakatani}, R., {Fialkov}, A., \& {Yoshida}, N. 2020, \apj, 905, 151,
\dodoi{10.3847/1538-4357/abc5b4}

[Nakatani} {et~al.}, 2018]
{Nakatani}, R., {Hosokawa}, T., {Yoshida}, N., {Nomura}, H., \& {Kuiper}, R.
2018{\natexlab{a}}, \apj, 857, 57, \dodoi{10.3847/1538-4357/aab70b}

[Nakatani} {et~al.}, 2018]
---. 2018{\natexlab{b}}, \apj, 865, 75, \dodoi{10.3847/1538-4357/aad9fd}

[Nakatani} \& {Takasao}, 2022]
{Nakatani}, R., \& {Takasao}, S. 2022, \apj, 930, 124,
\dodoi{10.3847/1538-4357/ac63a0}

[Nakatani} {et~al.}, 2023]
{Nakatani}, R., {Turner}, N.~J., {Hasegawa}, Y., {et~al.} 2023, \apjl, 959,
L28, \dodoi{10.3847/2041-8213/ad0ed8}

[Nakatani} \& {Yoshida}, 2019]
{Nakatani}, R., \& {Yoshida}, N. 2019, \apj, 883, 127,
\dodoi{10.3847/1538-4357/ab380a}

[Nomura} {et~al.}, 2007]
{Nomura}, H., {Aikawa}, Y., {Tsujimoto}, M., {Nakagawa}, Y., \& {Millar}, T.~J.
2007, \apj, 661, 334, \dodoi{10.1086/513419}

[Owen} \& {Alvarez}, 2016]
{Owen}, J.~E., \& {Alvarez}, M.~A. 2016, \apj, 816, 34,
\dodoi{10.3847/0004-637X/816/1/34}

[Owen} {et~al.}, 2010]
{Owen}, J.~E., {Ercolano}, B., {Clarke}, C.~J., \& {Alexander}, R.~D. 2010,
\mnras, 401, 1415, \dodoi{10.1111/j.1365-2966.2009.15771.x}

[Pascucci} {et~al.}, 2023]
{Pascucci}, I., {Cabrit}, S., {Edwards}, S., {et~al.} 2023, in Astronomical
Society of the Pacific Conference Series, Vol. 534, Protostars and Planets
VII, ed. S.~{Inutsuka}, Y.~{Aikawa}, T.~{Muto}, K.~{Tomida}, \& M.~{Tamura},
567, \dodoi{10.48550/arXiv.2203.10068}

[Pascucci} {et~al.}, 2012]
{Pascucci}, I., {Gorti}, U., \& {Hollenbach}, D. 2012, \apjl, 751, L42,
\dodoi{10.1088/2041-8205/751/2/L42}

[Pascucci} {et~al.}, 2014]
{Pascucci}, I., {Ricci}, L., {Gorti}, U., {et~al.} 2014, \apj, 795, 1,
\dodoi{10.1088/0004-637X/795/1/1}

[Pascucci} {et~al.}, 2020]
{Pascucci}, I., {Banzatti}, A., {Gorti}, U., {et~al.} 2020, \apj, 903, 78,
\dodoi{10.3847/1538-4357/abba3c}

[Picogna} {et~al.}, 2019]
{Picogna}, G., {Ercolano}, B., {Owen}, J.~E., \& {Weber}, M.~L. 2019, \mnras,
487, 691, \dodoi{10.1093/mnras/stz1166}

[Ricci} {et~al.}, 2021]
{Ricci}, L., {Harter}, S.~K., {Ercolano}, B., \& {Weber}, M. 2021, \apj, 913,
122, \dodoi{10.3847/1538-4357/abf5d8}

[Sellek} {et~al.}, 2022]
{Sellek}, A.~D., {Clarke}, C.~J., \& {Ercolano}, B. 2022, \mnras, 514, 535,
\dodoi{10.1093/mnras/stac1148}

[Sellek} {et~al.}, 2024]
{Sellek}, A.~D., {Bajaj}, N.~S., {Pascucci}, I., {et~al.} 2024, \aj, 167, 223,
\dodoi{10.3847/1538-3881/ad34ae}

[Shapiro} {et~al.}, 2004]
{Shapiro}, P.~R., {Iliev}, I.~T., \& {Raga}, A.~C. 2004, \mnras, 348, 753,
\dodoi{10.1111/j.1365-2966.2004.07364.x}

[Shoda} \& {Takasao}, 2021]
{Shoda}, M., \& {Takasao}, S. 2021, \aap, 656, A111,
\dodoi{10.1051/0004-6361/202141563}

[Shu} {et~al.}, 1993]
{Shu}, F.~H., {Johnstone}, D., \& {Hollenbach}, D. 1993, \icarus, 106, 92,
\dodoi{10.1006/icar.1993.1160}

[Suzuki} {et~al.}, 2016]
{Suzuki}, T.~K., {Ogihara}, M., {Morbidelli}, A., {Crida}, A., \& {Guillot}, T.
2016, \aap, 596, A74, \dodoi{10.1051/0004-6361/201628955}

[Takasao} {et~al.}, 2018]
{Takasao}, S., {Tomida}, K., {Iwasaki}, K., \& {Suzuki}, T.~K. 2018, \apj, 857,
4, \dodoi{10.3847/1538-4357/aab5b3}

[Takasao} {et~al.}, 2022]
---. 2022, \apj, 941, 73, \dodoi{10.3847/1538-4357/ac9eb1}

[Tanaka} {et~al.}, 2013]
{Tanaka}, K.~E.~I., {Nakamoto}, T., \& {Omukai}, K. 2013, \apj, 773, 155,
\dodoi{10.1088/0004-637X/773/2/155}

[Turner} {et~al.}, 2010]
{Turner}, N.~J., {Carballido}, A., \& {Sano}, T. 2010, \apj, 708, 188,
\dodoi{10.1088/0004-637X/708/1/188}

[Turner} {et~al.}, 2014]
{Turner}, N.~J., {Fromang}, S., {Gammie}, C., {et~al.} 2014, in Protostars and
Planets VI, ed. H.~{Beuther}, R.~S. {Klessen}, C.~P. {Dullemond}, \&
T.~{Henning}, 411--432, \dodoi{10.2458/azu_uapress_9780816531240-ch018}

[Vicente} {et~al.}, 2013]
{Vicente}, S., {Bern{\'e}}, O., {Tielens}, A.~G.~G.~M., {et~al.} 2013, \apjl,
765, L38, \dodoi{10.1088/2041-8205/765/2/L38}

[Virtanen {et~al.}, 2020]
Virtanen, P., Gommers, R., Oliphant, T.~E., {et~al.} 2020, Nature Methods, 17,
261, \dodoi{10.1038/s41592-019-0686-2}

[Wang} \& {Goodman}, 2017]
{Wang}, L., \& {Goodman}, J. 2017, \apj, 847, 11,
\dodoi{10.3847/1538-4357/aa8726}

[Wilms} {et~al.}, 2000]
{Wilms}, J., {Allen}, A., \& {McCray}, R. 2000, \apj, 542, 914,
\dodoi{10.1086/317016}

[Woods} {et~al.}, 1996]
{Woods}, D.~T., {Klein}, R.~I., {Castor}, J.~I., {McKee}, C.~F., \& {Bell},
J.~B. 1996, \apj, 461, 767, \dodoi{10.1086/177101}

\end{thebibliography}

\begin{appendix}
\section{Basic notes for the analytical model}
\label{sec:list_of_variables}
We list symbols used in our phenomenological model in \tref{tab:variables}.
\begin{table*}[]
\centering
\begin{tabular}{L l r}
\text{Expression}
\\
\hline
\alpha  &   Case~B recombination rate coefficient   &   \eqnref{eq:recombination_coeff}  \\
\alpha_\mathrm{eq}  &   $\alpha$ at the equilibrium temperature ($\approx 2\e{-13}\cm^3\sec^{-1}$)&   \eqnref{eq:recombination_coeff}  \\
\beta   &   Power index for the temperature profile of $\alpha$ &   \eqnref{eq:recombination_coeff}  \\
C       &   Empirical dimensionless factor for base density   &   \eqnref{eq:base_density}   \\
c_\mathrm{s}    &   Typical isothermal sound speed of a wind    &   \\
c_\mathrm{eq}   &   Isothermal sound speed at the equilibrium temperature ($\approx 10\kms$)    &   Beginning of \secref{sec:critical_radius}   \\
c_\mathrm{ch}   &   Characteristic sound speed  &   \eqnref{eq:c_ch_def} \\
c_\mathrm{ch}^\mathrm{I}    &   Characteristic sound speed for Type~I winds &   \eqnref{eq:cchI_def} \\
c_\mathrm{ch}^\mathrm{II}    &   Characteristic sound speed for Type~II winds &   \eqnref{eq:chII_original_def} \\
c_p     &   Specific heat at constant pressure  &   \eqnref{eq:critical_radius}\\
c_1     &   Approximated $c_\mathrm{ch}^\mathrm{II}$ in the slow-recombination limit    &   \eqnref{eq:c2_in_slow_recombination} \\
c_2     &   Approximated $c_\mathrm{ch}^\mathrm{II}$ in the rapid-recombination limit   &   \eqnref{eq:cchII_def}    \\
\Gamma_\mathrm{EUV} &   Specific EUV photoheating rate  &   \eqnref{eq:gamma_EUV} and \eqnref{eq:gamma_EUV:2}\\
\delta  &   Ratio of attenuated flux to the unattenuated flux   &   \eqnref{eq:gamma_EUV}    \\
\langle\Delta E\rangle_i    &   Frequency-average deposited energy per photoionization  &   \eqnref{eq:gamma_EUV}    \\
\aveenergy   &   $\langle\Delta E\rangle_i $ at $N_\mathrm{H} = 0$   &   \eqnref{eq:gamma_EUV:2}\\
\varepsilon &   Normalized $\bar{E}_0$ by $mc_p \ceq^2$  &   \eqnref{eq:tildes}   \\
F_\nu       &   Specific photon number flux     &   \secref{sec:euv_heating_rate}   \\
F_0     &   Unattenuated total photon number flux   &   \secref{sec:euv_heating_rate}   \\
f_\mathrm{eq}   &   Ratio of the secondary critical EUV emission rate to $\Phi_\mathrm{c}$      &   \eqnref{eq:secondary_critical_luminosity}\\
k_\mathrm{ioni} &   Photoionization rate coefficient    &   \eqnref{eq:tioni} \\
M_*     &   Mass of the radiation source    &   \secref{sec:critical_radius}    \\
m       &   Gas mass per hydrogen nucleus   &   \eqnref{eq:gamma_EUV}\\
\mach   &   Typical Mach number of a wind   &   \eqnref{eq:slow_recombination_condition} \\
\nh     &   Hydrogen nuclei number density  &   \secref{sec:euv_heating_rate}    \\
N_\mathrm{HI}   &   \ion{H}{1} column density & \secref{sec:euv_heating_rate}   \\
\Phi_\mathrm{EUV}   &   EUV emission rate of the source     &   \secref{sec:euv_heating_rate}   \\
\Phi_\mathrm{c} &   Critical EUV emission rate &    \eqnref{eq:critical_luminosity}  \\
\Phi_\mathrm{c}^\mathrm{II} &   Secondary critical EUV emission rate, $f_\mathrm{eq}\Phi_\mathrm{c}$ &  \eqnref{eq:secondary_critical_luminosity}    \\
\Phi_\mathrm{c,rec}^\mathrm{II} &   Secondary critical EUV emission rate for recombination-dominated Type~II winds &  \eqnref{eq:cchII_def}    \\
\varphi &   Normalized EUV emission rate, $\Phi_\mathrm{EUV}/\Phi_\mathrm{c}$    &   \eqnref{eq:cchI_def} \\
\varphi_\mathrm{II} &   EUV emission rate normalized by $\Phi_\mathrm{c,rec}^\mathrm{II}$   &   \eqnref{eq:cchII_def}    \\
\varphi_\mathrm{g}  &  Boundary of Regions~B-II and C at which $c_\mathrm{ch}^\mathrm{II} = v_\mathrm{g}$   &   \eqnref{eq:varphi_g}\\
\psi    &   Photoevaporation parameter in Bertoldi}, 1989 &   \eqnref{eq:trec_tcross}  \\
q   &   Spectrum hardness parameter for recombination   &   \eqnref{eq:q}  \\
r               &   Spherical radius from the source    &   \secref{sec:euv_heating_rate}   \\
R               &   Cylindrical distance from the source  &   \secref{sec:critical_radius}  \\
R_\mathrm{c}    &   Critical radius     &   \eqnref{eq:critical_radius}  \\
\langle\sigma\rangle    &   Frequency-average cross-section &   \eqnref{eq:gamma_EUV}    \\
\bar{\sigma}_0  &  $\langle\sigma\rangle$ at $N_\mathrm{H} = 0$   &   \eqnref{eq:gamma_EUV:2}\\
t_\mathrm{g}    &   Gravitational timescale & \eqnref{eq:gravitational_timescale}\\
t_\mathrm{ioni} &   Photoionization timescale   &   \eqnref{eq:tioni}\\
t_\mathrm{rec}  &   Recombination timescale     &   \eqnref{eq:trec_base_density}    \\
v_\mathrm{g}    &   Gravitational velocity  &   \eqnref{eq:vg_ceq}   \\
x       &   Normalized distance, $R/R_\mathrm{c}$    &   \secref{sec:critical_radius}\\
\chi_\mathrm{e} &   Attenuation factor due to absorption for photoheating & \eqnref{eq:gamma_EUV:2}  \\
\chi_\mathrm{i} &   Attenuation factor due to absorption for photoionization    &   \eqnref{eq:tioni}    \\
y_\mathrm{HI}   &   Atomic hydrogen abundance   &   \eqnref{eq:gamma_EUV}
\end{tabular}
\caption{Symbols for the fundamental variables used in our phenomenological model (Sections~\ref{sec:analytic_theory}--\ref{sec:cs}). We omit those that appear in later sections.}
\label{tab:variables}
\end{table*}

Throughout this paper,
we use ``region'' to denote spatial extents, such as the \ion{H}{2} region, region of $x > 1$, or region where $s < R$.
The word ``phase diagram'' refers to the $x$--$\varphi$ map.
In contrast, the $\varepsilon$--$q$ plane is referred to as the ``spectral-hardness  diagram.''
We use ``parameter space'' to indicate a certain area of the phase diagrams and the spectral-hardness diagram.
The term ``regime'' is employed to describe the parameter spaces where the same qualitative properties are shared, such as Regime~A, Type~I regime, among others.
Therefore, ``parameter space'' and ``regime'' are utilized in a different context from ``region'' consistently.

% \input{subsections/energy_equation}
% \input{subsections/pi_param}
% \section{Mathematical Proof for $c_\mathrm{ch}^\mathrm{II} > c_\mathrm{eq}$ in H-class}
\section{Mathematical proof for Type~II characteristic sound speed exceeding the equilibrium sound speed in H-class}
\label{sec:math_proof}

In \secref{sec:typeII:regime:H}, we mention that $c_\mathrm{ch}^\mathrm{II} > c_\mathrm{eq}$ generally holds for the H-class ($\left\{(\epsilon,q)|~\epsilon > 1 ~\cup~ \epsilon q > 1\right\}$), and it can be mathematically proven from \eqnref{eq:chII_original_def}.
We demonstrate this proof in this Appendix.

When $\epsilon > 1$, it is straightforward to see that $c_\mathrm{ch}^\mathrm{II} > c_\mathrm{eq}$, as it necessitates
\[
\braket{\frac{c_\mathrm{ch}^\mathrm{II}}{c_\mathrm{eq}}}^2 > \epsilon
\]
according to \eqnref{eq:chII_original_def}.
On the other hand, when $\epsilon < 1$ and $\epsilon q > 1$, the parameter space of interest is $\varphi > \mathrm{max}\braket{1, q} \epsilon x = \epsilon q x$.
In this range,
\[
\begin{split}
\braket{\frac{c_\mathrm{ch}^\mathrm{II}}{c_\mathrm{eq}}}^2
& = \epsilon
+ \braket{\frac{\varphi \epsilon  q}{ x }}^{1/2}
\braket{\frac{c_\mathrm{ch}^\mathrm{II}}{c_\mathrm{eq}}}^{-(\beta+1)}
\braket{1 - \varphi^{-1}\epsilon  x \frac{c_\mathrm{ch}^\mathrm{II}}{c_\mathrm{eq}}}
\\ &
> \epsilon + \epsilon q \braket{\frac{c_\mathrm{ch}^\mathrm{II}}{c_\mathrm{eq}}}^{-(\beta+1)}
\braket{1 - \frac{1}{q}\frac{c_\mathrm{ch}^\mathrm{II}}{c_\mathrm{eq}}}
\end{split}
\]
is required from \eqnref{eq:chII_original_def}.
The RHS of this inequity is a monotonically decreasing function of $c_\mathrm{ch}^\mathrm{II}$ that takes a value of $\epsilon q$ at $c_\mathrm{ch}^\mathrm{II}/c_\mathrm{eq} = 1$,
while the left-hand side (LHS) takes unity at the same $c_\mathrm{ch}^\mathrm{II}$.
This indicates that the intersection between the functions of the LHS and RHS is present at $c_\mathrm{ch}^\mathrm{II}/c_\mathrm{eq} > 1$.
Denoting $c^\prime $ as the coordinate of the intersection,
the above inequity is rewritten as
\[
\frac{c_\mathrm{ch}^\mathrm{II}}{c_\mathrm{eq}} > \frac{c^\prime}{c_\mathrm{eq}} > 1.
\]
Thus, $c_\mathrm{ch}^\mathrm{II} > c_\mathrm{eq}$ is proven for Type~II winds in the H-class.
% \input{subsections/appendix}
% \input{subsections/cross_section}
% \input{subsections/on_absorption}
% \input{trash/previous_sclass_note}

\end{appendix}

\end{document}