1 INTRODUCTION
Photoevaporation is a pressure-driven wind produced by high energy stellar radiation that heats and/or ionises gas located in the incident surface layers of protoplanetary atmospheres (Hollenbach et al. Reference Hollenbach, Johnstone, Lizano and Shu1994). If the thermal energy of the heated gas exceeds the gravitational binding energy of the central gravitating body, the gas is unbound and can escape in a slow, often centrifugally launched, wind. These winds are similar in nature to the familiar pressure-driven Parker winds in stars (Parker Reference Parker1958), but are made complicated by rotation, disc geometry, and/or off-axis radiation sources. For example, the flow solution for photoionised disc winds cannot rigorously be solved analytically because the solution depends on knowing a priori the exact streamline trajectories (or divergence; see Begelman, McKee, & Shields Reference Begelman, McKee and Shields1983). While trivial for spherically symmetric winds, the extension to discs can only be approximated (e.g. Waters & Proga Reference Waters and Proga2012).
The analytic solution for isothermal Parker winds has typically been used as a numerical test for hydrodynamic simulations involving astrophysical winds (e.g. Keppens & Goedbloed Reference Keppens and Goedbloed1999; Font et al. Reference Font, McCarthy, Johnstone and Ballantyne2004). However, apart from sharing a similar transonic wind structure, stellar winds and photoevaporation in discs are physically quite different (e.g. geometry, gravity, temperature, density). If one is only interested in photoevaporating discs, the numerical overhead of setting up alternate conditions necessary to produce stellar winds can be inconvenient. In such cases, it would be ideal to have an analytic solution to a problem that uses the same numerical setup and physical parameters as a real disc.
An analytic wind solution for photoevaporation in a disc-like environment can be derived using a non-rotating, stratified, plane-parallel atmosphere. On local scales, the vertical structure of protoplanetary discs is approximately plane-parallel so the physical parameters and numerical setup can be made to be almost identical to that of a disc at any given radius. The resulting wind’s simple 1-D geometry makes the solution analytically tractable and straight forward to use as an alternative test to the isothermal Parker wind—its utility has motivated this study.
2 ANALYTIC FLOW SOLUTION
The relevant equations describing a steady-state, pressure-driven, isothermal Parker wind come from setting ∂/∂t = 0 in the fluid equations:
$$\begin{equation}
\nabla \cdot (\rho \mathbf {v}) = 0,
\end{equation}$$\\
$$\begin{equation}
\rho \left( \mathbf {v} \cdot \nabla \mathbf {v} \right) = - \nabla P + \rho \mathbf {g},
\end{equation}$$\\
$$\begin{equation}
P = \rho \mathcal {R} T,
\end{equation}$$\\
$$\begin{equation}
T = T_0,
\end{equation}$$
where g is the gravitational force,
$\mathcal {R}$
is the gas constant, and ρ, v, P, and T are the gas density, velocity, pressure, and temperature, respectively. In anticipation of applying this test to photoevaporating circumstellar discs, we define the gravity g to be the vertical field produced by a massive central object,
$$\begin{equation}
\mathbf {g} = -\frac{\mathcal {G} M z}{(R^2+z^2)^{3/2}}\mathbf {\hat{z}},
\end{equation}$$
to ensure a disc-like density and temperature structure in the atmosphere. Here,
$\mathcal {G}$
is the gravitational constant, M is the mass of the central star, and R is the cylindrical distance from the central source to our local patch of atmosphere. Without loss of generality, we restrict the variables to be functions of z only. To close the set of equations, we adopt the equation of state of an ideal gas (
$pV = n\mathcal {R}T=$
constant, where n is the number of moles of the gas) such that the sound speed of the wind is constant and can be written as
$c_{\text{s}}^{2} = P/\rho$
.
Integrating Equation (1) gives
$\rho v = \text{constant}$
, or written in terms of an accretion rate (Bondi Reference Bondi1952),
$$\begin{equation}
\dot{M} = A \rho v,
\end{equation}$$
where A is a problem dependent characteristic surface area. Meanwhile, using the sound speed relationship to replace P, Equation (2) can be rewritten as
$$\begin{equation}
v \frac{\text{d} v}{\text{d} z} = - \frac{c_{\text{s}}^2 }{\rho } \frac{\text{d}\rho }{\text{d} z} - \frac{\mathcal {G} M z}{(R^2+z^2)^{3/2}}.
\end{equation}$$
The dependence here on ρ can be removed by taking the derivative of Equation (6). After some manipulation, we obtain
$$\begin{equation}
-\frac{1}{\rho } \frac{\text{d}\rho }{\text{d}z} = \frac{1}{v} \frac{\text{d}v}{\text{d}z},
\end{equation}$$
which can immediately be substituted back into Equation (7) to obtain
$$\begin{equation}
v \frac{\text{d} v}{\text{d} z} = \frac{c_{\text{s}}^2 }{v} \frac{\text{d} v}{\text{d} z} - \frac{\mathcal {G} M z}{(R^2+z^2)^{3/2}}.
\end{equation}$$
Collecting the derivatives on v and using the following relation,
$$\begin{equation}
v \frac{\text{d} v }{\text{d}z} = \frac{ c_{\text{s}}^2}{2} \frac{\text{d} (v^2/c_{\text{s}}^2)}{\text{d}z},
\end{equation}$$
we obtain a separable, ordinary differential equation for
$v^2/c_{\text{s}}^2$
:
$$\begin{equation}
\left( 1-\frac{c_{\text{s}}^2}{v^2} \right) \frac{\text{d} (v^2/c_{\text{s}}^2)}{\text{d}z} = - \frac{2 \mathcal {G} M z}{c_{\text{s}}^2 \, (R^2+z^2)^{3/2}}.
\end{equation}$$
Non-dimensionalising Equation (11) using
$\bar{v}^2 \equiv v^2/c_{\text{s}}^2$
,
$\bar{z} \equiv z/R$
, the Keplerian Mach number
$\mathcal {M} \equiv v_K/c_{\text{s}}$
, and
$v_K = \sqrt{\mathcal {G} M/R}$
, we obtain
$$\begin{equation}
\left( 1-\frac{1}{\bar{v}^2} \right) \frac{\text{d} \left(\bar{v}^2\right)}{\text{d}\bar{z}} = - \frac{2 \mathcal {M}^2 \bar{z}}{\left(1+\bar{z}^2\right)^{3/2}}.
\end{equation}$$
Note the presence of a critical point located at the sonic point,
$\bar{v} = 1$
, on the left-hand side of the equation. From inspection of the right-hand side, the corresponding position must be at
${|\bar{z}|} \rightarrow \infty$
. For comparison, the spherically symmetric isothermal Parker-wind solution is transonic with the sonic point located at
$r_{\text{s}} \equiv \mathcal {G}M/2c_{\text{s}}^2$
.
Integrating Equation (12), we obtain a transcendental equation for the outflow velocity as a function of
$\bar{z}$
,
$$\begin{equation}
\bar{v}^2 - \ln { \bar{v}^2 } = \frac{2\mathcal {M}^2}{\sqrt{1+\bar{z}^2}} + C,
\end{equation}$$
where C is an integration constant. Following Cranmer (Reference Cranmer2004), we can write the solution for the velocity in closed form using the Lambert W function (Corless et al. Reference Corless, Gonnet, Hare, Jeffrey and Knuth1996; Veberič Reference Veberič2012):
$$\begin{equation}
\bar{v}^2 = - \mathrm{W}_k \! \! \left[ -\exp { \left( -\frac{2\mathcal {M}^2}{\sqrt{1+\bar{z}^2}} -C \right) } \right],
\end{equation}$$
where
$$\begin{equation}
k = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}0, & \quad \text{if } \bar{v} \le 1 \\
-1, & \quad \text{if } \bar{v} > 1. \end{array}\right.
\end{equation}$$
For comparison, the velocity for the spherically symmetric Parker wind is
$$\begin{equation}
\bar{v}^2 = - \mathrm{W}_k \! \! \left[ - \frac{1}{\bar{r}^{4}} \exp { \left( -\frac{4}{\bar{r}} -C \right) } \right],
\end{equation}$$
where
$\bar{r} \equiv r/r_{\text{s}}$
and r is the spherical radius measured from the centre of mass M. Figure 1 contrasts the two solutions above. As the plane-parallel wind cannot support a finite sonic point without diverging streamlines (Begelman et al. Reference Begelman, McKee and Shields1983), it looks similar to an isothermal Parker wind with its sonic point remapped to infinity. Consequently, the plane-parallel ‘breeze’ solutions (always subsonic) are not hydrostatic at infinity. Another minor difference is that the plane-parallel solutions remain finite at z = 0 due to having a finite gravitational potential at the midplane of the disc. As a final point of interest, the rate of convergence of
$v \rightarrow c_{\text{s}}$
in the asymptotically transonic solution (C = 1) can more conveniently be expressed by expanding Equation (14) in a Taylor series in the limit
$|\bar{z}| \rightarrow \infty$
,
$$\begin{equation}
\bar{v} \approx 1-\frac{\mathcal {M}}{\sqrt{|\bar{z}|}} + \mathcal {O}\left( \frac{1}{z}\right),
\end{equation}$$
which, due to the
$\sqrt{\bar{z}}$
dependence, makes convergence very slow.
Comparison between the plane-parallel wind at R = 5AU (left) derived in this paper and the more familiar spherically symmetric Parker wind (right). The different contour levels are determined by the value of C in Equation (14). The sonic point for the plane-parallel case is only asymptotically crossed as z → ∞, whereas the Parker-wind model is transonic at
$r_{\text{s}} = \mathcal {G}M/2c_{\text{s}}^2$
.
The plane-parallel wind has only three possible classes of solutions:
-
(i) C < 1: v(z) is double-valued on
${z_i \le z\le z_{\text{max}}}$
. -
(ii) C = 1: v(z) is asymptotically transonic and monotonically increasing for outflow (k = 0) or decreasing for inflow (k = −1).
-
(iii) C > 1: v(z) is not transonic and monotonically increasing/decreasing.
Physically, the solution should be locally mono-valued for stability while continuity and symmetry of the disc indicate that the velocity should be stationary at z = 0. Meanwhile, studies of the Parker wind show its breeze solutions to be unstable (Velli Reference Velli1994), a result which holds in the zero-curvature limit. The only admissible solution remaining is C = 1 with k = 0, i.e.,
$$\begin{equation}
\bar{v} = \sqrt{- \mathrm{W}_0 \! \! \left[ -\exp{ \left( -\frac{2\mathcal{M}^2}{\sqrt{1+\bar{z}^2}} -1 \right) } \right]},
\end{equation}$$
but this too is only marginally stable to Velli’s global stability criterion (Velli Reference Velli1994; Grappin, Cavillier, & Velli Reference Grappin, Cavillier and Velli1997; Del Zanna, Velli, & Londrillo Reference Del Zanna, Velli and Londrillo1998). We must therefore turn to numerical simulations to verify the stability of the solution, as suggested by Waters & Proga (Reference Waters and Proga2012).
3 NUMERICAL STABILITY
Using Equations (6) and (18) and the equation of state, all of the fluid parameters are uniquely determined. A practical setup for our proposed test can be achieved in three steps:
-
(i) In a 2- or 3-D box with periodic horizontal boundary conditions, set up a vertically isothermal atmosphere using the thin-disc approximation and the gravitational force given in Equation (5).
-
(ii) Instantaneously heat any fluid that falls below some density threshold to a high temperature (e.g. T = 10000K to mimic ultraviolet photoevaporation; see Alexander, Clarke, & Pringle Reference Alexander, Clarke and Pringle2006).
-
(iii) Create a steady-state flow using a vertical boundary condition appropriate for the numerical method of choice. In smoothed particle hydrodynamics (SPH), this is done with a dynamic vertical boundary condition that is constrained to move at the prescribed analytic velocity from Equation (18). The benefit of this method is that steady-state solution is obtained almost immediately. Grid based codes, on the other hand, will typically converge to the steady-state solution using fixed outflow boundary conditions. However, if convergence is too slow, dynamically forcing a small section of the outflow until it exits the computational domain will help precipitate steady-state flow.
Implementing the setup and SPH boundary conditions described above, we perform the photoevaporation test using our SPH code gdphoto (Hutchison et al. Reference Hutchison, Price, Laibe and Maddison2016). Gdphoto has been benchmarked against the test suite described in Laibe & Price (Reference Laibe and Price2012) and achieves accuracies comparable to commonly used SPH codes. Using 200028 particles, we create a 2-D disc in isothermal hydrostatic equilibrium with the following physical parameters:
$M = 1\,\text{M}_\odot$
, R = 5AU, and ρ0 = 10−11g cm−3. We then initiate photoevaporation by ionising all particles with densities that are five orders of magnitude below the disc midplane density. Ionised particles are held isothermally at T = 10000K such that
$c_{\text{s}} \approx 10\,$
km s−1. Figure 2 shows the results after 100yr plotted together with the analytic solution from Equation (18). The L
2 errors for the velocity and density, computed using splash (Price Reference Price2007), are ~ 2 and ≲ 1%, respectively.
Velocity (top) and density (bottom) for 200028 SPH gas particles plotted against their respective analytic solutions (red solid and dashed lines) in a 2-D plane-parallel disc wind. The green points make up a semi-uniform lattice of gas particles that form a moving boundary condition constrained to move at the velocity prescribed by the analytic solution in Equation (18). The blue points are the unrestrained gas particles. The L 2 errors between the analytic and the numerical solutions are < 2%, consistently with the second-order SPH scheme used.
4 DISCUSSION
The plane-parallel flow described in this paper is comparable to the flow derived by Adams (Reference Adams2011) for magnetically controlled outflows from hot Jupiters when the stellar magnetic field completely dominates over the magnetic field produced by the planet (i.e. their β → ∞ limit). The main difference between our analytic solutions stems from assuming different gravitational potentials; however, to apply their solution to photoevaporating discs would require readers to rederive the equations themselves. The solution in this paper is significantly more transparent and its closed form makes it especially easy to implement as a numerical test.
Although we focus on using our plane-parallel model as a numerical test, it may have use in wider applications as well. For example, the model’s simple geometry, accurate approximation of winds close to the disc, and closed form wind solution could make it a perfect springboard for developing an analytic or semi-analytic model for coupled two-phase photoevaporation. To date, few studies have focused on dust dynamics in winds and a simple two-phase model would be very valuable.
ACKNOWLEDGEMENTS
We thank Daniel Price and Sarah Maddison for useful discussions and the anonymous referee whose careful report significantly improved this paper and pointed out the explicit form of the solution. M.H. acknowledges funding from a Swinburne University Postgraduate Research Award (SUPRA). G. Laibe acknowledges funding from the European Research Council for the FP7 ERC advanced grant project ECOGAL. The numerical calculations were performed on the SwinSTAR supercomputer at Swinburne University of Technology and the subsequent visualisation was made using splash (Price Reference Price2007).
