© ESO, 2011

1. Introduction

Because the measurement of the radii of close-in transiting planets continues to gain in accuracy, providing stringent constraints on exoplanet theoretical models, any source of errors in the radius determination must be determined with precision. Current ground and space-based photometric observations of the host stars of transiting planets enable us to address new problems. The first direct detection with Spitzer of the light emitted by the planet (Deming et al. 2007) opened a new path for probing the physical properties of the surface and the atmosphere of transiting exoplanets. Among the first results of the Kepler mission, the detection of ellipsoidal variations of the host star induced by tidal interaction with a low mass companion has been claimed Welsh et al. (2010). More recently, Carter & Winn (2010a,b) have shown that light curve analysis can put direct constraints on the actual shape of transiting planets. They also investigated the impact of the precession of an oblate object with a non zero obliquity around the orbital axis on the shape and timing of the transit signal.

These observations motivate us to investigate the deformation of the planet with respect to a spherical body and caused by tidal or rotational forces. While previous studies have focused on the detectability of the oblateness of a flattened body, we address the more general problem in the present paper, namely the determination of the general shape of a planet (or star) distorted by both a tidal and a centrifugal potential, and its impact on the transit depth, hence on the determination of its correct radius. In order to compute the ellipsoidal shape (flattening and triaxiality) of a gaseous body, we derive in Sect. 2 a simple analytical model of the internal structure of the object based on the polytropic equation of state (Lai et al. 1994). The polytropic indexes are calibrated in Sect. 4 by comparing with numerical models describing the structure of strongly irradiated gaseous planets (Leconte et al. 2009). In Sect. 3, we present directly usable analytical expressions giving the shape (oblateness and triaxiality) of a distorted planet (or star) as a function of its mass and polytropic index and compare our estimates with the numerical method outlined in Appendix A and with the measured values for the major planets of our solar system. As a by-product, our model yields analytical expression for the first gravitational moment J₂ and the love number k₂ of a self gravitating fluid body. Finally, Sect. 5 quantifies the effect of the non-sphericity of the planet on the transit depth.

We find that as the planet transits across the stellar disc, we only see the smallest cross section of its actual ellipsoidal shape so that the depth of the transit is decreased with respect to the expected signal for a spherical object, as discussed by Li et al. (2010) in the case of WASP-12 b. This implies that the radius inferred from the light curve analysis, derived under the assumption of spherical planet and star, underestimates the real equilibrium radius of the object. This bias needs to be corrected for a proper comparison of the structure and evolution of extrasolar planets with theoretical 1D numerical simulations and enhances the actual discrepancy between theory and observation for the so-called “bloated” planets.

2. Variational method for compressible ellipsoids

In this section, we briefly describe the energy variational method developed by Lai et al. (1993) and Lai et al. (1994) (hereafter LRS1 and LRS2) to construct general Darwin-Riemann equilibrium models. In Sect. 2.1 we briefly summarize the basic assumptions and the equilibrium relations are derived in Sect. 2.2. More details about the method in general, as well as the applications to compact objects, can be found in LRS1 and LRS2, and references to equations in these papers are denoted with numbers preceded by “I” and “II”, respectively, in the present paper. Solutions to first order in the deformation are derived for tidal and rotational deformations in Sects. 3.2 and 3.1, respectively.

2.1. Model description

Consider an isolated, self-gravitating fluid system in steady state. The system is characterized by conserved global quantities such as its total mass M and total angular momentum J. The basic idea in our method is to model our self gravitating system by a limited number of parameters x₁,   x₂,... and in such a way that the total energy can be written as (1)An equilibrium configuration is then determined by extremizing the energy according to (2)An expression like Eq. (1) can be written down for the total energy of a binary system (with components of mass M and M′). We adopt a polytropic equation of state between the pressure P and the mass density ρ, (3)This defines the polytropic index n and the “entropy” K – both are constant within the object and sufficient (with M) to describe the mechanical structure of a given object. Under the combined effects of centrifugal and tidal forces, the objects (stars or planets) achieve nonspherical shapes. We model these shapes as triaxial ellipsoids of principal axes (a₁,a₂,a₃) and , respectively. Throughout this paper, unprimed quantities refer to the component of mass M while primed quantities refer to the component of mass M′. The three directions along which our principal axes are measured are, respectively, the line connecting the center of mass of the two components, its normal contained in the orbital plane, and the direction of the orbital angular momentum vector. In the simple case of coplanar and synchronous rotation, a₃ is simply the polar radius and a₁ and a₂ are the equatorial radii of the component measured toward its companion and in the orthogonal direction, respectively.

Specifically, we assume that the surfaces of constant density within each object can be modeled as self-similar ellipsoids. The geometry is then completely specified by the three principle axes of the outer surface. Furthermore, we assume that the density profile ρ(m) inside each component, where m is the mass interior to an isodensity surface, is identical to that of a spherical polytrope with the same volume. The velocity field, v, of the fluid is modeled as either uniform rotation (corresponding to the case of a synchronized binary system) or uniform vorticity, ∇ × v, (for nonsynchronized systems). The vorticity vector is assumed to be parallel everywhere to the orbital rotation axis.

For an isolated rotating gaseous sphere, these assumptions are satisfied exactly when the fluid is incompressible (polytropic index n = 0), in which case the true equilibrium configuration is a homogeneous ellipsoid (Chandrasekhar 1969). For a binary system, our assumptions are strictly valid in the incompressible limit only if we truncate the tidal interaction at the quadrupole order. We adopt this quadrupole-order truncation of the interaction potential in this paper.

After adding up the orbital separation, r, our set of unknowns is or equivalently where ρ_c is the central density and λ₁ ≡ (a₃/a₁)^2/3,λ₂ ≡ (a₃/a₂)^2/3 (when no ambiguity exists or otherwise stated, primed quantities are defined in the same manner as unprimed ones by simply making the transformation x ⇄ x′ throughout the equations). The total energy can be written (4)where (5)is the internal energy of component 1 (cf. Eq. (I.3.1)), and (6)is the self-gravitational energy of component 1 (cf. Eq. (I.4.6)) with, and The dimensionless radius and density, ξ and θ (with a subscript 1 when taken at the surface), are the classical variables used to describe polytropic gaseous spheres (see Chandrasekhar 1939) and G is the gravitational constant. The kinetic energy in the inertial frame reads (11)where the spin kinetic energy of body 1 (T_s) is given by (cf. Eq. (I.5.6)) (12)with Ω being the rotational orbital velocity, Λ is a measure of the internal rotation rate in the co-rotating frame, (13)is a dimensionless coefficient measuring the inertia of the body, and (14)is the moment of inertia with respect to the rotation axis. Tables giving values of the polytropic constants , , κ_n, and ξ₁ as a function of n can be found in LRS1 and in Chandrasekhar (1939). For non synchronous rotation (Λ ≠ 0), our gaseous body is not in the state of solid-body rotation. A rotation rate can thus not have the usual meaning. To have a sense of the angular velocity, one can take the half of the vorticity (, where v is the fluid velocity vector in the inertial frame) as a proxy¹. ω is related to Λ by (15)The orbital kinetic energy T₀ is simply (16)Finally, the gravitational interaction energy W_i reads (17)with (18)

2.2. Equilibrium relations

We can now derive the set of equilibrium relations yielding seven algebraic equations for our seven unknowns . The details of the transformations necessary to express the total energy as a function of the unknowns and conserved quantities alone and to be able to carry out the differentiation are explained in Sect. 2.2.1 of LRS2 and just add technical details not needed here. We thus give only the results.

Differentiation with respect to r simply yields the modified Kepler’s law for the orbital mean motion Ω (19)with (20)Differentiation with respect to the central density ρ_c yields the virial relation, (21)with the mean radius R = (a₁a₂a₃)^1/3 and (22)Using expressions for U and W, we get the equilibrium mean radius (23)where R₀ is the radius of the unperturbed spherical polytrope given by (Chandrasekhar 1939) (24)Finally, the differentiation with respect to λ₁ and λ₂ yields after some algebra (cf. Eqs. (I.8.4)–(I.8.6)): where

μ_R ≡ GM′/r³, , and is the mean density of the ellipsoid.

3. The shape of gaseous bodies

In this section, we derive analytical expressions for the deformation – induced either by centrifugal or tidal potential – of a gaseous body. To test the validity of our assumptions, we compared our predictions with the measured values for the major planets of our solar system in Sect. 3.1. Since we do not make any assumption about the masses of the two components, these equations can be used indifferently to compute the shape of the star or of the planet by choosing M = M_p and M′ = M _⋆  when considering the planet, and vice versa when focusing on the star.

In general, the set of equations described in the previous section must be solved numerically, but we study here the first order development of these equations at large orbital separation. This approximation corresponds to neglecting terms of order , which is consistent with our truncation of the gravitational potential at the quadrupole order and is appropriate to address close-in transiting planetary systems. In practice, this is done by setting Δ = Δ′ = 0 and (27)in Eqs. (19), (23), (25), (26) and their primed equivalent and by expanding these equations to first order in α_i (). First we derive some general formulae by expanding the integrand in the definition (10): which yields To first order, (28)and (29)The principal moment of inertia of the body can also be computed and reads (30)The other moments of inertia can be computed by replacing 1 and 2 by the appropriate indices. The dimensionless moment of inertia κ_n for different planetary masses, age and stellar irradiation can be found in Tables B.2 and B.3.

3.1. Rotational deformation: Maclaurin spheroids

Our set of equations also allows us to compute the effect of the centrifugal force alone on a slowly rotationg fluid object. To do so, one just has to take the M′ → 0 limit in Eqs. (19), (23), (25) and (26). Therefore, Ω is a free parameter (the rotation rate of our body) and there is a degeneracy between Ω and Λ that allows us to choose Λ = 0.

We introduce the dimensionless angular velocity (31)as a small parameter of order in all expansions. The volume expansion factor can be calculated using (32)We get (33)To the same order of approximation, the two remaining equations are given by Eqs. (25) and (26) which yield (34)Combinations of Eqs. (33) and (34) give the three figure functions (cf. Eq. (A12) of LRS2) (35)For this configuration, the usual variables are the oblateness (36)and the dimensionless quadrupole moment of the gravitational field J₂ given by the theory of figures to first order (Zharkov et al. 1973) (37)For practical purposes, R₀ can be computed to first order by using Eq. (27) and reads (38)in our geometry, where R_eq and R_pol denote the usual equatorial and polar radii. For an incompressible body (n = 0) we retrieve the usual solution of the theory of planetary figures (Zharkov & Trubitsyn 1980).

Attempts have been made to constrain the oblateness and thus the rotation period of transiting planets by using the solar system planets as test cases (Carter & Winn 2010a,b). Because of the wide variety of exoplanets, it is important to have the ability to predict the flattening of fluid planets for a wider range of parameters than encountered in the solar system. Figure 1 shows the predicted oblateness for various planet masses as a function of the rotational period .

Fig. 1 Oblateness given by Eq. (36) as a function of the rotation period (in days) at 1 Gyr for planets of mass: 0.3 MJ (dotted), 0.5 MJ (dashed), 1 MJ (dash-dotted), 3 MJ (long dashed), 15 MJ(solid). The oblateness decreases when the mass of the planet increases because massive objects are more compressible (see Sect. 4), have a more intense self-gravity field and are thus less subject to perturbations.

3.2. Tidal deformation: determination of the Love number (k₂)

To compute the shape induced by the tidal force alone, we consider a non-rotating configuration (ω = 0). From Eq. (15), this is achieved if . Then Thus from Eq. (23) we see that there is no change of volume to lowest order,

Since Ω² = μ_R(1 + p) with p = M/M′ and , only the zeroth order must be taken in the left hand side of Eqs. (25) and (26), which yields (with help of Eq. (29)) (39)Thus (cf. Eq. (A25) of LRS2) (40)As long as the hydrostatic equilibrium holds, this equation can be used to compute the shape of the planet and its host star at each point of the orbit. We recover the usual dependence of the tidal deformation in , with a factor of order unity, q_n, which encompasses all the structural properties of the gaseous configuration.

Since we are in the linear approximation with a gravitational potential restricted to quadrupolar order, the shape of our body can be described with the usual Love number of second order, k₂ (Love 1909, which is twice the apsidal motion constant often called k₂ in the stellar binary literature). Indeed, once k₂ (and h₂ = 1 + k₂ for a body in hydrostatic equilibrium) is known, the external potential and the shape that a body will exhibit in response to any perturbing potential can be computed as detailed in Appendix A. To derive k₂, we compute the quadrupolar term of the gravitational potential energy of the system formed by our compressible ellipsoid and a point mass, by introducing Eq. (40) in the linearized version of Eq. (17), and identify this term to the potential energy due to tides given by (Darwin 1908) (41)This yields (42)As expected, in the n = 0 limit, we retrieve the Love number of an incompressible ideal fluid planet k₂ = 3/2. We can also see that k₂ is linked to the square of the dimensionless moment of inertia κ_n. This is because level surfaces are self-similar in our model and that the love number encompass both the deformation of the body ( ∝ κ_n) and the gravitational potential created by the deformation ( ∝ κ_n). The Love number for different planetary masses, age and stellar irradiation can be found in Tables B.2 and B.3.

We can see that the value of the Love number tends to decrease with mass above 1 M_J. This is due to the fact that more massive objects are more compressible and thus more centrally condensed (see Sect. 4). At constant mass, enrichment in heavy elements toward the center (possibly in a core) acts to decrease the value of k₂. In general, redistributing mass from the external to the internal layers, which are less sensitive to the disturbing potential, decreases the response of the body to an exciting potential, which translates into a lower k₂.

Our model predicts k₂ values in the range 0.3 − 0.6. As discussed by Ragozzine & Wolf (2009) such values of the Love number could be inferred by the measurement of the precession rate of very Hot Jupiters on eccentric orbits. Such measurements could be carried out by Kepler for WASP-12 b analogs with an eccentricity  > 3 × 10^-4 (most favorable case) or Tres-3 b analogs with an eccentricity  > 2 × 10^-3 (for k₂ ≈ 0.3) and lower eccentricities for higher Love number values. Such measurements would indeed be extremely valuable as they would put direct constraints on the central enrichment in heavy elements inside close Hot Jupiters, like the measurements of the gravitational moments of the solar system planets.

3.3. Synchronized planets

For values of the tidal dissipation factors inferred for Jupiter (Goldreich & Soter 1966; Leconte et al. 2010), the timescale of pseudo synchronization of close-in giant planets is less than about a million years. The planet is thus in a state of pseudo synchronization, with a rotation rate given by (Hut 1981; Leconte et al. 2010) (43)in the weak friction theory, with e being the eccentricity of the orbit. For the simple case of a circular orbit, the spin is thus synchronized and, either solving Eqs. (19), (23), (25) and (26) in the synchronized case (Λ = 0) or simply adding the results of Eqs. (40) and (35) (there is no cross correlation terms to first order) yields (44)and (45)

3.4. Model validation

There are two major assumptions in the present calculations:

• the absence of a central core. The aim of such an approximation isto avoid to introduce any free parameter in the model. In any case,the core mass and the global enrichment in giant extrasolarplanets are yet weakly constrained (Guillot 2005; Leconteet al. 2009). We show that thisapproximation introduces an uncertainty of  ≈10% on the derived shape;

• the polytropic assumption. This allows us to derive a completely analytical model. Comparison with a more detailed numerical integration (see Appendix A) shows that the deviation between the results of the two models (analytical vs. numerical) is smaller than the uncertainty due to the no-core approximation.

Since the oblateness (f), J₂, k₂, the mean radius and rotation rate are known for the major planets of our solar system, we can test our theory on these objects. The details of the calculation of the chosen polytropic indexes are presented in Sect. 4. The results are summarized in Table 1, which shows the actual values of the relevant parameters for the two major planets taken from Guillot (2005), the values of the oblateness and of J₂ calculated with our model and, for comparison, with the assumption of an incompressible body (n = 0). We see that, whereas the values of f derived from the incompessible model differ from the true values by almost a factor of 2, our polytropic model predicts the f-values to within 12%. Note that higher-order terms (of order ) are not completely negligible for rapidly rotating bodies such as Jupiter and Saturn (see Zharkov et al. 1973; Chabrier et al. 1992). The polytropic model, however, yields J₂-values that differ from the measured values by 30% (for Jupiter) and 59% (for Saturn). These discrepancies are mostly due to the large metal enrichment in these planet interiors, probably with the presence of a large dense core as detailed three layers models can reproduce exactly the measured moments (Chabrier et al. 1992). Note that this no-core approximation has less relative impact on the distortion of the shape predicted by the model than on the gravitational moments because these effects scale as h₂ = 1 + k₂ (see Appendix A) and k₂ respectively, k₂ being  ≲ 0.6 in the situations of interest. Such discontinuities in the density profile (and its derivatives) could be addressed more precisely with two different polytropes, but this would add extra free parameters and would not serve the actual purpose of the present paper.

While we decided to use a polytropic assumption to infer a fully analytical theory, the figures of a body in hydrostatic equilibrium can be derived without this assumption. As shown in Sterne (1939) (see also Zharkov et al. 1973; Zharkov & Trubitsyn 1980; Chabrier et al. 1992, for more detailed applications to the giant planets case) and outlined in Appendix A, this theory, however, requires a numerical integration even to first order. For an ideal n = 1 polytropic sphere, Eq. (42) agrees with the numerical results of Sterne (1939) with less that 1% error. To compare these methods in our context, we derive the values of k₂ using our analytical model (Eq. (42)), and by numerical integration of Eqs. (A.6) and (A.11) for our best representative models of Jupiter and Saturn in our grid (although without cores). For Jupiter, our Eq. (42) gives k₂ = 0.55, the numerical integration gives k₂ = 0.57 to be compared with the measured value of k_2,J = 0.49. For Saturn, these values are 0.60, 0.66 and k_2,S = 0.32 respectively. Both models predict k₂ values than are higher than the measured ones. A direct consequence of the presence of heavy elements inside our giant planets. Comparing the models, our Eq. (42) yields slightly smaller k₂ values than the numerical integration which tends to mimic a central over-density (see Sect. 3.2). As discussed above a more precise modeling requires the addition of central enrichment in heavy elements whose mass fraction would be a free parameter. Without better knowledge of the internal composition of giant exoplanets, we think that the two methods yield similar results up to the sought level of accuracy. For sake of completeness, the values of k₂ computed with both methods are presented in Tables B.2 and B.3.

4. Polytropic index in gaseous irradiated planets

To readily use the results of Sect. 3, one only needs to have a proper value for the polytropic index n to be used. In this section, we derive realistic polytropic indices from numerical models of gaseous irradiated planets. All the other polytropic functions (κ_n, q_n, ...) can be derived by integrating the Lane-Emden equation and are tabulated in Chandrasekhar (1939) and LRS1. They are given for different planetary masses, age and stellar irradiation in Tables B.2 and B.3. We focus on the polytropic index in the planet because, in the context of transiting exoplanets, both the stellar rotation and the stellar tides have a negligible impact on the transit depth, as will be discussed in Sect. 5. The main physics inputs (equations of state, internal composition, irradiated atmosphere models, boundary conditions) used in the present calculations have been described in detail in previous papers devoted to the evolution of extrasolar giant planets (Baraffe et al. 2003; Chabrier et al. 2004; Leconte et al. 2009), and will not be repeated here.

We computed a grid of evolution models of gaseous giant planets with solar composition for various masses M_p ∈  [0.35 M_J,   20 M_J]  and incoming stellar flux F _⋆  ∈  [0,   4.18 × 10⁶   W   m^-2] . Low irradiation model can be used to infer the oblateness of long period rotating planets (such as Jupiter), while strongly irradiated models can be used to infer the shape and its impact on the transit of close-in planets. For the non irradiated case, the grid extends to 75 M_J. As the effect of irradiation on the internal structure decreases with the effective temperature of the object, these models computed with non irradiated boundary conditions should give a fair description of massive brown dwarfs in the range of irradiation considered. The pressure-density profile of each model is then fitted by a polytropic equation of state (Eq. (3)) at each time step and an example of the result of such a fit is shown in Fig. 2. Note that the disagreement between the actual P − ρ profile and the polytrope in the lower left area of Fig. 2 is both expected and needed: this low-density region (the first 5% in mass below the atmospheric boundary surface) has a different effective polytropic index than the planetary interior. In order to capture the bulk mechnical property of the planet, we weight each shell in the internal structure profile by its mass during the fitting procedure.

Fig. 2 The internal pressure-density profile of an irradiated 1.8 MJ planet (solid line). The dashed line represents the best-fit polytropic equation of state. The pressure-density range covered in the inner part of the body (95% in mass) is represented by the thicker part of the solid curve, which is well modeled by a polytropic EOS. As the thin part of the P − ρ curve represents only 5% in mass of the body it is disregarded by the fit.

This provides us with a grid tabulating the polytropic index of the planet, n_p ≡ n_p(M_p,F _⋆ ,t), and its spherical equilibrium radius, R_0,p ≡ R_0,p(M_p,F _⋆ ,t), where t is the age of the object. These functions, along with other quantities (T  _eff, ...), are tabulated in Tables B.2 and B.3². Figures 3 and 4 show the variation of n with the mass for different ages with F _⋆  = 0 and 4.18 ×    10⁶   W   m^-2, respectively.

Fig. 3 Polytropic index for non-irradiated planets as a function of the planet’s mass Mp at 100 Myr (dotted), 1 Gyr (dashed) and 5 Gyr (solid). The shaded area represents the uncertainty on the polytropic index for the 5  Gyr case (see text).

Fig. 4 Polytropic index for strongly-irradiated planets as a function of Mp at 100 Myr (dotted), 1 Gyr (dashed) and 5 Gyr (solid). As the irradiated atmosphere impedes the radiative cooling of the objects, it retards its contraction. Therefore, the non-monotonous behavior observed at the early ages in the non-irradiated case (Fig. 3) is enhanced, even at a later epoch. The bump at the high mass end of the 100 Myr curve is caused by deuterium burning (see text). The shaded area represents the uncertainty on the polytropic index for the 5  Gyr case (see text).

As shown in Fig. 3, in the non-irradiated case, we recover qualitatively the results of Chabrier et al. (2009): except for the early stages of the evolution, the (dimensionless) isothermal compressibility of the hydrogen/helium mixture is a monotonically increasing function of the polytropic index, , and thus of the mass of the object. In the high mass regime, n slowly increases as the relative importance of ionic Coulomb effects compared with the degenerate electron pressure decreases, and approaches the n = 3/2 limit, the expected value for a fully degenerate electron gas, when M_p approaches the hydrogen burning minimum mass (≈70 M_J) as can be seen in Table B.2. In the low mass regime, the compressibility decreases with the mass because the repulsive Coulomb potential between the ions, and thus the ionic electrostatic energy becomes dominant. Ultimately, electrostatic effects dominate, leading eventually to for solid, terrestrial planets.

A new feature highlighted by the present calculations is the non-monotonic behavior occurring between 1−3 M_J at early ages. This occurs when the central regions of the planet, of pressure P_c and temperature T_c, previously in the atomic/molecular regime, become pressure-ionized, above 1−3 Mbar and 5000−10 000 K (Saumon et al. 1995; Chabrier et al. 1992; Saumon et al. 1992), and the electrons become degenerate. An effect more consequential for the lowest mass objects, whose interiors encompass a larger molecular region. This stems from the fact that (Chandrasekhar 1939) (46)Older (with smaller R_p) and more massive (M_p ≳ 2 M_J) objects have P_c > 10   P_ionization and the ionization extends all the way up to the outermost layers of the gaseous envelope, which then contains a small enough mass fraction of molecular hydrogen to significantly affect the value of the polytropic index. This contrasts with younger objects around 1−3 M_J, whose external molecular hydrogen envelope contains a significant fraction of the planet’s mass, leading to a larger value of the polytropic index, as molecular hydrogen is more compressible than ionized hydrogen (see e.g. Fig. 21 of Saumon et al. 1995). Once again, for these latter objects, the interior structure would be better described by using two different polytropes, but such a significant complication of the calculations is not needed at the presently sought level of accuracy.

As seen in Fig. 4, a strong irradiation enhances the aforementioned feature: the evolution is delayed because the irradiated atmosphere impedes the release of the internal gravothermal energy. This yields a slower contraction, thus a lower central pressure (and lower central temperature) for a longer period so that the object enters the ionization regime at a later epoch. The bump at the high mass end of the 100 Myr isochrone is due to deuterium burning which also occurs later for a given mass, because of the cooler central temperature (see above). At 100 Myr, the 20 M_J has already burned a significant amount of its deuterium content and starts contracting again, whereas lower mass planets are still burning some deuterium supply, leading to a less compact and thus less ionized structure. This leads to the non-monotonic behavior on the high-mass part of the n − M diagram at 100 Myr, which reflects a similar behavior in the mass-radius relationship.

To evaluate the uncertainty in our determination of the polytropic index, we use an alternative derivation of n. As shown by Chandrasekhar (1939), the knowledge of M, K and n is sufficient to infer the radius of the polytrope, with the help of Eq. (24) and the central density, using (47)Since our numerical simulations provide both the radius, R_0,p(M_p,F _⋆ ,t), and the central density of the object, ρ  _c,p(M_p,F _⋆ ,t), we can invert Eqs. (24) and (47) to compute K  _p and n_p. This new determination of the polytropic index is compared with the previous one, obtained by fitting the P − ρ profile, in Figs. 3 and 4 for the 5 Gyr case: the new n_p value corresponds to the upper envelope of the shaded area. Figure 3 shows that the two approaches yield very similar results in the non-irradiated case. For the irradiated case, the average uncertainty on our determination of n_p lies between about 5% and 15% for the low mass planets.

5. Implications for transit measurements

When limb darkening is ignored, the depth of a transit is given by the ratio of the planetary and stellar projected areas. When both bodies are spherical, this simply reduces to δL _⋆ /L _⋆  ∝ (R_p/R _⋆ )². For close-in planet-star systems, however, both tidal and rotational deformations yield a departure from sphericity, so that what is measured is no longer the mean radius but an effective “transit radius” defined such that the cross section of the planet is equal to and similarly for the star. Thus the transit depth δ reads (48)

5.1. Impact on transit depth

In general, the projected area of an ellipsoid can be computed for any orientation and then at each point of the orbit, as explained in Appendix B. Figure 5 shows the projected area of the planet () as a function of its anomaly (φ) normalized to the spherical case (). When the planet is seen from its “side” (φ/π = 0.5), the observer sees a bigger planet because the rotation of the latter on itself tends to increase its volume, as has been mentioned by Li et al. (2010) for WASP-12 b. The possibility to measure these effects from the light curve is discussed in Ragozzine & Wolf (2009) and Carter & Winn (2010a).

Fig. 5 Normalized projected area of the planet as a function of its anomaly (φ) for inclinations of the orbit going from i = 90° (lowest curve) to i = 65° (highest curve) by steps of 5° for a WASP-12 b analog on a circular orbit. Top: for the full orbit. Bottom: zoom on the (primary or secondary) transit. The ordinates of the dotted, solid and dashed horizontal lines are respectively (face-on orbit), and .

For the simple case of an edge-on orbit at mid transit (φ = 0), since the observer, the planet and the star are aligned with the long axis of the tidally deformed ellipsoid^3,4, and . Therefore, (49)where R_0,p and R_0, ⋆  are the respective radii the planet and the star would have in spherical equilibrium and ηis by definition the variation of the transit depth induced by the ellipsoidal shape of the components relative to the transit depth in the spherical case. To first order in the deformation, this is given by (50)The choice of the expression to be taken for the α_i depends on the physical context. In the general case, one can use a linear combination of Eqs. (35) and (40) and get a general expression which depends on r, ω_p and ω _⋆ . However, most of the planet hosting stars have a low rotation rate compared to the orbital mean motion. This entails that the rotational deformation is negligible compared to the tidal one and can generally be neglected. As mentioned in Sect. 3.3, hot Jupiters should be pseudo synchronized early in their evolution. Therefore, we assume such an approximation in our calculations in order not to introduce any other free parameter. The impact of the rotation alone is described Sect. 3.1. Under such an approximation, (51)where the parameter p now denotes the mass ratio M_p/M _⋆ , and q_p and q _⋆  are equal to q_n for n = n_p and n = n _⋆ , respectively. The first line in the above equation represents the contribution of the planet, which is always negative (for reasonable values of n). Our line of sight follows the long axis of the tidal bulge and we see the minimal cross section of the ellipsoid.

The contribution of the star is positive and, in most cases, negligible compared the planet’s contribution because

for a typical system (10^-3 for a Jupiter-Sun like system). As a consequence, the results presented hereafter do not depend on q _⋆  as long as realistic values of n _⋆  ∈  [1.5,   3]  are taken.

Figure 6 portrays the relative transit depth variation computed with Eq. (51) for several planet masses as a function of the orbital distance, for a Sun-like parent star. While all the curves are calculated at an age of 1 Gyr, they do not change much for older ages because both the radii and the polytropic indices remain nearly unchanged after 1 Gyr (see Fig. 4). Given the accuracy of the radius determination achieved by the latest observations (1 to 10%), the transit depth variation is significant for Saturn mass objects (M_p ≈ M_J/3) closer than 0.04 AU and Jupiter mass objects closer than 0.02−0.03 AU. Because we derived the equations to first order, the value of η derived from our model should be taken with caution when η ≳ 0.1−0.3 (and are clearly not meaningful for η ≳ 1). In this regime, corresponding to the upper left region of Fig. 6, one should use the theory of planetary figures to higher order, but then numerical calculations become necessary, loosing the advantage of our simple analytical expressions. Figure 6 also displays the transit depth variation computed for the most distorted known transiting exoplanets, with the observationally measured parameters. The error bars reflect the uncertainties in the model and in the measured data.

Fig. 6 Relative transit depth variation η computed with Eq. (51) as a function of the semimajor axis at 1 Gyr for planets of mass: 0.3 MJ (dotted), 0.5 MJ (dashed), 1 MJ (dash-dotted), 3 MJ (long dashed), 15 MJ (solid). The shaded area shows the zone where higher order terms become non-negligible. The decrease in the transit depth due to tidal interactions is smaller when the mass of the planet increases because massive objects are more compressible (see Sect. 4) and thus less subject to nonspherical deformations.

5.2. Which radius?

Before going further, it is important to summarize the differences between the various radii that we have defined above. In the literature, the term “radius” is used loosely, even for nonspherical objects. Note that this can lead to discrepant normalizations throughout different studies and published values of transit radius measurements when, for example, radii are shown in units of Jupiter radii (R_J) without precisely defining the latter.

One can define a₁, a₂ and a₃ as the distances between the center and some isobar surface along the three principal axes of inertia. For any distorted object, we can define the mean radius (R) as the radius of the sphere that would enclose the same volume as the described surface. In our case of a general ellipsoid, we have R = (a₁a₂a₃)^1/3. If axial symmetry holds (e.g. for a rotating fluid body), we have a₁ = a₂ ≡ R_eq, defining the equatorial radius, and a₃ ≡ R_pol the polar radius. Finally, R₀ is the radius of the spherical shape that the fluid body would assume if it was isolated and at rest in an inertial frame (the limiting case for which all the mentioned radii would be equal). This latter is the radius computed in usual 1D numerical evolution calculations. Note that in general R ≠ R₀ because the centrifugal force has a net outward component that increases the volume of the object, as can be seen from Eq. (33).

One must be aware that only a₁, a₂ and a₃ (reducing to R_eq and R_pol for solar system gaseous bodies) can be measured directly and are not model dependent. This is why we define R_J as the equatorial radius of Jupiter at the 1 bar level (R_J ≡ R_eq,J = 7.1492 × 10⁷m, Guillot 2005, and reference therein).

Unfortunately, transit measurements only give access to the projected opaque cross section of the planet ( ≡ ) defining a “transit radius” which depends on the shape of the planet, its orientation during the observation and the wavelength used. To convert this transit radius inferred from the observations (R_tr,p) to the spherical radius (R_0,p) – that can be compared to 1D numerical models – one must eliminate δ from Eqs. (48) and (49). As shown above, the stellar impact on η is negligible compared to the planet’s contribution (R_tr, ⋆  ≈ R_0, ⋆ ). Then, using the first term in Eq. (51) and expanding the expression giving the definition of R_tr,p, one gets (52)For the most distorted known planets, the relative variation between the transit radius and the equilibrium radius

is positive and amounts to 3.00% for WASP-12 b, 2.72% for WASP-19 b, 1.21% for WASP-4 b, 1.20% for CoRot-1 b, 0.89% and OGLE-TR-56 b.

Of course, since (53)Eq. (52) is an implicit equation on R_0,p. To obtain R_0,p to the sought accuracy, a perturbative development in powers of η_tr = η(R_0,p = R_tr,p) can be obtained using recursively Eq. (52) (54)However, terms of order are of the same order than the second order corrections to the shape that we have neglected throughout.

6. Conclusion

Because of the large variety of exoplanetary systems presently discovered, with many more expected in the near future, and the increasing accuracy of the observations, it is important to take into account the corrections arising from the nonspherical deformation of the planet or the star, due to rotational and/or tidal forces, as such a deformation yields a decrease in the transit depth. In order to do so, it is extremely useful to be able to compute analytically the shape of planets and stars in any configuration from the knowledge of only their mass, orbital separation and one single parameter describing their internal structure, namely the polytropic index, n. Such formulae are derived in Sect. 3, and can be easily used to determine the impact of the shape of the planet on its phase curve and on the shape of the transit light curve itself (Carter & Winn 2010a). They can also be used to model ellipsoidal variations of the stellar flux that are now detected in the CoRoT and Kepler light curves (Welsh et al. 2010). These formulae also give good approximations for various parameters describing the mass redistribution in the body’s interior and the response to a perturbing gravitational field, i.e. the moment of inertia, I, and the Love number of second degree k₂.

Another major implication of the present work is to show that departure from sphericity of the transiting planets produces a bias in the determination of the radius. For the closest planets detected so far (≲0.05 AU), the effect on the transit depth is of the order of 1 to 10% (see Fig. 4), by no means a negligible effect. The equilibrium radius of these strongly distorted objects can thus be larger than the measured radius, inferred from the area of the (smaller) cross section presented to the observer by the planet during the transit. The analytical formulae derived in the present paper, and the characteristic polytropic index values derived for various gaseous planet masses and ages, make possible to easily take such a correction into account. Interestingly, since this equilibrium radius is the one computed with the 1D structure models available in the literature, the bias reported here still enhances the magnitude of the puzzling radius anomaly (see Fig. 6 of Leconte et al. 2010) exhibited by the so-called bloated planets.

Acknowledgments

The authors acknowledge the hospitality of the Kavli Institute for Theoretical Physics at UCSB (funded by the NSF through Grant PHY05-51164), where this work started. This work has been supported in part by NASA Grant NNX07AG81G and NSF grants AST 0707628. We also acknowledge funding from the European Community via the P7/2007-2013 Grant Agreement No. 247060. The authors are grateful to the anonymous referee for his/her sharp and enlightening comments.

References

Appendix A: Theory of planetary figures: numerical methods

Here, we briefly outline the method described in Sterne (1939) to compute numerically the response of a body in hydrostatic equilibrium to a perturbing potential⁵ and derive additional formulae. To lowest order (which is consistent with the order of approximation used throughout the present paper) the body response is linear and the total deformation is the sum of the response to each term of the decomposition of the perturbing potential. Let us consider a term of the decomposition of the potential of the form (A.1)where the are tesseral harmonics defined by (A.2)The cos (sin) corresponds to positive (negative) values of m and are the usual associated Legendre polynomials. The reference axis defining θ and ψ may change from one term to the other. For example, the rotation axis is best suited to treat rotational distortion and the line connecting the center of mass of each body is better to describe the tidal distortion. It is shown by Sterne (1939) that to first order, the shape of the distorted level surface of mean radius s (see Zharkov & Trubitsyn 1980, for a detailed definition) takes the form (A.3)where r is the distance between the center and the level surface as a function of θ and ψ and a figure function yet to be calculated. Sterne (1939) shown that, ignoring terms of order , verifies the following differential equation (A.4)with (A.5)Using the variable , this rewrites (A.6)Then, η_l(R) (R being the external mean radius of the object) can be obtained by numerical integration (with η_l(0) = l − 2) and the shape and external potential () of the body are given by, respectively (A.7)and (A.8)In order to compare this numerical model with others, we can compute several observable quantities. By definition, the potential Love number of degree 2 (k₂) is given by (A.9)which yields (A.10)The level Love number (h₂) is given by (A.11)where g is the surface gravity acceleration, and (A.12)as expected for a body in hydrostatic equilibrium.

A.1. Axi-symmetric case

Thus the first gravitational moment (J₂) defined by (A.13)is given by (A.14)No distinction is made between R and R_eq (the equatorial radius) when comparing Eqs. (A.8) and (A.13) because this would only add higher order corrections to J₂ which is already a first order quantity.

For the rotational distortion of the body whose angular velocity is ω, and (A.15)If one is concerned with the external shape, the oblateness (f, see Eq. (36)) of a rotating body is given by (A.16)By extension, one can define J₂ for a tidal perturbation by a secondary of mass M′ at a distance r′, leading to (A.17)but the reference axis is the line connecting the two center of mass and not the rotational axis.

A.2. Triaxial case

While it is tempting to add Eqs. (A.15) and (A.17) to obtain the total J₂ of a body in a close binary, we must remember that the tidal and rotational deformations do not have the same axis of symmetry in general. Taking θ as the colatitude and ψ as the longitude of the body considered, the external gravitational field of the latter reads (A.18)The quadrupole moment in the linear approximation, is given at the surface by (A.19)where θ′ is the angle between the current point and the line connecting the two center of mass. For the coplanar case where the tides raising object orbits in the equatorial plane of the distorted body, cosθ′ = sinθcosψ and thus (A.20)Thus (A.21)and (A.22)All the other moments are equal to 0. Similar decompositions can be used to infer the precise shape of the surface from a sum of perturbing fields. This gives (A.23)which directly translates into a₁, a₂ and a₃ (once R is known) by setting (θ,ψ) equal to (π/2,0), (π/2,π/2) and (0,0), respectively. Translating this into α₁, α₂ and α₃ is a little more complicated because one needs to account for the fact that R > R₀ due to the centrifugal potential. This can be taken into account either numerically – by including the centrifugal force when solving the hydrostatic equilibrium – or analytically using Eq. (32) or (45).

Appendix B: Projected area of a triaxial ellipsoid

B.1. General case

Let us define two coordinate systems. The first one (, , ) is defined by the three main axes of the ellipsoid. In this frame, the equation of the surface of the ellipsoid is (B.1)To compute the projected area of this ellipsoid as it will be seen by the observer, it is easier to put ourselves in another coordinate system defined by the line connecting the center of mass of the system and the observer (toward the observer; ), the projection of the orbital angular momentum on the sky plane (ẑ) and a third axis in the sky plane chosen so that follows the right-hand vector sense. The current position vector (r = (x,   y,   z)) expressed in this frame is thus related to the one expressed in the first coordinate system by a rotation matrix ℛ such as (B.2)With ℛ^trℛ = 1. The equation of the ellipsoid in the new system thus writes (B.3)The exact value of the matrix will depend on the rotation needed and on the angles chosen to represent it. This can be worked out in each specific case. To keep some generality, we take of the form (B.4)The symmetry is ensured by the fact that both of our coordinate systems are orthonormal. The equation of the contour of the projected shadow is given by the fact that the normal to the ellipsoid is normal to the line of sight () there. This assumes a completely opaque body below the isobar chosen to be the surface. The complete calculation of the level at which optical rays that are grazing, at the terminator, have an optical depth close to unity (Hubbard et al. 2001; Burrows et al. 2003; Guillot 2010) – in the present geometry – should give rise to subtle effects but of smaller importance. This reads (B.5)This shows that these points are located on a plane whose equation is (since a ≠ 0) (B.6)Substituting x in Eq. (B.3) by Eq. (B.6) we see that the cross section is an ellipse following the equation (B.7)It is thus possible to find the rotation in the sky plane needed to reduce the ellipse and find its principal axes (p₁,   p₂). If only the cross section (πp₁p₂) is needed, we can use the fact that the determinant of a matrix is independent of the coordinate system so that (B.8)with (B.9)In the case of an edge-on orbit at mid transit, no rotation is needed, ℛ is the identity and thus , , and d = e = f = 0. We retrieve (B.10)

B.2. Coplanar case

If the planet equator and the orbital plane are coplanar, the unit vectors of first coordinate system defined above coincides with the unit vectors defined by the line connecting the two center of mass (from the secondary to the object under consideration; ), its normal in the orbital plane (in the direction of motion; ) and the rotation axis of the body (). If i is the inclination of the orbit with respect to the sky plane and φ the true anomaly defined to be 0 at mid transit, the rotation matrix defined by Eq. (B.2) reads (B.11)The matrix can be computed thanks to Eq. (B.3) giving the a, b, ..., f coefficients and thus Det(ℬ). This gives the project area of the planet or the star at any given point of the orbit as shown on Fig. 5.

Online material

All Tables

All Figures

Fig. 1 Oblateness given by Eq. (36) as a function of the rotation period (in days) at 1 Gyr for planets of mass: 0.3 MJ (dotted), 0.5 MJ (dashed), 1 MJ (dash-dotted), 3 MJ (long dashed), 15 MJ(solid). The oblateness decreases when the mass of the planet increases because massive objects are more compressible (see Sect. 4), have a more intense self-gravity field and are thus less subject to perturbations.

In the text

Fig. 2 The internal pressure-density profile of an irradiated 1.8 MJ planet (solid line). The dashed line represents the best-fit polytropic equation of state. The pressure-density range covered in the inner part of the body (95% in mass) is represented by the thicker part of the solid curve, which is well modeled by a polytropic EOS. As the thin part of the P − ρ curve represents only 5% in mass of the body it is disregarded by the fit.

In the text

	Fig. 3 Polytropic index for non-irradiated planets as a function of the planet’s mass Mp at 100 Myr (dotted), 1 Gyr (dashed) and 5 Gyr (solid). The shaded area represents the uncertainty on the polytropic index for the 5  Gyr case (see text).
In the text

Fig. 4 Polytropic index for strongly-irradiated planets as a function of Mp at 100 Myr (dotted), 1 Gyr (dashed) and 5 Gyr (solid). As the irradiated atmosphere impedes the radiative cooling of the objects, it retards its contraction. Therefore, the non-monotonous behavior observed at the early ages in the non-irradiated case (Fig. 3) is enhanced, even at a later epoch. The bump at the high mass end of the 100 Myr curve is caused by deuterium burning (see text). The shaded area represents the uncertainty on the polytropic index for the 5  Gyr case (see text).

In the text

Fig. 5 Normalized projected area of the planet as a function of its anomaly (φ) for inclinations of the orbit going from i = 90° (lowest curve) to i = 65° (highest curve) by steps of 5° for a WASP-12 b analog on a circular orbit. Top: for the full orbit. Bottom: zoom on the (primary or secondary) transit. The ordinates of the dotted, solid and dashed horizontal lines are respectively (face-on orbit), and .

In the text

Fig. 6 Relative transit depth variation η computed with Eq. (51) as a function of the semimajor axis at 1 Gyr for planets of mass: 0.3 MJ (dotted), 0.5 MJ (dashed), 1 MJ (dash-dotted), 3 MJ (long dashed), 15 MJ (solid). The shaded area shows the zone where higher order terms become non-negligible. The decrease in the transit depth due to tidal interactions is smaller when the mass of the planet increases because massive objects are more compressible (see Sect. 4) and thus less subject to nonspherical deformations.

In the text