# Absolute densities in exoplanetary systems: photodynamical modelling of Kepler-138

J. M. Almenara, R. F. Díaz, C. Dorn, X. Bonfils, S. Udry

► To cite this version:

J. M. Almenara, R. F. Díaz, C. Dorn, X. Bonfils, S. Udry. Absolute densities in exoplanetary systems: photodynamical modelling of Kepler-138. Monthly Notices of the Royal Astronomical Society, 2018, 478 (1), pp.460-486. 10.1093/mnras/sty1050. hal-02393195

HAL Id: hal-02393195

https://hal.science/hal-02393195

Submitted on 11 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Absolute densities in exoplanetary systems: photodynamical modelling of Kepler-138

J. M. Almenara,<fnref n="1" />* R. F. Díaz,<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> C. Dorn,<fnref n="4" />* X. Bonfils<fnref n="5" /> and S. Udry<fnref n="1" />

Accepted 2018 April 20. Received 2018 April 17; in original form 2017 December 20

## ABSTRACT

In favourable conditions, the density of transiting planets in multiple systems can be determined from photometry data alone. Dynamical information can be extracted from light curves, providing modelling is done self-consistently, i.e. using a photodynamical model, which simulates the individual photometric observations instead of the more generally used transit times. We apply this methodology to the Kepler-138 planetary system. The derived planetary bulk densities are a factor of 2 more precise than previous determinations, and we find a discrepancy in the stellar bulk density with respect to a previous study. This leads, in turn, to a discrepancy in the determination of masses and radii of the star and the planets. In particular, we find that interior planet, Kepler-138b, has a size in between Mars and the Earth. Given our mass and density estimates, we characterize the planetary interiors using a generalized Bayesian inference model. This model allows us to quantify for interior degeneracy and calculate confidence regions of interior parameters such as thicknesses of the core, the mantle, and ocean and gas layers. We find that Kepler-138b and Kepler-138d have significantly thick volatile layers and that the gas layer of Kepler-138b is likely enriched. On the other hand, Kepler-138c can be purely rocky.

**Key words:** planets and satellites: interiors – stars: fundamental parameters – planetary systems.

## 1 INTRODUCTION

The photometric eclipses that occur when an extrasolar planet moves across the face of its host star provide information on both the orbiting object and its parent star. For example, the size of the planet relative to the star can be measured from the depth of the eclipse, and an approximate bulk stellar density can also be inferred through the use of Kepler's third law of planetary motion (Seager & Mallén-Ornelas 2003). The radius ratio $R_p/R_*$ obtained from a transit light curve is the first step towards a measurement of the planetary bulk density, which contains a wealth of information on the planet composition and also allows us to study the diversity of these objects. In many cases, the radius ratio is coupled with a dynamical measurement of the planet-to-star mass and with stellar evolutionary models and spectroscopic analyses to achieve a measurement of the bulk planetary density (e.g. Pepe et al. 2013; Christiansen et al. 2017).

To date, the mass measurement has been secured for hundreds of transiting giant planets using precise Doppler velocimetry measurements. The majority of these planets were first discovered by ground-based wide-field surveys such as SuperWASP and HATnet, and then followed-up with high-precision radial velocity (RV) facilities to secure their nature and measure their masses. On the contrary, only a handful of transiting planets smaller than Neptune have had their masses characterized by precise RV measurements, and this at the price of long observing campaigns requiring a large amount of telescope time, usually spanning many seasons (e.g. Queloz et al. 2009; Batalha et al. 2011; Pepe et al. 2013). The *Kepler* mission (Borucki et al. 2010) discovered thousands of Neptune-sized or smaller candidates, but their host stars are faint, hindering the high-precision RV measurements needed to measure the small reflex motions induced by these planets on their host stars. For the few objects whose masses were determined by RVs, the precision in the planetary densities are usually better than 20 per cent, but the determination is still dependent on the accuracy of stellar models, whose systematic errors are not yet fully known or understood.

<footnote n="1">¹Observatoire de Genève, Département d'Astronomie, Université de Genève, Chemin des Maillettes 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="2">²Universidad de Buenos Aires, Facultad de Ciencias Exactas y Naturales, Buenos Aires, 1428, Argentina</footnote>

<footnote n="3">³CONICET - Universidad de Buenos Aires, Instituto de Astronomía y Física del Espacio (IAFE), Buenos Aires, 1428, Argentina</footnote>

<footnote n="4">⁴University of Zürich, Institut of computational sciences, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland</footnote>

<footnote n="5">⁵University Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote>* E-mail: jose.almenara@unige.ch (JMA) rodrigo@iafe.uba.ar (RFD) cdorn@physik.uzh.ch (CD)</footnote>
In multiple transiting systems, the gravitational interactions be-
tween the planets can be exploited to obtain mass ratios with re-
spect to the stellar mass without resorting to RVs. These interactions
make the objects' trajectories deviate from Keplerian motion, which
in turn leads to the presence of transit timing variations (TTVs),
i.e. departures from strict periodicity in the transit times. Modelling
TTVs has been in recent years a very fruitful alternative technique
to reach the planetary physical properties (e.g. Hadden & Lithwick
2017). However, these analyses must also rely on atmospheric and
evolutionary stellar models to measure the planetary masses and
radii. Furthermore, the bulk densities obtained by combining the
planetary masses and radii measured in this way are determined in
most cases with large uncertainties.<fnref n="1" />

One such system is Kepler-138, composed of three planets orbiting an M1-type star, with periods 10.3, 13.8, and 23.1 d (Kepler-138b, c, and d, respectively). Of particular interest in this system, Kepler-138b is a Mars-sized planet. Planets b and c are close to the 4:3 first-order mean-motion resonance, whereas planets c and d are close to second-order resonance 5:3. The system is dynamically packed with period ratios smaller than two. All planets exhibit timing variations (Ford et al. 2011; Mazeh et al. 2013; Kipping et al. 2014; Jontof-Hutter et al. 2015). The system was studied by Kipping et al. (2014), who confirmed the planetary nature of the two external planets, and measured their masses by modelling the dynamic effects between them. Despite having similar radii, the masses of Kepler-138c and Kepler-138d differ significantly, with planet d's composition probably being dominated by a gaseous envelope, while planet c is probably rocky. The mass of the interior Mars-size planet was first measured by Jontof-Hutter et al. (2015, hereafter JH15) from the TTVs exhibited by the planets.<fnref n="2" /> This is a remarkable result as the RV signal expected for this planet is few cm s⁻<fnref n="1" />, unattainable by current facilities, even for bright stars. However, planetary masses and radii are determined to a precision of 62–68 and 6 per cent, respectively, and therefore, the derived densities are reported with a precision between 62 and 65 per cent. This hindered distinguishing between various possible compositions for the Kepler-138b.

In this article, we obtain model-independent bulk densities for
Kepler-138 and its three planets by fully exploiting the information
available in the light curve of Kepler-138. This is possible thanks to
the presence of dynamical interactions in the system and the nature
of the gravitational force. The measured density gives us corre-
lated samples from the posterior distributions of masses and radii,
which are used as input for a detailed interior characterization of
the planets. We use the generalized model of Dorn, Hinkel & Ven-
turini (2017) to calculate confidence regions of interior parameters,
while accounting for the generally large degeneracy of interior pa-
rameters. The range of possible interior scenarios are key to better
understand the planet's possible formation and evolution.

In Section 2, we expose the foundation of the method to derive
absolute densities. Then we show the application of this method-
ology to Kepler-138; in Section 3, we describe the data used; in
Section 4, we detail the photodynamical modelling; and in Section

5, we present the results including a detailed interior characteriza-
tion of the planets. Finally, we discuss the results of our work in
Section 6 and present our conclusions in Section 7.

## 2 THE METHOD

It is usually stated that to obtain a precise estimation of the density of planetary bodies one needs to measure their masses and radii. Although this is true for most single planets,<fnref n="3" /> transiting multi-planet systems offer the possibility to measure the bulk density of the objects in the system without measuring masses and radii independently (Lissauer et al. 2011; Carter et al. 2012; Sanchis-Ojeda et al. 2012; Huber et al. 2013; Almenara et al. 2015).

This can be seen by simply doing the ratios of the planet to star
mean densities ($\rho_p$, and $\rho_*$, respectively), assuming sphericity for
all objects:

$$
\rho_p = \rho_* \left( \frac{M_p}{M_*} \right) \left( \frac{R_p}{R_*} \right)^{-3},
$$

where *M* and *R* represent masses and radii, and the subindexes
*p* and *⋆* refer to the planet and the star, respectively (Southworth
2009; Jontof-Hutter et al. 2014; Weiss & Marcy 2014). The radius
ratio *R*<sub>*p*</sub>*R*<sub>*⋆*</sub> is constrained by the transit shape, and the mass ratio *M*<sub>*p*</sub>*M*<sub>*⋆*</sub>
is constrained by the gravitational interactions. It is well known
that the mean density of a star can be estimated from the planetary
transit shape using the Kepler's third law, if the orbital eccentricity
and argument of pericentre are known (Seager & Mallén-Ornelas
2003; Kipping 2014a). Constraints on the orbital geometry and ori-
entation can be obtained, for example, from RV measurements or
directly from the light-curve modelling of the gravitational interac-
tions between the planets.

Therefore, if the gravitational interactions between the planets
of a multitransiting system are detected, the mean densities can
be inferred from the light-curve data alone, without relying on
stellar models, apart from the assumption of sphericity, and a limb-
darkening law. Under these assumptions, the determination of the
bulk densities can still be biased by undetected blended companions
and unaccounted for stellar spots (Kipping 2014a).

Interestingly, bulk density is the only non-dimensionless phys-
ical magnitude obtainable from the analysis of a light curve. The
ultimate reason for this is that the gravitational constant *G* has units
of inverse density over a squared time. The relative flux being a di-
mensionless quantity, the light curve only provides a mapping of the
motion of the bodies in time. In other words, the light-curve model,
including the gravitational interactions, is invariant to scaling the
lengths by a constant factor and the masses by the same factor at
cubic exponent. This is called the Newtonian *MR*⁻<fnref n="3" /> degeneracy. If
the light-time or relativistic effects are measured<fnref n="4" /> or radial veloci-
ties of the star are obtained, this invariance is broken and absolute
masses and radii are measurable (Agol et al. 2005; Almenara et al.
2015).

The photodynamical analysis consistently models the light curve
and permits deriving absolute bulk densities, taking into account all
the correlations between parameters. The photodynamical analysis

<footnote n="1">¹Notable exceptions are Kepler-36 (Carter et al. 2012) and Kepler-56 (Huber et al. 2013) for which analyses very similar to the one presented here were performed.</footnote>

<footnote n="2">²Note that the planet labels are not the same between Kipping et al. (2014) and JH15. Planet b in Kipping et al. (2014) (KOI-314b) is Kepler-138c in JH15, and KOI-314c is Kepler-138d. Here, we follow the notation of JH15: Kepler-138b (planet b) = KOI-314.03; Kepler-138c (planet c) = KOI-314.01 = KOI-314b; Kepler-138d (planet d) = KOI-314.02 = KOI-314c.</footnote>

<footnote n="3">³For short-period planets, if the ellipsoidal modulation is detected, it is possible to determine $M_p/M_*$.</footnote>

<footnote n="4">⁴Because *Kepler* measurements do not permit detecting the light-time effect in single-star systems (Ragozzine & Wolf 2009), nor the relativistic effects in Kepler-138, we hereafter avoid making this clarification each time we refer to the light-curve data.</footnote>
couples an *N*-body simulation that determines the movement of the bodies in time with a photometric model that computes the relative loss of light during the planetary transits. The output is a light curve that can be compared directly with the observed one. The assumptions and a more detailed description are given in Almenara et al. (2015).

## 3 DATA

*Kepler* observed Kepler-138 during the 4 yr of its prime mission. The *Kepler* light curves of all Quarters (Q0–Q17) were retrieved from the Mikulski Archive for Space Telescopes (MAST).<fnref n="5" /> We preferred short-cadence data (about one point per minute; quarters Q6–Q17) whenever available. We used the simple aperture photometry (SAP) light curve, which we corrected for the flux contamination<fnref n="6" /> (between 0 and 2 per cent depending on the quarter) using the value estimated by the *Kepler* team.<fnref n="7" /> Only the data spanning three transit durations around each transit were modelled, after normalization with a parabola, accounting for the variability of the out-of-transit light curve (Czesla et al. 2009), although this is a minor correction in this case (see Section 6.2.1). The observed transits are presented in Figs A1–A3, for Kepler-138b, Kepler-138c, and Kepler-138d, respectively.

## 4 PHOTODYNAMICAL MODELLING

Our model is parametrized by the stellar density and flux distribution across the disc, the planetary-to-star mass and radius ratios, and a set of orbital parameters at a given reference time $t_{\text{ref}}$. As the *N*-body integrator relies on the bodies' masses and their positions and velocities at $t_{\text{ref}}$, and we are dealing only with light-curve data, the input stellar density has to be converted to a mass using an arbitrary fixed value for the radius. This radius value is also employed to normalize the input semimajor axis, $a/R_*$.$^8$ The same value is then used to normalize the output positions of the integrator to use as input for the photometric model. Finally, the light is modelled using the projected centre-to-centre distance between star and planets. For the model to be valid and depend on the masses only through the bulk densities, light-time and relativistic effects must be negligible. For Kepler-138, the amplitude of the light-time effect is estimated at 3 ms, whereas the periastron advance due to relativistic effects is less than 1 ms for individual transits (Heyl & Gladman 2007). For comparison, the absolute accuracy of *Kepler* times is 7 s to 97.5 per cent confidence (*Kepler Data Characteristics Handbook*).

In this analysis, we use the REBOUND code (Rein & Liu 2012) with the WHFast integrator (Rein & Tamayo 2015), with a temporal resolution of 0.01 d (864 s). The positions of the objects at the times of the short-cadence light curve were interpolated between the integration points using a cubic spline. To model the long-cadence data (about one point every 30 min; quarters Q0–Q5), we obtained the position of the system bodies at 30 evenly spaced points around each observation date.

The light curve was computed using the analytic description of Mandel & Agol (2002), with the implicit assumption of spherical shape for the star and planets, which are assumed to be non-emitting bodies. The flux distribution across the stellar disc was modelled using a quadratic limb-darkening law (Manduca, Bell & Gustafsson 1977) with the Kipping (2013a) parametrization to consider only physical values. The model also includes additional white noise terms for the *Kepler* long and short-cadence data, and a free normalization factor for each data set, corresponding to the out-of-transit flux. The model at the times of the short-cadence light curve is compared directly to the observations. The model for the long-cadence light curve is obtained by binning the 30×oversampled model light curve back to the cadence of the observations, to avoid an artificial deformation of the signal (Kipping 2010). The resulting maximum-a-posteriori (MAP) model light curves are presented in Figs A1–A3.

To quantify the effect of the finite temporal resolution of the *N*-body integration on the model photometric error of the model, we calculated a light curve using a resolution of 8.64 s, i.e. a hundred times finer than the resolution used in the analysis, and oversampling of 10 for the short-cadence light curve. We chose a set of model parameters corresponding to the MAP estimate obtained in Section 5. This light curve was considered as the ground-truth data to which we compared a light-curve model obtained with the nominal *N*-body temporal resolution of 864 s and no oversampling. The maximum difference was found to be 5 ppm, which is much smaller than the typical uncertainty of the photometric measurements.

Overall, the model has a total of 30 free parameters: five orbital elements and a mass and radius ratio per planet, two relative longitudes of the ascending nodes,<fnref n="8" /> two limb-darkening coefficients, the amplitude of the two additional noise terms described above, and the out-of-transit flux levels for short and long-cadence data.

In an attempt to reduce the correlations between the orbital parameters, which hinder sampling efficiently from the posterior distribution (see below), we chose the parameters listed in Table 1 as jump parameters. To sample from the posterior distributions of the parameter models, we used the emcee algorithm (Goodman & Weare 2010; Foreman-Mackey et al. 2013).

The model is innately multimodal, as different configurations of the orbital inclinations produce similar, although not identical, results. As most Markov chain Monte Carlo (MCMC) algorithms are unable to correctly sample multimodal distributions (Foreman-Mackey et al. 2013), we used four different sets of chains in parameter space, each started at a different inclination configuration<fnref n="9" /> (Fig. 1).

Initializing the MCMC algorithm in an adequate point in parameter space is crucial to reduce burn-in time. This is particularly important for the photodynamical modelling. Model transits associated with poor parameter values may fall between the observed transits. As we only keep the data around observed transits, the algorithm can be very inefficient in going out of this region of pa-

<footnote n="5">⁵http://archive.stsci.edu. We used the data release 25.</footnote>

<footnote n="6">⁶We used the keyword 'CROWDSAP' on the header of the fits files containing the light curves. The correction was done as specified on *Kepler Data Processing Handbook* (KSCI-19081-001), Section 7.3.3.</footnote>

<footnote n="7">⁷Wang et al. (2015) found no companion in adaptive optics images of this star, reducing the probability that there is a contaminant star in the *Kepler* aperture whose contamination was not taken into account by the Kepler team. Besides, this contamination estimate seem to be self-consistent, as no differences in transit depth where detected between seasons.</footnote>

<footnote n="8">⁸For a spherical star, the model is independent of the individual longitudes of the ascending node. Here, we chose to fix $\Omega_c$ at $t_{\text{ref}}$ to 180°, and fit $\Delta\Omega_{ic} = \Omega_i - \Omega_c$ where the subindex $i = \{b, d\}$ refers to planets Kepler-138b and Kepler-138d.</footnote>

<footnote n="9">⁹The number of possible hemisphere configurations are $\frac{2^N}{2}$, with *N* the number of planets in the system. The fraction 2 comes from the symmetry for spherical stars and planets that permits to fix the transiting hemisphere for one of the planets. Here, we limit the possible inclinations of Kepler-138d below 90°.</footnote>
<table><caption>**Table 1.** Inferred system parameters from the light curve only: MAP, posterior median, and 68.3 per cent CI for the photodynamical analysis. The astrocentric orbital elements are given for the reference time $t_{\text{ref}} = 2454955$ BJD$_{\text{TDB}}$ (the orbital parameters can change significatively even during the time span of Kepler observations, see for example the eccentricity of Kepler-138b in Fig. A6).</caption><thead><tr><th>Parameter</th><th></th><th>MAP</th><th>Median and 68.3 per cent CI</th></tr></thead><tbody><tr><td colspan="4"><strong>Kepler-138</strong></td></tr><tr><td>Stellar mean density, ρ<sub>⋆</sub><sup>a</sup></td><td>(g cm<sup>-3</sup>)</td><td>3.735</td><td>3.92<sup>+0.81</sup><sub>-0.66</sub></td></tr><tr><td>q<sub>1</sub><sup>a, b</sup></td><td></td><td>0.3890</td><td>0.425<sup>+0.084</sup><sub>-0.062</sub></td></tr><tr><td>q<sub>2</sub><sup>a, b</sup></td><td></td><td>0.313</td><td>0.24<sup>+0.20</sup><sub>-0.15</sub></td></tr><tr><td>Linear limb darkening, u<sub>a</sub></td><td></td><td>0.391</td><td>0.32<sup>+0.23</sup><sub>-0.19</sub></td></tr><tr><td>Quadratic darkening, u<sub>b</sub></td><td></td><td>0.233</td><td>0.34±0.26</td></tr><tr><td colspan="4"><strong>Kepler-138b</strong></td></tr><tr><td>Scaled semimajor axis, a/R<sub>⋆</sub></td><td></td><td>27.59</td><td>28.0±1.8</td></tr><tr><td>Eccentricity, e</td><td></td><td>0.00469</td><td>0.0060<sup>+0.0082</sup><sub>-0.0044</sub></td></tr><tr><td>Inclination, i<sup>a</sup></td><td>(°)</td><td>88.43</td><td>88.45<sup>+0.17</sup><sub>-0.15</sub>, 91.51<sup>+0.21</sup><sub>-0.18</sub> c</td></tr><tr><td>Argument of pericentre, ω</td><td>(°)</td><td>65.4</td><td>23<sup>+47</sup><sub>-65</sub></td></tr><tr><td>Longitude of the ascending node, Ω<sup>a</sup></td><td>(°)</td><td>179.45</td><td>180.1±1.1</td></tr><tr><td>Mean anomaly, M<sub>0</sub></td><td>(°)</td><td>341.2</td><td>23<sup>+65</sup><sub>-46</sub></td></tr><tr><td>√e cos ω<sup>a</sup></td><td></td><td>0.0286</td><td>0.052±0.048</td></tr><tr><td>√e sin ω<sup>a</sup></td><td></td><td>0.0623</td><td>0.022<sup>+0.048</sup><sub>-0.060</sub></td></tr><tr><td>T'<sub>0</sub><sup>a</sup></td><td>(BJD<sub>TDB</sub>)</td><td>2454956.23845</td><td>2454956.2380±0.0034</td></tr><tr><td>P'<sub>a</sub></td><td>(d)</td><td>10.313737</td><td>10.31368<sup>+0.00018</sup><sub>-0.00022</sub></td></tr><tr><td>Radius ratio, R<sub>p</sub>/R<sub>⋆</sub><sup>a</sup></td><td></td><td>0.011207</td><td>0.01104±0.00031</td></tr><tr><td>Mass ratio, M<sub>p</sub>/M<sub>⋆</sub><sup>a</sup></td><td></td><td>1.187 × 10<sup>-6</sup></td><td>(1.04 ± 0.24) × 10<sup>-6</sup></td></tr><tr><td>Planet mean density, ρ<sub>p</sub></td><td>(g cm<sup>-3</sup>)</td><td>3.149</td><td>3.02<sup>+1.0</sup><sub>-0.87</sub></td></tr><tr><td colspan="4"><strong>Kepler-138c</strong></td></tr><tr><td>Scaled semimajor axis, a/R<sub>⋆</sub></td><td></td><td>33.47</td><td>34.0±2.2</td></tr><tr><td>Eccentricity, e</td><td></td><td>0.00884</td><td>0.0110<sup>+0.0069</sup><sub>-0.0036</sub></td></tr><tr><td>Inclination, i<sup>a</sup></td><td>(°)</td><td>88.78</td><td>88.81<sup>+0.15</sup><sub>-0.13</sub>, 91.16<sup>+0.17</sup><sub>-0.15</sub> c</td></tr><tr><td>Argument of pericentre, ω</td><td>(°)</td><td>353.9</td><td>341<sup>+27</sup><sub>-18</sub></td></tr><tr><td>Longitude of the ascending node, Ω</td><td>(°)</td><td>180<sup>d</sup></td><td></td></tr><tr><td>Mean anomaly, M<sub>0</sub></td><td>(°)</td><td>76.09</td><td>89<sup>+18</sup><sub>-27</sub></td></tr><tr><td>√e cos ω<sup>a</sup></td><td></td><td>0.0935</td><td>0.096<sup>+0.030</sup><sub>-0.023</sub></td></tr><tr><td>√e sin ω<sup>a</sup></td><td></td><td>-0.0100</td><td>-0.034<sup>+0.049</sup><sub>-0.027</sub></td></tr><tr><td>T'<sub>0</sub><sup>a</sup></td><td>(BJD<sub>TDB</sub>)</td><td>2454955.727995</td><td>2454955.72831±0.00060</td></tr><tr><td>P'<sub>a</sub></td><td>(d)</td><td>13.7815991</td><td>13.781564<sup>+8.2 × 10<sup>-5</sup></sup><sub>-9.5 × 10<sup>-5</sup></sub></td></tr><tr><td>Radius ratio, R<sub>p</sub>/R<sub>⋆</sub><sup>a</sup></td><td></td><td>0.026429</td><td>0.02628<sup>+0.00048</sup><sub>-0.00043</sub></td></tr><tr><td>Mass ratio, M<sub>p</sub>/M<sub>⋆</sub><sup>a</sup></td><td></td><td>3.069 × 10<sup>-5</sup></td><td>(2.90<sup>+0.45</sup><sub>-0.60</sub> × 10<sup>-5</sup>)</td></tr><tr><td>Planet mean density, ρ<sub>p</sub></td><td>(g cm<sup>-3</sup>)</td><td>6.21</td><td>6.1<sup>+1.9</sup><sub>-1.5</sub></td></tr><tr><td colspan="4"><strong>Kepler-138d</strong></td></tr><tr><td>Scaled semimajor axis, a/R<sub>⋆</sub></td><td></td><td>47.22</td><td>48.0±3.1</td></tr><tr><td>Eccentricity, e</td><td></td><td>0.02403</td><td>0.0270±0.0050</td></tr><tr><td>Inclination, i<sup>a</sup></td><td>(°)</td><td>88.9281</td><td>88.952±0.082</td></tr><tr><td>Argument of pericentre, ω</td><td>(°)</td><td>241.36</td><td>246.1<sup>+10</sup><sub>-5.6</sub></td></tr><tr><td>Longitude of the ascending node, Ω<sup>a</sup></td><td>(°)</td><td>179.997</td><td>180.21±0.42</td></tr><tr><td>Mean anomaly, M<sub>0</sub></td><td>(°)</td><td>165.98</td><td>161.2<sup>+6.1</sup><sub>-11</sub></td></tr><tr><td>√e cos ω<sup>a</sup></td><td></td><td>-0.0743</td><td>-0.066<sup>+0.029</sup><sub>-0.018</sub></td></tr><tr><td>√e sin ω<sup>a</sup></td><td></td><td>-0.1360</td><td>-0.150<sup>+0.014</sup><sub>-0.012</sub></td></tr><tr><td>T'<sub>0</sub><sup>a</sup></td><td>(BJD<sub>TDB</sub>)</td><td>2454957.822471</td><td>2454957.82216±0.00073</td></tr><tr><td>P'<sub>a</sub></td><td>(d)</td><td>23.093353</td><td>23.09302<sup>+0.00069</sup><sub>-0.00092</sub></td></tr><tr><td>Radius ratio, R<sub>p</sub>/R<sub>⋆</sub><sup>a</sup></td><td></td><td>0.026566</td><td>0.02643±0.00052</td></tr><tr><td>Mass ratio, M<sub>p</sub>/M<sub>⋆</sub><sup>a</sup></td><td></td><td>7.04 × 10<sup>-6</sup></td><td>(6.5<sup>+1.3</sup><sub>-1.5</sub>) × 10<sup>-6</sup></td></tr><tr><td>Planet mean density, ρ<sub>p</sub></td><td>(g cm<sup>-3</sup>)</td><td>1.403</td><td>1.36<sup>+0.44</sup><sub>-0.35</sub></td></tr><tr><td colspan="4"><strong>Data</strong></td></tr><tr><td>Kepler long-cadence jitter<sup>a</sup></td><td></td><td>1.0363</td><td>1.042±0.021</td></tr><tr><td>Kepler short-cadence jitter<sup>a</sup></td><td></td><td>1.00952</td><td>1.0089±0.0023</td></tr></tbody></table>

$$T'_{0} \equiv t_{\text{ref}} - \frac{p'}{2\pi} (M_{0} - E + e \sin E) \text{ with } E = 2 \arctan \left\{ \sqrt{\frac{1-e}{1+e}} \tan \left[ \frac{1}{2} \left( \frac{\pi}{2} - \omega \right) \right] \right\}, P' \equiv \sqrt{\frac{3\pi}{G\rho_{\star}} \left(\frac{a}{R_{\star}}\right)^{3}}$$

CODATA 2014: $G = 6.67408 \times 10^{-11}$ m$^3$ kg$^{-1}$ s$^{-2}$

IAU 2012: au = 149 597 870 700 m. IAU 2015: $(GM)_{\odot} = 1.327 124.4 \times 10^{20}$ m$^3$ s$^{-2}$

$k^2 = (GM)_{\odot} (86 400 \text{ s})^2 \text{ au}^{-3}$

<footnote><sup>a</sup> emcee jump parameter. <sup>b</sup> Kipping (2013b) parametrization for the limb-darkening coefficients to consider only physical values. <sup>c</sup> For configurations A and D respectively. <sup>d</sup> Fixed at $t_{\text{ref}}$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Inclination configurations (A, B, C, and D) for the planet (b, c, and d) orbits. The grey horizontal line represents the line of sight. The relative sizes of the star and the planets are to scale. The symbols of the planets are plotted at the measured impact parameters.</figcaption></figure>

parameter space. To insure a correct initialization, we performed a preliminary analysis using only the transit times from JH15. In this case, we only have information on the times of inferior conjunction of each planet. All the information on the stellar density and relative planetary sizes is lost. Therefore, the model is no longer dependent on the planetary radius ratios nor on the stellar density. The normalized semimajor axes for each planet, $a/R_*$, are replaced by $a^3/M_*$, which can be seen as the inverse mass density inside the planet orbit. For our analysis of the JH15 transit times, we assumed coplanar orbits, a condition we later relaxed for the final analysis. Samples from posterior distribution were used as starting points for the MCMC algorithm for the (complete) photodynamical model.

We used non-informative uniform priors for all emcee jump parameters. The walkers were started on each inclination configuration (Fig. 1) independently, and we noticed walkers in configurations B and C moved away to configurations A and D, respectively. This means that the latter configurations are preferred by the data. We finally ran 100 emcee walkers for $1.2 \times 10^6$ steps started on configurations A and D only. The last 200 000 steps of each configurations were used for the final inference assuming equal probability for each mode.<fnref n="10" /> We verified that the obtained Markov chains did not show signs of non-convergence by analysing the marginal samples for each parameters over different parts of the chain. The individual walkers have long autocorrelation lengths,<fnref n="11" /> ranging from 9000 to 150000 steps, depending on the model parameter and walker. The average over walkers is between 38000 and 64000 steps. This implies that on average three independent samples are produced by each walker, which means that the effective number of independent samples used for the inference is 300.

## 5 RESULTS

The sample obtained with the MCMC algorithm (Fig. A4) was used to measure the MAP of the model parameters, their median, and credible intervals (CI). In this case, as the prior distribution is uniform, the MAP is equal to the maximum likelihood. The results are listed in Table 1. The MAP model is shown in Figs A1–A3. Figs 2 and 3 present two views of the system and the inferred orbits of the planets. In Fig. 4, the posterior TTVs measured with the

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Observer view of the Kepler-138 system, projected trajectories of 1000 random MCMC steps for Kepler-138b (violet), Kepler-138c (green), and Kepler-138d (blue). The MAP model estimate is shown in black. Star and planet sizes are to scale.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top view of the Kepler-138 system. Orbit trajectories of 1000 random MCMC steps for Kepler-138b (violet), Kepler-138c (green), and Kepler-138d (blue). The MAP model is plotted as a black curve. Distances are to scale, but the size of the planets is multiplied by a factor of 100, and are shown at the position of $t_{\text{ref}}$. The orbit of the star is also represented but is not discernable at this scale. Positive Z-axis points towards the observer; planets move clockwise.</figcaption></figure>

photodynamical model are compared to the individual timings of JH15.

In configurations A and D Kepler-138b and Kepler-138c transit, the star on the same hemisphere (Fig. 2). This is reasonable, as these planets have a period ratio of 1.33 (Fig. 3). In both configurations, planets b and c have low mutual inclinations (Fig. A5).

The timing of – at least – the last two transits of Kepler-138d are not correctly modelled (Figs 4 and A3; see also panel c of fig. 1 in JH15). This may indicate the presence of an additional

<footnote n="10">¹⁰The posterior is available at https://zenodo.org/record/1227263.</footnote>

<footnote n="11">¹¹We considered the autocorrelation length as the lag needed for the autocorrelation to fall to $1/e$.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Posterior TTVs of Kepler-138b (top panel, violet), Kepler-138c (middle panel, green), and Kepler-138d (bottom panel, blue) from the photodynamical modelling. For comparison the TTVs of JH15 measured on individual transits are shown as empty circles with grey errorbars.</figcaption></figure>

non-transiting companion to the system. Because only the transits of Kepler-138d seem to be affected, we assume this would be an outer companion, but this is not necessarily true. We discuss the implication such a planet would have on our results in Section 6.3.1.

The planetary densities are measured with a precision of about 30 per cent. This is around two times more precise than the JH15 determination and independent of stellar models. As previously noted by Kipping et al. (2014), planets c and d have similar radii ($\sim 1.7$ R⊕) but very different densities, $\rho_c = 6.1^{+1.9}_{-1.5} \text{ g cm}^{-3}$ and $\rho_d = 1.36^{+0.44}_{-0.35} \text{ g cm}^{-3}$. This is unlike the Kepler-36 system (Carter et al. 2012), in which the planets have very different densities but one is more than twice the size of the other. We compare our results to previous literature in Section 6.1.

## 5.1 Masses and radii

Having only the bulk densities is not enough to constrain the nature of the planets (Hatzes & Rauer 2015). To obtain the planet masses and radii the Newtonian $MR^{-3}$ degeneracy (Section 2) must be broken. In this case, the expected RV amplitudes of the planets ranges from a few cm s⁻¹ to 2 m s⁻¹. These amplitudes are not easily detectable on such a faint star using current instrumentation. Therefore, we must rely on information on the stellar host.

On the other hand, it is well known that evolution models of M-type stars have issues to reproduce the observed masses and radii (Berger et al. 2006; Morales et al. 2009; Torres 2010; Boyajian et al. 2012; Torres 2013). Even invoking the effect of cool starspots is insufficient to explain most discrepancies (López-Morales 2007; Morales, Ribas & Jordi 2008; Morales et al. 2010). However, luminosity depends on the rate of nuclear reactions in the stellar core and is less affected by stellar activity. Therefore, to obtain the mass of Kepler-138, we used the empirical mass–luminosity relation (MLR) from Peretti et al. (in preparation), which is based on the work of Delfosse et al. (2000). We decided to use the $K_s$-band magnitude. Near-infrared (NIR) bands exhibit smaller dispersion than visible bands, which is interpreted as a lesser effect of stellar activity and metallicity (Bonfils et al. 2005). In particular, the $K_s$ is the one exhibiting the smallest dispersion among the NIR bands. The absolute magnitude in the $K_s$ band was determined spectroscopically by Pineda, Bottom & Johnson (2013), $M_{K_s} = 5.39 \pm 0.25$, and Térien et al. (2015a), $M_{K_s} = 5.42 \pm 0.18$, using empirical relations calibrated on nearby stars for which the parallax is known. Combining their values assuming normal errors for each, we obtain $M_{K_s} = 5.40 \pm 0.22$. This magnitude leads to a mass of $M_* = 0.551 \pm 0.068$ M☉ from the MLR. Our mass determination is compatible with the one employed by JH15 and Kipping et al. (2014).

With this mass measurement and the absolute density from the photodynamical model, we measured a stellar radius of $R_* = 0.582 \pm 0.045$ R☉. The corresponding planetary masses and radii are listed in Table 2 and plotted in Fig. 5. The determined planetary radii are significantly larger than the ones measured by JH15 but in agreement with Kipping et al. (2014). We discuss these differences in detail in Section 6.1.

Both masses and radii determined with the MLR and density from the photodynamical model are in agreement with the result using the $K_s$-band magnitude relation in Mann et al. (2015): $M_* = 0.565 \pm 0.039$ M☉ and $R_* = 0.541 \pm 0.041$ R☉.

## 5.2 Internal structure

### 5.2.1 Interior characterization: method

For a detailed interior characterization, we use the probabilistic analysis of Dorn et al. (2017) that calculates the full range of possible interiors, given a set of observational constraints. The data that we use as constraints are planetary mass and bulk density, as well as stellar abundance estimated by Souto et al. (2017) and their respective uncertainties. By using planetary mass and bulk density, we naturally account for the correlation between mass and radius. Stellar abundances of refractory elements are candidates for placing constraints on the relative abundance of rock-forming elements (i.e. Mg, Si, Fe, Ca, Na, Al) in the planet bulk as discussed and applied in, e.g. Sotin, Grasset & Mocquet (2007) and Dorn et al. (2015). For interior characterization, the above derived masses and densities represent part of the data ($d_{char}$). The complete data $d_{char}$ comprise:
<table><caption>**Table 2.** Derived system parameters using the results form the photodynamical modelling and the stellar mass from the MLR.</caption><thead><tr><th>Parameter</th><th></th><th>Median and 68.3 per cent CI</th></tr></thead><tbody><tr><td><strong>Star</strong></td><td></td><td></td></tr><tr><td>Mass, M<sub>*</sub></td><td>(M<sub>☉</sub>)</td><td>0.551±0.068</td></tr><tr><td>Radius, R<sub>*</sub></td><td>(R<sub>☉</sub><sup>N</sup>)</td><td>0.582±0.045</td></tr><tr><td>Surface gravity, log g</td><td>(cgs)</td><td>4.647±0.059</td></tr><tr><td><strong>Kepler-138b</strong></td><td></td><td></td></tr><tr><td>Semimajor axis, a</td><td>(au)</td><td>0.0760±0.0033</td></tr><tr><td>K'</td><td>(m s<sup>-1</sup>)</td><td>0.083±0.020</td></tr><tr><td>Mass, M<sub>p</sub></td><td>(M<sub>⊕</sub>)</td><td>0.187±0.050</td></tr><tr><td>Radius, R<sub>p</sub></td><td>(R<sub>e⊕</sub><sup>N</sup>)</td><td>0.701±0.066</td></tr><tr><td>Surface gravity, logg<sub>p</sub></td><td>(cgs)</td><td>2.58<sup>+0.11</sup><sub>-0.13</sub></td></tr><tr><td><strong>Kepler-138c</strong></td><td></td><td></td></tr><tr><td>Semimajor axis, a</td><td>(au)</td><td>0.0922±0.0040</td></tr><tr><td>K'</td><td>(m s<sup>-1</sup>)</td><td>2.10<sup>+0.34</sup><sub>-0.44</sub></td></tr><tr><td>Mass, M<sub>p</sub></td><td>(M<sub>⊕</sub>)</td><td>5.2±1.2</td></tr><tr><td>Radius, R<sub>p</sub></td><td>(R<sub>e⊕</sub><sup>N</sup>)</td><td>1.67±0.15</td></tr><tr><td>Surface gravity, logg<sub>p</sub></td><td>(cgs)</td><td>3.264<sup>+0.091</sup><sub>-0.11</sub></td></tr><tr><td><strong>Kepler-138 d</strong></td><td></td><td></td></tr><tr><td>Semimajor axis, a</td><td>(au)</td><td>0.1301±0.0056</td></tr><tr><td>K'</td><td>(m s<sup>-1</sup>)</td><td>0.395<sup>+0.082</sup><sub>-0.092</sub></td></tr><tr><td>Mass, M<sub>p</sub></td><td>(M<sub>⊕</sub>)</td><td>1.17±0.30</td></tr><tr><td>Radius, R<sub>p</sub></td><td>(R<sub>e⊕</sub><sup>N</sup>)</td><td>1.68±0.15</td></tr><tr><td>Surface gravity, logg<sub>p</sub></td><td>(cgs)</td><td>2.614<sup>+0.094</sup><sub>-0.12</sub></td></tr></tbody></table>

$$K' = \frac{M_p \sin i}{M_*^{2/3} \sqrt{1-e^2}} \left(\frac{2\pi G}{P^4}\right)^{1/3}$$

CODATA 2014: $G = 6.67408 \times 10^{-11} \text{ m}^3 \text{ kg}^{-1} \text{ s}^{-2}$

IAU 2012: au = 149 597 870 700 m

IAU 2015: $R_{⊙}^{N} = 6.957 \times 10^{8} \text{ m}$, $(GM_{⊙}^{N} = 1.327 124 4 \times 10^{20} \text{ m}^3 \text{ s}^{-2}$, $R_{e⊕}^{N} = 6.378 1 \times 10^{6} \text{ m}$, $(GM_{e⊕}^{N} = 3.986 004 \times 10^{14} \text{ m}^3 \text{ s}^{-2}$

$$M_{⊙} = (GM_{⊙}^{N}/G, M_{⊕} = (GM_{e⊕}^{N}/G)$$

**Figure 5.** Radius versus mass diagram. Colour points with errorbars represent the 68.3 per cent credible intervals of masses and radii of Kepler-138b (magenta), Kepler-138c (green), and Kepler-138 d (blue), based on the MLR. The colour contours represent the corresponding two-parameter 39.3 per cent joint credible regions. The bands represent the 68.3 per cent credible intervals for the densities obtained from the photodynamical modelling. Coloured solid lines represent theoretical models for different compositions (Zeng & Sasselov 2013), labelled in the upper left. Grey points with errorbars are the known exoplanets in this range compiled by Jontof-Hutter et al. (2016). The black dots represent Solar system planets, by order of increasing mass: Mercury, Mars, Venus, and Earth. Open circles are largest Solar system moons, by order of increasing mass: Moon, Io, Callisto, Titan, and Ganymede. Solar system bodies data from NASA.

<table><caption>**Table 3.** Prior ranges for interior parameters. $m_{\text{env, max}}$ refers to the maximum gas mass fraction based on the scaling law of Ginzburg, Schlichting & Sari (2016, equation 18).</caption><thead><tr><th>Parameter</th><th>Prior range</th><th>Distribution</th></tr></thead><tbody><tr><td>Core radius, r<sub>core</sub></td><td>(0.01–1) r<sub>core + mantle</sub></td><td>Uniform in r<sub>core</sub><sup>3</sup></td></tr><tr><td>Fe/Si<sub>mantle</sub></td><td>0–Fe/Si<sub>star</sub></td><td>Uniform</td></tr><tr><td>Mg/Si<sub>mantle</sub></td><td>Mg/Si<sub>star</sub></td><td>Gaussian</td></tr><tr><td>Size of rocky interior, r<sub>core + mantle</sub></td><td>(0.01–1) R</td><td>Uniform in r<sub>core + mantle</sub><sup>3</sup></td></tr><tr><td>Water mass fraction, m<sub>water</sub></td><td>0–0.98</td><td>Uniform</td></tr><tr><td>Gas mass fraction, m<sub>env</sub></td><td>0–m<sub>env, max</sub></td><td>Uniform in log-scale</td></tr><tr><td>Planet luminosity, L<sub>env</sub></td><td>10<sup>18</sup>–10<sup>23</sup> erg s<sup>-1</sup></td><td>Uniform in log-scale</td></tr><tr><td>Gas metallicity, Z<sub>env</sub></td><td>0–1</td><td>Uniform</td></tr></tbody></table>

(i) planet masses and densities (Table 2);

(ii) planet effective temperature<sup>12</sup> (Souto et al. 2017);

(iii) relative stellar abundances of Fe, Si, Mg, Na, Ca, Al (Souto et al. 2017).

Our assumptions for the interior model are similar to those in Dorn et al. (2017) and are summarized in the following. We consider planets being made of iron-rich cores, silicate mantles, layers of water ice and oceans, and gas. We define a set of interior parameters that we constrain given prior considerations and data $d_{\text{char}}$. The interior parameters comprise:

(i) core size;

(ii) mantle size;

(iii) mantle composition (i.e. Fe/Si<sub>mantle</sub>, Mg/Si<sub>mantle</sub>);

(iv) mass fraction of water;

(v) gas mass fraction;

(vi) intrinsic luminosity;

(vii) envelope metallicity.

The prior distributions of the interior parameters are stated in Table 3.

For the interior model, we use a self-consistent thermodynamic model (Dorn et al. 2017). The model calculates the density profile for any given set of interior parameters. This allows us to calculate the respective mass, bulk density, and bulk abundances and compare them to the actual observed data $d_{\text{char}}$. The thermodynamic model comprises the equation of state (EoS) of iron by Bouchet et al. (2013), the silicate-mantle model by Connolly (2009) to compute equilibrium mineralogy for general mantle compositions and density profiles. For the water layers, we follow Vazan et al. (2013) using a quotidian equation of state (QEOS) and above a pressure of 44.3 GPa, we use the tabulated EoS from Seager et al. (2007). We assume an adiabatic temperature profile within core, mantle, and water layers. Compared to previous work of Dorn et al. (2017), we impose the additional condition that in case a water layer is present, there must be an atmosphere on top. Specifically, the atmosphere must impose an atmospheric pressure of at least the vapour pressure of water. Thereby, we exclude vapour or supercritical vapour phases in the water layer.

For the atmosphere, we solve the equations of hydrostatic equilibrium, mass conservation, and energy transport. For the EoS of elemental compositions of H, He, C, and O, we employ the CEA (Chemical Equilibrium with Applications) package (Gordon & McBride 1994), which performs chemical equilibrium calculations for an arbitrary gaseous mixture, including dissociation and

<footnote n="12"><sup>12</sup>Computed from the stellar radius, stellar temperature, and semimajor axis, assuming zero albedo.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Two-dimensional (2D) marginalized posterior of interior parameters for (a–b) Kepler-138b, (c–d) Kepler-138c, (e–h) Kepler-138d. Shown interior parameters are the radius fractions of (a, c, e) atmosphere $r_{env}/R$ and water layer $r_{water}/R$ and (b, d, f) rocky interior $r_{core + mantle}/R$ and core $r_{core}/R$. The contours correspond to 2D joint probability levels of 68 and 95 per cent of the posterior distribution. An Earth-like interior is shown for reference.</figcaption></figure>

ionization and assuming ideal gas behaviour. The envelope metallicity $Z_{env}$ is the mass fraction of C and O in the gas, which can range from 0 to 1. An irradiated atmosphere is assumed at the top of the gaseous envelope, for which we assume a semigrey, analytic, global temperature averaged profile (Guillot 2010). The boundary between the irradiated atmosphere and the envelope is defined where the optical depth in visible wavelength is $100/\sqrt{3}$ (Jin et al. 2014). Within the envelope, the usual Schwarzschild criterion is used to distinguish between convective and radiative layers. The planet radius is defined where the chord optical depth becomes 2/3 (Mihalas 1978).

### 5.2.2 Interior characterization: results and discussion

Fig. 6 shows selected projections of posterior samples of the interiors of the Kepler-138 planets. The rocky interior (core and mantle) for Kepler-138b, c, and d can range from 0.66 to 0.83 $R_p$, 0.69 to 0.91 $R_p$, 0.49 to 0.68 $R_p$ (within $1\sigma$ of the posterior distribution), respectively. Kepler-138c is the only planet that may be dominated by a rocky interior. However, all three planets, including Kepler-138c, can have massive layers of volatiles. Whether or not these layers are in form of water or gas is impossible to determine given the available data. However, substantial radius fractions of gas layers are very likely for Kepler-138b and Kepler-138d. This is because we a priori exclude models where no gas layers are on top of water layers.

Kepler-138b has a small mass of nearly twice the mass of Mars and is highly irradiated. Our results indicate that planet b must have a significant thick envelope on top of the rocky interior. This thick envelope can be comprised of an enriched atmosphere

($Z_{env} > 0.3$) or contain layers of condensed volatiles (e.g. water). Both scenarios suggest that planet b must have accreted material from outside the snow-line. We find that a H-dominated envelope is very unlikely. Despite having a similar radius, Kepler-138c and Kepler-138d are significantly different in terms of bulk density. The lower density of Kepler-138d implies a higher volatile content. Possible formation scenarios include the accretion of primordial gas (H–He dominated) and subsequent partial retention of the gas or an efficient accretion of volatiles from outside the snow-line. Further considerations that account for the evolution of these planets and the possibility of atmospheric erosion are required to gain a better understanding of their formation and evolution.

# 6 DISCUSSION

## 6.1 Comparison with previous results

The results based on the photodynamical modelling of the Kepler light curve differ from those reported by JH15. Despite of using the same photometric data and similar hypotheses, there exists a significant difference between the system parameters reported by both studies (see Table 4).

### 6.1.1 Density

The main reason is the difference in the stellar parameters, and in particular in the stellar bulk density obtained from the transit light curve. Our analysis yields a density of $3.92_{-0.66}^{+0.81} \text{ g cm}^{-3}$. This is 42 per cent smaller than the value reported by JH15 (9.5 $\pm$ 2.2 g cm$^{-3}$). The high density determined by these authors is problematic.

First, this value is not compatible with theoretical models for their effective temperature $T_{eff} = 3841 \pm 49$ K and metallicity value $[Fe/H] = -0.280 \pm 0.099$. JH15 use the Dartmouth stellar evolution models (Dotter et al. 2008) to determine the absolute dimensions of the star and planets in the system. However, the Dartmouth tracks corresponding to their metallicity determination have lower densities (Fig. 7). Even a metallicity as low as $[Fe/H] = -0.5$ is not enough to solve the discrepancy. Furthermore, a recent study (Souto et al. 2017) has established that Kepler-138 has close-to-solar metallicity ($[Fe/H] = -0.09 \pm 0.09$). Invoking an unrealistically low metallicity to explain the high density seems, therefore, impossible. The discrepancy is further augmented by the fact that stellar models overestimate the effect of metallicity on the stellar radii (Boyajian et al. 2012).

Secondly, the mean densities measured on other stars of similar effective temperature hosting transiting companions are systematically smaller than the value reported by JH15 (see e.g. Pepper et al. 2017, fig. 11).

In Kipping et al. (2014), the stellar density is determined to be $2.75_{-0.47}^{+0.70} \text{ g cm}^{-3}$ from a fit to the two outer planets in the system. This value is closer to our determination, but it is not in agreement with the spectroscopic measurements. The discrepancy can, however, be explained by invoking a 2-min TTV not accounted for in the model (Kipping et al. 2014).

### 6.1.2 Radius

The larger density and lower metallicity of JH15 leads to a stellar radius which is 24 per cent smaller than the one determined here:
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Full Dartmouth tracks (lines) in the density versus $T_{\text{eff}}$ plane, wider lines represent times from the beginning of the main sequence to 13.82 Gyr, the shaded grey area is the envelope of the latter. Tracks corresponding to two values of the metallicity are represented: [Fe/H] = -0.5 (left, with the stellar mass annotated at the beginning of the main sequence), and [Fe/H] = 0.0 (right). The black point and errorbars represent the density determination from the photodynamical modelling, whereas the grey point and errorbars represent that from JH15. For both points, we considered the $T_{\text{eff}}$ that comes from Muirhead et al. (2012), which is used in JH15.</figcaption></figure>

$0.582 \pm 0.045 R_{\odot}^N$ versus $0.442 \pm 0.024 R_{\odot}$. Our value is in agreement with the one used by Kipping et al. (2014), $0.54 \pm 0.05 R_{\odot}$, from Pineda et al. (2013). The radius determined by JH15 is systematically smaller than previous determinations of Kepler-138 based either on theoretical evolution models or interferometric measurements (see Fig. 8). The average radius produced by these studies is $0.52 R_{\odot}$. The planet radii computed here are therefore also significantly larger than the values determined by JH15. This changes the inferred nature of Kepler-138d, as discussed below.

JH15 noted the differences in stellar parameters compared to Pineda et al. (2013), but argued that the calibration in this paper excludes active stars, and should therefore not be applicable to Kepler-138. However, the cut-off used by Pineda et al. (2013) excludes stars with rotational periods shorter than 5 d. Stars as active as Kepler-138 are therefore included in the calibration, as its rotational period is around 20 d (Section 6.2.1).

### 6.1.3 Mass

The stellar mass, on the other hand, is compatible along different studies (see Fig. 8). However, there exists an important difference between the mass ratios obtained here and those determined by JH15. We find mass ratios which are significantly larger than those in JH15. As a consequence, the derived planetary masses are larger in our analysis.

What is the origin of these discrepancies? A clear difference between this work and JH15 is the way in which the light curve is analysed. While the photodynamical analysis consistently models the observed light curve accounting for the gravitational interactions between the system bodies, the traditional TTV analysis artificially constructs two data sets based on the light curve. On the one hand, the mid-transit times measured on each individual transit are used to model the dynamical effects. On the other hand, the individual transit curves are stacked together to measure the transit shape parameters and stellar density.

In this kind of TTV analyses, the information of the transit shape is not taken into account when modelling the dynamics of the system. As a consequence, solutions which are not compatible with the observed planetary transit light curves are considered valid.<fnref n="13" /> This leads, to excessively wide eccentricity distributions, including values which are incompatible with the observed transit duration, let alone its variations. JH15 deals with this by imposing informative priors on the eccentricity distribution. A clear improvement of the method is therefore obtained by including the information of the transit duration in the dynamical modelling, as done, for example, by Kipping et al. (2014) and Nesvorný et al. (2013). In the case of Kepler-138, failing to account for the transit duration allows for larger eccentricity values. Because of the anticorrelation between eccentricity and mass ratios (see Extended Data fig. 5 of JH15) this can explain, at least partially, that our measured masses are larger than those reported by JH15.

Besides, to stack together the transit light curves, a single value of the transit times must be used. The transit times measured individually, however, are known with finite – usually poor – precision. This produces two effects: on the one hand, the uncertainties of the model parameters are underestimated, and on the other hand, the shape of the transit is distorted systematically towards a more V-shaped one. This will bias the determination of the stellar density, orbital inclination, and radius ratio. Even small differences in the TTV can produce observable effects. Here, the bias is in the opposite direction than the one studied by Kipping (2014b): the density of the JH15 analysis is smaller than the one measured with a photodynamical model. However, in Kipping (2014b) the bias is studied for an analysis that assumes no TTVs and circular orbits. Here, instead, the bias is for an analysis fixing the deviation of the transit time from a linear ephemeris. Furthermore, often the individual transit curves are shifted to agree with a linear ephemeris, artificially imposing an orbital period on the data. However, on interacting systems the period cannot be defined exactly.

Another difference is the treatment of limb darkening. JH15 fixed the values of the limb darkening parameters based on models. However, it has been shown that in general a better practice is to allow these parameters to vary (Csizmadia et al. 2013; Espinoza & Jordán 2015), although other researchers (Müller et al. 2013) show that for impact parameters larger than 0.8, transit light curves no longer constrain the limb darkening coefficients. However, their conclusion that the parameters should be fixed is not justified. Note that limb darkening coefficients are correlated with transit parameters $\frac{R_g}{R_\infty}$, inclination, and density. Fixing the limb-darkening law parameters effectively reduces the posterior size of these parameters and may augment the discrepancy produced by the biased analysis described above. Moreover, in the analysis by JH15 the limb-darkening coefficients are fixed before the transit fit, based on the spectroscopically determined $\log g$, but they are not updated once the transit stellar density is determined. The limb-darkening parameters are therefore inconsistent with their final stellar parameters.

We believe these effects may be at the root of the discrepancy in the stellar density measured by JH15 and by ourselves. However, we also explored other possibilities. To test if the discrepant result are not a result of inadequate exploration of parameter space, we initialized the MCMC algorithm at the solution reported by JH15.

<footnote n="13">¹³For example, it is possible to produce systems with very different eccentricities and periastron arguments that exhibit similar mid-transit times, although very different transit duration.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Stellar parameters comparison. Masses and radii derived using evolutionary stellar models are plotted with empty circles (Rowe et al. 2014, use the Yonsei-Yale models; the remaining works use Dartmouth tracks). Otherwise the parameters are determined using calibrations based on MLRs (Henry & McCarthy 1993; Delfosse et al. 2000; masses determined dynamically, luminosity with parallaxes) and radius from interferometric angular diameter measurements (Boyajian et al. 2012). Kipping et al. (2014) used the stellar parameters of Pineda et al. (2013). The values are listed in Table A1.</figcaption></figure>

<table><caption>**Table 4.** Comparison with previous works.</caption><thead><tr><th>Parameter</th><th></th><th>Kipping et al. (2014)</th><th>Jontof-Hutter et al. (2015)</th><th>This work</th></tr></thead><tbody><tr><td colspan="5"><strong>Kepler-138</strong></td></tr><tr><td>Stellar mean density, ρ<sub>★</sub></td><td>(g cm<sup>-3</sup>)</td><td>2.75<sup>+0.70a</sup><sub>-0.47</sub>, 5.3<sup>+2.1b</sup><sub>-1.9</sub></td><td>9.5±2.2</td><td>3.92<sup>+0.81</sup><sub>-0.66</sub></td></tr><tr><td>Mass, M<sub>★</sub></td><td>(M<sub>☉</sub>)</td><td>0.57±0.05</td><td>0.521±0.055</td><td>0.551±0.068</td></tr><tr><td>Radius, R<sub>★</sub></td><td>(R<sub>☉</sub>)</td><td>0.54±0.05</td><td>0.442±0.024</td><td>0.582±0.045</td></tr><tr><td>Surface gravity, log g</td><td>(cgs)</td><td>4.73<sup>+0.09</sup><sub>-0.09</sub></td><td>4.886±0.055</td><td>4.647±0.059</td></tr><tr><td colspan="5"><strong>Kepler-138 b</strong></td></tr><tr><td>Planet mean density, ρ<sub>p</sub></td><td>(g cm<sup>-3</sup>)</td><td></td><td>2.6<sup>+2.4</sup><sub>-1.5</sub></td><td>3.02<sup>+1.0</sup><sub>-0.87</sub></td></tr><tr><td>Mass, M<sub>p</sub></td><td>(M<sub>☉</sub>)</td><td></td><td>0.066<sup>+0.059</sup><sub>-0.037</sub></td><td>0.187±0.050</td></tr><tr><td>Radius, R<sub>p</sub></td><td>(R<sub>☉</sub>)</td><td>0.446<sup>+0.062</sup><sub>-0.050</sub></td><td>0.522<sup>+0.032</sup><sub>-0.032</sub></td><td>0.701±0.066</td></tr><tr><td>Radius ratio, R<sub>p</sub>/R<sub>★</sub></td><td></td><td>0.00753<sup>+0.00078</sup><sub>-0.00050</sub></td><td>0.0108<sup>+0.0003</sup><sub>-0.0003</sub></td><td>0.01104±0.00031<br>(1.04 ± 0.24) × 10<sup>-6</sup></td></tr><tr><td>Mass ratio, M<sub>p</sub>/M<sub>★</sub></td><td></td><td></td><td></td><td></td></tr><tr><td colspan="5"><strong>Kepler-138 c</strong></td></tr><tr><td>Planet mean density, ρ<sub>p</sub></td><td>(g cm<sup>-3</sup>)</td><td>5.0<sup>+3.0</sup><sub>-2.0</sub></td><td>6.2<sup>+5.8</sup><sub>-3.4</sub></td><td>6.1<sup>+1.9</sup><sub>-1.5</sub></td></tr><tr><td>Mass, M<sub>p</sub></td><td>(M<sub>☉</sub>)</td><td>3.83<sup>+1.51</sup><sub>-1.26</sub></td><td>1.970<sup>+1.912</sup><sub>-1.120</sub></td><td>5.2±1.2</td></tr><tr><td>Radius, R<sub>p</sub></td><td>(R<sub>☉</sub>)</td><td>1.61<sup>+0.16</sup><sub>-0.15</sub></td><td>1.197<sup>+0.070</sup><sub>-0.070</sub></td><td>1.67±0.15</td></tr><tr><td>Radius ratio, R<sub>p</sub>/R<sub>★</sub></td><td></td><td>0.02730<sup>+0.00087</sup><sub>-0.00070</sub></td><td>0.0247<sup>+0.0005</sup><sub>-0.0005</sub></td><td>0.02628<sup>+0.00048</sup><sub>-0.00043</sub></td></tr><tr><td>Mass ratio, M<sub>p</sub>/M<sub>★</sub></td><td></td><td>(2.03<sup>+0.76</sup><sub>-0.65</sub>) × 10<sup>-5</sup></td><td></td><td>(2.90<sup>+0.45</sup><sub>-0.60</sub>) × 10<sup>-5</sup></td></tr><tr><td colspan="5"><strong>Kepler-138 d</strong></td></tr><tr><td>Planet mean density, ρ<sub>p</sub></td><td>(g cm<sup>-3</sup>)</td><td>1.31<sup>+0.82</sup><sub>-0.54</sub></td><td>2.1<sup>+2.2</sup><sub>-1.2</sub></td><td>1.36<sup>+0.44</sup><sub>-0.35</sub></td></tr><tr><td>Mass, M<sub>p</sub></td><td>(M<sub>☉</sub>)</td><td>1.01<sup>+0.42</sup><sub>-0.34</sub></td><td>0.640<sup>+0.674</sup><sub>-0.387</sub></td><td>1.17±0.30</td></tr><tr><td>Radius, R<sub>p</sub></td><td>(R<sub>☉</sub>)</td><td>1.61<sup>+0.16</sup><sub>-0.15</sub></td><td>1.212<sup>+0.075</sup><sub>-0.075</sub></td><td>1.68±0.15</td></tr><tr><td>Radius ratio, R<sub>p</sub>/R<sub>★</sub></td><td></td><td>0.02731<sup>+0.00085</sup><sub>-0.00072</sub></td><td>0.0251<sup>+0.0007</sup><sub>-0.0007</sub></td><td>0.02643±0.00052</td></tr><tr><td>Mass ratio, M<sub>p</sub>/M<sub>★</sub></td><td></td><td>(0.53<sup>+0.21</sup><sub>-0.18</sub>) × 10<sup>-5</sup></td><td></td><td>(6.5<sup>+1.3</sup><sub>-1.5</sub>) × 10<sup>-6</sup></td></tr></tbody></table>

The walkers quickly moved out of this region in parameter space and approached our reported solution.

Moreover, in JH15 the transit shapes were modelled separately for each planet, and the planets were not even assumed to transit the same star. This provides an independent stellar density from each planet, which can be used to check if the planets orbit the same star and exclude scenarios of false-positives. However, once this hypothesis is assumed, then a combined consistent analysis should be preferred, and will produce a better precision in the system parameter determinations. It may even have the potential to reduce the biases produced by stellar activity or other effects.

## 6.2 Stellar characterization

### 6.2.1 Spot modelling

Kepler-138 is considered an active star based on the variability of its light curve (Reinhold, Reiners & Basri 2013, who use only Q3 data). The amplitude (around 1 per cent at most) is around two times larger than the arbitrary cut-off considered by the authors, with which 25 per cent of *Kepler* targets are in the active sample. With this criteria, the Sun is not in the active sample, even at its maximum activity level. The variability of young, active stars seems to be dominated by cool spots (Radick et al. 1998; Lock-

<footnote>a From a combined analysis of Kepler-138c and Kepler-138d assuming a common star.</footnote>

<footnote>b From an analysis of only the transits of Kepler-138b.</footnote>
<table><caption>**Table 5.** Modelling of the spectral energy distribution: Parameter, prior, posterior median, and 68.3 per cent CI. Prior on $T_{\text{eff}}$ and [Fe/H] are from Souto et al. (2017), and prior on log g and $R_*$ are from Table 2.</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Posterior median and 68.3 per cent CI</th></tr></thead><tbody><tr><td>Effective temperature, $T_{\text{eff}}$ (K)</td><td>$N(3835.64)$</td><td>3933 ± 37</td></tr><tr><td>Surface gravity, log g (cgs)</td><td>$N(4.647,0.059)$</td><td>4.662 ± 0.059</td></tr><tr><td>Metallicity, [Fe/H] (dex)</td><td>$N(-0.09,0.09)$</td><td>-0.127 ± 0.086</td></tr><tr><td>Radius, $R_*$ (R<sub>☉</sub>)</td><td>$N(0.582,0.045)$</td><td>0.586 ± 0.044</td></tr><tr><td>Distance (pc)</td><td>$U(0,1000)$</td><td>74.3 ± 5.8</td></tr><tr><td>$E_{(B-V)}$ (mag)</td><td>$U(0,3)$</td><td>0.0127<sup>+0.018</sup><sub>-0.0093</sub></td></tr></tbody></table>

$N(\mu, \sigma)$: Normal distribution prior with mean $\mu$ and standard deviation $\sigma$.

$U(1,u)$: Uniform distribution prior in the range [1, u].

wood et al. (2007) rotating with the stellar surface. One can worry about the influence of the stellar activity in the transit modelling and the determination of stellar parameters. In particular, spot crossing during transit could induce an underestimation of the stellar density (Léger et al. 2009; Barros et al. 2014). To gain insight on the variability of Kepler-138, we modelled the *Kepler* light using the *macula* code (Kipping 2012), with 39 spots. Details on the model and the procedure to choose the number of spots are given in Appendix B.

We found that the inclination of the rotational axis of the star with respect to the line of sight, $i_*$, must be close to 90°, and a rotational period in agreement with some of the previous determinations (see Table A2). Additionally, we detect a significant differential rotation of $P_{\text{POLE}}-P_{\text{EQ}} = 1.72^{+0.10}_{-0.17}$ d, where $P_{\text{POLE}}$ and $P_{\text{EQ}}$ are the rotational periods at the poles and at the equator, respectively. With the measured rotational period, we estimate a $\log_{10}(R'_{\text{HK}}) = -4.72 \pm 0.07$ (Suárez Mascareño et al. 2015) and $-4.48 \pm 0.02$ (Astudillo-Defru et al. 2017). The typical modelled spot size is $\alpha_{\text{max}} \sim 6^\circ$ but can be up to 14°, and the typical spot-to-star flux ratio is $f_{\text{spot}} \sim 0.69$, which corresponds to a photosphere to spot temperature difference of ~240 K (Silva-Valio et al. 2010, for an effective observation wavelength of 680 nm, based on the *Kepler* response function and the star spectral energy distribution). Spots have lifetimes of up to 2.5 yr, and are preferably located around latitude $\Phi \sim 27^\circ$ in each hemisphere, where they have longer lifetimes. We obtain a stellar surface spot coverage (spot-filling factor) between 0.3 and 3 per cent, although this should be considered a lower limit, as we do not model small spots nor large spots that do not produce significant variability in the light curve (like polar spots). Besides, we do not observe part of the stellar surface due to $i_*$. When the spot-to-star flux ratio is taken into account we obtain an effective spot coverage – that is, the equivalent covered area for zero-temperature spots, $f_{\text{spot}} = 0$– between 0.1 and 1.0 per cent.

For comparison, in the Sun, the typical spot size is ~5°, and the mean umbral core intensity is $f_{\text{umbra}} \sim 0.75$ (de Toma et al. 2013). A sunspot-photosphere difference temperature accounting for the umbra/penumbra is around 540 K (Lanza et al. 2009). The maximum spotted area observed in the Sun is about 0.2 per cent of the surface (Lanza et al. 2009), and the spots appear on belts ~35° wide in each hemisphere right above and below the equator.

In principle, the amount of spot coverage inferred in Kepler-138 could not be responsible for a considerable influence in the stellar radius or luminosity (Chabrier & Baraffe 2007; Morales et al. 2010; Jackson & Jeffries 2013).

The equivalent spot angular radius of the planets, projected on to the centre of the stellar disc, is ~0.7° for Kepler-138b, and ~1.6° for planets c and d.

The typical spot size (~6°) at the centre of the stellar disc corresponds to a $R_p/R_* \sim 0.1$, i.e. much larger than the planets. This kind of spot-modelling analysis can therefore be used to obtain information on the stellar flux distribution across the disc and the limb-darkening parameters more efficiently than using the transit light curves of much smaller planets.

The MAP normalized light curve (Fig. B1) is used to correct the transits (Czesla et al. 2009) prior to the photodynamical analysis (Section 3). Figs B2–B4 show the stellar surface at the time of the transits for the MAP model assuming a projected spin-orbit angle, $\lambda = 0$ (see e.g. Benomar et al. 2014), which seems a reasonable assumption in the light of the results by Albrecht et al. (2013).

### 6.2.2 Gyrochronology

The rotational period of an M-star is a more robust age indicator than proxies based on activity measurements like $Ca\bar{i}$, $H\alpha$ or $L_X$, because the latter can be affected by starspots, plages, activity cycles, and flaring (Engle & Guinan 2011). Besides, it can be complicated to determine the isochronal age of M stars because of their slow evolution. With the equatorial rotational period (Section 6.2.1) and the stellar mass (Section 5.1), we derived a gyrochronological age of $1.08^{+0.29}_{-0.11}$ Gyr (Barnes 2010; Barnes & Kim 2010, using a $P_0$ between 0.12 and 3.4 d), where we have added a systematic 10 per cent error to the statistical one (Meibom et al. 2015). A value in agreement with $1.13 \pm 0.23$ Gyr based only on the rotational period using the relation in Engle & Guinan (2011), obtained from a sample of M-stars. The small mass planets in this system, with orbital periods above 10 d, should not have affected significantly the rotational evolution of this M-star by tidal interactions (Lanza 2010). However, other factors can affect this gyrochronological age estimation (Epstein & Pinsonneault 2014), particularly for low-mass stars. This age determination must therefore be taken with caution.

### 6.2.3 Spectral energy distribution

The spectral energy distribution of Kepler-138 constructed using magnitudes from AAVSO (American Association of Variable Star Observers) Photometric All-Sky Survey (APASS),<fnref n="14" /> 2-Micron All-Sky Survey (2MASS; Skrutskie et al. 2006), and *Wide-field Infrared Survey Explorer* (WISE; Wright et al. 2010) is shown in Fig. 9. The measurements are listed in Table A3. We modelled the data using the PHOENIX/BT-Settl synthetic spectral library (Allard, Homeier & Freytag 2012) and the procedure described in Díaz et al. (2014), with the priors listed in Table 5. The results are reported in Table 5 as well and plotted in Fig. 9. We obtained a distance of $74.3 \pm 5.8$ pc. Using the stellar radius determination of JH15 as prior, as well as their spectroscopic parameters (Teff, log g, and [Fe/H]), we found a distance of $55.5 \pm 3.1$ pc. The corresponding discrepancy in the parallaxes is well within the measurement capability of the Gaia satellite (de Bruijne, Rygl & Antoja 2014).

## 6.3 Advantages and caveats of the photodynamical modelling

The photodynamical model exploits the available data in a more thorough manner than the traditional TTV analysis. In the first place, each transit constrains the timing of all other transits in the light curve, as discussed in Almenara et al. (2015). This leads to a

<footnote n="14">¹⁴http://aavso.org/apass</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Spectral energy distribution of Kepler-138. The solid line is the MAP PHOENIX/BT-Settll interpolated synthetic spectrum, red circles are the absolute photometric observations, and grey open circles are the result of integrating the synthetic spectrum in the observed bandpasses.</figcaption></figure>

much better precision in the transit timings, specially for low signal-to-noise transits. In Kepler-138, transit times are determined with a precision 9, 10, and 6 times better than the individual measurements of JH15 (Fig. 4). As a consequence, system parameters are also determined more precisely. For example, for the Kepler-138 system, densities are determined with a factor of 2 improvement in the precision with respect to JH15. Secondly, the information on the transit shapes is also included naturally in the model. This leads, for example, to eccentricity distributions that naturally exclude values that produce transit shapes incompatible with observations. As multiplanet systems tend to have low eccentricities (Van Eylen & Albrecht 2015), this means that large eccentricities often appearing in TTV analyses are usually excluded. Because there exists a negative correlation between orbital eccentricities and mass ratios, the masses obtained via a TTV analysis tend to be biased towards smaller values (Weiss & Marcy 2014).

In the case of Kepler-117 (Almenara et al. 2015), where the transits have a large signal-to-noise ratio, these features are clearly demonstrated. The density of the stellar host Kepler-117 is determined with a precision of 2 per cent, and that of the planets with a precision of 3 per cent. These are among the most precise determinations obtained to date, and are independent of stellar models. However, Kepler-117 is one of the most favourable cases in which to use this technique. For Kepler-138, on the other hand, densities are over an order of magnitude less precise (between 28 and 33 per cent for the planets, and 19 per cent for the star). This is mainly due to the low signal-to-noise ratio of the Kepler-138 transits.

### 6.3.1 Caveats

We mention in this section some caveats concerning the photodynamical modelling of multiplanet systems in general, and of our analysis of the Kepler-138 system in particular.

One important difficulty when performing statistical inference on the model parameters appears from the involved parameter space that arises when using photodynamical models. Thoroughly exploring the space of parameters associated with models of multiplanet systems is a challenging task. They contain a number of disconnected local maxima (Ragozzine & Holman 2010) that hinder correct exploration and makes it impossible to guarantee, in most cases, that the dominant mode in parameter space has been found (see e.g. Gillon et al. 2017). Therefore, the solution presented in this paper could just be one maximum among many others. The problem grows in complexity rapidly with the number of planets included in the model, and is aggravated by the strong non-linear correlations between certain parameters. One such correlation exists between the stellar density, planet orbital inclination, and radius ratios. If the impact parameters of the transits is lower, the star is denser and the planetary-to-star radius ratios are smaller.

Results obtained with a photodynamical modelling are also sensitive, in principle, to the filtering and detrending techniques employed on the light curve before it is analysed. We are currently unable to include these steps of the data reduction process in the iterative exploration of parameter space, as the number of parameters as well as the computation resources needed would become prohibitively large. However, we have tried detrending the individual transits using linear, quadratic and cubic polynomials. The results are independent of the degree of the polynomial. In particular, the posterior distributions of the stellar bulk density obtained with the three data sets are almost identical. Nevertheless, we cannot guarantee that our solution is completely independent of the chosen normalization method.

Similarly, unaccounted for stellar activity can bias the results by affecting the transit shapes systematically. Barros et al. (2013) showed that transit timings can be severely affected by planets crossing dark stellar spots. This should be particularly important for small-planet transits, where individual spot crossings cannot be detected. On the other hand, the use of a photodynamical model that relies on the entire set of transits to infer the timing of each individual event should be less affected than the traditional two-step approach, in which timings are measured on individual transits. Using Gaussian process regression on the light curve to model unwanted signals produced by the star should allow us to reach more robust parameter determinations, at expenses of a more complex model, increasing the computing time needed, which is already large.

The measured densities can also be affected by the presence of undetected blend sources, either physically bound or fortuitously aligned with the target star in the line of sight from Earth. The effect is known as the photoblend effect (Kipping 2014b), and causes a systematic underestimation of the measured bulk densities. In the case of Kepler-138, however, high-resolution imaging discards most chance-aligned scenarios (Wang et al. 2015), therefore reducing the probability that our results are affected by a blend.

Last, but not least, the analysis is performed conditional to the assumption on the number of planets in the system. Here, we have chosen to include only the three transiting planets. However, under this model at least the two last transits of the exterior Kepler-138d seem to be incorrectly reproduced by the model, which hints to the presence of additional bodies in the system. Numerical simulations suggest that massive companion in outer orbits may be a common feature accompanying tightly packed low-mass planetary systems (Hands & Alexander 2016). On the other hand, it is not expected that low-mass stars such as Kepler-138 form massive companions efficiently. We have tried to model the data assuming four planets in the system. However, the difficulties listed above and the limited time span of the observations make it hard to precisely determine the parameters under this hypothesis. In the solution we found, the fourth exterior planet would be on an unrealistic high-eccentricity orbit, precisely aligned to produce a close encounter with Kepler-138d at the moment of the last transits of this planet. It is therefore hard to establish if a fourth planet is actually present in the system and how it would precisely affect our results. However, we can safely assume that the effect on the inner planet pair – b and c
– would be smaller than on the outer planet. First, these planets are close to the 4:3 resonance, and therefore the secular effect of an outer companion in a larger order resonance would be smaller. Secondly, short-term effects should also be dominated by the close companions, as the perturber would be much farther away. Radial velocity follow-up of this system could place constrains on the properties of a potential outer companion. However, the relative faintness of the star would probably require a large number of observations to obtain meaningful constrains.

## 7 CONCLUSIONS

The photodynamical analysis models consistently a transit light curve accounting for the gravitational interactions between the system bodies. Contrary to TTV-only analysis, it fully exploits the information contained in the transit shape, including transit duration variations, transit depth variations, but more generally, transit shape variations. This leads to an increased precision in all model parameters with respect to the traditional TTV analysis based on measurements of the transit times of individual transits. In fact, the photodynamical analysis allows measuring the transit times, conditional to the model hypotheses, with a precision 6–10 times better than the individual measurements (Fig. 4). Therefore, even if gravitational interactions are not detected based on individual TTV measurements, it is in principle worth performing a complete photodynamical analysis, which may reveal subtle signs of interactions not detectable with the traditional methods. In particular, the detection of gravitational interactions between the transiting companions can work as a confirmation of their planetary nature. Even if no significant interactions are detected, this methods will provide the most stringent constrains to the system parameters.

A further advantage of this way of analysing photometric time series is the independence on stellar evolution models. As a consequence, only certain physical parameters are measurable, unless additional data sets or information are provided that break the Newtonian degeneracy. One prime example of such measurements are stellar radial velocity variations induced by the planetary companions. This measurement is challenging, and in some cases, such as the Kepler-138 system, the required precision is unattainable by current instrumentation. Alternatively, a direct measurement of the radius can be obtained from interferometry (Ligi et al. 2016), although a slight dependence with the limb-darkening exists if the star is not completely resolved (Mérand et al. 2010). This technique is not accessible for *Kepler* targets, but will be determinant for TESS and PLATO targets. A very promising option in the light of the Gaia mission is the inclusion of distance measurements, which also break the Newtonian degeneracy.

The method is susceptible to a series of assumptions related to the number of planets in the system, the level of flux contamination caused by unseen stellar companions, and limitations, linked mostly to the involved parameter space associated with the photodynamical model. However, under the hypothesis of a three-planet system, and assuming that parameter space was thoroughly explored and that we have found the global posterior maximum, the photodynamical modelling of the available *Kepler* photometry of Kepler-138 revealed a system not in agreement with previous studies JH15. The main discrepancy lies in the stellar density, which in turn affects the determination of the stellar radius and therefore the planetary sizes. We inferred a system composed of a rocky interior planet with a mass of $0.187 \pm 0.050 M_\oplus$, and a radius of $0.701 \pm 0.066 R_{e\oplus}^N$. From this, we deduced the presence of a substantial gaseous envelope constituting around 15 per cent of the planetary radius. The two

outer planets, with periods 13.8 and 23.1 d, have similar radius, but very different densities. While Kepler-138c may be purely rocky, Kepler-138d requires an outer envelope to explain the observations (see also Hadden & Lithwick 2017). The results for the exterior planets are in agreement with the study by Kipping et al. (2014).

Kepler-138 planets are subject to low incident fluxes ($11.5 \pm 1.7 F_\oplus$, $7.8 \pm 1.1 F_\oplus$, and $3.91 \pm 0.58 F_\oplus$, for Kepler-138 b, c, and d, respectively) and therefore the proposed evaporation mechanism (Lopez 2017; Fulton et al. 2017) should not be efficient. Thus, the inferred variability in bulk compositions may be primordial.

We warn that the understanding of the system may not be complete, as an exterior perturber planet may be present. Interestingly some studies suggest that exterior giant planets may determine the architecture of tightly packed multiplanet systems like Kepler-138 (Hands & Alexander 2016). However, this would affect more strongly, in principle, the parameters of Kepler-138d, the interior pair (Kepler-138b and Kepler-138c) being in a lower order resonance.

We found that the radius of Kepler-138b is more than 26 per cent larger and its mass is 2.8 times larger than previously determined. The mass and radius of Mars lie in the boundary of the 97 per cent credible region of a bivariate normal distribution in the mass–radius plane centred in the posterior means of these parameters for Kepler-138b (Table 2) and with a covariance matrix equal to the empirical sample covariance. This is one of the smallest planets known to date with a mass measurement. With a size and mass between those of Mars and Venus, this planet does not have an analogue in the Solar system.

## ACKNOWLEDGEMENTS

We thank the anonymous referee from a previous failed paper who made us realized the underlying symmetry in Newtonian gravity, C. Damiani, J. Stadel, and R. Mardling for discussions about dynamics, S. Peretti for kindly sharing with us the updated MLR prior to its publication, J. Rowe for sending us the short-cadence *Kepler* data prior to data release 24, Y. Revaz for his assistance with the computing cluster used in this work, and L. Kreidberg for her Mandel & Agol code. We wish to thank the International Space Science Institute (ISSI) and the members of the international team lead by J. Cabrera. This paper includes data collected by the *Kepler* mission. Funding for the *Kepler* mission is provided by the NASA Science Mission directorate. Data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. The original description of the VizieR service was published in A&AS 143, 23. This research was made possible through the use of the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet
Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. Simulations in this paper made use of the REBOUND code which can be downloaded freely at http://github.com/hannorein/rebound. These simulations have been run on the Regor cluster kindly provided by the Observatoire de Genève. JMA and XB acknowledges funding from the European Research Council under the ERC Grant Agreement n. 337591-ExTra. XB acknowledges the support of the French Agence Nationale de la Recherche (ANR), under the programme ANR-12-BS05-0012 Exo-atmos. This work has been carried out within the framework of the National Centre for Compe- tence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF.

## REFERENCES

Agol E., Steffen J., Sari R., Clarkson W., 2005, MNRAS, 359, 567

Albrecht S., Winn J. N., Marcy G. W., Howard A. W., Isaacson H., Johnson J. A., 2013, ApJ, 771, 11

Allard F., Homeier D., Freytag B., 2012, Phil. Trans. R. Soc. A, 370, 2765

Almenara J. M., Díaz R. F., Mardling R., Barros S. C. C., Damiani C., Bruno G., Bonfils X., Deleuil M., 2015, MNRAS, 453, 2644

Astudillo-Defru N., Delfosse X., Bonfils X., Forveille T., Lovis C., Rameau J., 2017, A&A, 600, A13

Barnes S. A., 2010, ApJ, 722, 222

Barnes S. A., Kim Y.-C., 2010, ApJ, 721, 675

Barros S. C. C., Boué G., Gibson N. P., Pollacco D. L., Santerne A., Keenan F. P., Skillen I., Street R. A., 2013, MNRAS, 430, 3032

Barros S. C. C. et al., 2014, A&A, 569, A74

Batalha N. M. et al., 2011, ApJ, 729, 27

Benomar O., Masuda K., Shibahashi H., Suto Y., 2014, PASJ, 66, 94

Berger D. H. et al., 2006, ApJ, 644, 475

Bonfils X., Delfosse X., Udry S., Santos N. C., Forveille T., Ségransan D., 2005, A&A, 442, 635

Borucki W. J. et al., 2010, Science, 327, 977

Bouchet J., Mazevet S., Morard G., Guyot F., Musella R., 2013, Phys. Rev. B, 87, 094102

Boyajian T. S. et al., 2012, ApJ, 757, 112

Carter J. A. et al., 2012, Science, 337, 556

Chabrier G., Baraffe I., 2007, ApJ, 661, L81

Christiansen J. L. et al., 2017, AJ, 154, 122

Connolly J. A. D., 2009, Geochemi. Geophys. Geosyst., 10, Q10014

Csizmadia S., Pasternacki T., Dreyer C., Cabrera J., Erikson A., Rauer H., 2013, A&A, 549, A9

Czesla S., Huber K. F., Wolter U., Schröter S., Schmitt J. H. M. M., 2009, A&A, 505, 1277

de Bruijne J. H. J., Rygl K. L. J., Antoja T., 2014, in Nicholas W., eds proceedings of the GREAT-ITN conference "The Milky Way Unravelled by Gaia: GREAT Science from the Gaia Data Releases", 1-5 December 2014, EAS Publ. Ser. Vol. 67, University of Barcelona, Spain, EAS Publications Series, Francesca Figueras, and Caroline Soubiran, p. 23

Delfosse X., Forveille T., Ségransan D., Beuzit J.-L., Udry S., Perrier C., Mayor M., 2000, A&A, 364, 217

de Toma G., Chapman G. A., Cookson A. M., Preminger D., 2013, ApJ, 771, L22

Dorn C., Hinkel N. R., Venturini J., 2017, A&A, 597, A38

Dorn C., Khan A., Heng K., Connolly J. A., Alibert Y., Benz W., Tackley P., 2015, A&A, 577, A83

Dorn C., Venturini J., Khan A., Heng K., Alibert Y., Helled R., Rivoldini A., Benz W., 2017, A&A, 597, A37

Dotter A., Chaboyer B., Jevremović D., Kostov V., Baron E., Ferguson J. W., 2008, ApJS, 178, 89

Díaz R. F., Almenara J. M., Santerne A., Moutou C., Lethuillier A., Deleuil M., 2014, MNRAS, 441, 983

Engle S. G., Guinan E. F., 2011, in Qain S., Leung K., Zhu L., Kwok S., eds, ASP Conf. Ser. Vol. 451, 9th Pacific Rim Conf. Stellar Astrophys. Astron. Soc. Pac., San Francisco, p. 285

Epstein C. R., Pinsonneault M. H., 2014, ApJ, 780, 159

Espinoza N., Jordán A., 2015, MNRAS, 450, 1879

Ford E. B. et al., 2011, ApJS, 197, 2

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Fulton B. J. et al., 2017, AJ, 154, 109

Gillon M. et al., 2017, Nature, 542, 456

Ginzburg S., Schlichting H. E., Sari R., 2016, ApJ, 825, 29

Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Gordon S., McBride B. J., 1994, Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications, Vol. 1. National Aeronautics and Space Administration, Office of Management, Scientific and Technical Information Program, NASA. NASA Reference Publication 1311, Washington, DC

Guillot T., 2010, A&A, 520, A27

Hadden S., Lithwick Y., 2017, AJ, 154, 5

Hands T. O., Alexander R. D., 2016, MNRAS, 456, 4121

Hatzes A. P., Rauer H., 2015, ApJ, 810, L25

Henry T. J., McCarthy D. W., Jr., 1993, AJ, 106, 773

Heyl J. S., Gladman B. J., 2007, MNRAS, 377, 1511

Huber D. et al., 2013, Science, 342, 331

Jackson R. J., Jeffries R. D., 2013, MNRAS, 431, 1883

Jin S., Mordasini C., Parmentier V., Van Boekel R., Henning T., Ji J., 2014, ApJ, 795, 65

Jontof-Hutter D., Lissauer J. J., Rowe J. F., Fabrycky D. C., 2014, ApJ, 785, 15

Jontof-Hutter D., Rowe J. F., Lissauer J. J., Fabrycky D. C., Ford E. B., 2015, Nature, 522, 321 (JH15)

Jontof-Hutter D. et al., 2016, ApJ, 820, 39

Kipping D. M., 2010, MNRAS, 408, 1758

Kipping D. M., 2012, MNRAS, 427, 2487

Kipping D. M., 2013a, MNRAS, 435, 2152

Kipping D. M., 2013b, MNRAS, 435, 2152

Kipping D. M., 2014a, MNRAS, 440, 2164

Kipping D. M., 2014b, MNRAS, 440, 2164

Kipping D. M., Nesvorný D., Buchhave L. A., Hartman J., Bakos G. Á., Schmitt A. R., 2014, ApJ, 784, 28

Lanza A. F., 2010, A&A, 512, A77

Lanza A. F. et al., 2009, A&A, 493, 193

Ligi R., et al., 2016, A&A, 586, 94

Lissauer J. J. et al., 2011, Nature, 470, 53

Lockwood G. W., Skiff B. A., Henry G. W., Henry S., Radick R. R., Baliunas S. L., Donahue R. A., Soon W., 2007, ApJS, 171, 260

Lopez E. D., 2017, MNRAS, 472, 245

Léger A. et al., 2009, A&A, 506, 287

López-Morales M., 2007, ApJ, 660, 732

Mandel K., Agol E., 2002, ApJ, 580, L171

Manduca A., Bell R. A., Gustafsson B., 1977, A&A, 61, 809

Mann A. W., Feiden G. A., Gaidos E., Boyajian T., von Braun K., 2015, ApJ, 804, 64

Mann A. W., Gaidos E., Ansdell M., 2013, ApJ, 779, 188

Mann A. W., Gaidos E., Kraus A., Hilton E. J., 2013, ApJ, 770, 43

Mazeh T., Perets H. B., McQuillan A., Goldstein E. S., 2015, ApJ, 801, 3

Mazeh T. et al., 2013, ApJS, 208, 16

McQuillan A., Aigrain S., Mazeh T., 2013, MNRAS, 432, 1203

McQuillan A., Mazeh T., Aigrain S., 2013, ApJ, 775, L11

Meibom S., Barnes S. A., Platais I., Gilliland R. L., Latham D. W., Mathieu R. D., 2015, Nature, 517, 589

Mihalas D., 1978, Stellar atmospheres, WH Freeman and Co., San Francisco, p. 650

Morales J. C., Gallardo J., Ribas I., Jordi C., Baraffe I., Chabrier G., 2010, ApJ, 718

Morales J. C., Ribas I., Jordi C., 2008, A&A, 478, 507
Morales J. C. et al., 2009, ApJ, 691, 1400

Muirhead P. S., Hamren K., Schlawin E., Rojas-Ayala B., Covey K. R.,
Lloyd J. P., 2012, ApJ, 750, L37

Muirhead P. S. et al., 2014, ApJS, 213, 5

Mérand A. et al., 2010, A&A, 517, A64

Müller H. M., Huber K. F., Czesla S., Wolter U., Schmitt J. H. M. M., 2013, A&A, 560, A112

Nesvorný D., Kipping D., Terrell D., Hartman J., Bakos G. Á., Buchhave L. A., 2013, ApJ, 777, 3

Newton E. R., Charbonneau D., Irwin J., Mann A. W., 2015, ApJ, 800, 85

Paz-Chinchón F. et al., 2015, ApJ, 803, 69

Pepe F. et al., 2013, Nature, 503, 377

Pepper J. et al., 2017, AJ, 153, 177

Pineda J. S., Bottom M., Johnson J. A., 2013, ApJ, 767, 28

Queloz D. et al., 2009, A&A, 506, 303

Radick R. R., Lockwood G. W., Skiff B. A., Baliunas S. L., 1998, ApJS, 118, 239

Ragozzine D., Holman M. J., 2010, pre print (arXiv: 1006.3727)

Ragozzine D., Wolf A. S., 2009, ApJ, 698, 1778

Rein H., Liu S.-F., 2012, A&A, 537, A128

Rein H., Tamayo D., 2015, MNRAS, 452, 376

Reinhold T., Reiners A., Basri G., 2013, A&A, 560, A4

Rowe J. F. et al., 2014, ApJ, 784, 45

Sanchis-Ojeda R. et al., 2012, Nature, 487, 449

Seager S., Kuchner M., Hier-Majumder C. A., Militzer B., 2007, ApJ, 669, 1279

Seager S., Mallén-Ornelas G., 2003, ApJ, 585, 1038

Silva-Valio A., Lanza A. F., Alonso R., Barge P., 2010, A&A, 510, A25

Skrutskie M. F., et al., 2006, AJ, 131, 1163

Sotin C., Grasset O., Mocquet A., 2007, Icarus, 191, 337

Southworth J., 2009, MNRAS, 394, 272

Souto D. et al., 2017, ApJ, 835, 239

Suárez Mascareño A., Rebolo R., González Hernández J. I., Esposito M., 2015, MNRAS, 452, 2745

Terrien R. C., Mahadevan S., Bender C. F., Deshpande R., Robertson P., 2015a, ApJ, 802, L10

Terrien R. C., Mahadevan S., Deshpande R., Bender C. F., 2015b, ApJS, 220, 16

Torres G., 2010, AJ, 140, 1158

Torres G., 2013, Astron. Nachr., 334, 4

Van Eylen V., Albrecht S., 2015, ApJ, 808, 126

Vazan A., Kovetz A., Podolak M., Helled R., 2013, MNRAS, 434, 3283

Walkowicz L. M., Basri G. S., 2013, MNRAS, 436, 1883

Wang J., Fischer D. A., Xie J.-W., Ciardi D. R., 2015, ApJ, 813, 130

Weiss L. M., Marcy G. W., 2014, ApJ, 783, L6

Wright E. L. et al., 2010, AJ, 140, 1868

Zeng L., Sasselov D., 2013, PASP, 125, 227
APPENDIX A: OTHER FIGURES AND TABLES

See Tables A1–A3 and Figs A1–A6.

<table><caption>**Table A1.** Parameters of Kepler-138 from the literature.</caption><thead><tr><th>Work</th><th>T<sub>eff</sub> (K)</th><th>log g</th><th>[Fe/H]</th><th>[M/H]</th><th>M<sub>*</sub> (M<sub>☉</sub>)</th><th>R<sub>*</sub> (R<sub>☉</sub>)</th><th>ρ<sub>*</sub> (g cm<sup>-3</sup>)</th><th>L<sub>*</sub> (L<sub>☉</sub>)</th><th>d (pc)</th><th>SpT</th></tr></thead><tbody><tr><td>Muirhead et al. (2012)</td><td>3841<sup>+50</sup><sub>-51</sub></td><td></td><td></td><td>-0.18 ± 0.10</td><td>0.52±0.06</td><td>0.50±0.06</td><td></td><td></td><td></td><td>M1V</td></tr><tr><td>Pineda et al. (2013)</td><td></td><td></td><td>-0.28±0.10</td><td></td><td>0.57±0.05</td><td>0.54±0.05</td><td></td><td></td><td>66.5±7.3</td><td></td></tr><tr><td>Mann et al. (2013)</td><td></td><td></td><td>-0.07±0.08</td><td></td><td></td><td>0.52±0.04</td><td></td><td></td><td></td><td></td></tr><tr><td>Mann, Gaidos & Ansdell (2013)</td><td>3871±58</td><td></td><td></td><td></td><td>0.593±0.065</td><td>0.553±0.034</td><td></td><td>0.060±0.008</td><td></td><td></td></tr><tr><td>Kipping et al. (2014)</td><td></td><td>4.73±0.09</td><td></td><td></td><td></td><td></td><td>5.3<sup>+2.1</sup><sub>-1.9</sub><br>2.75<sup>+0.70a</sup><sub>-0.47</sub></td><td></td><td></td><td></td></tr><tr><td>Rowe et al. (2014)</td><td>4079±238</td><td>4.747±0.150</td><td>-0.20±0.10</td><td></td><td></td><td>0.528±0.036</td><td>5.436±0.645</td><td></td><td></td><td></td></tr><tr><td>Muirhead et al. (2014)</td><td>3847<sup>+39</sup><sub>-60</sub></td><td></td><td>-0.25±0.12</td><td>-0.18±0.12</td><td>0.52±0.03</td><td>0.49±0.03</td><td></td><td></td><td>20</td><td>M1V</td></tr><tr><td>Newton et al. (2015)</td><td>3841±73</td><td></td><td></td><td></td><td></td><td>0.512±0.027</td><td></td><td></td><td></td><td></td></tr><tr><td>Jontof-Hutter et al. (2015)</td><td>3841±49</td><td>4.886±0.055</td><td>-0.280±0.099</td><td></td><td>0.521±0.055</td><td>0.442±0.024</td><td>9.5±2.2</td><td>10<sup>-1.33</sup> ± 0.05</td><td></td><td></td></tr><tr><td>Terrien et al. (2015b)</td><td></td><td></td><td>-0.21 ± 0.11</td><td>-0.16±0.10</td><td></td><td></td><td></td><td></td><td></td><td>M1.0</td></tr><tr><td>Souto et al. (2017)</td><td>3835±64</td><td>4.64±0.10</td><td>-0.09 ± 0.09</td><td></td><td>0.59±0.06</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>This Work</td><td></td><td>4.647±0.059</td><td></td><td></td><td>0.551±0.068</td><td>0.582±0.045</td><td>3.92<sup>+0.81</sup><sub>-0.66</sub></td><td></td><td>74.3±5.8</td><td></td></tr></tbody></table>

<table><caption>**Table A2.** Rotational periods of Kepler-138 from the literature.</caption><thead><tr><th>Work</th><th></th><th>P<sub>rot</sub> (d)</th></tr></thead><tbody><tr><td>McQuillan, Aigrain & Mazeh (2013)</td><td></td><td>18.860±0.111</td></tr><tr><td>McQuillan, Mazeh & Aigrain (2013)</td><td></td><td>19.394±0.013</td></tr><tr><td>Reinhold et al. (2013)</td><td></td><td>18.937±0.093</td></tr><tr><td>Walkowicz & Basri (2013)</td><td></td><td>9.77±7.98</td></tr><tr><td>Mazeh et al. (2015)</td><td></td><td>19.34</td></tr><tr><td>Paz-Chinchón et al. (2015)</td><td></td><td>9.70±0.17</td></tr></tbody></table>

<table><caption>**Table A3.** Photometric measurements used for the SED analysis of Kepler-138.</caption><thead><tr><th>Filter</th><th>Magnitude</th><th>±1σ</th><th>Source</th></tr></thead><tbody><tr><td>Johnson-B</td><td>14.621</td><td>0.052</td><td>APASS DR9</td></tr><tr><td>Johnson-V</td><td>13.168</td><td>0.033</td><td>APASS DR9</td></tr><tr><td>SDSS-G</td><td>13.870</td><td>0.056</td><td>APASS DR9</td></tr><tr><td>SDSS-R</td><td>12.529</td><td>0.055</td><td>APASS DR9</td></tr><tr><td>SDSS-I</td><td>11.943</td><td>0.036</td><td>APASS DR9</td></tr><tr><td>2MASS-J</td><td>10.293</td><td>0.022</td><td>2MASS</td></tr><tr><td>2MASS-H</td><td>9.680</td><td>0.018</td><td>2MASS</td></tr><tr><td>2MASS-Ks</td><td>9.506</td><td>0.011</td><td>2MASS</td></tr><tr><td>WISE-W1</td><td>9.378</td><td>0.023</td><td>WISE</td></tr><tr><td>WISE-W2</td><td>9.355</td><td>0.019</td><td>WISE</td></tr><tr><td>WISE-W3</td><td>9.261</td><td>0.030</td><td>WISE</td></tr></tbody></table>

<footnote><sup>a</sup> The first value is estimated from Kepler-138b, while the second is from a combined analysis of planets c and d.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure A1.** Transits of Kepler-138b observed by *Kepler*. Dots represent the individual short-cadence observations and larger circles are 30-min averaged values. In those panels without short-cadence points, the circles represent the long-cadence data. Each panel is labelled with the transit epoch, and centred relative to a linear ephemeris. The model distribution is constructed from 1000 random MCMC steps. The black line is the median model, and 1, 2, and 3σ confidence intervals are shown in three different grey-scales. In the lower part of each panel, the residuals after subtracting the mean model are shown.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure A2. Idem Fig. A1 for Kepler-138c.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure A3. Idem Fig. A1 for Kepler-138d.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure A4.** 2D projections of the joint posterior samples obtained with the MCMC algorithm. The 39.3, 86.5, and 98.9 per cent 2D joint confidence regions (in the case of a Gaussian posterior, these regions project on to the one-dimensional (1D) 1, 2, and 3σ intervals) are denoted by three different grey levels. The 1D histogram of each parameter is shown at the top of each column, except for the parameter on the last line that is shown at the end of the line. Units are the same as in Table 1.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure A5.** Mutual inclination between planet pairs at $t_{\text{ref}}$. Configuration A (top) and D (bottom).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure A6.** Evolution of the eccentricity of Kepler-138b during *Kepler* observations from the photodynamical modelling. The 68.3, 95.5, and 99.7 per cent Bayesian credible intervals are plotted in different shades. The red curve marks the median of the posterior distribution. The grey curve correspond to the model based on the MAP values.</figcaption></figure>
APPENDIX B: SPOT MODELLING

In this section, we give further detail on the spot modelling described in Section 6.2.1. For this we used the PDCSAP version of the long-
cadence *Kepler* data, which is designed to preserve the astrophysical signal, such as rotational modulations due to starspots. We removed the
transits and we averaged the data in 4.9 h bins. In some cases, artificial variability is seen in the light curve after gaps in the data, which are
due to the stabilization of the satellite after a change in pointing. These parts were not considered for the modelling.

The code macula models circular, non-overlapping, small (≈10°) starspots, with linear size evolution. The star is described by seven free
parameters: the equatorial rotational period (P<sub>EQ</sub>), a quadratic differential rotation coefficient (κ),<fnref n="15" /> the inclination of the star rotation axis
with respect to the line of sight (i<sub>*</sub>), and two quadratic law limb-darkening coefficients for each the spotted and non-spotted stellar surface
(u<sub>a</sub>, u<sub>b</sub>, u<sub>a,spot</sub>, u<sub>b,spot</sub>). Each spot is modelled with eight parameters: the spot maximum angular radius (α<sub>max</sub>), its time span at α<sub>max</sub> (t<sub>max</sub>), the
mid time of the spot at α<sub>max</sub> (T<sub>max</sub>), ingress duration (t<sub>ingress</sub>, i.e. the time it takes the spot to grow from zero size to α<sub>max</sub>), egress duration
(t<sub>egress</sub>, i.e. time between spot α<sub>max</sub> to zero size), spot-to-star flux ratio (f<sub>spot</sub>), latitude (Φ), and longitude at T<sub>max</sub> (longitude zero is defined at
the stellar surface intersection with the line of sight).

The number of spots was chosen so that the periodic variability in the light curve is explained by the model. The process is based on
plot on the upper panel of Fig. B1, where periodic flux oscillations can be seen with increasing and decreasing amplitude. We placed spots
in an iterative manner. We start by placing a spot to represent the highest amplitude oscillation and roughly adjust its parameters. Then we
computed the residual plot and repeated the operation until no obvious periodic oscillation was seen in the residuals. Of course, this is a
strongly subjective method that should be automatized for a more rigorous study. We found that most of the light-curve variability can be
reproduced using 39 spots. At this first stage, we observed that, to reproduce the shape of the variability, i* must be close to 90°, and that
different spots had different periodicity pointing to a differential rotation stellar surface.

Also, as each quarter has a different flux level, one needs to use a physical evolution of the spots to normalize all quarter light curves to
a common flux level. The normalized light curve is presented in the top panel of Fig. B1. At each iteration of the modelling, each chunk
(defined between time gaps larger than 0.5 d) is normalized with a line to be compared with the model.

With 39 spots, the model has 320 free parameters including an additional white noise term for the data (jitter). To find the optimal model
parameters, we run one-thousand MCMC chains (the algorithm is described in Díaz et al. 2014) randomly started around the values estimated
in the previous step. At each MCMC step, there is a 10 per cent probability that a spot swap latitude sign, as this parameter is hemisphere
degenerated for i* = 90°. The chains evolved independently to different local maxima in parameter space. Then, 650 emcee walkers were
started around the global maximum found in the previous step. Here, we did not allowed for swapping of the spot latitude any more. We
ran 1.9 × 10⁶ steps of the emcee algorithm and considered only the last 20 000 steps for the final inference. Our results are summarized in
Table B1 and the posteriors of the parameters are shown in the lower panels of Fig. B1.

The amplitude of the residuals are up to one-tenth of the full variability and eight times the data uncertainty. This can be explained by a
number of reasons. The individual spots modelled can in fact represent groups of spots, with averaged umbra and penumbra. Not considering
umbra and penumbra separately, in addition to the assumption of circular spots and linear spot size evolution, may limit the ability of the
model to reproduce the light curve variability. Besides, in the complex parameter spaced defined by the model, the convergence of the Markov
chains cannot be guaranteed. In spite of all this, the current spot modelling allowed us to get information about the rotating stellar surface,
too complex to be understood by direct inference of the observed light curve.

<footnote n="15"><sup>15</sup>The rotational period at a latitude Φ is P<sub>Φ</sub> = P<sub>EQ</sub> / (1 - κ sin<sup>2</sup> Φ).</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure B1.** Spot modelling. From top to bottom, and left to right: MAP model of the spots (in red) compared to *Kepler* data (in black), the individual contribution of each spot in the model is shown in the upper part of the main panel (with different colours and numbered at the time of maximum spot-size), the residuals are shown in the lower panel. Spot coverage during the *Kepler* observations for the MAP model in grey, whereas in black the equivalent covered area for zero-temperature spots is shown. Posteriors distributions for spot latitude, spot-to-star flux ratio, maximum angular radius, lifetime, time span of the spots at maximum size, ingress, and egress, area, and equivalent area for zero-temperature spots. Correlations for the MAP model (the spot number is shown close to the points): spot-to-star flux ratio versus maximum angular radius, spot-to-star flux ratio versus lifetime, maximum angular radius versus lifetime, spot lifetime versus latitude, maximum angular radius versus latitude, and equivalent area for zero-temperature versus latitude.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure B2.** For each observed transit of Kepler-138b (Fig. A1), the spots on the visible surface of the star are plotted based on the MAP model of the analysis in Section 6.2.1, assuming a zero-projected spin-orbit angle. The lines in different colours represent the transit path. The size of the planet to scale is shown in the middle of the transit path.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure B3. Idem Fig. B2 but for Kepler-138c.</figcaption></figure>
<table><caption>**Table B1.** Parameters from spot modelling, with median and 68.3 per cent CI for the stellar parameters, and 99 per cent Highest Density Interval (HDI) for the spot parameters (to obtain this range, we merged all spots for each spot parameter).</caption><thead><tr><th>Parameter</th><th></th><th>68.3 per cent CI, 99 per cent HDI</th></tr></thead><tbody><tr><td>P<sub>EQ</sub></td><td>(d)</td><td>18.984± 0.050</td></tr><tr><td>P<sub>POLE</sub></td><td>(d)</td><td>20.707± 0.091</td></tr><tr><td>κ</td><td></td><td>0.0832<sup>+0.0043</sup><sub>-0.0080</sub></td></tr><tr><td>i<sub>*</sub></td><td>(°)</td><td>101.86± 0.95</td></tr><tr><td>u<sub>a</sub></td><td></td><td>0.401± 0.053</td></tr><tr><td>u<sub>b</sub></td><td></td><td>-0.157 ± 0.050</td></tr><tr><td>u<sub>a, spot</sub></td><td></td><td>0.042<sup>+0.088</sup><sub>-0.031</sub></td></tr><tr><td>u<sub>b, spot</sub></td><td></td><td>0.043± 0.055</td></tr><tr><td>α<sub>max</sub></td><td>(°)</td><td>[2.6, 13]</td></tr><tr><td>f<sub>spot</sub></td><td></td><td>[0, 0.88]</td></tr><tr><td>t<sub>ingress</sub></td><td>(d)</td><td>[0, 360]</td></tr><tr><td>t<sub>max</sub></td><td>(d)</td><td>[0, 120]</td></tr><tr><td>t<sub>egress</sub></td><td>(d)</td><td>[0, 590]</td></tr><tr><td>Lifetime<sup>a</sup></td><td>(d)</td><td>[19, 920]</td></tr><tr><td>jitter</td><td></td><td>7.87± 0.14</td></tr></tbody></table>

Notes. <sup>a</sup> t<sub>ingress</sub>+t<sub>max</sub>+t<sub>egress</sub>
<figure><img src="image_20.png" /><figcaption>Figure B4. Idem Fig. B2 but for Kepler-138 d.</figcaption></figure>

This paper has been typeset from a *TEX/LATEX* file prepared by the author.