Tidally distorted barytropes and their Roche limits, with application to WASP-12b

The hot Jupiter WASP-12b has been found to be on a decaying 1.09-day orbit. The mean density of the planet inferred from transit and radial-velocity data is near its Roche limit; just how near depends on the planet's uncertain internal structure. There is also spectroscopic evidence of mass loss. We accurately calculate the Roche density on the assumption of a synchronously rotating $n=1$ polytrope, and find this to be only $15-20\%$ below the observational estimates for the mean density. We discuss the implied shape of the planet, its lifetime before complete disruption, and its current rate of mass loss based on our improved tidal model. The present mass-loss rate is at least as sensitive to the opacity and temperature profiles of the planet's atmosphere as to its internal structure, however.


INTRODUCTION
This work is motivated by the hot Jupiter WASP-12b, which uniquely among hot Jupiters has been found to be on a decaying orbit (Yee et al. 2020), and which shows spectroscopic indications of possible mass loss (Fossati et al. 2010). As discussed below, the photometric and radial-velocity data suggest that the planet is close to contact with its Roche lobe, but this depends upon the planet's radial mass profile: two limits in which Roche contact can be estimated analytically give contrasting answers.
The mean density of the host star can be estimated directly from the planet's orbital period (P ) and transit light curve (Seager & Mallén-Ornelas 2003). The light curve also provides the ratio R p /R * of planetary and stellar radii. Given the mass ratio q ≡ M p /M * from radial-velocity data, one also has the mean density of the planet,ρ p = 3M p /4πR 3 p . For a given planetary structure, there is a critical value of the dimensionless group Gρ p P 2 below which Roche overflow occurs, presuming that the planet's orbit is circular and its rotation tidally locked. Collins et al. (2017) find ρ p = 0.266 ± 0.015 g cm −3 for WASP-12b. This is on the low side of earlier values from the literature that they tabulate in their Table 4. A weighted average of all of these estimates based on their quoted errors is ρ p = 0.283 ± 0.012 g cm −3 . We emphasize that through-out this paper,ρ p is defined in terms of the transit radius. Because of the planet's tidal distortion, its volume must be somewhat larger than 4πR 3 p /3, and its true mean density averaged over that volume somewhat less thanρ p .
If one estimates the Roche lobe in the usual way as if all of the mass of the planet were concentrated at its center, then the projected area of the lobe at mid-transit would be approximately 1.296R 2 H , where R H ≡ (q/3) 1/3 a is the Hill radius and a the orbital semimajor axis. This would correspond to a critical transit radius R p,crit ≈ 0.642R H . Given the observed orbital period P = 1.09142 d (Yee et al. 2020), this model for the planet's structure predicts a minimum mean densitȳ ρ p ≈ 2.71G −1 (2π/P ) 2 ≈ 0.180 g cm −3 : significantly below the observed value. Thus on this model, the planet would be safely within its Roche lobe.
If instead the planet had uniform density, then the limiting equipotential at Roche contact would be ellipsoidal with semiaxes b 1 : b 2 : b 3 = 1 : 0.511 : 0.483, and b 1 b 2 b 3 = 0.06757qa 3 (e.g. Chandrasekhar 1987). The corresponding mean density based on the transit radius √ b 2 b 3 would then be 0.473 g cm −3 : well above what is inferred from the data. On this model, the planet would already have been tidally disrupted. The minimum densityρ p,crit for a realistic planetary structure lies somewhere between these limiting cases.
To a first approximation, the radial density profile of a jovian planet can be described as an n = 1 Emden polytrope (Stevenson 2020). One would have thought that the tidal limit of such polytropes would have been worked out long ago. Chandrasekhar (1933) has calculated weak tidal distortions for such bodies by perturbative methods, while Rasio & Shapiro (1992, 1994 have obtained numerical results via smoothparticle hydrodynamics for mass ratios q ≡ M 2 /M 1 near unity, as they were motivated by neutron-star mergers. Taniguchi & Gourgoulhon (2002) have made calculations for n = 1 (adiabatic index γ = 2) and mass ratios as small as q = 0.1. Numerical difficulties prevented them from making models in full Roche contact, but by double extrapolation of their tabulated results in both q and tidal strength, we estimate axis ratios b 1 : b 2 : b 3 ≈ 1 : 0.62 : 0.59 for such a case when q = 10 −3 . At WASP-12b's orbital period, this would implyρ p,crit ≈ 0.24 g cm −3 , 10-15% below the observational estimates quoted above; in other words, the planet would not yet fill its Roche lobe. In view of this slender margin and the importance of the WASP-12 system, a more careful calculation would seem to be called for. That is the goal of the present paper.
Section 2 describes our implementation of the Ostriker-Mark (OM) method for tidally distorted polytropes and some tests against analytic results. Section 3 and Table 1 present Roche-limited models for Emden indices n ∈ {1, 3/2, 3} models and mass ratios q = M p /M * ∈ {10 −4 , 10 −3 , 10 −2 }. Section 4 and the Appendix discuss mass-loss rates for models that best fit the observedρ p ; these do not quite fill the Roche lobe but we endow them with realistic finite-temperature atmospheres, which we borrow from recent work in spherical symmetry. We also estimate from the observed orbital-decay rate (Yee et al. 2020) and our best-fitting polytropic models how long it will be before the main body of the planet (i.e., its roughly isentropic core) comes into Roche contact. Section 5 sums up and points to directions for future work.

METHODOLOGY
Our general problem is determining the shape of a small fluid body (jovian planet) in circular orbit around a larger body (the star). WASP-12b's orbit is indeed consistent with circular (Yee et al. 2020). Presumably this is enforced by tidal dissipation within the planet, in which case the planet's (unobserved) spin is expected to be synchronous apart from thermally or convectively driven zonal flows. Even with these simplifications to the orbital dynamics, a very accurate calculation of the equilibrium shape would require accurate knowledge not only of the equation of state, but also of the planet's internal stratification of entropy and composition.
These things being uncertain even for Jupiter, let alone WASP-12b, we represent the planet as a polytrope, i.e. pressure and density are locally related by p(ρ) = Kρ (n+1)/n for some constant K. For the most part, we take n = 1. This is believed to be a rather good first approximation to Jupiter's equation of state (Stevenson 2020). Furthermore, the Juno mission has measured Jupiter's quadrupolar Love number k 22 = 0.565 ± 0.006 (Durante et al. 2020); for a hypothetical non-rotating planet with Jupiter's radial density profile, this would be reduced by about 9%, i.e. to k 2 ≈ 0.51. An n = 1 polytrope has k 2 = 15/π 2 ≈ 0.520, whereas k 2 = 3/2 for n = 0, k 2 ≈ 0.285 for n = 3/2 and of course k 2 ≈ 0 for a strongly centrally concentrated body. While Jupiter's internal structure is clearly more complicated than a polytrope, Love numbers are a direct measure of tidal response, albeit for infinitesimal tidal fields, and depend only on the radial density profile. On the other hand, WASP-12b is substantially inflated, with a mean density less than a quarter of Jupiter's, so that its internal structure (d ln P/d ln ρ) may be rather different from Jupiter's. Therefore we have considered other values of the polytropic index, especially when estimating the planet's mass-loss rate (Appendix A).
We use the Ostriker-Mark method for finding the self-consistent density and potential for the planet in the presence of the tidal field exerted by the star (Ostriker & Mark 1968). Some minor additions to the method are needed to cope with spherical harmonics of odd degree, which were not needed in the original application to rotating stars. For completeness, we briefly review the method here, with emphasis on our additions to it.
Since we assume the planet to rotate synchronously at angular velocity Ω, hydrostatic equilibrium in the corotating frame takes the form ∇p/ρ + ∇Ψ = 0 (1) in terms of an effective potential in which is the gravitational potential of the planet, Φ * is that of the star, represented here by a point mass M * , and r 0 is the barycenter of the orbit. It follows from eq. (1) that the pressure is constant on equipotentials Ψ = constant; furthermore from the curl of the equation, ∇ρ ×∇ p = 0, so that the density is also constant on equipotentials. Hence p and ρ may be regarded as functions of Ψ and therefore of one another. One may introduce a formal enthalpy function in terms of which hydrostatic equilibrium can be reexpressed as Ψ s being the value of Ψ at the surface of the planet where its pressure vanishes. Note that eqs. (4)-(5) are consequences of hydrostatic equilibrium regardless of the equation of state, which in general involves more than one independent thermodynamic variable (e.g. density, entropy, and chemical composition). To use the OM method, however, one needs an explicit formula or tabulation for p(ρ). A polytropic relation p = Kρ (n+1)/n leads to It's convenient to choose the arbitrary constant in eq. (3) so that Ψ s = 0 and so Ψ < 0 within the planet. For the reasons given earlier, the choice n = 1 appears to be most relevant to jovian planets, but we also test our code against analytic results for n = 0, and we compute models for n = 3/2 and n = 3 as well as n = 1. We put the origin of r at the center of mass of the planet. The distance from there to the barycenter is r 0 = a/(1 + q), with a the separation between the centers of the planet and the star, and q = M p /M * the mass ratio. We ignore the oblateness and tidal distortion of the star (WASP-12 appears to be slowly rotating) and expand its potential in zonal harmonics about the planet: Notice that the term at l = 0, which exerts no force on the planet, is omitted from the sum, so that Φ * = 0 at r = 0. The polar angle θ is measured from the line between the planet and the star. Following the conventions of the restricted three-body problem, we make this the x axis, except that x = y = z = 0 at the center of the planet rather than the barycenter, and the star lies at (x, y, z) = (a, 0, 0); y increases in the direction of the planet's orbital motion, and z in the direction of Ω.
With these conventions, the centrifugal part of the effective potential (2) becomes upon dropping a constant term − 1 2 Ω 2 r 2 0 that exerts no force. Together with eq. (7), this ensures that Ψ(0) = Φ p (0), i.e. the effective potential at the center of the planet is due to solely to its self-gravity. The next term, Ω 2 r 0 x, is almost balanced by the l = 1 part of the stellar potential (7), but not exactly, because of the finite size and nonspherical shape of the planet: Let us define the l th mass moment of the planet by so that µ 0 = M p ; also µ 1 = 0 because we expand around the planet's center of mass. The (l + 1) st term in the expansion (7) for the stellar potential Φ * exerts (through ∂Ψ * /∂x) a net attractive force on µ l . Thus the balance of centrifugal and gravitational forces on the planet is The left side is due to the linear part of the centrifugal potential (8); the quadratic part exerts no net force because µ 1 = 0. Since r 0 = a/(1 + q), eq. (10) can be rewritten as with Ω K being the "Keplerian" orbital frequency (mean motion) that would obtain if the planet were a point mass or a sphere. The notation δ for the correction follows LRS93. For a planet of mean radiusR p ≪ a, δ scales as (R p /a) 2 , which becomes O(q 2/3 ) at Roche contact. Since we are interested in mass ratios q ≈ 10 −3 , the correction might be ∼ 1%; in fact it is an order of magnitude smaller because of the central concentration of the (n ≈ 1) planet. Our OM iterations are carried out in Emden units: 4πG = (n + 1)K = ρ c = 1. So for example, eq. (6) becomes simply ρ(Ψ) = [max(−Ψ, 0)] n . To emphasize that a physical quantity is expressed in Emden units, we sometimes mark it with a tilde. In particular, is our dimensionless measure of the strength of the tide, corresponding to (1 + q)μ/4 in the notation of LRS93. Our OM procedure is standard overall: 1. At the start of the (k + 1) st iteration, we have an estimate Ψ (k) for the effective potential. Using this, we estimate the planet's density ρ (k+1) on a uniform grid in each coordinate (r, θ, φ) using eq. (6) with Ψ → Ψ (k) .
2. At each radial grid point r i , we transform the density samples on the (θ, φ) grid to coefficients of spherical harmonics Y l,m (θ, φ) up to some maximum degree l max using the publicly available Python package Pyshtools (Wieczorek & Meschede 2018). 1 In this way we get a discrete approximation to the series representation 3. The corresponding radial functions Φ l,m (r) in the expansion of the potential Φ p (r, θ, φ) are obtained from the usual radial integrals over ρ l,m (r), except that the boundary condition on the We approximate these radial integrals via spline interpolation and Simpson's Rule. We now have preliminary updates to the planet's potential Φ p and mass moments µ l .
4. In general the center of mass of the planet is no longer exactly at the origin. We alter the radial functions Φ l,0 (0) in the expansion of Φ p slightly to correct this, as described below.
5. Using the new estimate µ (k+1) 0 for the planet's mass M p , we re-estimate the separation of the planet and star by solvingΩ 2 K = M p (1+1/q)/4πa 3 for a (recall G → 1/4π in Emden units). Then we update Ω 2 and δ using eq. (11). 6. We add terms to the coefficients to represent the stellar potential (7) and the centrifugal potential (8). The latter contributes at l ∈ {0, 1, 2}; the contribution to the radial function for l = 0 is quadratic in r rather than constant, because y 2 + z 2 averages to (2/3)r 2 over a sphere. With these tidal and centrifugal additions, we have a preliminary representation for the next iterate Ψ (k+1) of the effective potential in the form of a series of spherical harmonics with coefficients depending on radius.
We elaborate here on step 4. For q ≪ 1, the tidalplus-centrifugal potential ∆Ψ ≡ Ψ−Φ p is predominantly quadrupolar (l = 2), apart from the centrifugal contribution at l = 0 noted above. If we ignore all odd values of l in ∆Ψ, then step 4 in not needed. However, the oddl components of ∆Ψ contribute to the l = 1 parts of the planet's density distribution. This is true even for n = 1 where the relation between ρ and Ψ is linear within the planet, but nevertheless is globally nonlinear because of the truncation of the density at Ψ = 0. Physically, the distortion of the planet is stronger on the side facing the star than on its opposite side. Thus, each iteration generally perturbs the center of mass. 2 We cannot simply ignore the (l, m) = (1, 0) part of Φ p , because it should be nonzero within the planet. Thus at step 4, we find the displacement ∆x = µ 1 /µ 0 of the center of mass, and correct for it via the first-order approximation Since ∂/∂x = cos θ∂/∂r + r −1 ∂/∂ cos θ, the corrections to the radial functions are (to first order in ∆x) in which the coefficients are if l ≥ |m|, else 0. (15) Without step 4, the method usually does not converge when the odd parts of the tide are included.
Hachisu's Self-Consistent Field method (Hachisu & Eriguchi 1984), while otherwise similar to Ostriker-Mark, has the advantage in this respect that the endpoints of the bodies on the x axis are directly constrained; however, Hachisu's method does not directly constrain the mass ratio (unless q = 1), which has to be discovered through iteration.

Tests of the method
We have satisfied ourselves that the output parameters of our models presented in Table 1 [columns (3)-(10)] converge as we refine the radial grid and increase the maximum spherical-harmonic degree (l max ). In fact, values shown there are obtained by Aitken acceleration with respect to both ∆r ∈ {0.005, 0.0025, 0.00125} (in Emden units) and l max ∈ {8, 12, 16}.
We have also tested our code against analytic results. Our code accurately reproduces the properties of spherical Emden polytropes. In particular, for n = 0 and n = 1, the calculated masses and radii agree with the exact values to at least 5 places (in fact 7 places except for the mass of the n = 0 model); this is for a radial grid spacing ∆r = 0.0025 with rational grid points that do not coincide with the exact radii (ξ max,n=0 = √ 6 and ξ max,n=1 = π): ξ max is estimated by linear interpolation between neighboring grid points.
To test the accuracy of tidal distortions, we calculate the Love numbers is the corresponding component of the distorted planet's potential, and r max is the radius of the spherical planet. In particular, the quadrupolar Love number k 2 = 15/π 2 − 1 ≈ 0.5198 for n = 1, whereas our code yields 0.5208 with ∆r = 0.00125 and a tidal strength Ω 2 K = 10 −8 . This is an error of 0.2%, and it is linear in ∆r: Aitken extrapolation based on the results for ∆r ∈ {0.005, 0.0025, 0.00125} yields the correct value to 5 places. The Love numbers for n = 0 are also known analytically but are much more difficult for our code. During the OM iterations, we evaluate the density from the potential on a fixed grid in (r, θ, φ), making no attempt to interpolate between grid points. For n = 0, the density depends only on the sign of the potential (vanishing for Ψ > 0); therefore, if we apply a very small tidal potential that should move the bounding equipotential a distance < ∆r, there is no change to the density on the grid at all. When we apply a finite tidal distortionΩ 2 K = 10 −3 ≈ 0.044Ω 2 K,Roche , we obtain the estimate k 2 = 1.547, which is 3% larger than the exact value k 2 = 3/2 for an infinitesimal tide.
To test the ability of the code to measure finite distortions, we search for the Roche limit when the Emden index n = 0 and the mass ratio q → 0 + . In this case, the exact result is known via elliptic integrals (e.g. Chandrasekhar 1987):Ω 2 K,Roche ≈ 0.02252325. To estimate this with our code, we set q = 10 −18 , l max = 12, ∆r = 0.00125, and gradually increaseΩ 2 K until the code no longer converges. The last converging model isΩ 2 K ≈ 0.02266, or 0.6% more than the exact value. The size of the error is disappointing, but as explained above, n = 0 is more stressful for our code than n ≥ 1.

RESULTS
We define the critical model for a positive polytropic index n and mass ratio q to be that for which the surface gravity vanishes at the pole facing the star: i.e., ∂Ψ/∂x = 0 at the point where surface Ψ = 0 inter-sects the positive x axis ("L1"). 3 With this definition, the critical model is at the brink of mass transfer onto its companion even in the approximation that the density and temperature vanish at its surface, and therefore can be considered to be in contact with its Roche lobe. (In Appendix A, following Jackson et al. (2017) and earlier authors, we estimate the mass-loss rate for slightly weaker tides when a finite atmospheric temperature is allowed for.) LRS93 have shown within their ellipsoidal approximation that secular or dynamical instability typically sets in before the Roche limit, but the differences in the tidal parameter (12) are very small for q ≪ 1. Table 1 displays the calculated dimensionless parameters of the critical model for three values of n and q that seem likely to be relevant for gaseous exoplanets. The parameters are estimated by Aitken extrapolation in ∆r and ℓ max and should be accurate to the number of significant digits shown.
We also calculate models that match the mean densities estimated from observations of WASP 12b. For n ≥ 1, such models lie within their Roche lobes: the surface gravity at the inner Lagrange point is positive, though small compared to that of a spherical planet of the same mass and mean density. For n = 1, the dimensionless tidal parameters (12) areΩ 2 K = 0.00760 and 0.00691 for mean densities of 0.266 and 0.283 g cm −3 , respectively.

DISCUSSION
By modeling WASP-12b as a polytrope, we pretend that the density and temperature vanish on its limiting equipotential. This is not physical; the surface of the planet must have a finite temperature, which will be rather high on the side facing the star: 3000 − 4000 K. The atmosphere therefore has no sharp edge, but attenuates roughly exponentially with scale height H = c 2 s /g, g being the surface gravity and c s = (k b T /µm p ) 1/2 the thermal speed. In fact the atmosphere can be expected to extend to the Lagrange points, 4 where it will pass through sonic points and escape the planet. The massloss rate from WASP-12b has been estimated on this basis by several groups, notably Li et al. (2010), Lai et al. (2010), Bisikalo et al. (2013), Jackson et al. (2017), and Dwivedi et al. (2019). Our methods are detailed in Appendix A. We use the formula of Jackson et al. (2017)   but improve upon it by replacing their assumption of a constant thermal speed with one based on recent models for the day-side atmospheric temperature profile by Lothringer et al. (2018) and Arcangeli et al. (2021). Our result for the mass-loss timescale |M p /Ṁ p | ≈ 5 Myr is within a factor of 2 of that of Li et al. (2010), although this seems to be partly coincidence: their estimate for the density at the L1 point is more than two orders of magnitude less than ours, but they assume that the mass loss occurs through the entire Roche lobe, while we estimate the width of the "nozzle" around L1 following Jackson et al. (2017). At the same time, ouṙ M p is 50 times larger than the latter authors'. A factor 10 is due to the differences in potential (our n = 1 model vs. their point-mass approximation), the rest being due to differences in the temperature structure and to the assumed opacity at the transit radius. These things remain uncertain. We also conclude that the mass-loss rate on the opposite side of the planet through the L2 point should be negligible, in contrast to the conclusions of Dwivedi et al. (2019), whose total mass-loss rate is ∼ 10 −4 times ours because they allow for heating of the gas only by extreme ultraviolet radiation that carries only ∼ 10 −6 of the star's bolometric luminosity.

SUMMARY AND CONCLUSION
Using a version of the Ostriker-Mark method, we have computed polytropic models of "planets" in Roche contact. We are motivated by WASP 12b, because it is very close to this condition, is on a decaying orbit (so that if it is not yet in Roche contact, as we find, then it soon will be), and because it shows spectroscopic evidence for mass loss. These methods, however, could be applied to other very-short-period hot Jupiters that have well-constrained densities, such as some of those in Table 2 of Jackson et al. (2017).
Because real jovian planets have atmospheres that are not on the same adiabat as their interiors, the notion of "Roche contact" is ambiguous. A more precise statement is that the equipotential and isobaric surfaces (these are not the same) of WASP 12b appear to be detached from the corresponding surfaces that pass through the inner saddle point of the effective potential (L1). Our improved tidal models have allowed us to make more precise estimates of just how far WASP 12b is detached from its Roche lobe in this sense. We have also used these models to revisit the thermally-driven atmospheric mass-loss rate of the atmosphere, which we find to be comparable to the decay time of the orbit as measured by Yee et al. (2020)-a few Myr. However, the planet will achieve Roche contact in only a few hundred thousand years. We have not attempted to analyze what may happen after that: the planet may disrupt, or the orbit may re-expand, as in some cataclysmic variables (Warner 2003).
While an improvement over previous work as regards the tidal distortion and gravitational potential of the planet, our models remain highly idealized. We have worked with polytropes rather than realistic equations of state and profiles of entropy and composition. These remain quite uncertain, especially for highly "inflated" (low-density) cases such as WASP 12b. The atmospheric structures that we have borrowed from previous work to estimate the mass-loss rate were computed without accounting for the tide; ideally one would treat the atmosphere and the tide self-consistently. But it is not clear that the circulation of the atmosphere is well enough understood at present to justify the effort.
Here ρ s is the mass density at the "surface" of the planet where Ψ ≡ 0, Ψ L1 is the value of the effective potential at the inner Lagrange point, and c s = k b T /µm amu is the thermal speed for temperature T and molecular weight µ; m amu is the atomic mass unit. This formula is based on the following approximations: 1. The gas temperature is presumed constant above the surface.
2. The atmosphere is approximately hydrostatic below the sonic point, which lies near L1.
The first two assumptions give the exponential in eq. (A1). The last factor is the estimated cross section of the "nozzle" through which the flow passes, the second partial derivatives being evaluated at L1, where ∇Ψ = 0 and ∂ 2 Ψ/∂x 2 < 0. Clearly, because of the exponential it is essential to have accurate models both for the effective potential and for the thermal state of the atmosphere. Jackson et al. (2017) evaluate the former in the usual restricted three-body approximation where the gravitational potential of the planet is that of a point mass. We improve on this with our tidally-distorted polytropic models, and find that it makes a significant difference. The indicated second partial derivatives in eq. (A1)are quite similar to those derived by Jackson et al. (2017) from the restricted 3-body problem, however. For the temperature, Jackson et al. (2017) use T eq ≡ (F * /2σ sb ) 1/4 , where F * ≡ L * /4πa 2 is the insolation at the substellar point, as if the insolation were re-radiated uniformly over the "day" side without redistribution to the other hemisphere. With the parameters of Collins et al. (2017), F * ≈ 10 10 erg cm −2 s −1 , so that T eq ≈ 3000 K. We improve on this using recent models for the temperature structure and circulation of the atmosphere by Lothringer et al. (2018) and Arcangeli et al. (2021). 5 Since these atmospheres are not perfectly isothermal vertically, we replace the exponential factor in eq. (A1) (except for the e −1/2 , which represents the transition from hydrostatic to sonic flow) with eqs. (4)-(5), recast as with the understanding that the integral is to be taken vertically, perpendicular to the constant-pressure surfaces: since the entropy varies along the equipotentials, it is not possible for the atmosphere to be in exact equilibrium, but one assumes that hydrostatic equilibrium is better approximated vertically than horizontally because the atmospheric scale height H ≡ c 2 s /g ≪ R. The function T (P ) is taken from Lothringer et al. (2018)'s fiducial model (their Figs. 2 & 14), which is hotter than T eq -closer to 4000 K-above the 10 −4 bar level, presumably because of the wavelength dependence of the opacity. The molecular weight µ(P, T ) varies with height because of dissociation of molecular hydrogen. For this we assume LTE and the partition function of Irwin (1987) for H 2 . The gas is predominantly molecular at the transit radius, but quickly becomes atomic with decreasing pressure. We take (X, Y, Z) = (0.71, 0.27, 0.02) for the abundances of hydrogen, helium, and metals, withμ Z ≈ 16.
Following Jackson et al. (2017) and Howe & Burrows (2012), we identify the pressure level at the transit radius R T from the condition that the optical depth along the transit chord at the planet's limb should be 0.56, i.e. τ T ≈ κ λ ρ √ 2πHR T = 0.56, where κ λ , ρ, and H ≡ c 2 s /g are evaluated at the tangent point. We take R T = 1.9 R J and M p = 1.47 M J from Collins et al. (2017). By scaling our n = 1 model to these values, we find that the surface gravity at the tangent point is g = 0.913GM p /R 2 T = 922 cm/s 2 (whereas the gravity at the pole facing the star is 410 cm/s 2 ). For the temperature at this point, we use T = 2500 K ≈ (F * /4σ sb ) 1/4 , which is roughly consistent with the plots in Arcangeli et al. (2021). With these choices, hydrogen turns out to be overwhelmingly molecular at the tangent point, so µ ≈ 2.36. For the wavelength-dependent opacity κ λ , we consider only the contribution of H − , calculating its volume mixing ratio from the formula provided by Lothringer et al. (2018), and its optical cross section from Ohmura & Ohmura (1960), evaluating this at λ = 0.62 µm, the effective wavelength of the SDSS r ′ band (Collins et al. 2017). We find κ 0.62 µm (H − ) = 0.002 cm 2 /g and P ≈ 9 mbar at the tangent point. This is quite important: from Lothringer et al. (2018), H − appears to dominate the opacity at this wavelength, but a larger opacity would predict a smaller pressure and density at the tangent point, and ultimately a smaller mass-loss rate, as discussed below.
The circulation models of Arcangeli et al. (2021) indicate that the temperature is almost constant with longitude at P 1 bar. Therefore, by setting P 2 = 1 bar in eq. (A2), we find the value of P 1 on the day side that corresponds to (i.e., lies on the same equipotential as) the pressure at the transit point by using the day-side temperature profile instead of the constant 2500 K that we assumed for the terminator; we also take into account the variation in molecular weight as described above. This gives P 1,day = 20 mbar; presumably this would be slightly higher at the substellar point, where the temperatures will be higher than the average for the day side. This equipotential is considered the "surface" of the planet in our tidal model, Ψ(P 1,day ) = 0. The potential at the Lagrange pointΨ L1 = 0.0988 in Emden units, in which also the planet massM p = 42.1195, transit radiusR T = 2.97396, and squared orbital frequency (1 + δ)Ω 2 = 0.0076075. By scaling the latter three to the observed values from Collins et al. (2017), we find that the predicted potential difference between the planet surface and the Lagrange point is Ψ L1 ≈ (10.69 km/s) 2 . Using eq. (A2) yet again, and including the factor e −1/2 from eq. (A1), we have P L1 ≈ 0.15 mbar, ρ L1 ≈ 5.7 × 10 −10 g cm −3 , and c s,L1 ≈ 5.16 km/s (the hydrogen is fully atomic there). In this way we obtain, finally,Ṁ ≈ 1.7 × 10 16 g s −1 , corresponding to a mass-loss timescale M p /Ṁ ≈ 5 Myr. A fraction (Ṁ Ψ L1 )(L * R 2 p /a 2 ) ≈ 0.003 of the insolation is needed to drive this outflow.
To test the sensitivity ofṀ to our assumptions about the internal structure of the planet, we have repeated the procedure described above for an n = 3 rather than n = 1 polytrope (scaled to the same mean density, etc.). An n = 3 model is of course much more centrally concentrated, and its potential is much better described in the restricted-threebody approximation. In this case, we find that Ψ L1 ≈ (13.28 km s −1 ) 2 , about 50% larger than for the n = 1 model. As a result, the density at L1 is about one order of magnitude smaller than before, and so is the predicted mass-loss rate:Ṁ n=3 ≈ 2 × 10 15 g s −1 .
These two models (i.e. n = 1 and n = 3) probably bracket the truth as far as the dependence ofṀ on the structure of the planet's interior is concerned. However,Ṁ is also sensitive to the structure of the atmosphere. In particular, the pressure at the transit point (P T ) is inversely proportional to the assumed opacity (including scattering as well as absorption). Sing et al. (2013) argue from transmission spectroscopy for the presence of Al 2 O 3 aerosols, and therefore P T between 0.02 and 0.5 mbar depending upon particle size. Also, Sing et al. (2016) classify WASP-12b as "hazy" based on the difference between optical and infrared transit radii. If we put P T = 0.02 mbar rather than 9 mbar as we found from H − opacity alone, the day-side pressure on the same equipotential becomes 0.73 mbar rather than 20 mbar, andṀ → 1.2 × 10 15 g s −1 (with our n = 1 model) rather than 1.7 × 10 16 g s −1 . Notice thatṀ varies sublinearly with P T because of the larger pressure scale height on the day side. By measuring the secondary eclipse of WASP-12 in the optical with HST, however, Bell et al. (2017) constrain the day-side geometric albedo of the planet to be less than 0.064 at 97.5% confidence, indicating very little scattering, and they reject the aerosol models of Sing et al. (2013).
We have used the same methods to estimate the mass-loss rate of the n = 1 model through the L2 point, i.e. on the side of the planet opposite the star. This is much smaller thanṀ L1 . The difference in potential between the Lagrange points and the "surface" of the planet is 33% larger for L2 than for L1. If the night side had the same temperature structure as the day side, this would lowerṀ L2 by a factor of about 4. But the temperature structure is not the same, of course. While we could not find a detailed model of the night-side temperature profile, Fig. 12 of Arcangeli et al. (2021) offers color contour plots of the estimated temperature at five discrete pressure levels over the whole surface of the planet; the temperature at the anti-stellar point appears to decrease from ≈ 2500 K at 1 bar to ≈ 1100 K at 1 mbar. If T = 1100 K at all P ≤ 1 mbar, thenṀ L2 ≈ 10 2 g s −1 ∼ 10 −14Ṁ L1 . This enormous contrast is due to the exponential sensitivity of eq. (A1) to the thermal speed. (In addition to the direct effect of the temperature on the thermal speed, there is an indirect effect via the molecular weight: at lower temperatures, the hydrogen remains molecular to much lower pressures, assuming LTE.) At sufficiently low pressures, however, the cooling time of the gas will exceed the time that it takes to be advected from the day side, so we expect that the temperature should rise well above 1100 K at such low pressures. Based on our own rough estimates (involving, e.g., cooling by rovibrational lines of CO), we are confident that the mass-loss rate through L2 should be at least two orders of magnitude smaller than through L1 insofar as eqs. (A1)-(A2) apply.