© ESO 2018

1 Introduction

A knowledge of the mass and radius of an exoplanet allows determining its mean density, which is the most basic indicator of its composition. Using an example from the solar system, the similarity of the densities of the Earth and Mercury along with the diversity of their sizes allows us to infer, under the assumption that they are both composed of rocks and metals, that the metallic component in Mercury is larger than in the Earth (Ash et al. 1971). Planets are roughly spherical objects because self-gravitation overcomes material strength for bodies larger than about a few hundred kilometers. However, a spherical shape does not guarantee a differentiated interior. An integral measure of the concentration of mass – and thus, indirectly of differentiation and interior structure – is provided by the normalized moment of inertia (MoI), defined for a spherical body of mass M and radius R as $$\text{MoI}=\frac{1}{M{R}^2}{\displaystyle {\int }_0^V\rho }(r)r_{\perp }^2\text{d}V,$$(1)

where V is the volume, ρ (r) is the density as a function of the radius r, and r_⊥ is the distance from an axis passing through the center of mass of the body. Planets are not perfectly spherical, and the value of the MoI will in general depend on the chosen axis. However, there are only three independent moments of inertia, usually indicated with A, B, and C, and the MoI as defined in Eq. (1) can be taken as representing the normalized mean moment of inertia given by (A + B + C)∕(3MR²). A homogeneous spherical body has a value of MoI of 0.4, while in a gravitationally stable body where density increases with depth, 0 ≤ MoI < 0.4, with 0 representing the value for a point mass.

The most common causes of departure from sphericity at large spatial scales for planets are rotation and gravitational interactions with other bodies: parent star, other planets, and moons. In general it is possible to express these perturbations in terms of a potential. In the case of a tide-inducing perturber of mass M_p at a distance d, the tide-generating potential W within the planet can be expanded in spherical harmonics as (e.g., Kaula 1966) $$W\left(r,\psi \right)=\frac{G{M}_{\text{p}}}{d}{\sum }_{n=2}^{\infty }{\left(\frac{r}{d}\right)}^n{P}_n\left(\cos \psi \right)={\sum }_{n=2}^{\infty }{W}_n,$$(2)

where r is the coordinate of a point within the planet, ψ is the angle with respect to the center of mass between the point and the perturber, and P_n is the Legendre polynomial of degree n. Similarly, a rotational potential can be written as a harmonic term of degree 2 (e.g., Murray & Dermott 1999). Under the assumption that the response of the planet to these perturbing potentials is linear, each harmonic degree of the perturbing potential W will generate a response with the same degree, $V_n^{\text{p}}$. At the surface, $$V_n^{\text{p}}=k_n{W}_n,$$(3)

where k_n is the gravitational Love number, introduced by A. E. H. Love (although not with this name, Love 1911), and can be regarded as a useful way of condensing in a single parameter the many unknowns controlling the gravitational response of the planet to the perturbation. In addition to k_n, the Love number h_n describes the radial displacement of the surface that results from the presence of the perturbing potential. The equipotential surface is defined by the external potential W_n and the additional potential k_nW_n, corresponding to (1 + k_n)W_n∕g₀, with g₀ the gravitational acceleration at the surface. For a fluid planet the location of the surface, which is an equipotential surface, corresponds to the radial displacement h_n W_n∕g₀, thus resulting in the simple relation h_n = 1 + k_n (Munk & MacDonald 1960).

In general, the forward calculation of the Love numbers for a given planetary interior structure requires the knowledge of the density, viscosity, and rigidity profiles within the planet, in addition to the timescale of the perturbing potential (e.g., Alterman et al. 1959). However, if the planet is in hydrostatic equilibrium, that is, if it responds as a fluid, only the density profile is required (e.g., Sterne 1939; Gavrilov et al. 1975). Thus, both the moment of inertia and the fluid Love numbers only depend on the distribution of matter in the interior of the planet, and there exists an equation, known as the Darwin–Radau equation, which provides a link between the MoI and k₂ (or h₂). However, this relation is only an approximation (Murray & Dermott 1999; Kramm et al. 2011; Helled et al. 2011), and in this paper, the MoI and the Love numbers are calculated independently, using Eq. (1) and the method presented in Sect. 3, respectively.

The expressions we derive for the Love numbers k_n are based on the matrix-propagator method, which seeks the solution to a system of differential equations through the use of matrices, and traces back to the ideas of Thomson (1950) and Haskell (1953), which in turn are part of the theoretical framework developed by Volterra (1887), as illustrated by Gilbert & Backus (1966). Here, we apply a procedure similar to that developed by Wolf (1994). The matrix-propagator method requires the knowledge of the radial density profile, discretized at a given number of internal interfaces, without requiring the knowledge (or numerical calculation) of the local derivative of the density profile.

The paper is structured as follows. In Sect. 2 we illustrate the calculation of the Love numbers k_n for the simple case of a homogeneous planet or a planet with two constant density layers, which we then generalize to the case of a planet defined by any number of layers in Sect. 3. As an illustration of the method, in Sect. 4.1 we apply the theory to the case of GJ 436b. The outer gaseous envelope of this planet might harbor a variety of structures in the deep interior (e.g., Adams et al. 2008). We also apply the theory to two models of Jupiter, whose core could have a well-defined outer surface or be diluted due to erosion (e.g., Guillot et al. 2004). We conclude with a discussion of the observability and meaningfulness of using forthcoming measurements of k₂ in the investigation of the interior structure of extrasolar planets.

2 Potential method for calculating the fluid Love numbr k_n

The computation of the fluid Love numbers requires the solution of the equation (e.g., Gavrilov et al. 1975) $${T^{″}}_n\left(r\right)+\frac{2}{r}{T^{\prime }}_n\left(r\right)+\left[\frac{4\pi G{\rho }^{\prime }\left(r\right)}{V^{\prime }\left(r\right)}-\frac{n\left(n+1\right)}{r^2}\right]T_n\left(r\right)=0,$$(4)

where r, V, and ρ are the radial coordinate, the gravitational potential, and the density of the unperturbed body (i.e., spherically symmetric). A (double) prime indicates (double) derivation with respect to r, and V′ (r) = −g(r) is the gravitational acceleration. The function T has the dimensions of a potential and describes the total change in potential according to $$V^{\text{p}}+W=\frac{1}{R{g}_0}T{\left(\frac{R}{r}\right)}^{n}W,$$(5)

where the perturbation-inducing potential W generates the additional potential V^p within the body. The gravitational acceleration at the surface is g₀. The perturbed potential is proportional to the perturbation-inducing potential and at the surface of the deformed body (r = R), the relation is $$V_n^{\text{p}}\left(R\right)=k_n{W}_n\left(R\right),$$(6)

where the subscript n indicates the degree in the harmonic expansion of W. From the two equations above, the Love number k_n is obtained as $$k_n=\frac{T_n\left(R\right)}{R{g}_0}-1.$$(7)

2.1 Boundary and interface conditions

At the center, the solution for T must be finite. At internal density discontinuities of radius r_i, both T and its derivative must be continuous. This requirement corresponds to $$\begin{array}{lll}{T}_n\left(r_i^{-}\right)\hfill & =\hfill & T_n\left(r_i^{+}\right),\hfill \\ {T^{\prime }}_n\left(r_i^{-}\right)\hfill & =\hfill & {T^{\prime }}_n\left(r_i^{+}\right)+\frac{4\pi G}{V^{\prime }\left(r_i\right)}\left[\rho \left(r_i^{-}\right)-\rho \left(r_i^{+}\right)\right]T_n\left(r_i\right),\hfill \end{array}$$

where a plus (minus) indicates that the variable is evaluated right above (below) the corresponding radius. At the surface $\left(r=R\right)$, the continuity condition is (Zharkov et al. 1985) $$T^{\prime }\left(R^{-}\right)=-\frac{\left(n+1\right)}{R}T\left(R^{-}\right)+4\pi G\rho \left(R\right)H\left(R\right)+\left(2n+1\right)g_0,$$(10)

where the function H is in relation with T and the fluid radial displacement Love number h_n through $$\begin{array}{lll}H\hfill & =\hfill & -\frac{T}{V^{\prime }},\hfill \\ h_n\hfill & =\hfill & \frac{H\left(R\right)}{R}.\hfill \end{array}$$

2.2 Solution for a homogeneous sphere

If the body has constant density, ρ′ = 0, and Eq. (4) reduces to an Euler-Cauchy equation. Using the trial solution T(r) = r^c results in the characteristic equation $$c^2+c-n(n+1)=0.$$(13)

Since the discriminant is positive (n is an integerlarger than 0) and the product of the two solutions is negative, there are two real solutions with opposite signs for the exponent c, namely, c₁ = n and c₂ = −(n + 1). The general solution is then the linear combination $$T\left(r\right)=A{r}^{c_1}+B{r}^{c_2},$$(14)

where A and B are unknown constants of integration. The requirement that the function T is finite at the center implies that B = 0. To determine the second integration constant A, the application of the surface boundary condition from Eq. (10) gives $$T_n(r)=\frac{R(2n+1)g_0}{(2n-2)}{\left(\frac{r}{R}\right)}^n.$$(15)

The Love number k_n of order n is, from Eq. (7), $$k_n=\frac{3}{2n-2}.$$(16)

The value ofk_n as a function of n is shown in Fig. 1. For n = 2 the value of k_n is 1.5, whichrepresents the limit corresponding to a normalized moment of inertia of 0.4.

Fig. 1 Fluid Love number k as a function of the degree n for a homogeneous sphere.

2.3 Surface deformation

With the value of the fluid Love number k_n, the shape ofthe equipotential surface, which corresponds to the physical surface in the fluid limit, can be evaluated from the value of h_n= 1 + k_n and the knowledge of the perturbing potential (Sect. 1). From the degree-2 term of the tide-inducing potential W in Eq. (2), the tidally induced radial deformation at the surface, $u_R^{\text{Tidal}}$, can be written as $$\begin{array}{lll}{u}_R^{\text{Tidal}}\left(\psi \right)\hfill & =\hfill & \frac{h_2}{g_0}{W}_2\left(R,\psi \right),\hfill \\ \hfill & =\hfill & h_{2}R\left(\frac{M_{\text{S}}}{M}\right){\left(\frac{R}{d}\right)}^3{P}_2\left(\cos \psi \right),\hfill \end{array}$$

where M_S, M, and d are the stellar (i.e., perturber) mass, the planetary mass, and the distance of the perturber. The shape is rotationally symmetric with respect to the line connecting the center of mass of the body with the perturber. Equation (18) shows that the larger, less massive, and the closer to the perturber the body, the larger its surface radial deformation (since $u_R^{\text{Tidal}}\propto R^4/M/d^3$). We note that in general the distance to the perturber d varies, being constant and equal to the semimajor axis of the orbit only for circular orbits. Thus, the shape of the planet, which is assumed to respond as fluid, continually evolves as the distance d varies.

The rotational potential can be written as a degree-2 harmonic term as (e.g., Murray & Dermott 1999) $$Z\left(r,\theta \right)=\frac{1}{3}{\omega }^2{r}^2\left[P_2\left(\cos \theta \right)-1\right],$$(19)

where ω is the angular rotational rate and θ is the colatitude measured from the rotation axis. The rotational potential is symmetrical with respect to the axis of rotation. Asfor the tidally induced radial deformation, the degree-2 rotationally induced surface radial deformation is $$\begin{array}{lll}{u}_R^{\text{Rotational}}\left(\theta \right)\hfill & =\hfill & \frac{h_2}{g_0}Z\left(R,\theta \right),\hfill \\ \hfill & =\hfill & \frac{h_2}{3G}\frac{R^4{\omega }^2}{M}\left[P_2\left(\cos \theta \right)-1\right].\hfill \end{array}$$

Equation (21) shows that the larger and less massive the body and the faster it rotates, the larger its deformation (since $u_R^{\text{Rotational}}\propto R^4{\omega }^2/M$).

In general, the spin axis does not coincide with the axis pointing at the tidally inducing perturber, and, assuming that the two perturbations can be added linearly, their combined effect is obtained by expressing them in the same frame of reference, that is, by applying the addition theorem of spherical harmonics. In the frame of reference of the rotational potential, $$u_R\left(\theta ,\varphi \right)=u_R^{\text{Rotational}}\left(\theta \right)+u_R^{\text{Tidal}}\left(\theta ,\varphi \right)\propto \frac{R^4{h}_2}{M},$$(22)

where the dependence on the longitude ϕ appears through the rotation of the reference system where the tidal perturbation is evaluated. The perturbation thus depends on the fundamental properties ofthe planet R, M, and h₂. For degree 2, the combination of rotation and tidal distortion corresponds to the sum of a rotationally induced oblate spheroid (Eq. (21)) with a tidally induced prolate spheroid (Eq. (18)).

2.4 Planet with two constant density layers

This model represents the simplest approximation for a differentiated terrestrial (gaseous) planet, where a constant-density mantle (gaseous envelope) overlies a constant-density core. We indicate with ρ_c and ρ_m the inner and outer layer densities and with α the ratio of the inner layer radius r_c to the radius R. The basic relations for the mean density ρ and the mean moment of inertia MoI are $$\begin{array}{l}\rho ={\rho }_{\text{c}}{\alpha }^3+{\rho }_{\text{m}}\left(1-{\alpha }^3\right),\\ \text{MoI}=\frac{2}{5}\left[\frac{{\rho }_{\text{c}}}{\rho }{\alpha }^5+\frac{{\rho }_{\text{m}}}{\rho }\left(1-{\alpha }^5\right)\right].\end{array}$$

Figure 2 illustrates the interplay among the densities of the two layers (normalized to the mean density of the object) and the radius of the core (normalized to the planetary radius). The normalized core density is plotted on a log scale, which includes the case of a small dense core overlaid by a thick, almost massless, envelope. For the Love numbers k_n, in each layer the solution is the same as in the homogeneous case, Eq. (14). Continuity across the two layers must be enforced through Eqs. (8) and (9). It is then possible to obtain a closed-form solution for k_n.

The model with two constant density layers allows identifying the basic dependences among the parameters. For a planet with a mass and radius similar to GJ 1214b (M = 6.55 M_⊕, R = 2.68 M_⊕, Charbonneau et al. 2009), Fig. 3 shows the values of k₂, the normalized moment of inertia MoI, and the maximum surface tidal deformation, which is reached at the sub-stellar and anti-stellar points where P₂ (cosθ) = 1 in Eq. (18), as a function of the core radius, which corresponds, through Eq. (23), to a core density. For the ratio between the mantle density and the average density, we use either 0.73, a value similar to the Earth case, or 0.01, corresponding to a model where mass is mostly concentrated in a high-density core. From an inspection of the figure, we note the following:

• 1.

  The smaller (and thus, the denser) the core, the smaller both k₂ and MoI, consistent with the interpretation of the fluid Love number as a measure of the concentration of mass (e.g., Kramm et al. 2011).

• 2.

  A body where most of the mass is concentrated in the core, as the case for ρ_m∕ρ = 0.01 illustrates, has a value of k₂ that rapidly approaches 0 as the core decreases in size (specifically, for r_C∕R = 0.08, k₂ = 6 × 10⁻³).

• 3.

  When k₂ approaches 0, this does not imply that the body is not deformed, since the deformation depends on h₂ = k₂ + 1. The limiting case of k₂ = 0 (all the mass in a point core overlaid by a massless envelope of radius R) simply corresponds to h₂ = 1 in Eq. (18).

• 4.

  A comparison of the two cases in Fig. 3, which only differ in the distribution of the mass in the interior, shows that the three quantities plotted, k₂, MoI, and the deformation, decrease with increasing concentration of the mass.

Fig. 2 Model with two constant density layers. Normalized density of the mantle vs. density of the core. Colors indicate the normalized radius of the core. The density of the core is plotted in logarithmic scale.

Fig. 3 Model of a planet with two constant density layers with M = 6.55 M⊕ and R = 2.68 R⊕. Values of k2 (first row), normalized moment of inertia MoI (second row), and maximum surface deformation (third row, assuming MS = 0.157 M⊙ and a = 0.0143 AU in Eq. (18)), for ρm∕ρ = 0.73 (left column), and for ρm∕ρ = 0.01 (right column). Colors indicate the normalized density of the core. The color bar is logarithmic.

3 Matrix-propagator approach to the calculation of k_n

With ${\rho }^{\prime }\left(r\right)=0,$ Eq. (4) can be recast as a system of two first-order differential equations by introducing the function P = dT∕dr. Before writing the system in matrix form, the variables T and P are non-dimensionalized as follows (e.g., Wolf 1994): $$\begin{array}{lll}T\hfill & =\hfill & \frac{GM}{R}{y}_1,\hfill \\ P\hfill & =\hfill & \frac{dT}{dr}=\frac{GM}{R^2}\left(\frac{R}{r}\right)y_2,\hfill \end{array}$$

where y_i are non-dimensional variables. The derivatives of T and P with respect to r are $$\frac{\text{d}T}{\text{d}r}=\frac{GM}{R^2}\frac{d{y}_1}{ds};\frac{dP}{dr}=\frac{GM}{R^3}\left(\frac{1}{s}\frac{d{y}_2}{ds}-\frac{y_2}{s^2}\right),$$(27)

where s = r∕R is the non-dimensional radius. The non-dimensional system of equations is $$\begin{array}{lll}\frac{\text{d}{y}_1}{\text{d}s}\hfill & =\hfill & \frac{y_2}{s},\hfill \\ \frac{\text{d}{y}_2}{\text{d}s}\hfill & =\hfill & \frac{n\left(n+1\right)}{s}{y}_1-\frac{1}{s}{y}_2.\hfill \end{array}$$

The system can be written in matrix form as $$\frac{\text{d}}{\text{d}s}\left[y\right]=\frac{\text{d}}{\text{d}s}\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]=\left[\begin{array}{cc}0& \frac{1}{s}\\ \frac{n(n+1)}{s}& -\frac{1}{s}\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right].$$(30)

Assuming power solutions of the form $y_i=s^{k^{(j)}}{\text{y}}_i^{(j)}$, with ${\text{y}}_i^{(j)}$ a constant, the system (Eq. (30)) becomes $$\left[\begin{array}{cc}{k}^{(j)}& -1\\ -n(n+1)& 1+k^{(j)}\end{array}\right]\left[\begin{array}{c}{\text{y}}_{\text{1}}{}^{(j)}\\ {\text{y}}_{\text{2}}{}^{(j)}\end{array}\right]=0,$$(31)

which admits non-trivial solutions only if the determinant vanishes, that is, if $$k^{(j)}{}^2+k^{(j)}-n\left(n+1\right)=0.$$(32)

The last equation corresponds to the characteristic Eq. (13), with solutions $$k^{(1)}=n,k^{(2)}=-\left(n+1\right).$$(33)

Thus, the general solution for y_i is the linear combination $$y_i={\sum }_{j=1}^2{C}_{(j)}{s}^{k(j)}{\text{y}}_i^{(j)},$$(34)

where the solution vectors ${\text{y}}_i^{(j)}$ are $${\text{y}}_i^{(j)}=\left[\begin{array}{c}1\\ k^{(j)}\end{array}\right].$$(35)

The general solution in matrix form is then $$y=PC,$$(36)

where $$\begin{array}{lll}P\hfill & =\hfill & \left[\begin{array}{cc}{s}^n& s^{-(n+1)}\\ n{s}^n& -(n+1)s^{-(n+1)}\end{array}\right],\hfill \\ C\hfill & =\hfill & {\left[C_1{C}_2\right]}^{\text{T}}.\hfill \end{array}$$

P is the propagator matrix and C is the vector of constants of integration.

3.1 Internal boundary conditions in matrix form

Indicating with r_i the radius of a density discontinuity, the two internal boundary conditions, Eqs. (8) and (9), can be expressed in matrix form as $${\left[y^{+}\right]}_{r_i}={\left[\begin{array}{c}{y}_1^{+}\\ y_2^{+}\end{array}\right]}_{r_i}=\left[\begin{array}{cc}1& 0\\ -\frac{4\pi G{r}_i\Delta \rho }{g\left(r_i\right)}& 1\end{array}\right]{\left[\begin{array}{c}{y}_1^{-}\\ y_2^{-}\end{array}\right]}_{r_i},$$(39)

where a plus (minus) indicates that the variable is evaluated right above (below) the boundary. The density difference $\Delta \rho =\rho (r_i^{-})-\rho (r_i^{+})$ is positive for gravitationally stable planets, for which density increases downward.

3.2 Propagator

The planetary models used in this study are made of a series of constant property layers. Equation (36) provides the general solution in each layer. The interface conditions between layers in matrix form corresponds to the square matrix in Eq. (39): $$B\left(r_i,\Delta {\rho }_i\right)=\left[\begin{array}{cc}1& 0\\ -\frac{4\pi G{r}_i\Delta \rho }{g\left(r_i\right)}& 1\end{array}\right].$$(40)

With reference to the interface at r = r_i−1 in Fig. 4, the two solutions to be matched are $$\begin{array}{lll}\hfill & \hfill & y^{(i)}(s_{i-1})=P_i(s_{i-1})C^{(i)},\hfill \\ \hfill & \hfill & y^{(i-1)}(s_{i-1})=P_{i-1}(s_{i-1})C^{(i-1)}.\hfill \end{array}$$

The interface matrix B(r_i−1) provides the connection between the two solutions: $$y^{(i)}(s_{i-1})=B(r_{i-1},\Delta {\rho }_{i-1})y^{(i-1)}(s_{i-1}).$$(43)

With the last three expressions, the vector of constants C⁽ⁱ⁾ can be expressed in terms of the vector of constants C⁽ⁱ⁻¹⁾: $$C^{(i)}={\left[P_i(s_{i-1})\right]}^{-1}B(r_{i-1},\Delta {\rho }_{i-1})P_{i-1}(s_{i-1})C^{(i-1)}.$$(44)

From Eqs. (41) and (44), the solution for y⁽ⁱ⁾ at r = r_i is $$y^{(i)}(s_i)=P_i(s_i){\left[P_i(s_{i-1})\right]}^{-1}B(r_{i-1},\Delta {\rho }_{i-1})P_{i-1}(s_{i-1})C^{(i-1)},$$(45)

and the vector of integration constant C⁽ⁱ⁾ no longer appears. By extending this approach to each layer, the solution y^(N) at r_N can then be written as the product of a sequence of terms and the vector of constants at the center of the planet: $$\begin{array}{lll}{y}^{(N)}(s_N)\hfill & =\hfill & \left\{{\displaystyle \prod _{k=2}^N{P}_k}(s_k){\left[P_k(s_{k-1})\right]}^{-1}B(r_{k-1},\Delta {\rho }_{k-1})\right\}\hfill \\ \hfill & \hfill & \cdot P_1(s_1)C^{(1)}.\hfill \end{array}$$

Fig. 4 Schematic structure of the vector of solution y and radius r for a planetary model. The center of the planet is at r = 0, the surface at r = R. Each homogeneous layer is indicated by the index corresponding to the outer boundary of the layer. Similar indexing applies to the non-dimensional radius s.

3.3 Simplification of the propagator matrix product

To compute the product of P with its inverse in Eq. (47), it is convenient to write, using Eq. (37), $$P(r)=\left[\begin{array}{cc}1& 1\\ n& -(n+1)\end{array}\right]\left[\begin{array}{cc}{r}^n& 0\\ 0& r^{-(n+1)}\end{array}\right],$$(48)

so that $${\left[P(r)\right]}^{-1}=-\frac{r}{(2n+1)}\left[\begin{array}{cc}{r}^{-(n+1)}& 0\\ 0& r^n\end{array}\right]\left[\begin{array}{cc}-(n+1)& -1\\ -n& 1\end{array}\right].$$(49)

With the last two expressions $$P(r_j){\left[P(r_{j-1})\right]}^{-1}=D(n)X(n,r_j,r_{j-1}){\left[D(n)\right]}^{-1},$$(50)

where $$\begin{array}{lll}\hfill & \hfill & D(n)=\left[\begin{array}{cc}1& 1\\ n& -(n+1)\end{array}\right],\hfill \\ \hfill & \hfill & X(n,r_j,r_{j-1})=\left[\begin{array}{cc}{\left(\frac{r_j}{r_{j-1}}\right)}^n& 0\\ 0& {\left(\frac{r_j}{r_{j-1}}\right)}^{-(n+1)}\end{array}\right].\hfill \end{array}$$

3.4 Solution

With the results of Sects. 3.2 and 3.3, we can write the general solution at the surface as $$y(R)=My(0),$$(53)

where the matrix M is completely defined by the interior model, that is, by the density as a function of radius. Explicitly, $$\left[\begin{array}{c}{y}_1(R)\\ y_2(R)\end{array}\right]=\left[\begin{array}{cc}{M}_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right]\left[\begin{array}{c}{y}_1(0)\\ y_2(0)\end{array}\right],$$(54)

where y₂(0) = 0, according to the definition of P in Eq. (26). The surface boundary conditions, Eq. (10), in terms of y₁ and y₂ read $$y_2(R)=-(n+1)y_1(R)+\frac{4\pi G\rho (R)R}{g_0}{y}_1(R)+(2n+1).$$(55)

Using Eqs. (55) and (54) $$\left[\begin{array}{cc}{M}_{11}& -1\\ M_{21}& (n+1)-\frac{4\pi G\rho (R)R}{g_0}\end{array}\right]\left[\begin{array}{c}{y}_1(0)\\ y_1(R)\end{array}\right]=\left[\begin{array}{c}0\\ (2n+1)\end{array}\right],$$(56)

whose inversion provides the boundary conditions of the problem. From the value of y₁ (R), through Eqs. (25) and (7) the Love number for any n can be determined.

3.5 Sub-layering

Given the discretized nature of the matrix-propagator approach, the solution does depend on the assumed number of layers, and a suitable choice should be made to ensure an accurate result. To illustrate the dependence on the layering, we generated seven interior models where the density profile is obtained as the portion of a Gaussian curve with zero mean and standard deviation σ: $$\rho \left(r\right)=\exp \left[-\frac{r^2}{2{\sigma }^2}\right].$$(57)

By using only values of r between 0 and 1, these models cover a range of concentration of mass in the interior, from an almost homogeneous distribution (σ = 9.9) to a concentrated one (σ = 0.1). Their density profiles are plotted in Fig. 5. We computed the moment of inertia and the Love number k₂ as a function of the number of sub-layers used to discretize the analytical profile expressed in Eq. (57), and the results are shown in Fig. 6. Both k₂ and the MoI converge to a value that represents the “continuous” limit. This is the limit of an infinite number of sub-layers with zero thickness. In the right column of Fig. 6 we plot the relative variation of the values of k₂ and the MoI with respect to the limit value, which we take as the one corresponding to 10⁵ layers. The more homogeneous a planet, the smaller the density variations as a function of the radius, and correspondingly, the smaller the number of layers necessary to retrieve an accurate result, as the comparison of the blue and yellow curves indicates. In all cases, using a number of sub-layers in excess of about 10³ allows the retrieval of k₂ and the MoI with a precision better than 1%. The required accuracy of the modeling can guide the choice of the number of sub-layers to be used. Figure 7 shows that the computation of k₂ for a model with 10⁵ sub-layer requires approximately 1 s. When thousands of models are involved, Figs. 6 and 7 provide a guide to the tradeoff between precision and speed.

The results regarding the precision as a function of the number of sub-layers do depend on the assumed interior model, and the Gaussian curves used in this section represent an example that can be adapted on a case-by-case basis. However, in computing the Love numbers of the Earth (Sect. 3.6), of planet GJ 436b (Sect. 4.1), and of Jupiter (Sect. 4.2), we found that a number of sub-layers between 10³ and 10⁴ results in a precision of 1% or better.

Fig. 5 Set of normalized interior density profiles obtained using Gaussian curves with different standard deviations, as indicated in the legend. The yellow curve represents a homogeneous planet, and the blue curve shows a planet with a high concentration of mass in the interior.

Fig. 6 For the models of Fig. 5, the left column shows values of the Love number k2 and of the normalized moment of inertia (MoI) as a function of the number of layers used to discretize the density profile. For a given number of layers, the right column shows the relative variation with respect to the value for 105 layers.

3.6 Validation

We used three analytical density profiles for Jupiter, linear, quadratic, and polytropic, as proposed by Gavrilov et al. (1976). With these, we inferred values for the Love number k_n for n =2,..., 7. Gavrilov et al. (1976) solved the full Eq. (4), while we set the term proportional to ρ′ equal to zero and used 10³ sub-layers to reproduce the smooth variation of density with radius. Our results for the Love numbers match those of Gavrilov et al. (1976), showing that the matrix-propagator method returns the values obtained with the solution of the full equation, provided an appropriate number of sub-layers is chosen (Sect. 3.5).

In order to assess the precision of the method, we tested it using PREM, the preliminary reference Earth model (Dziewonski & Anderson 1981), as tabulated in the file PREM_1s.csv of the IRIS database (Trabant et al. 2012). For the crust and upper mantle, we used the PREM data directly (108 layers for r > 5701 km). For each of inner core, outer core, and lower mantle, we fit the PREM data with a third-degree polynomial, which was then used to generate 10^m sub-layers, with m = 2,...,5 (i.e., we built four Earth models with a total number of layers of 3 × 10^m + 108). The corresponding values of k₂ and momentof inertia are listed in Table 1. Even with a few hundred layers, the results are below one per thousand of the reference values.

Finally, we compared the results obtained with the matrix-propagator method with those computed with the method of the concentric Maclaurin spheroids (Hubbard 2013). Using the preferred density profile of Jupiter from Wahl et al. (2016), we obtain a value of k₂ that differs by less than 0.1% from the non-rotating Jupiter case (Table 4 in Wahl et al. 2016). The comparison with the rotating Jupiter case cannot be made in the framework of the Love numbers as defined here (Sect. 5.1).

4 Applications

4.1 GJ 436b

The extrasolar planet GJ 436b is a classical example of the degeneracy¹ of mass-radius relationships (e.g., Adams et al. 2008; Nettelmann et al. 2010; Kramm et al. 2011). Here, we followed the analysis and nomenclature of Adams et al. (2008) and considered three possible interior structures, where a gaseous helium or hydrogen/helium envelope surrounds an Earth-like, a rocky, or an ocean-like interior. Using a newly developed interior structure code (Baumeister et al. 2018), we recalculated the models of Adams et al. (2008) to make them compatible with the most recent values for the mass (21.4 M_⊕, Trifonov et al. 2018) and radius (4.191 R_⊕, Turner et al. 2016)of the planet. The profiles we obtained fit the observed mass and radius with a relative error smaller than 0.1%. To these profiles we applied the matrix-propagator method of Sect. 3 to compute the Love numbers k_n for n =2,..., 8. The density profiles and the corresponding Love numbers are plotted in Fig. 8. We note that the models considered here do not span the entire range of proposed interior structures of the planet (see, e.g., Nettelmann et al. 2010; Kramm et al. 2011), but are used as illustration.

The value of k_n decreases with increasing n, as in the case of the homogeneous sphere (Fig. 1). To first order, the absolute value of k₂ is controlled by the density difference between the center and the surface, which is a proxy for the degree of mass concentration. The Earth-like interior model, with a metallic core that reaches a density in excess of 30 g cm⁻³, has k₂ = 0.055, while the ocean-like interior model, with a central density of about 12 g cm⁻³, has k₂ = 0.160. The rocky interior model, with an intermediate central density, has a value of k₂ = 0.082. Unlike the moment of inertia, the value of k₂ tends to rapidly decrease toward zero as the mass concentration increases (compare the left panels of Fig. 6). This observation explains why these models, which have an extended light gaseous envelope and thus are very concentrated, have all similar low values for k₂.

In general, the higher n, the closer the region of the interior that contributes to the corresponding k_n (e.g., Gavrilov et al. 1976). Thus, it is expected that for increasing n the curves tend to converge, since these models have a similar gaseous envelope. This is true in particular for the Earth-like and rocky models, where the thickness of the envelope is similar.

The most likely parameter to be measured for exoplanets is k₂ (Ragozzine & Wolf 2009; Hellard et al. 2018), and Fig. 8 shows that an accurate determination of k₂ could more easily distinguish the ocean-like interior model from the other two. Although these three models are only illustrative of the degeneracy of a mass and radius determination (Adams et al. 2008), the additional modeling of the Love number k₂ shows how its potential measurement would help in breaking the mass-radius degeneracy, at least partly.

4.2 Jupiter-like hot Jupiter

Jupiter is a fast rotator, and the calculation of its Love numbers requires the use of the concentric Maclaurin spheroids method (Sect. 5.1). Here, we use it as a representative model of hot Jupiters, for which the linear theory developed above is accurate. The deep interior of this gas giant is enriched in heavy elements, which could be segregated into a compact, high-densitycore or diluted into a more extended, enriched central region (e.g., Guillot et al. 2004). We used two density profiles appropriate for these two scenarios (from Wahl et al. 2017b) and computed the Love numbers k_n for n =2,..., 8. In Fig. 9 we plot the profiles and the differences in the values of k_n, given that they are quite similarexcept for n = 2. The higher n, the shallower the region of the interior that mostly contributes to the corresponding k_n (e.g., Gavrilov et al. 1976, Fig. 2b). Thus, given the similarity of the two profiles for r ≳ 0.15, the Love numbers are almost identical except for k₂, which is larger for the diluted core (k₂ = 0.5378 versus k₂ = 0.5287), given its smoother density profile. Even with the high-quality data returned by the Juno mission, the degree of mixing of the core is still under investigation (Stevenson 2018).

5 Discussion

The possibility of improving our knowledge of the interior structure of exoplanets beyond the information provided by the mass and radius rests on the availability of additional data regarding atmospheric composition (e.g., Madhusudhan et al. 2016), stellar mass and composition (e.g., Dorn et al. 2015), and measurement of the fluid Love number k₂ (possibly of k_n for n > 2). While some or all of these parameters are expected to be measured in the future for an increasing number of extrasolar planets, k₂ is the most direct constraint on the interior structure, given its dependence on the density profile.

5.1 Nonlinear effects in the calculation of the Love numbers

The Love numbers as defined here depend only on the density profile, and thus, they represent intrinsic properties of a given planet, much like the case of the moment of inertia defined in Eq. (1). From the knowledge of the Love numbers, the response of the planet to a given perturbation can be determined. In the case of tidal perturbations, we would apply Eqs. (3) and (18) to estimate the tidally induced modification of the gravitational field and of the shape of the surface. A similar approach applies to the case of rotational perturbation (Sect. 2.3). Thus, under the assumption of linear response (Eq. (3)) and of linear combination of different perturbations (e.g., Eq. (22)), from the observation of the gravitational field and/or the shape, and by knowing the parameters of the perturbing potential, we may be able to invert for the Love numbers, and thus for the degree of concentration of mass in the interior. This simple strategy would not be accurate if the response were not linear, that is, if a given degree n of the perturbing potential induced a response in a degree q≠n, and if the perturbations induced by different processes, typically, rotation and tides, could not be added linearly. An example of a nonlinear response is provided by the Earth, where the degree-2 tidal perturbations of the Sun and the Moon induce a response in the degree-4 gravitational harmonics. However, the amplitude is more than three orders of magnitude smaller than the corresponding correction for n = 2 (Petit & Luzum 2010). When the rotational perturbation is much stronger than the tidal perturbation, nonlinear effects appear, which induce an increase of k₂ by about 10% for the fast rotators Jupiter and Saturn (Wahl et al. 2016, 2017a). These nonlinear effects can be accurately estimated with the method of the concentric Maclaurin spheroids (Hubbard 2013), an approach that requires more involved computations than the propagator-matrix developed here.

The appearance of nonlinear effects make the Love numbers dependent on both the planet interior structure and the dynamical environment of the planet. Thus, they are no longer a fundamental property of the planet. In this work we focused on the linear theory. This approach has the main advantage that the Love numbers do represent a measure of the internal concentration of mass (Fig. 3), they can be straightly compared among different objects (Figs. 8 and 9), and their computation is quite fast (Fig. 7). In addition, tidally locked hot Jupiters represent one of the best targets for the measurement of k₂ (Ragozzine & Wolf 2009; Batygin et al. 2009), and their rotational and tidal perturbations are comparable, thus making the linear theory applicable (Ragozzine & Wolf 2009; Wahl et al. 2017a).

Fig. 7 Speed of computation vs. number of layers in the calculation of the Love number k2 for a single-density profile. The curve is similar for n = 2,..., 8. A 2.5 GHz Intel Core i7 processor was used.

Fig. 8 Top panel: three possible density profiles for GJ 436b with the same mass and radius, 21.4 M⊕ and 4.19 R⊕. A low-density hydrogen or hydrogen/helium envelope surrounds an Earth-like (brown), rocky (orange), or ocean-like (blue) interior. Bottom panel: corresponding values of the Love numbers kn.

Fig. 9 Top panel: interior models of Jupiter with a compact high-density core (black) and a diluted and extended central region (red). Data are taken from Wahl et al. (2017b). Bottom panel: difference in the values of the Love numbers kn for the two models.

5.2 Observability and interpretation

The value of k_n can be obtained by inverting the transit light curve for the shape of the surface of the planet (Correia 2014; Hellard et al. 2018), which can be modeled with the Love numbers h_n and an expression for the radial surface deformation (e.g., Eq. (22)). Under the assumption of hydrostatic equilibrium, which assumes that the surface of the planet represents an equipotential surface corresponding to the body’s tidal and rotational potentials, there is a simple relation between the Love numbers h_n and the Love numbers k_n (Sect. 1 and, e.g., Munk & MacDonald 1960): $$h_n=1+k_n.$$(58)

Thus, the determination of h_n would simply translate into the determination of k_n, which provides the additional constraint on the interior structure.

The tidal and rotational potentials that modify the shape of the planet induce a related modification of the gravitational field, which in turn modifies the gravitational interaction of the planet with the parent star and with additional planets in the system, if any are present. In general, this modification results in orbital perturbations and evolution. The evolution is associated with dissipative processes, which depend on the Love number k₂ and the dissipation parameter Q of both the star and the planet (e.g., Goldreich & Soter 1966; Lainey et al. 2017). However, the variation of the orbital parameters occurs on a variety of timescales, and for some specific dynamical configurations, there exist fast components, like the apsidal precession, that do not depend on the dissipation parameters and are controlled by the value of k₂ of the planet (Mardling 2007; Ragozzine & Wolf 2009; Batygin et al. 2009). For such cases, even with the relatively short temporal baseline of extrasolar planet observations at the present time, there are some first successful attempts at constraining or placing bounds on the values of k₂ for some exoplanets (Batygin et al. 2009; Csizmadia et al. 2018).

Independent of the method used to infer k₂, the possibility of interpreting its value in terms of the interior structure rests on the assumption that the planet is relaxed. Since gases do not have shear strength, gaseous planets should satisfy the hypothesis of hydrostatic equilibrium. All the odd zonal harmonics of the gravitational field of the gas giant Jupiter, if in perfect hydrostatic equilibrium, would be equal to zero. However, its gravitational field is north-south asymmetric (Iess et al. 2018), and this non-hydrostatic component is informative of the wind dynamics (Kaspi et al. 2018). The non-hydrostatic component also modifies the even zonal harmonics, whose interpretation is thus affected by the depth of the wind dynamics (Guillot et al. 2018). Still, hydrostatic models represent the starting point for investigating its interior (Wahl et al. 2017b; Guillot et al. 2018). In the foreseeable future, there is no possibility of inferring the high-order spherical harmonics of the gravitational fields of extrasolar planets, thus making hydrostatic models the starting (and likely ending) point for the modeling of their interior structure.

Terrestrial planets are objects whose main constituents are metals, rocks, and ice. These materials, unlike gases, have finite shear strengths and do not respond instantaneously to a perturbing potential. Their response is a function of the material properties, which in general are strongly affected by the temperature, and of the timescale and history of the perturbation (e.g., Padovan et al. 2014). Thus, the assumption of hydrostatic equilibrium for this class of objects requires an assessment of their orbitalhistory and global internal evolution.

6 Conclusions

We provided a semi-analytical method for computing the fluid Love numbers k_n of any planet from the knowledge of its density profile (Sect. 3). This parameter depends on the radial distribution of mass, thus providing additional information on the interior structure beyond the mass and radius. The method has been benchmarked against several results available in the literature (Sect. 3.6). The computation is very fast (Fig. 7), and the code is freely available². We used a few simple cases to investigate the basic dependencies of the Love numbers on the interior structure of a planet (Fig. 3). In Fig. 8, we illustrate the application of the code to the planet GJ 436b, whose observed mass and radius are compatible with an ocean-like interior or an interior dominated by rocks and metals (Adams et al. 2008), and in Fig. 9 we apply the code to a Jupiter-like hot Jupiter, whose core might be diluted due to erosion during the age of the solar system (Wahl et al. 2017b). These basic applications show that measuring k₂ would improve our understanding of the interior of extrasolar planets, but of course even a perfect knowledge of its value would not completely remove the degeneracy of planetary interior models (e.g., Kramm et al. 2011).

In applying the calculation of the fluid Love numbers presented here, the following points are worth noting:

• 1.

  The fluid Love numbers obtained from solving Eq. (4) describe the tidal response of a fluid nonrotating planet. The theory can also be used to describe the deformation of a planet in a state of rotation synchronous with its circular orbit, that is, a planet for which the rotational frequency ω in Eq. (21) corresponds to the orbital frequency $n=\sqrt[]{G{M}_{\text{S}}/d^3}$. In this case, the response corresponds to the linear combination of the perturbing potentials W and Z, as in Eq. (22).

• 2.

  Planets that rotate more rapidly than the tidal perturbation orbital frequency and in which the rotational distortion is much larger than the tidal distortion (e.g., Jupiter) cannot be treated with a linear theory for two reasons. First, the theory presented here can only treat small perturbations to a spherical planet (as a reference, there is about a 6.5% variation between the polar and equatorial radii of Jupiter, compared to 0.3% in the case of the Earth). Second, a dynamical theory of tides may be required to fully treat the response when the tidal bulge moves rapidly in the planet’s corotating frame, depending on the proximity of the tidal frequency to the planet’s resonant frequencies (see, e.g., Wahl et al. 2016, and references therein).

• 3.

  The response of planets with solid layers depends in general on the timescale of the perturbation and on the density, rigidity, and viscosity of their interiors (e.g., Moore & Schubert 2000; Padovan et al. 2014). In general, the theory applied here cannot be used for these objects. However, it is possible that for a given perturbation the planet responds as fluid, as in the case of the rotational flattening of the Earth (Lambeck 1980), and in such case an observed k₂ may thus be interpreted with the theory applied here. However, it is not known a priori whether an observed response for extrasolar planets corresponds to a fluid response, and an assessment of the orbital and thermal history of each object would be required (Sect. 5.2).

• 4.

  Given the previous points, the theory developed here is applicable to close-in, tidally locked gaseous exoplanets. From an observational point of view, these objects are the first objects for which estimates of the Love number k₂ will become available (e.g., Hellard et al. 2018).

Acknowledgements

This work was supported by the DFG within the Research Unit FOR 2440 “Matter Under Planetary Interior Conditions”. P.B. and N.T. acknowledge the support of the DFG priority program SPP1992 “Exploring the diversity of extrasolar planets” (TO 704/3-1). Figures were created with the software described in Hunter (2007). S.P. thanks W.B. Moore and D.J. Stevenson for informative discussions. The feedback of an anonymous referee improved the clarity of the paper.

References

All Tables

All Figures

	Fig. 1 Fluid Love number k as a function of the degree n for a homogeneous sphere.
In the text

	Fig. 2 Model with two constant density layers. Normalized density of the mantle vs. density of the core. Colors indicate the normalized radius of the core. The density of the core is plotted in logarithmic scale.
In the text

Fig. 3 Model of a planet with two constant density layers with M = 6.55 M⊕ and R = 2.68 R⊕. Values of k2 (first row), normalized moment of inertia MoI (second row), and maximum surface deformation (third row, assuming MS = 0.157 M⊙ and a = 0.0143 AU in Eq. (18)), for ρm∕ρ = 0.73 (left column), and for ρm∕ρ = 0.01 (right column). Colors indicate the normalized density of the core. The color bar is logarithmic.

In the text

	Fig. 4 Schematic structure of the vector of solution y and radius r for a planetary model. The center of the planet is at r = 0, the surface at r = R. Each homogeneous layer is indicated by the index corresponding to the outer boundary of the layer. Similar indexing applies to the non-dimensional radius s.
In the text

	Fig. 5 Set of normalized interior density profiles obtained using Gaussian curves with different standard deviations, as indicated in the legend. The yellow curve represents a homogeneous planet, and the blue curve shows a planet with a high concentration of mass in the interior.
In the text

	Fig. 6 For the models of Fig. 5, the left column shows values of the Love number k2 and of the normalized moment of inertia (MoI) as a function of the number of layers used to discretize the density profile. For a given number of layers, the right column shows the relative variation with respect to the value for 105 layers.
In the text

	Fig. 7 Speed of computation vs. number of layers in the calculation of the Love number k2 for a single-density profile. The curve is similar for n = 2,..., 8. A 2.5 GHz Intel Core i7 processor was used.
In the text

	Fig. 8 Top panel: three possible density profiles for GJ 436b with the same mass and radius, 21.4 M⊕ and 4.19 R⊕. A low-density hydrogen or hydrogen/helium envelope surrounds an Earth-like (brown), rocky (orange), or ocean-like (blue) interior. Bottom panel: corresponding values of the Love numbers kn.
In the text

	Fig. 9 Top panel: interior models of Jupiter with a compact high-density core (black) and a diluted and extended central region (red). Data are taken from Wahl et al. (2017b). Bottom panel: difference in the values of the Love numbers kn for the two models.
In the text