The rate of thermal atmospheric escape. Andrei Gruzinov CCPP, Physics Department, New York University, 4 Washington Place, New York, NY 10003 1 1 0 2 ABSTRACT n a A formula is derived for the rate of thermal atmospheric escape, valid, and asymptotically J exact, at low Knudsen number. 5 1. The thermal escape problem ] P Consider a planetary atmosphere without stir- E ring and radiative heating or cooling above the . h nominalsurface. Theatmospherecannotbestatic p since a static atmosphere reaches thermodynamic - o equilibrium. The Boltzmann distribution would r then give a finite density at infinite distances, t s which is unphysical. a 0 [ The atmosphere must be in a state of perma- nent escape. We want to calculate the rate of es- 1 v cape and the resulting temperature and density 3 profiles. 0 The thermal escape problemhas a long history -5 1 (Jeans 1904, Parker 1958, Hunten 1982). But 1 in recent applications to planetary bodies, there . 1 hasbeendisagreementonwhichmodel–hydrody- 0 namic or free molecular flow – should be applied 1 1 (Johnson 2010). We believe that all issues have -10 : been finally settled by a direct molecular dynam- v ics simulations of Volkov et. al. (2010). Here we i X just show how the Volkov et. al. (2010) results r can be obtained analytically. a -15 It would seem that the escape rate calculation 0 10 20 30 40 is only possible via molecular dynamics or Boltz- mann numerical simulation, because the exobase Fig. 1.— Escape rate for Kn = 10−8, c = 5/2, region,wherethetransitionfromhydrodynamicto p α = 1. Je < 2.5 is the Parker regime, 2.5 < Je < free molecular flow occurs, cannot be treated an- 24 is the Fourier regime, Je > 24 is the Jeans alytically. But we show, for low Knudsen number regime. at the surface, that the exobase boundary condi- tions are either irrelevant or can be deduced due to large overlap between hydrodynamic and free molecular flow regions. This allows one to calcu- late the atmospheric escape using hydrodynamics and to derive a formula for the escape rate. For clarity, we approximate the molecular in- 1 teraction potential by a power law (this includes Formula (5) limits the Je range over which the hard-spheregasasalimitingcase). Then,byscal- Parker rate applies. The Parker regime is de- ing arguments, the escape rate must be given by scribed in 3. § M is the Jeans (1904) escape rate. It de- Jeans Ma=M (Je,Kn). (1) escape pends on Je only: Here we measurethe escape rate by the radialve- (Je+1)exp( Je) locity at the surface v, and the velocity is repre- MJeans = √2π − (7) sented by the isothermal Mach number Formula (5) extends the Kn range over which the v Jeans rate applies. The Jeans regime is described Ma . (2) ≡ T/m in 4. § p M is proportional to Kn, with the Je- Fourier The Jeans number is dependent coefficient of proportionality: GMm Je , (3) MFourier(Je,Kn)=F(Je)Kn. (8) ≡ TR Here F(Je) depends on Je, c , and p where G is the gravitational constant, M is the dlnκ planetmass,m isthe moleculemass, T isthe sur- α . (9) face temperature, R is the surface radius. The ≡ dlnT Knudsen number is Kn λ, where λ is the mean For example, for α=1, ∼ R free path. We define 1 F(Je)= (10) mκ Je c Kn , (4) − p ≡ Rρ T/m For uniformity of notation, we call it the Fourier p regime. The Fourier regime is described in 5. where κ is the thermal conductivity, ρ is the den- § To illustrate the formula, Fig. 1 shows the es- sity. cape rate as a function of Je, for Kn = 10−8, The escape rate formula is given in 2. The § α = 1, and c = 5/2. The different regimes were formula includes three branches. In 3-5, we give p § extended to the intersection points. Fig.1 is in the recipes for calculating the branches. In 6,we § qualitativeagreementwithmoleculardynamicsre- explain the origin of the recipes, describe numeri- sults of Volkov et. al. (2010). Quantitative com- calhydrodynamicswhichconfirmstheescaperate parison is done in 6, where we also show that formula,andcompareourresultstothemolecular § a simplified dissipative hydrodynamics model re- dynamics simulations of Volkov et. al. (2010). producesthemoleculardynamicsresultstoagood accuracy. 2. The escape rate formula In the next three sections we calculate the es- In the low-Knudsen limit, cape rate postulating the properties of the flow. The assumptions made can be shown to be self- consistent. But, in fact, what makes the results MParker(Je), Je<cp; trustworthy is that they do match the numerical Mescape = MFourier(Je,Kn), cp <Je.ln(Kn−1); hydrodynamics results described in 6, and they MJeans(Je), Je&ln(Kn−1). also match the molecular dynamics si§mulations of (5) Volkov et al (2010). Here c is the heat capacity. The formula is p The Parker regime is a high-velocity ideal flow asymptotically exact (Appendix). with negligible thermal conductivity. The Fourier MParker is just the Parker wind (Parker 1958) regime is a flow with algebraically small velocity escaperate. ItdependsonJeandcp. Forexample, and with thermal conduction balancing the ideal- for cp =5/2, hydrodynamic energy flux. The Jeans regime is a flow with exponentially small velocity, such that MParker = 5/3. (6) thermal conduction dominates. p 2 3. Parker regime “algebraic” equation which can be derived from eq.(12)2. For Je < c , Kn 1 one can use ideal hydro- p ≪ dynamics 1. To write equations in the simplest 4. Jeans regime form, we use scalings to set M =m=R=1, and also ρ = 1, T = 1, at the surface R = 1. Now Fix the Knudsen number and let the Jeans we need to calculate the outflow velocity at the numberincrease. Thentheescaperateapproaches surface v =M , for any Je=G<c . 0 Parker p (Je+1)exp( Je) We write stationary ideal hydrodynamics as MJeans = − (14) √2π flux conservation, energy flux conservation, and adiabaticity – the standard Jeans rate – density normal- izedfluxofpositive-energypositive-radial-velocity ρvr2 =v , 0 molecules, calculated for Maxwelliandistribution. ρv22=+TccppT−1−. Gr = v202 +cp−G, (11) maEtiqo.n(1(49)%isearnrorexaptonJeent=ial2ly.5,ac0c.u0r0a8t%e aeprrporroxait- Je = 10) of the true escape rate for collision- From eqs.(11) less molecular outflow with Maxwellian injection at the surface, which in our notations is Thivs22“+alcgpeb(cid:16)rrav2i0vc”(cid:17)cepq1−u1a=tiovn202fo+rcvp−mGust+hGrav.e s(o1l2u)- MJeans = R0∞0∞uudduuR−∞−∞∞∞ddvvffv. (15) tions for all r > 1, such that the corresponding R R Herevistheradialvelocity,uistheabsolutevalue density ρ=v /(r2v) goes to zero at large r. This 0 ofthetangentialvelocity,f θexp( (u2+v2)/2) requirement limits the range of possible values of ∝ − is the truncated Maxwellian. θ = 0 if simultane- v . 0 ously the energy is positive, (u2+v2)/2 Je>0, For c < 3, and also for c > 3, G < − p p and the radial velocity v is negative; otherwise 2c /(c 1), the smallest possible value is v = p p 0 θ =1. − c /(c 1), in which case the sonic point is p p − pat the surface. For cp > 3 and larger G, 5. Fourier regime G > 2c /(c 1), the sonic point moves out to p p − r >1, and v0 decreases. FixtheJeansnumberJe>cpandlettheKnud- It is the smallest possible v that we call sen number decrease. The escape rate then ap- 0 M . Time-dependent ideal hydrodynamics proaches Parker simulationsandsimulationsofVolkovetal(2010) MFourier =F(Je)Kn. (16) confirm that the flow does approach the smallest The coefficient F(Je) is calculated below. v allowed by stationary ideal hydrodynamics. 0 In the Fourier regime, the negative hydrody- The above results can be written as namic energy flux carried by the flow is almost exactlycompensatedbythethermalconductivity: c p M = , (13) Parker rcp−1 r2κT′ =Φ(c T G). (17) p − r for c <3, and also for c >3, Je<2c /(c 1). p p p p For c > 3, 2c /(c 1) < Je < c , M −de- Here Φ = ρvr2 =const is the flux per steradian, p p p p Parker creases, reaching zer−o at Je = c . The precise M = m = 1 was set by scaling, prime denotes p value of M is then given implicitly by an the r-derivative, kinetic energy and viscosity are Parker negligibly small. 1 More precisely, dissipative effects must become negligible Usingscalings,weputT =1atr=1,andthen after a possible boundary layer is passed. Whether the solveeq.(17)forr >1. WeadjustΦinsuchaway boundarylayerispresent,dependsonhowtheatmosphere isfed. Iftheboundarylayerispresent,wesimplymovethe nominalsurfaceup,beyondtheboundarylayer. 2 (cid:16)M2a2 +cp−Je(cid:17)cp−52 =(cid:16)cpc−p1(cid:17)cp−1(cp− 52)cp−254JMe2a. 3 as to get the temperature going to zero at large usestationarynon-idealhydrodynamics. Theflux radii (otherwise the flow either stalls or becomes conservation,theenergyfluxconservation,andthe singularat finite radii,neverreachinganexobase, momentum equations are 6). We want to solve the problem for arbitrary § Φ=ρvr2, α≡Fidrsltnκco/ndslindTer3j.ust the case α = 1. One can Qρvv=′ =Φ(v22(ρ+Tc)p′T −GρGr+)−4 134η(rr23vη((vv′′− vrv)))−′.r2κT′, check numerically or analytically, that there is − − r2 3r3 − r (25) only one solution of eq.(17) with correct asymp- Here Φ, Q are constant fluxes, η is the viscosity4. totic behavior at large r: The boundary conditions at r =1 are 1 T = (18) ρ=1,T =1,v =Φ (26) r GivenG=Jeandκ(T =1)=Kn,wemustfind Substitution into eq.(17) then gives a pair (Φ,Q) that gives the “true” flow 5. Which κ flows are “true” is decided by the boundary con- =(G c )Φ. (19) T − p ditions at the exobase, and it would appear that onecannotproceedwithoutamoleculardynamics Or, in different notations, simulation of the exobase region. 1 But it turns out that for Kn 1 the mere F(Je)= . (20) ≪ Je c requirement that the exobase exists allows to cal- p − culate the escape flow in the Parker and Fourier For arbitrary α < 1, numerical integration of regimes. In the Jeans regime, this trick no longer eq.(17) shows that there is still a unique choice of works– there are many outflows with anexobase. Φ which gives a non-singular zero-at-infinity tem- What allows to solve the problem in the Jeans perature. One can show that for 0<Je c 1, regime is a large overlap between hydrodynamic p − ≪ and free molecular flow regions. 1 F = . (21) Parker regime – Take G < cp, and, as always, Je c − p κ = Kn 1. If the initial velocity is large, ≪ v Φ 1,the non-idealhydrodynamics(25)de- while for Je 1, 0 ≡ ∼ ≫ generates into the ideal hydrodynamics (11), and 2 1 we get the Parkerescape. F = . (22) 1+αJe Only for small velocities, v . κ, dissipative 0 effects can be significant. But numerical integra- Joining the asymptotics as follows tion of (25) shows that for small v the exobase 0 1 is never achieved. Kn stays small at all radii, no F(Je)= (23) (Je c )f(Je) matter what energy flux Q one chooses. p − What sets the boundary of the Parker regime 1 α Je c f(Je) 1 − − p . (24) to Je = cp is the sign of the ideal hydrodynamic ≈ − 2 Je 0.55c p energy flux in the small-velocity limit − gives an error less than 0.3% for all relevant α, cp G Q =Φ(c T ). (27) and Je. h p − r 6. Hydrodynamics 4 Wehaveassumedthatthesecondviscosityvanishes. This is trueforthe hard-spheregas and allows precisecompar- To show that the escape rate formula (5) is ison to the results of Volkov et. al. For polyatomic gases, thesecondviscositymustbeincluded,butitturnsoutthat asymptotically exact at Kn 0, and to improve → the viscosity effects are negligible at low Knudsen, unless the accuracy in the non-asymptotic regime, we Jeisveryclosetocp. 5 Tosolve(25)wealsoneed v′(1). Butthis boundary con- 3α = 1 is a good approximation for real gases at relevant dition is irrelevant. Wrong v′(1) would resultin abound- temperatures; α = 0.5 for the hard-sphere gas, and we arylayer,andwehaveagreedtoplacethenominalsurface wanttoreproducetheresultsofVolkovet. al. (2010) abovetheboundarylayer. 4 For Je < c this is positive. Now we show that nearly constant temperature profile. Eq.(28) then p this means no exobase for small v . gives an exponential density. These temperature 0 Forsmallv ,thedensityisapproximatelygiven and density are also predicted by the free molec- 0 by the equilibrium condition ular outflow at high Je. We see that the two de- scriptions agree. Gρ (ρT)′ =0. (28) We can therefore use hydrodynamic descrip- − − r2 tion well above the point r = R where the e AtmoderateG,theonlywaytoreducethedensity density scale height becomes comparable to the enough to get an exobase, is by decreasing the mean free path, Kn(Re) = Je−1(Re). But for temperatureatlargedistances. Ifthetemperature Kn(r) Je−1(r), we can also use the free molec- ≫ remains finite at infinity, the density saturates as ular outflow results to calculate the fluxes (Φ,Q). Now, we use hydrodynamics to integrate back to G r = 1. One can show, or just check numerically, ρ exp( ). (29) ∝ Tr that this procedure, for Je lnKn−1, gives the ≫ Jeansrate(toleadingorderinJeintheprefactor). Atsmallv ,we candropthe v2 termsfromthe 0 Simplified hydrodynamics – Both the Fourier energy flux expression in (25), and get and the Jeans regimes can be described simulta- G neously, in a much simplified version of hydrody- r2κT′ =Φ(c T ) Q. (30) p − r − namics which we give below. One can then hope that even between the asymptotics this hydrody- For the temperature to decrease, we must take a namicdescriptionwillwork(itdoesworkforhard positive flux Q. But finite Q then dominates the spheres, see below). r.h.s. of eq.(30) at large radii, and we get ThemodelhasoneadjustableparameterKn 0 ∼ r2κT′ = Q. (31) 1. We set T = 1, ρ = 1, Kn= Kn0 at r = 1, and − use the free molecular flow to calculate Φ and Q This gives at r =1: − 1 T ∝r 1+α. (32) Φ= e−G(G+1), Q= e−G((cp− 32)G+cp− 21). This temperature decreases – but too slowly. For √2π √2π α > 0, eqs.(28, 32) give a density growing at in- (33) finity. We then integrate eq.(30) backwardsin radius, to In summary, small-velocity flows have no r < 1. Simultaneously we integrate eq.(28) and exobase, and large-velocity flows are ideal, giv- calculate Kn(r). The integration is terminated ing the Parker flow. The lack of collisions beyond when the desired value of Kn is reached. We ad- the exobase cannot change the escape rate of the just G so that the desired value of Je is reached highly supersonic outflow. simultaneouslywiththe desiredvalueofKn. This Fourier regime – At fixed Je > c , as Kn de- gives the escape rate for these Kn and Je. One p creases, the set of pairs (Φ,Q) which give solu- also gets correct temperature and density profiles tions with an exobase shrinks to a single point for G 1 and also for r 1 at any G. We have ≫ ≪ (M ,0). If (Φ,Q)=(M ,0), the flow ei- checked, by numerical integration of the full sys- Fourier Fourier 6 ther becomes singular at a finite radius or goes to tem (25), that neglecting viscosity and inertia is infinite radii, without ever reaching local Kn 1. justified, again for G 1 and also for r 1 at ∼ ≫ ≪ So in this case too, we can claim that M is any G. Fourier an exact asymptotic even without knowing what We have used the simplified hydrodynamics happens at the exobase. (30,33) to establish the logarithmic accuracy of Theonlywayto reachanexobaseistosetQ= the escape rate formula (5) in the Fourier-Jeans 0 in eq.(30), and then tune Φ, as described in 5. transition region (Appendix). § Jeans regime – Fix Kn and let Je increase. Comparison with Volkov et. al. (2010) – Now The outflow velocity becomes so small that ther- comparethepredictionsofourescaperateformula malconductionineq.(30) dominatesandwe geta eq.(5) to the molecular dynamics simulations of 5 Volkov et. al. (2010). We put c =5/2, α=1/2. the Parker wind model ( 3). Both temperature p § WealsorememberthatourdefinitionofKnisdif- and density decrease algebraically. In the Jeans ferent by a constant factor from the conventional regime, we have a nearly isothermal slowly escap- hard-sphere Knudsen Kn =0.34Kn. ing atmosphere, with exponential density. In the hs We find that for Je = 3,10,15, Kn = Fourierregime,thermalconductivity balancesthe hs 10−2,10−3,10−4 the simplified hydrodynamics ideal energy flow ( 5). The temperature and the § (30,33) (with Kn = 3) predicts the escape rates density decrease algebraically. 0 which agree with the results of Volkov et. al. (2010) to about 10% for all points but one. I thank Robert Johnson for useful discussions. Let us discuss this discrepancy. For Je = 3, the escape rate formula (5) gives for Kn = Appendix hs 10−2,10−3,10−4 (using eqs.(16, 23, 24)) The escape rate formula (5) is asymptotically exact in the following sense: Ma =6.4, 6.4, 6.4 (34) Knhs Je=const<c ,Kn 0: Mescape 1, (37) p → M → Parker The simplified hydrodynamics (33) gives M escape Ma Je=const>cp,Kn 0: 1, (38) =5.6, 6.0, 6.2 (35) → MFourier → Kn hs M escape Kn=const,Je : 1, (39) Fig.4 of Volkov et al (2010) gives →∞ M → Jeans ln(M ) Ma Je>c ,Kn 0: escape 1. =5.0, 6.8, 4.1 (36) p → ln(min(M ,M )) → Kn Fourier Jeans hs (40) The simplified hydrodynamics works as expected In words: (i) As Knudsen decreases at fixed – as we go deeper into the Fourier regime, pre- Jeans, the Parker and Fourier formulas approach dicted Ma approaches MFourier. The first two val- the true escape rate everywhere but near the ues of Volkov et. al. are also reasonable – the Parker-Fourier transition. As Knudsen decreases, error of MFourier decreases as we get deeper into the Jeans number width of the Parker-Fourier theFourierregime. Butthelastpointshouldhave transition region shrinks to zero. (ii) As Jeans been evencloserto the Fouriervalue of 6.4. Some increases at fixed Knudsen, the Jeans formula ap- 50% error is unexpected. proachesthe true escaperate. (iii) Inthe Fourier- We believe that it is the molecular dynamics Jeans transition interval, both formulas are valid, results at Je=3, Kn =10−4 which have a large but only with logarithmic accuracy. hs error. Wehavechecked,byintegratingthefullsys- tem (25) numerically, that a hydrodynamic flow, REFERENCES startingatKn =10−4,Ma=4.1 10−4,Je=3, hs never goes above Kn =1.5 10−×3 for any Q. It Jeans, J. H., 1904, “The dynamical theory of hs × gases”, Chapter XVII, University Press, Cam- seems unlikely that molecular dynamics can no- bridge. ticeably deviate from hydrodynamics at such low Knudsen numbers. Parker,E. N., 1958, ApJ, 128, 664 Hunten, D. M., 1982, Planet. Space Sci., 30, 773 7. Summary Johnson, R. E., 2010,ApJ, 716, 1573 Formula (5) gives the thermal escape rate at Volkov, A. N., et al, 2010, arXiv:1009.5110 low Knudsen. The escape occurs either in Parker, in Fourier, or in the Jeans regime. In the Parker regime, we get a supersonic ideal hydrodynamical wind. This2-columnpreprintwaspreparedwiththeAASLATEX The temperatureanddensity profilesaregivenby macrosv5.2. 6