Issue 
A&A
Volume 620, December 2018



Article Number  A178  
Number of page(s)  10  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201834181  
Published online  12 December 2018 
Matrixpropagator approach to compute fluid Love numbers and applicability to extrasolar planets^{★}
^{1}
Institute of Planetary Research,
German Aerospace Center (DLR),
Berlin,
Germany
email: sebastiano.padovan@dlr.de
^{2}
Department of Astrophysics,
Technische Universität,
Berlin,
Germany
Received:
3
September
2018
Accepted:
19
October
2018
Context. The mass and radius of a planet directly provide its bulk density, which can be interpreted in terms of its overall composition. Any measure of the radial mass distribution provides a first step in constraining the interior structure. The fluid Love number k_{2} provides such a measure, and estimates of k_{2} for extrasolar planets are expected to be available in the coming years thanks to improved observational facilities and the everextending temporal baseline of extrasolar planet observations.
Aims. We derive a method for calculating the Love numbers k_{n} of any object given its density profile, which is routinely calculated from interior structure codes.
Methods. We used the matrixpropagator technique, a method frequently used in the geophysical community.
Results. We detail the calculation and apply it to the case of GJ 436b, a classical example of the degeneracy of massradius relationships, to illustrate how measurements of k_{2} can improve our understanding of the interior structure of extrasolar planets. We implemented the method in a code that is fast, freely available, and easy to combine with preexisting interior structure codes. While the linear approach presented here for the calculation of the Love numbers cannot treat the presence of nonlinear effects that may arise under certain dynamical conditions, it is applicable to closein gaseous extrasolar planets like hot Jupiters, likely the first targets for which k_{2} will be measured.
Key words: planets and satellites: interiors / planets and satellites: fundamental parameters / planets and satellites: gaseous planets / planets and satellites: terrestrial planets / planets and satellites: individual: GJ 436b
A copy of the code is available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/620/A178 and at https://bitbucket.org/sebastianopadovan/planetaryfluidlovenumbers/src/master/
© 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 selfgravitation 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 (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^{2}). 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 tideinducing perturber of mass M_{p} at a distance d, the tidegenerating potential W within the planet can be expanded in spherical harmonics as (e.g., Kaula 1966) (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, . At the surface, (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_{n}W_{n}, corresponding to (1 + k_{n})W_{n}∕g_{0}, with g_{0} 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_{0}, 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_{2} (or h_{2}). 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 matrixpropagator 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 matrixpropagator 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 welldefined 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_{2} 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) (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 (5)
where the perturbationinducing potential W generates the additional potential V^{p} within the body. The gravitational acceleration at the surface is g_{0}. The perturbed potential is proportional to the perturbationinducing potential and at the surface of the deformed body (r = R), the relation is (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 (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
where a plus (minus) indicates that the variable is evaluated right above (below) the corresponding radius. At the surface , the continuity condition is (Zharkov et al. 1985) (10)
where the function H is in relation with T and the fluid radial displacement Love number h_{n} through
2.2 Solution for a homogeneous sphere
If the body has constant density, ρ′ = 0, and Eq. (4) reduces to an EulerCauchy equation. Using the trial solution T(r) = r^{c} results in the characteristic equation (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_{1} = n and c_{2} = −(n + 1). The general solution is then the linear combination (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 (15)
The Love number k_{n} of order n is, from Eq. (7), (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. 

Open with DEXTER 
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 degree2 term of the tideinducing potential W in Eq. (2), the tidally induced radial deformation at the surface, , can be written as
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 ). 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 degree2 harmonic term as (e.g., Murray & Dermott 1999) (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 degree2 rotationally induced surface radial deformation is
Equation (21) shows that the larger and less massive the body and the faster it rotates, the larger its deformation (since ).
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, (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_{2}. 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 constantdensity mantle (gaseous envelope) overlies a constantdensity 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
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 closedform 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_{2}, the normalized moment of inertia MoI, and the maximum surface tidal deformation, which is reached at the substellar and antistellar points where P_{2} (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 highdensity core. From an inspection of the figure, we note the following:
 1.
The smaller (and thus, the denser) the core, the smaller both k_{2} 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_{2} that rapidly approaches 0 as the core decreases in size (specifically, for r_{C}∕R = 0.08, k_{2} = 6 × 10^{−3}).
 3.
When k_{2} approaches 0, this does not imply that the body is not deformed, since the deformation depends on h_{2} = k_{2} + 1. The limiting case of k_{2} = 0 (all the mass in a point core overlaid by a massless envelope of radius R) simply corresponds to h_{2} = 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_{2}, 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. 

Open with DEXTER 
Fig. 3 Model of a planet with two constant density layers with M = 6.55 M_{⊕} and R = 2.68 R_{⊕}. Values of k_{2} (first row), normalized moment of inertia MoI (second row), and maximum surface deformation (third row, assuming M_{S} = 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. 

Open with DEXTER 
3 Matrixpropagator approach to the calculation of k_{n}
With Eq. (4) can be recast as a system of two firstorder differential equations by introducing the function P = dT∕dr. Before writing the system in matrix form, the variables T and P are nondimensionalized as follows (e.g., Wolf 1994):
where y_{i} are nondimensional variables. The derivatives of T and P with respect to r are (27)
where s = r∕R is the nondimensional radius. The nondimensional system of equations is
The system can be written in matrix form as (30)
Assuming power solutions of the form , with a constant, the system (Eq. (30)) becomes (31)
which admits nontrivial solutions only if the determinant vanishes, that is, if (32)
The last equation corresponds to the characteristic Eq. (13), with solutions (33)
Thus, the general solution for y_{i} is the linear combination (34)
where the solution vectors are (35)
The general solution in matrix form is then (36)
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 (39)
where a plus (minus) indicates that the variable is evaluated right above (below) the boundary. The density difference 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): (40)
With reference to the interface at r = r_{i−1} in Fig. 4, the two solutions to be matched are
The interface matrix B(r_{i−1}) provides the connection between the two solutions: (43)
With the last three expressions, the vector of constants C^{(i)} can be expressed in terms of the vector of constants C^{(i−1)}: (44)
From Eqs. (41) and (44), the solution for y^{(i)} at r = r_{i} is (45)
and the vector of integration constant C^{(i)} 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:
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 nondimensional radius s. 

Open with DEXTER 
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), (48)
With the last two expressions (50)
3.4 Solution
With the results of Sects. 3.2 and 3.3, we can write the general solution at the surface as (53)
where the matrix M is completely defined by the interior model, that is, by the density as a function of radius. Explicitly, (54)
where y_{2}(0) = 0, according to the definition of P in Eq. (26). The surface boundary conditions, Eq. (10), in terms of y_{1} and y_{2} read (55)
whose inversion provides the boundary conditions of the problem. From the value of y_{1} (R), through Eqs. (25) and (7) the Love number for any n can be determined.
3.5 Sublayering
Given the discretized nature of the matrixpropagator 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 σ: (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_{2} as a function of the number of sublayers used to discretize the analytical profile expressed in Eq. (57), and the results are shown in Fig. 6. Both k_{2} and the MoI converge to a value that represents the “continuous” limit. This is the limit of an infinite number of sublayers with zero thickness. In the right column of Fig. 6 we plot the relative variation of the values of k_{2} and the MoI with respect to the limit value, which we take as the one corresponding to 10^{5} 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 sublayers in excess of about 10^{3} allows the retrieval of k_{2} and the MoI with a precision better than 1%. The required accuracy of the modeling can guide the choice of the number of sublayers to be used. Figure 7 shows that the computation of k_{2} for a model with 10^{5} sublayer 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 sublayers do depend on the assumed interior model, and the Gaussian curves used in this section represent an example that can be adapted on a casebycase 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 sublayers between 10^{3} and 10^{4} 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. 

Open with DEXTER 
Fig. 6 For the models of Fig. 5, the left column shows values of the Love number k_{2} 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 10^{5} layers. 

Open with DEXTER 
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^{3} sublayers 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 matrixpropagator method returns the values obtained with the solution of the full equation, provided an appropriate number of sublayers 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 thirddegree polynomial, which was then used to generate 10^{m} sublayers, 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_{2} 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 matrixpropagator 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_{2} that differs by less than 0.1% from the nonrotating 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).
Fluid Love number k_{2} and normalized moment of inertia of the Earth.
4 Applications
4.1 GJ 436b
The extrasolar planet GJ 436b is a classical example of the degeneracy^{1} of massradius 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 Earthlike, a rocky, or an oceanlike 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 matrixpropagator 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_{2} is controlled by the density difference between the center and the surface, which is a proxy for the degree of mass concentration. The Earthlike interior model, with a metallic core that reaches a density in excess of 30 g cm^{−3}, has k_{2} = 0.055, while the oceanlike interior model, with a central density of about 12 g cm^{−3}, has k_{2} = 0.160. The rocky interior model, with an intermediate central density, has a value of k_{2} = 0.082. Unlike the moment of inertia, the value of k_{2} 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_{2}.
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 Earthlike and rocky models, where the thickness of the envelope is similar.
The most likely parameter to be measured for exoplanets is k_{2} (Ragozzine & Wolf 2009; Hellard et al. 2018), and Fig. 8 shows that an accurate determination of k_{2} could more easily distinguish the oceanlike 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_{2} shows how its potential measurement would help in breaking the massradius degeneracy, at least partly.
4.2 Jupiterlike 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, highdensitycore 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_{2}, which is larger for the diluted core (k_{2} = 0.5378 versus k_{2} = 0.5287), given its smoother density profile. Even with the highquality 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_{2} (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_{2} 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 degree2 tidal perturbations of the Sun and the Moon induce a response in the degree4 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_{2} 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 propagatormatrix 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_{2} (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 k_{2} for a singledensity profile. The curve is similar for n = 2,..., 8. A 2.5 GHz Intel Core i7 processor was used. 

Open with DEXTER 
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 lowdensity hydrogen or hydrogen/helium envelope surrounds an Earthlike (brown), rocky (orange), or oceanlike (blue) interior. Bottom panel: corresponding values of the Love numbers k_{n}. 

Open with DEXTER 
Fig. 9 Top panel: interior models of Jupiter with a compact highdensity 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 k_{n} for the two models. 

Open with DEXTER 
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): (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_{2} 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_{2} 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_{2} for some exoplanets (Batygin et al. 2009; Csizmadia et al. 2018).
Independent of the method used to infer k_{2}, 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 northsouth asymmetric (Iess et al. 2018), and this nonhydrostatic component is informative of the wind dynamics (Kaspi et al. 2018). The nonhydrostatic 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 highorder 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 semianalytical 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^{2}. 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 oceanlike interior or an interior dominated by rocks and metals (Adams et al. 2008), and in Fig. 9 we apply the code to a Jupiterlike 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_{2} 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 . 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_{2} 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 closein, tidally locked gaseous exoplanets. From an observational point of view, these objects are the first objects for which estimates of the Love number k_{2} 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/31). 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
 Adams, E. R., Seager, S., & ElkinsTanton, L. 2008, ApJ, 673, 1160 [NASA ADS] [CrossRef] [Google Scholar]
 Alterman, Z., Jarosch, H., & Pekeris, C. L. 1959, Proc. R. Soc. Lond., Ser. A, 252, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Ash, M. E., Shapiro, I. I., & Smith, W. B. 1971, Science, 174, 551 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., Bodenheimer, P., & Laughlin, G. 2009, ApJ, 704, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Baumeister, P., MacKenzie, J., Tosi, N., & Godolt, M. 2018, Eur. Planet Sci. Conf., 678 [Google Scholar]
 Charbonneau, D., Berta, Z. K., Irwin, J., et al. 2009, Nature, 462, 891 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Correia, A. C. M. 2014, Astron. Astrophys., 570, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Csizmadia, Sz., Hellard, H., & Smith, A. 2018, Eur. Planet Sci. Conf., 858 [Google Scholar]
 Dorn, C., Khan, A., Heng, K., et al. 2015, A&A, 577, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dziewonski, A. M., & Anderson, D. L. 1981, Phys. Earth Planet. Inter., 25, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Gavrilov, S. V., Zharkov, V. N., & Leontev, V. V. 1975, Astron. Zh., 52, 1021 [NASA ADS] [Google Scholar]
 Gavrilov, S. V., Zharkov, V. N., & Leontev, V. V. 1976, Sov. Astron., 19, 618 [NASA ADS] [Google Scholar]
 Gilbert, F., & Backus, G. E. 1966, Geophysics, 31, 326 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Soter, S. 1966, Icarus, 5, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Guillot, T., Stevenson, D. J., Hubbard, W. B., & Saumon, D. 2004, in Jupiter. The Planet, Satellites and Magnetosphere, eds. F. Bagenal, T. E. Dowling, & W. B. McKinnon (Cambridge: Cambridge University Press), 35 [Google Scholar]
 Guillot, T., Miguel, Y., Militzer, B., et al. 2018, Nature, 555, 227 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Haskell, N. 1953, Bull. Seismol. Soc. Am., 43, 17 [Google Scholar]
 Hellard, H., Csizmadia, Sz., Padovan, S., et al. 2018, Eur. Planet Sci. Conf., 310 [Google Scholar]
 Helled, R., Anderson, J. D., Schubert, G., & Stevenson, D. J. 2011, Icarus, 216, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Hubbard, W. B. 2013, ApJ, 768, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Iess, L., Folkner, W. M., Durante, D., et al. 2018, Nature, 555, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Kaspi, Y., Galanti, E., Hubbard, W. B., et al. 2018, Nature, 555, 223 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kaula, W. M. 1966, Theory of Satellite Geodesy Applications of Satellites to Geodesy (Waltham, Mass: Dover) [Google Scholar]
 Kramm, U., Nettelmann, N., Redmer, R., & Stevenson, D. J. 2011, A&A, 528, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Lambeck, K. 1980, The Earth’s Variable Rotation: Geophysical Causes and Consequences, Cambridge Monographs on Mechanics (Cambridge: Cambridge University Press) [Google Scholar]
 Love, A. E. H. 1911, Some Problems of Geodynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Madhusudhan, N., Agúndez, M., Moses, J. I., & Hu, Y. 2016, Space Sci. Rev., 205, 285 [NASA ADS] [CrossRef] [Google Scholar]
 Mardling, R. A. 2007, MNRAS, 382, 1768 [NASA ADS] [Google Scholar]
 Moore, W. B., & Schubert, G. 2000, Icarus, 147, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Munk, W. H. & MacDonald, G. J. F. 1960, The Rotation of the Earth: A Geophysical Discussion (Cambridge: Cambridge University Press.) [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Nettelmann, N., Kramm, U., Redmer, R., & Neuhäuser, R. 2010, A&A, 523, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovan, S., Margot, J.L., Hauck, S. A., Moore, W. B., & Solomon, S. C. 2014, J. Geophys. Res. Planets, 119, 850 [NASA ADS] [CrossRef] [Google Scholar]
 Petit, G. & 2010, IERS Conventions 2010, Tech. Rep., Frankfurt am Main: Verlag des Bundesamts für Kartographie und Geodäsie [Google Scholar]
 Ragozzine, D., & Wolf, A. S. 2009, ApJ, 698, 1778 [NASA ADS] [CrossRef] [Google Scholar]
 Sterne, T. E. 1939, MNRAS, 99, 451 [NASA ADS] [Google Scholar]
 Stevenson, D. J. 2018, Asia Oceania Geosciences Society Annual Meeting, PS07A010 [Google Scholar]
 Thomson, W. T. 1950, J. Appl. Phys., 21, 89 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Trabant, C., Hutko, A. R., Bahavar, M., et al. 2012, Seismol. Res. Lett., 83, 846 [CrossRef] [Google Scholar]
 Trifonov, T., Kürster, M., Zechmeister, M., et al. 2018, A&A, 609, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Turner, J. D., Pearson, K. A., Biddle, L. I., et al. 2016, MNRAS, 459, 789 [NASA ADS] [CrossRef] [Google Scholar]
 Volterra, V. 1887, Mem. Soc. It. Sci., 3, 1 [Google Scholar]
 Wahl, S. M., Hubbard, W. B., & Militzer, B. 2016, ApJ, 831, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Wahl, S. M., Hubbard, W. B., & Militzer, B. 2017a, Icarus, 282, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017b, Geophys. Res. Lett., 44, 4649 [NASA ADS] [CrossRef] [Google Scholar]
 Wolf, D. 1994, Geophys. J. Int., 116, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Zharkov, V. N., Leontjev, V. V., & Kozenko, V. A. 1985, Icarus, 61, 92 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Fluid Love number k as a function of the degree n for a homogeneous sphere. 

Open with DEXTER  
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. 

Open with DEXTER  
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 k_{2} (first row), normalized moment of inertia MoI (second row), and maximum surface deformation (third row, assuming M_{S} = 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. 

Open with DEXTER  
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 nondimensional radius s. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 6 For the models of Fig. 5, the left column shows values of the Love number k_{2} 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 10^{5} layers. 

Open with DEXTER  
In the text 
Fig. 7 Speed of computation vs. number of layers in the calculation of the Love number k_{2} for a singledensity profile. The curve is similar for n = 2,..., 8. A 2.5 GHz Intel Core i7 processor was used. 

Open with DEXTER  
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 lowdensity hydrogen or hydrogen/helium envelope surrounds an Earthlike (brown), rocky (orange), or oceanlike (blue) interior. Bottom panel: corresponding values of the Love numbers k_{n}. 

Open with DEXTER  
In the text 
Fig. 9 Top panel: interior models of Jupiter with a compact highdensity 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 k_{n} for the two models. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.