A Plane-Parallel Wind Solution For Testing Numerical Simulations of Photoevaporation
PPublications of the Astronomical Society of Australia (PASA)c (cid:13)
Astronomical Society of Australia 2018; published by Cambridge University Press.doi: 10.1017/pas.2018.xxx.
A Plane-Parallel Wind Solution For Testing NumericalSimulations of Photoevaporation
Mark A. Hutchison ∗ and Guillaume Laibe Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Hawthorn, VIC 3122, Australia School of Physics and Astronomy, University of St. Andrews, North Haugh, St. Andrews, Fife KY16 9SS, UK
Abstract
Here we derive a Parker-wind like solution for a stratified, plane-parallel atmosphere undergoing photoioni-sation. The difference compared to the standard Parker solar wind is that the sonic point is crossed only atinfinity. The simplicity of the analytic solution makes it a convenient test problem for numerical simulationsof photoevaporation in protoplanetary discs.
Keywords: protoplanetary discs – planets and satellites: atmospheres – circumstellar matter
Photoevaporation is a pressure-driven wind producedby high energy stellar radiation that heats and/orionises gas located in the incident surface layers of pro-toplanetary atmospheres (Hollenbach et al. 1994). If thethermal energy of the heated gas exceeds the gravita-tional binding energy of the central gravitating body,the gas is unbound and can escape in a slow, often cen-trifugally launched, wind. These winds are similar innature to the familiar pressure-driven Parker winds instars (Parker 1958), but are made complicated by rota-tion, disc geometry, and/or off-axis radiation sources.For example, the flow solution for photoionised discwinds cannot rigorously be solved analytically becausethe solution depends on knowing a priori the exactstreamline trajectories (or divergence; see Begelmanet al. 1983). While trivial for spherically symmetricwinds, the extension to discs can only be approximated(e.g. Waters & Proga 2012).The analytic solution for isothermal Parker winds hastypically been used as a numerical test for hydrody-namic simulations involving astrophysical winds (e.g.Keppens & Goedbloed 1999; Font et al. 2004). How-ever, apart from sharing a similar transonic wind struc-ture, stellar winds and photoevaporation in discs arephysically quite different (e.g. geometry, gravity, tem-perature, density). If one is only interested in photoe-vaporating discs, the numerical overhead of setting upalternate conditions necessary to produce stellar windscan be inconvenient. In such cases, it would be ideal to ∗ E-mail: [email protected] have an analytic solution to a problem that uses thesame numerical setup and physical parameters as a realdisc.An analytic wind solution for photoevaporation ina disc -like environment can be derived using a non-rotating, stratified, plane-parallel atmosphere. On localscales, the vertical structure of protoplanetary discs isapproximately plane-parallel so the physical parametersand numerical setup can be made to be almost identi-cal to that of a disc at any given radius. The resultingwind’s simple 1-D geometry makes the solution analyt-ically tractable and straight forward to use as an alter-native test to the isothermal Parker wind—its utilityhas motivated this study.
The relevant equations describing a steady-state,pressure-driven, isothermal Parker wind come from set-ting ∂/∂t = 0 in the fluid equations: ∇ · ( ρ v ) = 0 , (1) ρ ( v · ∇ v ) = −∇ P + ρ g , (2) P = ρ R T, (3) T = T , (4)where g is the gravitational force, R is the gas con-stant, and ρ , v , P , and T are the gas density, velocity,pressure, and temperature, respectively. In anticipationof applying this test to photoevaporating circumstellardiscs, we define the gravity g to be the vertical field1 a r X i v : . [ a s t r o - ph . S R ] M a r Hutchison & Laibe produced by a massive central object, g = − G M z ( R + z ) / ˆz , (5)to ensure a disc-like density and temperature structurein the atmosphere. Here G is the gravitational constant, M is the mass of the central star, and R is the cylindri-cal distance from the central source to our local patchof atmosphere. Without loss of generality, we restrictthe variables to be functions of z only. To close the setof equations, we adopt the equation of state of an idealgas ( pV = n R T = constant, where n is the number ofmoles of the gas) such that the sound speed of the windis constant and can be written as c = P/ρ .Integrating equation (1) gives ρv = constant, or writ-ten in terms of an accretion rate (Bondi 1952),˙ M = Aρv, (6)where A is a problem dependent characteristic surfacearea. Meanwhile, using the sound speed relationship toreplace P , equation (2) can be rewritten as, v d v d z = − c ρ d ρ d z − G M z ( R + z ) / . (7)The dependence here on ρ can be removed by takingthe derivative of equation (6). After some manipulationwe obtain, − ρ d ρ d z = 1 v d v d z , (8)which can immediately be substituted back into equa-tion (7) to obtain, v d v d z = c v d v d z − G M z ( R + z ) / . (9)Collecting the derivatives on v and using the followingrelation, v d v d z = c v /c )d z , (10)we obtain a separable, ordinary differential equation for v /c : (cid:18) − c v (cid:19) d( v /c )d z = − G M zc ( R + z ) / . (11)Nondimensionalising equation (11) using ¯ v ≡ v /c ,¯ z ≡ z/R , the Keplerian Mach number M ≡ v K /c s , and v K = (cid:112) G M/R , we obtain, (cid:18) − v (cid:19) d (cid:0) ¯ v (cid:1) d¯ z = − M ¯ z (1 + ¯ z ) / . (12)Note the presence of a critical point located at thesonic point, ¯ v = 1, on the left-hand side of the equa-tion. From inspection of the right-hand side, the cor-responding position must be at | ¯ z | → ∞ . For compari-son, the spherically symmetric isothermal Parker-wind solution is transonic with the sonic point located at r s ≡ G M/ c .Integrating equation (12) we obtain a transcendentalequation for the outflow velocity as a function of ¯ z ,¯ v − ln ¯ v = 2 M √ z + C, (13)where C is an integration constant. Following Cran-mer (2004), we can write the solution for the velocityin closed form using the Lambert W function (Corlesset al. 1996; Veberiˇc 2012):¯ v = − W k (cid:20) − exp (cid:18) − M √ z − C (cid:19)(cid:21) , (14)where k = (cid:40) , if ¯ v ≤ − , if ¯ v > . (15)For comparison, the velocity for the spherically-symmetric Parker wind is,¯ v = − W k (cid:20) − r exp (cid:18) − r − C (cid:19)(cid:21) , (16)where ¯ r ≡ r/r s and r is the spherical radius measuredfrom the centre of mass M . Figure 1 contrasts the twosolutions above. As the plane-parallel wind cannot sup-port a finite sonic point without diverging streamlines(Begelman et al. 1983), it looks similar to an isother-mal Parker wind with its sonic point remapped to infin-ity. Consequently, the plane-parallel “breeze” solutions(always subsonic) are not hydrostatic at infinity. An-other minor difference is that the plane-parallel solu-tions remain finite at z = 0 due to having a finite gravi-tational potential at the midplane of the disc. As a finalpoint of interest, the rate of convergence of v → c s inthe asymptotically transonic solution ( C = 1) can moreconveniently be expressed by expanding equation (14)in a Taylor series in the limit | ¯ z | → ∞ ,¯ v ≈ − M (cid:112) | ¯ z | + O (cid:18) z (cid:19) , (17)which, due to the √ ¯ z dependence, makes convergencevery slow.The plane-parallel wind has only three possibleclasses of solutions:(i) C < v ( z ) is double-valued on z i ≤ z ≤ z max .(ii) C = 1: v ( z ) is asymptotically transonic and mono-tonically increasing for outflow ( k = 0) ordecreasing for inflow ( k = − C > v ( z ) is not transonic and monotonically in-creasing/decreasing.Physically, the solution should be locally mono-valuedfor stability while continuity and symmetry of the disc PASA (2018)doi:10.1017/pas.2018.xxx
Plane-Parallel Wind Solution z = 0.Meanwhile, studies of the Parker wind show its breezesolutions to be unstable (Velli 1994), a result whichholds in the zero-curvature limit. The only admissiblesolution remaining is C = 1 with k = 0, i.e.¯ v = (cid:115) − W (cid:20) − exp (cid:18) − M √ z − (cid:19)(cid:21) , (18)but this too is only marginally stable to Velli’s globalstability criterion (Velli 1994; Grappin et al. 1997; DelZanna et al. 1998). We must therefore turn to numericalsimulations to verify the stability of the solution, assuggested by Waters & Proga (2012). Using equations (6) and (18) and the equation of state,all of the fluid parameters are uniquely determined. Apractical setup for our proposed test can be achieved inthree steps:(i) In a 2- or 3-D box with periodic horizontal bound-ary conditions, set up a vertically-isothermal atmo-sphere using the thin-disc approximation and the grav-itational force given in equation (5).(ii) Instantaneously heat any fluid that falls belowsome density threshold to a high temperature (e.g. T = 10 000 K to mimic ultraviolet photoevaporation;see Alexander et al. 2006).(iii) Create a steady-state flow using a vertical bound-ary condition appropriate for the numerical method ofchoice. In smoothed particle hydrodynamics (SPH), thisis done with a dynamic vertical boundary condition thatis constrained to move at the prescribed analytic veloc-ity from equation (18). The benefit of this method isthat steady-state solution is obtained almost immedi-ately. Grid based codes, on the other hand, will typi-cally converge to the steady-state solution using fixedoutflow boundary conditions. However, if convergence istoo slow, dynamically forcing a small section of the out-flow until it exits the computational domain will helpprecipitate steady-state flow.Implementing the setup and SPH boundary condi-tions described above, we perform the photoevapora-tion test using our SPH code gdphoto (Hutchisonet al. 2016). Gdphoto has been benchmarked againstthe test suite described in Laibe & Price (2012) andachieves accuracies comparable to commonly used SPHcodes. Using 200 028 particles we create a 2-D discin isothermal hydrostatic equilibrium with the follow-ing physical parameters: M = 1 M (cid:12) , R = 5 AU, and ρ = 10 − g/cm . We then initiate photoevaporationby ionising all particles with densities that are fiveorders of magnitude below the disc midplane density.Ionised particles are held isothermally at T = 10 000 Ksuch that c s ≈
10 km/s. Figure 2 shows the results after 100 yr plotted together with the analytic solution fromequation (18). The L errors for the velocity and den-sity, computed using splash (Price 2007), are ∼ (cid:46) The plane-parallel flow described in this paper is com-parable to the flow derived by Adams (2011) for mag-netically controlled outflows from hot Jupiters whenthe stellar magnetic field completely dominates over themagnetic field produced by the planet (i.e. their β → ∞ limit). The main difference between our analytic solu-tions stems from assuming different gravitational po-tentials; however, to apply their solution to photoevap-orating discs would require readers to rederive the equa-tions themselves. The solution in this paper is signifi-cantly more transparent and its closed form makes itespecially easy to implement as a numerical test.Although we focus on using our plane-parallel modelas a numerical test, it may have use in wider applica-tions as well. For example, the model’s simple geometry,accurate approximation of winds close to the disc, andclosed form wind solution could make it a perfect spring-board for developing an analytic or semi-analytic modelfor coupled two-phase photoevaporation. To date, fewstudies have focused on dust dynamics in winds and asimple two-phase model would be very valuable. Acknowledgements
We thank Daniel Price and Sarah Maddison for usefuldiscussions and the anonymous referee whose careful re-port significantly improved this paper and pointed outthe explicit form of the solution. M.H. acknowledgesfunding from a Swinburne University Postgraduate Re-search Award (SUPRA). G. Laibe acknowledges fund-ing from the European Research Council for the FP7ERC advanced grant project ECOGAL. The numericalcalculations were performed on the SwinSTAR super-computer at Swinburne University of Technology andthe subsequent visualisation was made using splash (Price 2007).
REFERENCES
Adams, F. C. 2011, ApJ, 730, 27Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006,MNRAS, 369, 216Begelman, M. C., McKee, C. F., & Shields, G. A. 1983,ApJ, 271, 70Bondi, H. 1952, MNRAS, 112, 195Corless, R., Gonnet, G., Hare, D., Jeffrey, D., & Knuth,D. 1996, Advances in Computational Mathematics, 5,329
PASA (2018)doi:10.1017/pas.2018.xxx
Hutchison & Laibe - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - r=r s - - - - z=R
0 50 100 150 v = c s Figure 1.
Comparison between the plane-parallel wind at R = 5 AU (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 theplane-parallel case is only asymptotically crossed as z → ∞ whereas the Parker Wind model is transonic at r s = G M/ c . v z [ k m / s ] -505 t=100 [yr] l og ρ [ g / c m ] z [AU]-100 -50 0 50 100-20-15-10 Figure 2.