Magnetohydrodynamic Simulations of Hot Jupiter Upper Atmospheres
aa r X i v : . [ a s t r o - ph . E P ] A p r Preprint typeset using L A TEX style emulateapj v. 5/2/11
MAGNETOHYDRODYNAMIC SIMULATIONS OF HOT JUPITER UPPER ATMOSPHERES
George B. Trammell, Zhi-Yun Li & Phil Arras
Department of Astronomy, University of Virginia, P.O. Box 400325, Charlottesville, VA 22904-4325
ABSTRACTTwo-dimensional simulations of hot Jupiter upper atmospheres including the planet’s magneticfield are presented. The goal is to explore magnetic effects on the layer of the atmosphere that isionized and heated by stellar EUV radiation, and the imprint of these effects on the Ly α transmissionspectrum. The simulations are axisymmetric, isothermal, and include both rotation and azimuth-averaged stellar tides. Mass density is converted to atomic hydrogen density through the assumptionof ionization equilibrium. The three-zone structure – polar dead zone, mid-latitude wind zone, andequatorial dead zone – found in previous analytic calculations is confirmed. For a magnetic fieldcomparable to that of Jupiter, the equatorial dead zone, which is confined by the magnetic field andcorotates with the planet, contributes at least half of the transit signal. For even stronger fields, thegas escaping in the mid-latitude wind zone is found to have a smaller contribution to the transit depththan the equatorial dead zone. Transmission spectra computed from the simulations are comparedto HST STIS and ACS data for HD 209458b and HD 189733b, and the range of model parametersconsistent with the data is found. The central result of this paper is that the transit depth increasesstrongly with magnetic field strength when the hydrogen ionization layer is magnetically dominated,for dipole magnetic field B &
10 G. Hence transit depth is sensitive to magnetic field strength, inaddition to standard quantities such as the ratio of thermal to gravitational binding energies. Anothereffect of the magnetic field is that the planet loses angular momentum orders of magnitude faster thanin the non-magnetic case, because the magnetic field greatly increases the lever arm for wind brakingof the planet’s rotation. Spin-down timescales for magnetized models of HD 209458b that agree withthe observed transit depth can be as short as ≃
30 Myr, much shorter than the age of the system.
Subject headings: (stars:) planetary systems - (magnetohydrodynamics:) MHD INTRODUCTION
Hot Jupiters are gas giants orbiting close to their par-ent stars. The large stellar EUV flux heats and ion-izes the upper atmosphere of these planets, increasingthe thermal energy to a value approaching the gravita-tional binding energy, leading to a region weakly boundto the planet. The resulting large gas scale heights andatmospheric escape form an extended upper atmospherearound the planet, which may be probed by transmissionspectroscopy using strong atomic resonance lines.The existence of an extended upper atmosphere hasbeen established through a variety of observations. Spec-troscopic UV observations of HD 209458b (Henry et al.2000) indicate a ∼
10% decrease in flux during transitat ∼
100 km s − from the center of the hydrogen Ly α line. This transit depth has been attributed to an atmo-sphere of neutral H extending to a radius ≃ . R p , where R p is the radius of the broadband photosphere of theplanet (see Vidal-Madjar et al. 2008). As this radius iscomparable to the Roche lobe radius, Vidal-Madjar et al.(2004) suggested that the planet is losing mass throughRoche lobe overflow.Additional observations of HD 209458b at transithave indicated absorption in other resonance lines, in-cluding NaI (Charbonneau et al. 2002; Sing et al. 2008),OI (Vidal-Madjar et al. 2004), CII (Vidal-Madjar et al.2004; France et al. 2011), and SiIII (Linsky et al. 2010;France et al. 2011). Follow-up observations and re-analysis of HST-ACS data, in comparison with HST-STIS low and medium-resolution spectra, confirmed the [email protected], [email protected],[email protected] reduction of Ly α flux (Ehrenreich et al. 2008).Transmission spectra of the hot Jupiter HD 189733bhave also revealed absorption due to HI (in both Ly α andH α ; Lecavelier Des Etangs et al. 2010 and Jensen et al.2012) and NaI (Redfield et al. 2008; Snellen et al. 2008).HST-COS observations by Linsky et al. (2010) have in-dicated absorption at up to ±
50 km s − from line centerin CII and SiIII that may be indicative of high veloc-ity absorbers in the upper atmosphere (although theseobservations probe deeper layers than those probed bythe HI observations). Multi-epoch spectra have alsorevealed significant changes in the Ly α transit depth,which are correlated with flares in ionizing radiationfrom the host star detected with HST and SWIFT(Lecavelier des Etangs et al. 2012).This paper will focus on the Ly α absorption ob-served in the upper atmospheres of HD 209458b and HD189733b. One interpretation of this absorption invokeshydrogen with thermal velocity ∼
10 km s − with sucha large column density that the damping wings of Ly α become optically thick (e.g. Yelle 2004). Alternatively,a much smaller column is required if hydrogen atomsat thermal velocities ∼
100 km s − , created by chargeexchange with stellar wind protons (Holmstr¨om et al.2008; Ekenb¨ack et al. 2010; Tremblin & Chiang 2013),produce a sufficiently broad line profile.This paper considers the former scenario in which thetransit depth is due to a layer of thermal hydrogen in theplanet’s atmosphere. A number of studies have alreadyexplored the properties of strongly irradiated exoplanetatmospheres, and the possibility of thermally-drivenhydrodynamic outflow (Yelle 2004, 2006; Tian et al.2005; Garc´ıa Mu˜noz 2007; Murray-Clay et al. 2009;Ehrenreich & D´esert 2011) and/or Roche Lobe over-flow (Gu et al. 2003; Li et al. 2010; Lai et al. 2010;Ehrenreich & D´esert 2011). The present study standsapart from the previous ones by including, through de-tailed magnetohydrodynamic (MHD) simulations, the ef-fect of the planetary magnetic field. It is a follow-up ofTrammell et al. (2011), which considered the magneticeffects semi-analytically.Trammell et al. (2011, hereafter, Paper I) showed thatthe addition of the planetary magnetic field leads to theformation of an equatorial “dead-zone” (DZ) — a staticregion where the wind ram pressure is insufficient to over-whelm magnetic stresses and open the field lines into anoutflow. This effect is well known in the classical MHDstellar wind theory (e.g. Mestel 1968). Paper I found asecond static region near the poles where the wind can beshut off by the increased gravitational potential barrierfrom the stellar tide. In the strong tide limit, a wind-zone(i.e., the outflow region; WZ) is then expected to existonly at intermediate latitudes. A goal of the presentpaper is to verify this analytically-obtained three-zonestructure with detailed numerical simulations.Another conclusion from Paper I was that observationsof Ly α absorption at the 5-10% level for HD 209458b maybe detecting neutral H which is collisionally coupled toionized gas confined to the equatorial DZ by the planet’smagnetic field. The bulk of the absorbing gas observedat transit thus may not be escaping, but rather is inthe static equatorial dead zone. This qualitative resultdiffers from the basic assumption in the hydrodynamicescape and Roche Lobe overflow models, that the transitobservations are probing gas in the act of escaping fromthe planet.A limitation of the analytic models in Paper I isthat they ignore magnetic forces, which means that thepoloidal magnetic geometry was assumed rather thancomputed self-consistently. Another limitation is thatthe fluid was assumed to corotate with the planet every-where, including the wind zone, where the corotation isexpected to break down at large distances. The treat-ment in this paper overcomes these limitations by per-forming MHD simulations, which compute the magneticfield structure and fluid rotation self-consistently. Thisallows a more accurate calculation of the mass and an-gular momentum loss rates, as well as the density andvelocity profiles required to compute transmission spec-tra.The plan of the paper is as follows. Section 2 de-scribes the simulation setup and model parameters, andSection 3 presents the simulation results. Section 4 de-scribes the method for computing model Ly α spectra.The frequency-dependent and frequency-integrated tran-sit depths for a range of simulation parameters are com-pared with observations of HD 209458b and HD 189733b.Findings are summarized in Section 5. The Appendixcontains a discussion of numerical effects at the shearlayer separating the dead and wind zones. SIMULATION SETUP
Consider a planet of mass M p and radius R p in a cir-cular orbit at a distance D from a star of mass M ⋆ . Theplanet’s rotation is synchronized to the orbit with angu-lar velocity Ω p = [ G ( M ⋆ + M p ) /D ] / , and the spin axis is aligned with the orbital angular momentum. Outsidethe planet, in the region modeled by the simulations, thegas is not required to corotate with the planet.Two-dimensional (2D), axisymmetric simulations inspherical coordinates are carried out with the publiclyavailable MHD code ZEUS-MP (see Hayes et al. 2006,and references therein), which solves the ideal-MHDequations: DρDt = − ρ ∇ · v (1) ρ D v Dt = −∇ P + 14 π ( ∇ × B ) × B − ρ ∇ U (2) ∂ B ∂t = ∇ × ( v × B ) , (3)where the comoving (Lagrangian) derivative is definedas DDt ≡ ∂∂t + v · ∇ . (4)Equations 1-3 are the mass continuity, momentum andinduction equation, respectively. The neglect of explicitfluid viscosity in Equation 2 and resistivity in Equation3 are discussed in Appendix B of Paper I. The quantity U is an effective potential to be defined below, and othersymbols have their standard meaning.Instead of solving the energy equation, an isothermalequation of state, P = ρa , is used, where a is the(constant) isothermal sound speed. This assumption isequivalent to adding energy to the flow to counter adi-abatic cooling. It gives rise to a transonic outflow (e.g.Lamers & Cassinelli 1999). The isothermal assumptionis convenient for the present study, where the focus is noton the initial launching of the wind, but rather on mag-netic effects. A more detailed study, beyond the scopeof this paper, would include heating and cooling effectsin an energy equation. Note, however, that since themagnetic field and rotation are included, the flow is alsoaccelerated in part by the “magneto-centrifugal” effect(Blandford & Payne 1982), as well as stellar tides.The computational grid extends from an inner radialboundary at r = R p to the outer boundary at r = 30 R p ,with R p set to the planet’s observed transit continuumradius (Southworth 2010), and from the north pole at θ = 0 to the south pole at θ = π . The radial box sizewas chosen through experimentation so that all MHDcritical points in the wind zone were contained withinthe computational domain for a wide range of model pa-rameters. The standard resolution is 272 × r = r i +1 − r i increasing outwardaccording to ∆ r i +1 / ∆ r i = 1 .
02; the ratio was chosento adequately resolve the wind acceleration region nearthe base. The θ grid is uniformly spaced. Surroundingthe active grid are two layers of ghost zones at each ofthe four boundaries; they are used to impose boundaryconditions.At time t = 0, the fluid is uniformly rotating withvelocity ( v r , v θ , v φ ) = (0 , , Ω p r sin θ ), and in hydrostaticbalance over most of the computational grid. The initialmagnetic field is assumed to be a potential field, so thatthe magnetic force is everywhere zero, an assumptionconsistent with the equation of hydrostatic balance. Ina reference frame corotating with the planet, and withthe origin comoving with the planet, hydrostatic balancetakes the form (Paper I)0 = − a ∇ ρ − ρ ∇ U rot . (5)The potential U rot includes contributions from the grav-ity of the planet, the stellar gravity, the dipole term aris-ing from the acceleration of the origin, and the centrifu-gal force, and takes the form U rot ( x ) = − GM p | x | − GM ⋆ | x − x ⋆ | + GM ⋆ x · x ⋆ | x ⋆ | − | Ω p × x | (6)The dipole term, which acts to accelerate the center ofmass of the planet, cancels off part of the stellar gravity,leaving only a tidal acceleration. Expressing the positionvector in spherical coordinates x = ( r, θ, φ ), the positionof the star as x ⋆ = ( D, π/ ,
0) = D e x , and making thetidal approximation, r ≪ D , gives U rot ( x ) ≃ − GM p r −
12 Ω p r (cid:0) f rot sin θ − (cid:1) , (7)where the longitude-dependent function f rot = 1 +3 cos φ . As the simulations are axisymmetric, we sub-stitute the azimuthal average cos φ → /
2, yielding f rot = 5 /
2. At the equator, the radial acceleration ∂U rot /∂r = 0 at the radius r H = (cid:18) GM p p (cid:19) / ≃ D (cid:18) M p M ⋆ (cid:19) / , (8)a factor of 2 / larger than the physically correct valuefor the Lagrange points, which are evaluated along thestar-planet line. In this paper r H will be called the Hillradius.Substituting Equation 7 into Equation 5, the initialdensity distribution over the inner part of the computa-tional domain takes the form ρ ( r, θ ) = ρ ss exp (cid:20) − (cid:18) U rot ( r, θ ) − U rot ( R p , π/ a (cid:19)(cid:21) , (9)where ρ ss is the density at ( r, θ ) = ( R p , π/ a and ρ ss as param-eters, in addition to the parameters of the planet, starand orbit. Well outside r H , the hydrostatic density pro-file rises steeply to large values, due to the net accelera-tion pointing outward. For the initial condition only, aceiling is placed on the density at ρ ≤ ρ ss to limit thisgrowth. In the outer regions of the computational gridwhere ρ = ρ ss , the pressure force is then zero and thetidal and centrifugal forces pull mass outward, initiatingthe outflow. Since the physical, steady-state solutionsexhibit ρ ≪ ρ ss in the outer regions, capping the densityis a device to allow the density to decrease to physicallevels more quickly.It is simpler to perform the simulation not in a corotat-ing frame, which would require the addition of Coriolisand centrifugal forces, but rather in a non-rotating frame.The origin still moves with the center of the planet. Inthis reference frame, the centrifugal term can be omitted from the potential, giving U ( x ) = − GM p | x | − GM ⋆ | x − x ⋆ | + GM ⋆ x · x ⋆ | x ⋆ | ≃ − GM p r −
12 Ω p r (cid:0) f sin θ − (cid:1) , (10)where now f = 1 + cos φ . The azimuthal average gives f = 3 / r and θ compo-nents of − ∇ U from Equation 10 are introduced into theZeus-MP code as a source term in the momentum equa-tions. Since the gas is not required to corotate, Equation8 may underestimate the radius at which the equatorialacceleration changes sign. An upper limit is found byignoring the centrifugal force. Using f = 3 / GM p / Ω p ) / for this radius, largerthan the expression in Equation 8 by a factor 3 / .The initial condition for the magnetic field is a dipolewith magnetic axis aligned with the rotation axis: B r = B (cid:18) rR p (cid:19) − cos θ (11) B θ = B (cid:18) rR p (cid:19) − sin θ (12) B φ = 0 , (13)where the field at the magnetic pole is B . The devel-opment of nonzero B φ at t > β at the innerradius: β = 8 πP ss ( B / (14)where P ss = a ρ ss is the base pressure at the equator, and B / ρ at the innerradial boundary. The densities in the inner radial ghostzones and the first radial active zones are kept at theirinitial values at all times. Even though the densitiesin the first active zones are updated at each time step,the updated values are discarded and replaced by theirinitial values. This guarantees that the base density isheld fixed at the prescribed value, even in the wind zone.For the inner radial boundary, v r = 0 is set at the in-ner face of the first active zone, as well as in the ghostzones. In other words, v r ( R p , θ ) = 0 and v r ( r − , θ )=0,where r − ( < R p ) denotes the inner radial ghost region.The reflection boundary condition is applied to v θ , sothat v θ ( r − , θ ) = v θ ( r + , θ ), where r + is the symmetrypoint (with respect to the r = R p surface) in the ac-tive domain of the location r − in the ghost region. Theboundary condition on the azimuthal velocity compo-nent, v φ , is v φ ( r − , θ ) = 2 v φ ( R p , θ ) − v φ ( r + , θ ) where v φ ( R p , θ ) = Ω p R p sin θ . That is, the average of the firstghost zone and the first active zone should equal the coro-tation velocity. It has been verified that, in the absenceof magnetic field, rotation and stellar tides, the innerhydro boundary conditions produce a thermally drivenwind that matches, in steady state, the well-known ana-lytic solution.The magnetic boundary conditions at the inner ra-dial boundary are more complicated to implement. Theyare enforced through the electromotive force (EMF) ǫ = v × B , as this will automatically preserve ∂ ( ∇· B ) /∂t = 0during the time evolution. In 2D (axisymmetric) geom-etry, only the r - and θ -components of ǫ affect B φ : ∂B φ ∂t = 1 r (cid:20) ∂∂r ( rǫ θ ) − ∂ǫ r ∂θ (cid:21) . (15)Although B φ is assumed to be zero initially in our sim-ulation, it can grow with time, particularly in the out-flow region. The boundary conditions on ǫ r and ǫ θ aredesigned to enable B φ in the ghost zones to grow at thesame rate as in the active zones. Specifically, we demand ǫ r ( r − ) r − = ǫ r ( r + ) r + , (16)and (cid:20) r ∂∂r ( rǫ θ ) (cid:21) r − = (cid:20) r ∂∂r ( rǫ θ ) (cid:21) r + . (17)The value of ǫ θ in the ghost zone is determined by equa-tion 17 together with the condition ǫ θ ( R p , θ ) = Ω p R p sin θB r ( R p , θ ) , (18)which ensures that the footpoints of the magnetic fieldlines corotate with the planet.For ǫ φ = v r B θ − v θ B r , the boundary condition ǫ φ ( R p , θ ) = 0 enforces poloidal velocity parallel topoloidal magnetic field at r = R p . The use of this“flux freezing” condition for the magnetic field is jus-tified in Appendix B of Paper I. It also guarantees that B r ( r = R p ) remains unchanged, i.e., the footpoints of thefield lines are firmly anchored on the rotating inner radialboundary. For ǫ φ in the ghost zone, ǫ φ ( r − ) = − ǫ φ ( r + ) isenforced so that the radial gradient of ǫ φ , which controlsthe evolution of B θ , is continuous across the inner radialboundary. This set of magnetic boundary conditions issimilar to that used successfully by Krasnopolsky et al.(1999, 2003) to simulate disk-driven magnetocentrifugalwinds.The standard “outflow” boundary condition imple-mented in ZeusMP is used at the outer radial boundary,with all hydrodynamic variables and the three compo-nents of the EMF projected to zero slope. In addition, ∂B φ /∂t = 0 is set for the outer radial ghost zones, whichwas found to prevent the growth of unphysically largecurrents that sometimes develop near the outer bound-ary. At the θ = 0 and θ = π boundaries, the standard“axial” boundary condition as implemented in ZeusMPis used, which enforces reflection symmetry for the r -component of the velocity and magnetic field; the θ - and φ -components are reflected with a change of sign.The simulations were evolved until a steady-state so-lution was achieved. A summary of the main model pa-rameters is shown in Tables 1 and 2. TABLE 1Primary Model Parameters
Parameter Range Description R p R Jup planet radius ρ ss ∼ − − − g cm − substellar point mass density a − isothermal sound speed B M p M Jup planet mass M ⋆ M ⊙ host stellar mass D Note . — Description and range of the model parameters usedin the simulations.
TABLE 2Parameters for HD 209458b Runs
Run D (AU) P ss ( µ bar) a (km/s) B (G) β Model 1 0.047 0.05 10.0 10.0 0.051Model 2 0.047 0.05 10.0 1.0 5.1Model 3 0.047 0.05 10.0 50.0 0.002Model 4 0.047 0.05 9.0 10.0 0.041Model 5 0.047 0.05 11.0 10.0 0.061Model 6 0.047 0.005 10.0 10.0 0.0051Model 7 0.047 0.5 10.0 10.0 0.51Model 8 0.035 0.05 10.0 10.0 0.051Model 9 0.06 0.05 10.0 10.0 0.051Model 10 0.047 0.05 10.0 100.0 0.0005
Note . — This table contains simulation parametersvarying a single parameter ( B , D, a, P ss ) relative to the fidu-cial case (Model 1). The planetary radius and mass are fixedto R p = 1 . R J , M p = 0 . M J . The values of β can becompared to the value β , Jup = 0 .
069 using Jupiter’s mag-netic field and a base pressure P ss = 0 . µ bar.3. SIMULATION RESULTS
A fiducial model is chosen with HD 209458b’s param-eters and B = 10 G (Model 1 in Table 2). The othersimulations listed in Table 2 vary the model parameterslisted in Table 1. A qualitative discussion of the simu-lation results is given in § § § Qualitative Results: Magnetic Field and TidalStrength
One of the most important qualitative results of thispaper is the confirmation of the three-zone structure ofthe magnetosphere predicted analytically in Paper I. Thethree distinct regions are clearly visible in Figures 1 and2 — (1) an equatorial dead-zone (DZ) containing staticgas confined by the magnetic field, (2) a wind-zone (WZ)where an outflow is driven along open magnetic fieldlines, and (3) a second polar DZ where the stellar tidehas shut off the outflow (see also Fig. 7 of Paper I). Therange of magnetic field and stellar tide over which theequatorial and polar DZ’s exist has been discussed inPaper I. Roughly, the equatorial DZ requires β .
1, i.e.the magnetic pressure dominates gas pressure at the baseof the atmosphere (the hydrogen ionization zone) at theequator. The existence of the polar dead zone, and theinability to drive a transonic outflow there, occurs in-side a critical orbital separation. Roughly, this criterion
Fig. 1.—
Contours of total gas density ρ ( r, θ ) in units of g cm − as viewed during transit, illustrating the effect of the magnetic field(white lines) on the size of the equatorial DZ. Left : Model 1 (fiducial model with B = 10 G) Center : Model 2 ( B = 1 G) Right : Model3 ( B = 50 G). For these three models, the remaining parameters given in Table 2 are otherwise identical. White arrows indicate thedirection and magnitude of the poloidal fluid velocity. Only the inner 8 R p portion of the grid is shown (the magnetic field lines near theplanet are not drawn for clarity). Fig. 2.—
Similar to Figure 1, but for contours of poloidal velocity v p in units of the isothermal sound speed a . Left : Model 1 (fiducialmodel)
Center : Model 2
Right : Model 3 for parameters given in Table 2. White-dashed contours trace the (slow magneto-)sonic points,which trace the vicinity of the shear layers separating the static DZs (i.e., darkest contours) from the transonic WZs. translates into the rotation velocity Ω p at the fiducialsonic point radius r s , = GM p / a must be supersonic,Ω p r s , & a (cf. Equation 36 of Paper I).The DZ-WZ boundaries in the simulation contain ashear layer separating the outflowing gas in the WZ fromthe static gas in the DZ. In addition, the magnetic fieldchanges rapidly in this boundary layer, implying a cur-rent sheet. The origin of this current sheet is that, foridentical Bernoulli constant at the inner boundary, the WZ has smaller density compared to the neighboring DZby a factor ∼ exp( − v p /a ) (Mestel and Spruit 1987; Pa-per I), where v p is the poloidal wind speed. Since thetotal pressure, gas plus magnetic, must be continuousacross the boundary, the decrease in gas pressure impliesan increase in magnetic pressure, and hence a currentsheet. Numerical issues related to the shear in velocityand magnetic field will be discussed further in the Ap-pendix. Fig. 3.—
Similar to Figure 1, illustrating how the stellar tide influences the size of the equatorial/polar DZs. All panels have the sameparameters as Model 1 of Table 2, except for the varying orbital distance D . Left : the fiducial Model 1 for HD 209458b with D = 0 . Center : Model 8 ( D = 0 .
035 AU), which illustrates the effect of increasing the stellar tide by shrinking the planet’s orbit.
Right :Model 9 ( D = 0 .
06 AU). Note that the polar DZ size is larger and equatorial DZ smaller for the stronger tide case (middle panel), as canbe seen by the range of angles with zero-length velocity vectors (white arrows). As in Figure 1, the innermost magnetic field lines havebeen suppressed for clarity.
Fig. 4.—
Similar to Figure 3, showing the 2D structure of the poloidal velocity v p field, in units of the isothermal speed a for the samethree models shown in Figure 3. Higher stellar tide pushes the slow magnetosonic point inward at mid-latitudes and leads to a strongeroutflow ram pressure, which can open a larger region of the planetary magnetic field lines and decrease the size of the equatorial DZ. Figure 1 illustrates the effect of the magnetic fieldon the density profile in the magnetosphere. The pa-rameters for the runs in each panel are identical exceptfor the magnetic field, with B = 10 , θ , the magnetic pres-sure at the looptop at the equator decreases outward as B ∝ (1 − θ/
4) sin θ / sin θ , and so field linesnearer the pole, with smaller θ , suffer a larger decreasein magnetic pressure from pole to equator. The larger DZsize for larger B then reflects the inability of ram pres-sure to overcome magnetic pressure, except in a smallerregion near the pole where the field decreases outwardmore rapidly.The observational implication of the increase of DZ sizewith magnetic field is that more of the circum-planetarymaterial is expected to be confined within the static deadzone, which should make this region easier to probe withtransit observations (see Section 4). A weaker magneticfield B ≪ B , the structure of the mag-netosphere is also influenced by the stellar tide. The tidaleffects are illustrated in Figures 3 and 4. All parametersexcept D are held fixed, even though temperature wouldlikely increase as the planet is moved nearer the star.The left panel is the fiducial Model 1. In the middlepanel (Model 8), the orbital distance has been decreasedto D = 0 .
035 AU, so that the stellar tide is stronger thanthat for HD 209458b ( D = 0 .
047 AU). As predicted inPaper I, the stronger stellar tide increases the outwardacceleration of the mid-latitude outflow by moving thesonic point inward. It results in an equatorial DZ that isslightly smaller in size but denser at the same distancefrom the planet relative to Model 1. Figure 4 shows thatthe polar DZ size is also larger for the stronger tide case,as can be seen by the range of angles occupied by largelysubsonic gas with small poloidal velocities, again in broadagreement with the analytic results of Paper I.
Quantitative Analysis: Density and VelocityProfiles
To examine the numerical simulations in more detail,Figure 5 shows the run of density and poloidal velocityalong three different co-latitudes for the fiducial Model 1.The θ values are chosen to highlight the separate polarDZ, WZ and equatorial DZ regions, respectively.Along θ ≈
0, near the pole, the density decreasesrapidly with r because the downward gravity of theplanet and star must be balanced by pressure gradient inhydrostatic equilibrium. For this region, the flow speedremains well below the sound speed, in agreement withPaper I, which predicts the absence of a transonic solu-tion in the polar region. The θ ≈ π/ r = 8 R p . The densitydrops with distance accordingly.The density distribution along the equator at θ ≈ π/ r ≃ R p , and thenresumes a slow decline. Such a “bump” in the density profile was predicted in Paper I, for the gas outside theHill radius, yet still confined inside the static magneto-sphere. This is a consequence of the outward pointinggravity outside the Hill radius in Equation 8, causingthe density to increase outward instead of inward. How-ever, the outward increase in figure 5 occurs well inside r H ≃ . R p ! Hence it cannot be due to the change in thesign of gravity. In the Appendix the origin of this den-sity increase is explored, and seems to be due to viscousstresses associated with numerical effects near the equa-torial DZ/WZ boundary. As the numerical resolution isincreased, the density bump in figure 5 decreases in size.In the Appendix it is shown that increasing the resolu-tion has a .
1% effect on the integrated transit depth,even as the density bump decreases. This gives confi-dence that resolution-dependent effects are not leadingto large errors in the transit depth.
Fig. 5.—
1D radial profiles of density ρ ( r ) ( top ) and poloidalvelocity v p ( r ) ( bottom ) in units of the sound speed a for the fiducialModel 1 at three representative angles. The three angles shown ineach panel are θ ≈ θ ≈ π/ θ ≈ π/ r S . Mass and Angular Momentum Loss Rates
The MHD simulations presented in this paper allowa more accurate determination of the rates of massand angular momentum losses ( ˙ M and ˙ J ) as comparedto the semi-analytic solutions from Paper I, since herethe magnetic field geometry and fluid velocity are self-consistently computed. These quantities are computedas a function of r by integrals over θ :˙ M ( r ) = 2 πr Z π dθ sin θρ ( r, θ ) v r ( r, θ ) (19)and ˙ J ( r ) = 2 πr Z π dθ sin θ × (cid:20) ρ ( r, θ ) v r ( r, θ ) v φ ( r, θ ) − B r ( r, θ ) B φ ( r, θ )4 π (cid:21) . (20)Typically ˙ M ( r ) and ˙ J ( r ) are constant with radius to bet-ter than 1%, which provides a check on the accuracy ofthe numerical solutions. Table 3 summarizes the resultsfor ˙ M and ˙ J for Models 1-10.The planet’s magnetic field affects the dynamics in sev-eral ways. A stronger magnetic field increases the size ofthe equatorial DZ, restricting the WZ to a smaller rangeof latitudes. Therefore, one might expect that the mass-loss rate will decrease for a stronger magnetic field. Thisexpectation is born out in the ˙ M values presented inTable 3, where the mass loss decreases by ∼
35% for afactor of 5 increase in B from Model 1 to Model 3. De-spite the reduction in ˙ M , the total angular momentumloss rate ˙ J increased by a factor of ∼ .
5, implying anincrease in loss of specific angular momentum, ˙ J/ ˙ M , dueto a longer magnetic lever arm for the torque. The effectof the magnetic field and tides on the specific angularmomentum loss is most clearly displayed in column 4 ofTable 3. The quantity ˙ J/ ( ˙ M Ω p R ) has the value 2 / M islarger due to the larger range of latitudes in the WZ (seethe center panel of Figure 2).Stronger tides result in a slightly smaller equatorialDZ, because the outward tidal force can open more mag-netic field lines, but a larger polar DZ, due to the in-creased potential barrier. Stronger tide also moves thesonic point inward, which tends to increase ˙ M . Forexample, ˙ M of the stronger tide Model 8 is increasedslightly, by a factor of ≃ M comes from varying thebase pressure P ss (Models 6 and 7) or the isothermalsound speed a (Models 4 and 5). For example, when P ss increases by a factor of 10, from 0.05 to 0.5 µ bar,˙ M rises by a factor of 16.9. When a increases by 10%,from 10 to 11 km/s, ˙ M shoots up by a factor of 4.16! Inthe more heavily mass-loaded winds, the field lines bendbackward significantly in the azimuthal direction rela-tively close to the planet, forcing the fluid to rotate sub-stantially below the corotation speed. The self-consistenttreatment of the deviation from corotation here is an im-provement over the analytic solutions of Paper I. Con- TABLE 3Mass/Ang. Mom. Loss Rates (HD 209458b)
Run ˙ M ˙ J ˙ J/ ( ˙ M Ω p R p ) δF/F Model 1 3.29 10.20 162.24 0.100Model 2 6.12 1.17 9.97 0.125Model 3 2.11 26.40 655.64 0.157Model 4 0.48 3.74 404.68 0.048Model 5 13.70 19.05 72.81 0.209Model 6 0.25 2.20 455.36 0.028Model 7 55.50 37.59 35.46 0.470Model 8 3.89 15.24 131.92 0.154Model 9 3.21 7.04 165.70 0.088Model 10 2.87 33.40 610.59 0.253
Note . — Mass-loss rates ˙ M [10 g s − ]and angular momentum loss rates ˙ J [10 gcm s − ] for the 9 models with parametersspecified in Table 2, along with the correspond-ing integrated Ly α transit depth δF/F (seeEquation 26) from -200 to +200 km s − fromline center. versely, a smaller P ss or a leads to a lower ˙ M , and awind that is dominated by the magnetic field out to alarger distance. It is interesting to note that the ratio˙ J/ ( ˙ M Ω p R p ) has rather large values of 404.68 and 455.36for Model 4 ( a = 9 km/s) and 6 ( P ss = 0 . µ bar),respectively. They are very different from the purely hy-dro winds from the planet, where the ratio is 2 /
3. Therelatively low mass loss rate in these cases allows themagnetic field to effectively enforce corotation up to adistance of ∼ R p .The large spin-down torques found in the stronglymagnetized models may torque the planet away from syn-chronous rotation (Paper I). Defining Γ = ˙ J/ ( ˙ M Ω p R )and J = αM R Ω p , the spindown timescale is J ˙ J = α Γ (cid:18) M ˙ M (cid:19) ≃ × yr (cid:16) α . (cid:17) (cid:18) (cid:19) (cid:18) M . M J (cid:19) (cid:18) × g s − ˙ M (cid:19) , (21)for Model 3 parameters. In torque equilibrium be-tween magnetic spin-down torques and gravitationaltidal torques, a steady-state asynchronous spin ratewould occur, with associated steady-state gravitationaltide heating. However, deviations from synchronous ro-tation depend on strength of the planet’s tidal dissipa-tion, which is uncertain, but likely to give synchroniza-tion timescales orders of magnitude shorter than Equa-tion 21 (e.g. Wu & Murray 2003). For HD 209458b,the heating rate can be estimated to be far smaller thanJupiter’s luminosity (3 × erg) for gravitational tidesynchronization timescales shorter than 1 Myr. Hencethe magnetic spin-down torque and associated asyn-chronous rotation are not likely to give rise to a heat-ing rate large enough to effect the thermal history of theplanet significantly. TRANSIT DEPTHS IN LY α Section 3 described numerical solutions for the MHDvariables ρ and v for different model parameters. In thissection the mass density ρ is converted into atomic hy-drogen number density n H , and the transmission spectrafor the models in Table 2 are discussed.As a point of departure when considering transmis-sion spectra of the MHD simulation results, the simplemodel of Lecavelier Des Etangs et al. (2008) is first sum-marized. They consider a plane parallel, isothermal at-mosphere with base radius R b and altitude z = r − R b .The number density is then n ( z ) = n exp( − z/H ), where H = k b T / ( µm p g ) is the scale height, µ is the meanmolecular weight, and g = GM p /r . The path lengththrough the atmosphere is ℓ ≃ √ πR b H , giving an op-tical depth τ ν ( z ) = n σ ν ℓ exp( − z/H ), where σ ν is theLy α (1s → τ ν ( z ν ) = 1 givesthe altitude z ν ≃ H ln (cid:18) n σ ν √ πR b H (cid:19) (22)up to which the atmosphere is optically thick. The tran-sit depth is then R ( ν ) R ⋆ ≃ R + 2 R b z ν R ⋆ . (23)The altitude z ν ∝ H ∝ T /µg , so hot atmospheres oflow mean molecular weight gas around planets with lowgravity will have large scale heights and transit depths.Due to the steeply falling density, the transit depth hasonly a weak logarithmic dependence on σ ν .For the Ly α transit depths of the MHD models consid-ered here, the DZ is hydrostatic, but the tidal/rotationalforces are important, and so gravity is weaker than GM p /r . The corresponding larger scale heights makethe plane parallel limit inaccurate, and the density pro-file, even of isothermal models, tends not to fall as steeplyas it does deeper in the atmosphere. One consequenceof the large scale heights is that there can be a signif-icant contribution from gas with optical depth τ ≤ τ = 1 contour – may no longer beaccurate. Hence, for careful work numerical integrationsare required. However, the analytic model gives usefulintuition and is simple. Details of the Calculation
Stellar Ly α photons passing through the planet’s at-mosphere can be absorbed or scattered out of the line ofsight to the observer, causing a decrease in flux. In addi-tion, the interstellar medium (ISM) can absorb/scatterthe light, most prominently in the Doppler core of theline. The spectrum observed at Earth is the combina-tion of these two effects. If the in-transit flux is F ν andthe out-of-transit flux is F (0) ν , the fractional decrease influx, the transit depth, is ( F (0) ν − F ν ) /F (0) ν .The optical depth through the planet’s atmosphere isgiven by τ ν ( y, z ) = Z dx n H ( x, y, z ) σ ν ( x, y, z ) (24)where n H is the number density of the atomic hydrogenin the 1s state, x specifies the direction along the line ofsight to the star, y and z are the perpendicular coordi-nates on the sky. This line profile is taken to be a Voigt function (e.g. Rybicki & Lightman 1979) evaluated us-ing the isothermal temperature T , and bulk fluid velocityis included by transforming the photon frequency fromthe planet frame to the rest frame of the fluid.The transit depth will be expressed in terms of a fre-quency dependent planet radius, R p ( ν ), which is definedas the radius of an opaque disk that is required to pro-duce the same transit depth as the integral over themodel atmosphere: F (0) ν − F ν F (0) ν ≡ R ( ν ) R ⋆ = Z star dy dz h − e − τ ν ( y,z ) i , (25)where corrections due to limb darkening have been ig-nored for simplicity. The fractional decrease in flux inEquation 25 is independent of ISM absorption, and de-pends solely on the planetary atmosphere. The integra-tion over y and z extends over the stellar disk, where starhas radius R ⋆ .The frequency-integrated transit depth for the modelsis calculated as δFF = R dν I ( ⋆ ) ν (cid:16) R p ( ν ) R ⋆ (cid:17) e − τ (ISM) ν R dν I ( ⋆ ) ν e − τ (ISM) ν (26)where I ( ∗ ) ν = " (cid:12)(cid:12)(cid:12)(cid:12) ∆ v
67 km s − (cid:12)(cid:12)(cid:12)(cid:12) − . (27)is a fit to the shape of the Ly α intensity of the Sununder quiet solar conditions (Feldman et al. 1997). Inunits of velocity from line center at frequency ν , ∆ v = c ( ν − ν ) /ν . The limits of integration in Equation 26are −
200 km s − ≤ ∆ v ≤
200 km s − as in Ben-Jaffel(2008). The ISM optical depth τ (ISM) ν is computedusing the Voigt line profile evaluated with a temper-ature T ism = 8000 K and a neutral hydrogen column N H , ism = 10 . cm − (Wood et al. 2005). The Ly α lineis completely absorbed within ∆ v ≃ ±
50 km s − fromline center by the ISM.The HI number density n H ( x, y, z ) is computed by as-suming a balance between optically-thin photoionizationand radiative recombination (cf. Paper I, Section 8), J n H = α B n e n p , (28)where J ≈ (6 hr) − (0 .
047 AU /D ) is the ionizationrate for a Solar EUV spectrum (Paper I), and α B ( T ) ≃ . × − cm s − (10 K/ T ) . is the case B ra-diative recombination rate (Osterbrock & Ferland 2006).Assuming charge neutrality, n e = n p , and setting ρ = m p ( n p + n H ), Equation 28 has the analytic solution n H = " p J /α B + 4 ρ/m p − p J /α B (29)At a number density n eq = J /α B the gas at density ρ is 50% ionized with n H = n p . For n H & n eq , the gas ismostly neutral, and vice versa for n H . n eq . The use ofa constant J above simplifies the problem by requiringonly the local gas density ρ to evaluate n H .0 Fig. 6.—
Model 1 transit curve for comparison to the transitdepths for HD 209458b from Ben-Jaffel (2008).
Results for HD 209458b
Figure 6 compares the Ly α transit radius versus wave-length for the fiducial Model 1 to HST STIS data fromBen-Jaffel (2008). Points near line center are heavilycontaminated by ISM absorption and geocoronal emis-sion and are omitted. Model 1 was designed to agreewith the data through adjusting P ss and a (see Table3). The integrated transit depth, δF/F ≈
10% (see Ta-ble 3), is in good agreement with Ben-Jaffel (2008) andVidal-Madjar et al. (2008).Figure 7 shows R p ( ν ) versus wavelength for some of themodels from Table 2. The left (right) panel shows theeffect of changing B ( D ), holding all other parametersfixed. The model number for each line is given in thefigure caption. For clarity, Doppler shifts ∆ ν = v x ( ν − ν ) /c due to line of sight motion were ignored in σ ν inthe left panel, but are included in the right panel, toassess the role of the tidal force in accelerating the fluid.Bulk fluid motion is able to increase the cross sectionsignificantly at wavelengths on the steeply falling part ofthe Doppler core, roughly within ∆ v = ± B = 1 and 10 G models, but in the range B =10 −
100 G, the transit depth is observed to grow onthe wings of the line. Since bulk fluid velocity effectshave been omitted, the increase in transit depth mustbe due to an increase in hydrogen column over a largearea surrounding the planet. Relative to B = 10 G,there is an increase in δF/F of 50% for the B = 50 Gmodel and 250% for the B = 100 G model. This resultclearly shows that the planetary magnetic field can havean important effect on the transit depth.In the present paper, the base pressure and isothermaltemperature are parameters of the model, and the tran-sit depth is most sensitive to these two parameters. Therange of these parameters (Table 2) was based on thedetailed one-dimensional hydrostatic models, includingionization and heating/cooling balance, presented in Pa-per 1. More complete MHD simulations including heat- ing and cooling would determine these quantities self-consistently as part of the solution, and for a given stellarEUV heating rate they would no longer be parameters.In such more complete models, the magnetic field wouldstill be an essentially unconstrained parameter, as it isnot measured or constrained by any observation as yet.Figure 7 shows that, if B was the main uncertainty inthe model, an upper limit may be placed on magneticfield so that the transit depth is not too large comparedto observations. For the fiducial parameters adopted forHD 209458b, that upper limit would be B ≈
50 G. How-ever, we caution the reader that the large uncertainty inthermal structure due to uncertainty in stellar EUV andaccelerated particles fluxes likely limit the practical abil-ity to constrain the planetary magnetic field. Neverthe-less, for the parameters used in this paper, sufficientlystrong magnetic fields may in principle have a strong ef-fect on the transit depth.Next consider the effect of changing the rotation rateand tidal force, by changing D with all other parame-ters held fixed. Comparison of the Model 1 lines (solidblack line) in the left and right panels shows that Dopplershifts due to bulk velocity in the WZ are small for thefiducial model and the model in which the planet hasbeen moved outward. However, moving the planet in-ward by 25% to D = 0 .
035 AU has the effect of broaden-ing the wavelength range where R p ( ν ) is large (comparethe dashed orange and solid black lines). This is due tobulk fluid velocities Doppler shifting those wavelengthsto the Doppler core, where the cross section is large.Gas that has escaped from the planet may still bestrongly bound to the star, and may achieve high bulkvelocity due to the gravity of the star. In the presentcase where the tidal force has been axisymmetrized, theeffect is symmetric on either side of the line. In the 3Dcase, red-shifted absorption due to gas falling toward thestar may achieve even larger velocities. For the chosenbox-sizes r ∼ R p , the tidal force can accelerate fluid topoloidal velocities ≃ Ω p r ≃
60 km s − in the simulationbox. For larger box sizes, even higher velocities may beachieved. However, it is unclear from the present sim-ulations if bulk velocities &
100 km s − in the flow canaffect the line profile, since the steeply falling gas densitymay not be sufficiently large to give τ ν & τ ν ( y, z ) at 100 km/sfrom line center are shown in the y-z plane in Figures8 and 9. The contours are evenly spaced in log τ , andwhite dashed lines show the τ = 0 . , τ & τ = 0 . − × τ ν , and so may contributesignificantly if the increase in area can overcome the de-crease in optical depth. For this to occur, the densitymust not decrease too rapidly outward from the planet.Nearly the entire planetary upper atmosphere is opticallythick when observed near line center (∆ v = 0), but mov-ing away from line center the transit depth falls rapidly1 Fig. 7.— Ly α transit depth (Equation 25) versus wavelength in velocity units for selected models. The left panel compares modelswith different B (Model 1, B = 10 G, solid black line; Model 2, B = 1 G, dotted orange line; Model 3, B = 50 G, solid red line;Model 10, B = 100 G, solid blue line), while the right panel compares models with different D (Model 1, D = 0 .
047 AU, solid black line;Model 8, D = 0 .
035 AU, dashed red line; Model 9, D = 0 .
06 AU, dashed blue line). In each panel, only the one parameter is changed,with all others held fixed at Model 1 values. For clarity, bulk Doppler shifts are ignored in the left panel; they only affect the line profileat ∆ v .
50 km s − from line center. Bulk Doppler shifts are included in the right panel, as the tidal force may give rise to additionalacceleration. Vertical lines are placed at ∆ v = ±
100 km s − for reference. once the atmosphere becomes optically thin, which oc-curs at a different value of ∆ v for the range of modelsshown. Fig. 8.—
Contours of Ly α optical depth τ ( y, z ) at ∆ v =100 km/s from line center for the steady-state simulation resultsfrom Model 1. The inner (outer) dashed white line shows τ = 1(0.1). Figure 8 shows the steady-state result for Model 1( B = 10 G) at an illustrative frequency correspond-ing to 100 km/s from the line center. A large equatorialDZ with τ ∼ . − r = 4 R p , whilethe same contour only extends to r = 2 R p at the poles.The gas in the mid-latitude WZ has significantly smalleroptical depth compared to points in the neighboring po- Fig. 9.—
Same as Figure 8 for the Model 10. lar and equatorial DZ’s. The left-right asymmetry, mostapparent in the τ = 0 . B = 100 G). The τ = 0 . B ≃
10 Gseen in the left panel of Figure 7.Figures 8 and 9 clearly show the contribution to thetransit depth from the DZ and WZ at a single, illustra-tive photon frequency (100 km s − from line center). Itis of interest to know what contribution the DZ and WZmake at all other wavelengths, and which wavelengthscontribute most to the integrated transit depth in Equa-tion 26. A technical point is that in order to know ifa certain point is inside the DZ (WZ), one must tracealong the field line to determine if it is closed (open),and if the fluid velocity is everywhere small (or accel-erates to the sonic point). A simpler but approximateapproach, followed here, is to compare the transit depthdue to “slow” and “fast” material. The dividing line be-tween the two is set by a threshold v p , thr on the poloidalvelocity; slow material has v p ≤ v p , thr and vice versafor fast material. The slow material does not strictlytrace out the DZ, since the fluid velocity at the base ofthe WZ is also small. By examination of optical depthcontour plots using different velocity thresholds, we findthat v p , thr = 0 . a leads to only a small amount of slowmaterial at the base of the WZ. Given the optical depths τ slow , fast ( y, z ) for the slow and fast material, the inte-grand of Equation 26 can be computed. Note that whilethe optical depth is linear in the contribution from slowand fast material, the transit depth δF/F is not, since τ occurs in an exponent. Fig. 10.—
Integrand of the transit depth δF/F in Equation26, broken into the contribution from “slow” fluid with v p ≤ . a and “fast” fluid with v p > . a . Two magnetic field strengths areshown: Model 1 ( B = 10 G) and Model 10 ( B = 100 G). Figure 10 shows the integrand of Equation 26, sep-arated into slow and fast material, and computed fortwo different field strengths, Model 1 (10 G) and Model10 (100 G). First consider the B = 10 G lines. Thecontribution from the slow material shows the expectedpeak near 100 km s − , with small contribution at large∆ v due to small stellar flux, as well as near line center,due to ISM absorption. The fast material shows narrowpeaks just outside the region of ISM absorption. This isdue to poloidal fluid motions Doppler shifting the pho- tons from the Lorentzian wing back into the Dopplercore of the line, where the cross-section increases rapidlytoward line center. Even though the slow and fast ma-terials make comparable contribution to the frequency-integrated transit depth (i.e., the areas under the dashedand solid curves are comparable), the shapes of the rel-ative transit depth profiles are very different. This dif-ference can in principle provide a way to distinguish theabsorption due to slowly moving DZ material and fastmoving WZ material. Next, comparing the B = 10and 100 G lines shows the far larger transit depth forthe strong field case. This is due to the increased den-sity at large radii, the which is the result of the largerdead zone extending out to a region where the gravity − ∇ U rot becomes quite small, so that the scale height islarge and the density nearly constant. For the fast mate-rial, the higher field case shows absorption further fromline center due to the higher poloidal velocity as the fieldis increased (see Figure 2).Having investigated the transit spectra in detail, themass and angular momentum loss rates are now consid-ered. Comparing the values of ˙ M and ˙ J to δF/F inTable 3, the dominant effect is that higher P ss and a lead to both larger ˙ M and ˙ J , as well as δF/F . This isdue to the higher gas density. The role of magnetic fieldwell into the strongly magnetized regime is also clear, inthat larger B leads to larger δF/F and ˙ J , and slightlysmaller ˙ M .How well did the semi-analytic solutions for ρ and v inPaper I do at predicting the shape of the magnetospheresand the optical depth contours? Recall that in Paper I aninner dipole field was fitted to an outer monopole to rep-resent the transition between wind and dead zones. Thesize of the dead zone was found by stress balance at theequatorial cusp point. Figure 14 of Paper I shows the hy-drogen column for a set of models, and can be comparedto the optical depth contours in Figures 8 and 9. Model 1(6) of Paper I is similar to Model 1 (3) here. The 3-zonestructure is evident on both treatments, and many of thetrends, e.g. the growth of the equatorial DZ, are evidentin both. One key difference is that the equatorial DZlooks more “cuspy” in the simulations where the field isself-consistently calculated, while the semi-analytic cal-culations give round-shaped dead zones. The dead zonesizes in the two calculations are comparable. The shapesof the wind zones differ between the two calculationsdue to the different shapes of the field lines, since thepoloidal velocity is parallel to the poloidal field. Over-all, the semi-analytic methods of paper I give results inrough quantitative agreement with the simulations here,at least near the planet where the assumed field shapeof Paper I is approximately correct. Further from theplanet, the backward bending of field lines, and the de-viations from corotation become large, and are not takeninto account in Paper I. Results for HD 189733b
Lecavelier Des Etangs et al. (2010) used the HST ACSto measure an integrated transit depth of 5% in the Ly α line, smaller than the transit depth of HD 209458b bya factor ≃
2. In this section the exercise from Section4.2 is repeated for the planet HD 189733b, and the ef-fect of different M p , R p and D is discussed. The plan-3 Fig. 11.—
The fiducial model for HD 209458b in comparisonto three models for HD 189733b. The dot-dashed line is thefrequency-dependent transit depth for parameters matched to HD189733b, while the blue-dotted and dashed lines vary the base den-sity and sound speed, respectively. etary mass and radius are set to M p = 1 . M J and R p = 1 . R J , and the semi-major axis is D = 0 .
031 AU( http://exoplanet.eu/catalog ). The Jean’s parame-ter λ = GM p /R p a is a factor ≃ a .Figure 11 compares the transit depth for the fidu-cial Model 1 of HD 209458b to three models for HD189733b: a fiducial model with P ss , a and B identicalto that of Model 1 of HD 209458b; a model with pres-sure 4 times larger; and a model with sound speed 40%larger. The larger base pressure and sound speed are mo-tivated by the higher stellar EUV flux for HD 189733b(Sanz-Forcada et al. 2011), due to the closer orbital sep-aration, as well as higher stellar activity.The fiducial model for HD 189733b has smaller tran-sit depth ( δF/F = 0 .
04) than that of Model 1 for HD209458b ( δF/F = 0 . M p /R p forHD 189733b causing the density to decrease outwardfaster. The transit depth is surprisingly close to the ob-served value. Since the EUV flux is higher, both thebase pressure and temperature are expected to be higherthan the HD 209458b case. As an illustration, increas-ing the base pressure slightly has the effect of increasingthe density everywhere, leading to larger transit depth( δF/F = 0 .
05) comparable to the observed value forHD 189733b. Lastly, increasing the sound speed by 40%gives the model the same λ as Model 1 of HD 209458b.The higher transit depth ( δF/F = 0 . SUMMARY
The MHD wind model simulations presented in this pa-per demonstrate that, for an extended range of latitudes,the planet’s magnetic field can qualitatively change theproperties of a thermally-driven outflow. In Paper I,dipole field geometry and expected field strengths of hotJupiters were used to estimate the size of the closed fieldline regions, where an outflow is quenched by rigid mag- netic field lines. The inclusion of the stellar tide was alsoshown to quench the outflow in the polar region due tothe higher potential barrier. The estimated equatorialdead-zone sizes were ∼ − R p for the parameters ofinterest.MHD simulations permitted the relaxation of a pre-scribed field geometry, as hydrodynamic stresses andmagnetic stresses in the thermally-driven outflow fromthe hot inner boundary were computed self-consistentlythrough force balance both across and along field lines.The simulation results verify the main features of thesemi-analytic models of Paper I, including the depen-dence of their structure of the upper atmosphere on B (Figures 1 and 2) and stellar tide strength (Figures 3and 4). The MHD wind model also permitted a self-consistent calculation of the mass and angular momen-tum loss rates, which were presented in Section 3.3 fordifferent stellar tide strengths and field strengths, andthe integrated Ly α transit depth, δF/F (the observablequantity, see Section 4). The results are most consistentwith a pressure of 50-60 nbar and a temperature of ∼ K at the base of the thermosphere for HD 209458b, cor-responding to a mass loss rate ˙ M ≃ × g s − andan angular momentum loss rate ˙ J ≃ × g cm s − .A central result of this paper is that for sufficientlylarge magnetic field, the large resulting equatorial DZmay dominate the optically thick area which gives riseto the transit depth signal. In this strongly magnetizedregime, we find that the transit depth increases with themagnetic field. If the thermal structure were well known,this strong dependence on the magnetic field would allowan upper limit to be placed on the planet’s magnetic fieldso as not to over-predict the transit depth. However, dueto uncertainties in heating rates due to stellar EUV andaccelerated particles, such an exercise is likely not pos-sible. However, the parameter study in this paper doesmake it clear what magnetic field strength is requiredto increase the transit depth in the parametrized mod-els, which may be compared to the thermal structure inmore realistic models.A consequence of the existence of a large DZ in thestrong field case is that Ly α absorption occurring out tonear the Roche-lobe radius does not directly imply theabsorbing gas is escaping (e.g. Vidal-Madjar et al. 2003),as emphasized in Paper I. The MHD model (both ana-lytic and numerical) does exhibit gas in the mid-latitudeWZ which is escaping, however, the Ly α transmissionspectrum is less sensitive to the gas in this region as theoptical depths are lower (see Figure 9).Because of large observational uncertainties, the tran-sit depth as a function of wavelength cannot preciselyconstrain the pressure at the base of the warm H layer(see Figure 7). The integrated transit depths com-puted from the model Ly α spectra presented in Section 4provided another quantitative comparison with observa-tions. The high sensitivity of the integrated transit depthon the pressure at the base of the warm H layer suggeststhat this observable quantity can probe and constrain theconditions in the thermosphere of highly irradiated hotJupiters. At the same time, the numerical models pre-sented here provide complementary information aboutthe resulting expected mass and angular momentum lossrates, which are inaccessible by observations.4The authors thank Duncan Christie and David Sing forhelpful conversations. The authors also thank the refereefor a thoughtful and thorough report which improved this paper. This work was supported in part by NSF (AST-0908079) and NASA Origins (NNX10AH29G) grants. APPENDIX
SHEAR LAYER AND CURRENT SHEET NEAR THE EQUATORIAL DEAD-ZONE/WIND-ZONE BOUNDARY
At the boundary between the equatorial DZ and mid-latitude WZ, there are sudden changes in fluid velocity andmagnetic field over short distance. The origin of this shear layer and current sheet was discussed in Section 3.1. Asthe simulations presented in this paper do not include explicit viscous forces and Ohmic diffusion, it is the grid-scalenumerical effects contained with the ZEUS-MP code that control the behavior of the solutions at these discontinuities.The possible effects from numerical diffusion came to our attention due to the spurious bump in density at theequator shown in Figure 5, where a rise in density occurs inside the Hill radius. This behavior contradicts basicanalytic considerations. It was shown in the Appendix of Paper I that, in steady state, the Bernoulli constant W ≡ v + a ln ρ + U rot (A1)(where v and U rot are defined in a frame corotating with the planet) must be a constant along a given field line inthe dead-zone (see Equation A10 of Paper I). Since v = 0, our choice of base density, Equation 9, indicates that W should be a constant throughout the dead-zone. The force-balance equation A8 of Paper I implies immediately thatthe Lorentz force must vanish in the dead-zone. If this is the case, the density along the equator can only increasewith distance outside the Hill radius. This contradicts the numerical results in Figure 5.The spurious density bump seems to be due to numerical effects near the dead-zone/wind-zone boundary, althoughthe origin is rather subtle. In hindsight, it is not surprising that the boundary would be difficult to treat numerically,because there are discontinuities in all quantities: density and each of the three components of the velocity andmagnetic field. The discontinuity in magnetic field, in particular, is difficult to treat accurately (see, e.g., Fig. 16 ofStone & Norman 1992b). One way to increase the accuracy is to increase the spatial resolution, which we have donefor the fiducial Model 1. In Fig. 12, we show the equatorial density profile at four different resolutions (100 × × × × ×
800 it does not disappear completely.This density enhancement is caused by magnetic forces (see Fig. 13), which should vanish throughout the dead-zoneaccording to the analytic considerations mentioned earlier. At the heart of these considerations is the constancy ofthe Bernoulli constant W throughout the dead-zone. It breaks down in the numerical simulations, as illustrated inFigure 14, where we show the distribution of the Bernoulli constant and poloidal velocity for a selected region for thestandard resolution 272 × W is indeed very close to the expected value over most of the dead-zone(where the poloidal velocity is small, see the right panel), except in a layer near the DZ/WZ boundary, where deviationof order 1% is evident; this is also the region where the magnetic forces become appreciable, and the density startsto increase outward spuriously. As the resolution increases, the boundary layer shrinks in size. This should serveas a cautionary tale for future simulations of hot Jupiter magnetospheres, especially in 3D, where the resolution willnecessarily be coarser than in 2D. Nevertheless, the basic three-zone structure of the magnetosphere is robust.Although the density profile in the equatorial DZ is resolution dependent, the transit depth for the four modelsshown in Figure 12 varies only slightly, with values δF/F = 0 . , . , . , . REFERENCESBallester, G. E., Sing, D. K., & Herbert, F. 2007, Nature, 445, 511Ben-Jaffel, L. 2007, ApJ, 671, L61Ben-Jaffel, L. 2008, ApJ, 688, 1352Ben-Jaffel, L., & Sona Hosseini, S. 2010, ApJ, 709, 1284Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883Braginskii, S. I. 1965, Reviews of Plasma Physics , 1, 205Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L.2002, ApJ, 568, 377Christensen, U. R., Holzwarth, V., & Reiners, A. 2009, Nature,457, 167Dalgarno, A., & McCray, R. A. 1972, ARA&A, 10, 375Ehrenreich, D., & D´esert, J.-M. 2011, A&A, 529, A136Ehrenreich, D., Lecavelier Des Etangs, A., H´ebrard, G., D´esert,J.-M., Vidal-Madjar, A., McConnell, J. C., Parkinson, C. D.,Ballester, G. E., & Ferlet, R. 2008, A&A, 483, 933Ekenb¨ack, A., Holmstr¨om, M., Wurz, P., et al. 2010, ApJ, 709,670Feldman, U., Behring, W. E., Curdt, W., Schuehle, U., Wilhelm,K., Lemaire, P., & Moran, T. M. 1997, ApJS, 113, 195Fossati, L., et al. 2010, ApJ, 714, L222 France, K., Linsky, J. L., Yang, H., Stocke, J. T., & Froning,C. S. 2011, Ap&SS, 335, 25Garc´ıa Mu˜noz, A. 2007, Planet. Space Sci., 55, 1426Gold, T. 1959, J. Geophys. Res., 64, 1219Grießmeier, J.-M., Stadelmann, A., Penz, T., et al. 2004, A&A,425, 753Gu, P.-G., Lin, D. N. C., & Bodenheimer, P. H. 2003, ApJ, 588,509Hapke, B. 1993, Topics in Remote Sensing, Cambridge, UK:Cambridge University Press, —c1993,Hayes, J. C., Norman, M. L., Fiedler, R. A., Bordner, J. O., Li,P. S., Clark, S. E., ud-Doula, A., & Mac Low, M.-M. 2006,ApJS, 165, 188Heinemann, M., & Olbert, S. 1978, J. Geophys. Res., 83, 2457Henry, G. W., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2000,ApJ, 529, L41Heyvaerts, J., & Norman, C. 1989, ApJ, 347, 1055Holmstr¨om, M., Ekenb¨ack, A., Selsis, F., et al. 2008, Nature, 451,970Ip, W.-H., Kopp, A., & Hu, J.-H. 2004, ApJ, 602, L53 Fig. 12.—
The equatorial density profile for Model 1 at three different grid resolutions in steady-state, as well as an inset that shows amore detailed view of a region from 2 . < r/R p < Fig. 13.—