# Refining the Transit-timing and Photometric Analysis of TRAPPIST-1: Masses, Radii, Densities, Dynamics, and Ephemerides

Eric Agol, Caroline Dorn, Simon Grimm, Martin Turbet, Elsa Ducrot, Laetitia Delrez, Michaël Gillon, Brice-Olivier Demory, Artem Burdanov, Khalid Barkaoui, et al.

► To cite this version:

Eric Agol, Caroline Dorn, Simon Grimm, Martin Turbet, Elsa Ducrot, et al. Refining the Transit-timing and Photometric Analysis of TRAPPIST-1: Masses, Radii, Densities, Dynamics, and Ephemerides. The Planetary Science Journal, 2021, 2 (1), pp.1. 10.3847/PSJ/abd022. hal-03452484

HAL Id: hal-03452484

https://hal.science/hal-03452484v1

Submitted on 11 Jul 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Refining the Transit-timing and Photometric Analysis of TRAPPIST-1: Masses, Radii, Densities, Dynamics, and Ephemerides

Eric Agol¹, Caroline Dorn², Simon L. Grimm³, Martin Turbet⁴, Elsa Ducrot⁵, Laetitia Delrez⁴,⁵,⁶, Michaël Gillon⁵, Brice-Olivier Demory³, Artem Burdanov⁷, Khalid Barkaoui⁵,⁸, Zouhair Benkhaldoun⁸, Emeline Bolmont⁴, Adam Burgasser⁹, Sean Carey¹⁰, Julien de Wit⁷, Daniel Fabrycky¹¹, Daniel Foreman-Mackey¹², Jonas Haldemann¹³, David M. Hernandez¹⁴, James Ingalls¹⁰, Emmanuel Jehin⁶, Zachary Langford¹, Jérémy Leconte¹⁵, Susan M. Lederer¹⁶, Rodrigo Luger¹², Renu Malhotra¹⁷, Victoria S. Meadows¹, Brett M. Morris³, Francisco J. Pozuelos⁵,⁶, Didier Queloz¹⁸, Sean N. Raymond¹⁵, Franck Selsis¹⁵, Marko Sestovic³, Amaury H. M. J. Triaud¹⁹, and Valerie Van Grootel⁶

¹ Astronomy Department and Virtual Planetary Laboratory, University of Washington, Seattle, WA 98195, USA; agol@uw.edu

² University of Zurich, Institute of Computational Sciences, Winterthurerstrasse 190, CH-8057, Zurich, Switzerland

³ Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, CH-3012, Bern, Switzerland

⁴ Observatoire de Genève, Université de Genève, 51 Chemin des Maillettes, CH-1290 Sauverny, Switzerland

⁵ Astrobiology Research Unit, Université de Liège, Allée du 6 Août 19C, B-4000 Liège, Belgium

⁶ Space Sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, Allée du 6 Août 19C, B-4000 Liège, Belgium

⁷ Department of Earth, Atmospheric and Planetary Science, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA

⁸ Oukaimeden Observatory, High Energy Physics and Astrophysics Laboratory, Cadi Ayyad University, Marrakech, Morocco

⁹ Center for Astrophysics and Space Science, University of California San Diego, La Jolla, CA 92093, USA

¹⁰ IPAC, California Institute of Technology, 1200 E California Boulevard, Mail Code 314-6, Pasadena, CA 91125, USA

¹¹ Department of Astronomy & Astrophysics, University of Chicago, 5640 S. Ellis Avenue, Chicago, IL 60637, USA

¹² Center for Computational Astrophysics, Flatiron Institute, Simons Foundation, 162 5th Avenue, New York, NY 10010, USA

¹³ University of Bern, Gesellschaftsstrasse 6, CH-3012, Bern, Switzerland

¹⁴ Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, MS 51, Cambridge, MA 02138, USA

¹⁵ Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, Allée Geoffroy Saint-Hilaire, F-33615 Pessac, France

¹⁶ NASA Johnson Space Center, 2101 NASA Parkway, Houston TX 77058, USA

¹⁷ Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA

¹⁸ Cavendish Laboratory, JJ Thomson Avenue, Cambridge, CB3 0H3, UK

¹⁹ School of Physics & Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK

Received 2020 June 26; revised 2020 September 30; accepted 2020 November 3; published 2021 January 22

## Abstract

We have collected transit times for the TRAPPIST-1 system with the Spitzer Space Telescope over four years. We add to these ground-based, HST, and K2 transit-time measurements, and revisit an *N*-body dynamical analysis of the seven-planet system using our complete set of times from which we refine the mass ratios of the planets to the star. We next carry out a photodynamical analysis of the Spitzer light curves to derive the density of the host star and the planet densities. We find that all seven planets' densities may be described with a single rocky mass-radius relation which is depleted in iron relative to Earth, with Fe 21 wt% versus 32 wt% for Earth, and otherwise Earth-like in composition. Alternatively, the planets may have an Earth-like composition but enhanced in light elements, such as a surface water layer or a core-free structure with oxidized iron in the mantle. We measure planet masses to a precision of 3%-5%, equivalent to a radial-velocity (RV) precision of 2.5 cm s⁻¹, or two orders of magnitude more precise than current RV capabilities. We find the eccentricities of the planets are very small, the orbits are extremely coplanar, and the system is stable on 10 Myr timescales. We find evidence of infrequent timing outliers, which we cannot explain with an eighth planet; we instead account for the outliers using a robust likelihood function. We forecast JWST timing observations and speculate on possible implications of the planet densities for the formation, migration, and evolution of the planet system.

*Unified Astronomy Thesaurus concepts: Extrasolar rocky planets (511); Exoplanet dynamics (490); Infrared photometry (792); Habitable planets (695); Transit timing variation method (1710); Transit photometry (1709); Exoplanet astronomy (486); Planetary interior (1248); Fundamental parameters of stars (555); Markov chain Monte Carlo (1889); N-body simulations (1083); Few-body systems (531)*

*Supporting material: machine-readable tables*

## 1. Introduction

The TRAPPIST-1 planetary system took the exoplanet community by surprise thanks to the high multiplicity of small transiting planets orbiting a very-low-mass star (~0.09M<sub>☉</sub>;

Gillon et al. 2016, 2017; Luger et al. 2017b; Van Grootel et al. 2018). This unexpected nature stems from the fact that this system was found in a survey of only 50 nearby ultracool dwarf stars (Jehin et al. 2011; Gillon et al. 2013), suggesting either a high frequency of such systems around the latest of the M dwarfs (He et al. 2016), or perhaps, this discovery was fortuitous (Sagear et al. 2020; Sestovic & Demory 2020). The proximity of the host star (~12 pc) makes it brighter in the infrared (*K* = 10.3) than most ultracool dwarfs. Its small size

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
(~0.12R<sub>☉</sub>) means that its planets' masses and radii are large relative to those of the star, which enables precise characterization of the planets' properties. The system provides the first opportunity for a detailed study of potentially rocky, Earth-sized exoplanets with incident fluxes spanning the range of the terrestrial planets in our solar system. As such, it has galvanized the exoplanet community to study this system in detail, both observationally and theoretically, and has fueled hopes that atmospheric signatures (or even biosignatures) might be detected with the James Webb Space Telescope (JWST; Barstow & Irwin 2016; Morley et al. 2017; Batalha et al. 2018; Krissansen-Totton et al. 2018; Fauchez et al. 2019; Lustig-Yaeger et al. 2019; Wunderlich et al. 2019).

More conservatively, the system provides a potential laboratory for comparative planetology of terrestrial planets and may provide insight and constraints on the formation and evolution of terrestrial planets around the lowest-mass stars. In particular, transiting multiplanet systems afford an opportunity to constrain the interior compositions of exoplanets. Sizes from transit depths combined with masses from transit-timing variations (TTVs) yield the densities of the planets (e.g., Agol & Fabrycky 2018). In the case of rocky planets with a thin atmosphere, the bulk density can constrain the core-mass fraction (CMF) and/or Mg/Fe mass ratio (Valencia et al. 2007), although for any given planet, there is still a degeneracy between a larger CMF and a volatile envelope (Dorn et al. 2018). In a multiplanet system, the bulk density as a function of planet orbital distance may be used to partly break the compositional degeneracy by assuming a common refractory composition and a water composition that increases with orbital distance (Unterborn et al. 2018; Lichtenberg et al. 2019).

The TRAPPIST-1 system was initially found with a ground-
based pilot survey using a 60 cm telescope, revealing two short-
period transiting planets and two additional orphan transits
(Gillon et al. 2016; Burdanov et al. 2018). Subsequent ground-
based study of the system revealed several additional orphan
transits, leading to an incomplete picture of the number of
planets and the architecture of the system. A 20 day observation
campaign with the Spitzer Space Telescope (Werner et al. 2004)
resolved the confusion, revealing the periods of six of the seven
transiting planets (Gillon et al. 2017), while only a single transit
observed of the outermost planet left its orbit in question. A
subsequent observation campaign of the system with the K2
mission included four additional transits of the outer planet,
identifying its period, and revealed a series of generalized three-
body Laplace relations (GLRs)<fnref n="20" /> between adjacent triplets of
planets (Luger et al. 2017b). Additional observations with
Spitzer continued to monitor the transit times of the seven
planets at a higher precision than afforded by ground-based
observations. An initial analysis of the Spitzer data to
determine planetary radii and masses was presented in Delrez
et al. (2018b) and Grimm et al. (2018). In total, Spitzer
observed TRAPPIST-1 for more than 1075 hr (nearly 45 days),
and the resulting time-series photometry includes 188 transits
(Ducrot et al. 2020). In this paper, we complement and extend
the analysis of Ducrot et al. (2020) to include a transit-timing
and photodynamic analysis of the system.

Although the planets in the TRAPPIST-1 system have short
orbital periods, ranging from 1.5 to 19 days, the dynamical

interactions accumulate gradually with time, which requires
longer-timescale monitoring to accurately constrain the orbital
model. The GLRs also cause adjacent pairs of planets to reside
near mean-motion resonances, for which $jP_i^{-1} \approx (j+k)P_{i+1}^{-1}$
for integers $j$ and $k$ for the $i$th and $(i+1)$th planets. This
proximity causes a resonant timescale for $k=1$ given by

$$
P_{\text{TTV}} = \frac{1}{jP_i^{-1} - (j+1)P_{i+1}^{-1}} \quad (1)
$$

(Lithwick et al. 2012), which is the characteristic timescale of
the TTVs of the outer five planets. The period of the resonant
terms for each of these pairs of planets is $P_{\text{TTV}} \approx 491 \pm 5$
days (ranging from 485 to 500 days for each pair). This
timescale has two consequences for measuring the masses of
the planets from TTVs: (1) the transit times for each planet
need to be sampled on this timescale, preferably covering two
cycles so that the amplitude and phase of the cycles may be
distinguished from the planets’ orbital periods; (2) this resonant
period also sets the timescale for the amplitude variability of
“chopping” (short-timescale TTVs), which can help to break a
degeneracy between mass and eccentricity for the resonant
terms (Lithwick et al. 2012; Deck & Agol 2015). As a
consequence, we expect the measurements of the masses of the
system to require sampling on a timescale of
$t_{\text{min}} \approx 2P_{\text{TTV}} \approx 2.7$ yr. The current paper is the first with a
survey time, $t_{\text{survey}} = 4.114$ yr, such that $t_{\text{survey}} > t_{\text{min}}$ for the
TRAPPIST-1 system.

Prior studies used the data available at the time (Delrez et al. 2018b), with $t_{\text{survey}} < t_{\text{min}}$, causing ample degeneracy in the dynamical model and hence larger uncertainties in the masses of the planets (Gillon et al. 2017; Grimm et al. 2018). Even so, these papers were ground-breaking as they enable the first density determinations of temperate, Earth-sized planets exterior to the solar system. Both papers indicated densities for the planets that were lower than the value expected for an Earth-like composition (with the exception of planet e), indicating that these planets might have significant volatile content. However, these conclusions were subject to significant uncertainty in the planet masses, making the determination of the compositions less definitive as the uncertainties were still consistent with rocky bodies at the $1\sigma-2\sigma$ level. In addition, the masses of all of the planets are highly correlated due to the fact that the dynamical state of all of the planets needs to be solved together and their masses and radii are measured relative to the star, so model comparisons with individual planets are not independent.

In this paper, we revisit a transit-timing and photometric
analysis with the completed Spitzer program using the more
extensive transit data set we now have in hand. The goal of this
program is to provide a more precise understanding of the
masses, radii, and densities of the planets. These measurements
may be used for planetary science with the extrasolar planets in
the TRAPPIST-1 system, whose similarity to the sizes, masses,
and effective insolation range of the terrestrial planets in our
solar system is the closest match known. In addition, we refine
the dynamical state of the system, revisiting some of the
questions explored in Grimm et al. (2018). Our final goal is to
prepare for upcoming observations with the JWST (Gardner
et al. 2006). More precise constraints on the parameters of the
planets will not only improve the precision with which we can
schedule observations, but also provide the best possible

<footnote n="20">²⁰ This refers to the condition $pP_i^{-1} - (p+q)P_2^{-1} + qP_3^{-1} \approx 0$, which is a generalization of the Laplace resonance with $p=1$ and $q=2$ (Papaloizou 2014).</footnote>
predictions of the potential environmental characteristics that could be discriminated observationally. This work will therefore help to optimize both the acquisition and interpretation of observations of the TRAPPIST-1 system with JWST.

In Section 2, we summarize the observational data which are analyzed in this paper. In Section 3, we discuss the nature of transit-timing outliers and the robust likelihood function we use for characterizing the system. This is followed by a description of our *N*-body transit-timing analysis in Section 4. With the improved *N*-body model, we revisit the photometric fit to the Spitzer data using a photodynamical model in Section 5. The results of these two analyses are combined to obtain the planet bulk properties in Section 6. In Section 7, we derive revised parameters for the host star. In Section 8, we search for an eighth planet with transit timing. In Section 9, we interpret the mass–radius measurements for the planets in terms of the interior and atmospheric structure models. The discussion and conclusions are given in Sections 10 and 11.

We provide Julia, Python, and Matlab code for running the Markov chains, creating the figures, and creating the paper PDF in https://github.com/ericagol/TRAPPIST1_Spitzer. The 3.5 GB data/ directory in the repository may be found at doi:10.5281/zenodo.4060252. In each figure, we embed links to the code (<>/) that produced that figure.

## 2. New TRAPPIST-1 Observations

Since the work described in Grimm et al. (2018) based on 284 transits, we have added an additional 163 transit times from a combination of Spitzer (Section 2.1) and ground-based observations (Section 2.2) for a total of 447 transits. With preliminary transit-timing fits, we found evidence for outliers among the measured times (Section 3), which we account for with a robust likelihood model. Each transit time is measured as a Barycentric Julian Date (BJD<sub>TDB</sub>), correcting for the location of Earth/spacecraft relative to the solar system barycenter (Eastman et al. 2010) at the time of each transit observation. We next describe our data.

### 2.1. Spitzer Observations

The data set used in this work includes the entire photometry database of TRAPPIST-1 observations with Spitzer Space Telescope’s Infrared Array Camera (IRAC; Carey et al. 2004) since the discovery of its planetary system. This represents all time-series observations gathered within the DDT programs 12126 (PI: M. Gillon), 13175 (PI: L. Delrez) and 14223 (PI: E. Agol). These cover a total of 188 transits observed from 2016 February to 2019 October and include 64, 47, 23, 18, 16, 13, and 7 transits of planets b, c, d, e, f, g, and h, respectively (Ducrot et al. 2020). All of these data can be accessed through the online Spitzer Heritage Archive database.<fnref n="21" /> Spitzer IRAC Channels 1 (3.6 µm, 0.75 µm wide) and 2 (4.5 µm, 1.015 µm wide) were used during the Spitzer Warm Mission (Fazio et al. 2004; Storrie-Lombardi & Dodd 2010) with 61 and 127 transits observed in each band, respectively. Observations were obtained with IRAC in subarray mode (32 × 32 pixel windowing of the detector) with an exposure time of 1.92 s and a cadence of 2.02 s. In order to minimize the pixel-phase effect (Knutson et al. 2008), the peak-up mode was used (Ingalls et al. 2016) to fine-tune the positioning of the target on the detector

following the IRAC Instrument Handbook.<fnref n="22" /> Finally, calibration was performed using Spitzer pipeline S19.2.0 to output data as cubes of 64 subarray images of 32 × 32 pixels (the pixel scale being 1′′2). Each set of exposures was summed over a 2.15 minute cadence to allow for a tractable data volume for carrying out the photometric analysis, which is described in detail in Delrez et al. (2018b) and Ducrot et al. (2020).

### 2.2. Ground-based Observations

In addition to the new Spitzer times, 125 transits were observed by the SPECULOOS-South Observatory at Cerro Paranal, Chile (SSO; Burdanov et al. 2018; Delrez et al. 2018a; Gillon 2018; Jehin et al. 2018), TRAPPIST-South at La Silla Chile, (TS; Jehin et al. 2011; Gillon et al. 2011), and TRAPPIST-North at Oukaïmeden, Morocco, (TN; Barkaoui et al. 2019). These observations were carried out in an *I* + *z* filter with exposure times 23 s, 50 s, and 50 s, respectively; the characteristics of this filter are described in Murray et al. (2020). Observations were also performed with the Liverpool Telescope (LT; Steele et al. 2004) and the William-Herschel Telescope (WHT), both installed at the Roque de los Muchachos Observatory, La Palma. Only one transit of planet b and one of d were targeted with the WHT, whereas 15 transits of several planets were targeted with LT. For LT observations, the IO:O optical wide-field camera was used in the Sloan *z'* band with a 20 s exposure time. One transit of b was observed with the Himalayan Chandra Telescope (HCT). Finally, a total of 26 transits were observed in the near-IR (1.2–2.1 µm) with the WFCAM near-IR imager of the United Kingdom InfraRed Telescope (UKIRT; Casali et al. 2007), the IRIS2IR-imager installed on the Anglo-Australian Telescope (AAT; Tinney et al. 2004), and the HAWK-I cryogenic wide-field imager installed on Unit Telescope 4 (Yepun) of the ESO Very Large Telescope (VLT; Siebenmorgen et al. 2011). These observations are summarized in Table 1: 504 transit observations were collected with 57 duplicate (or triplicate) transits, which were observed by a second (or third) observatory simultaneously, for a total of 447 unique planetary transit times which are used in our analysis. Additional information may be found in Gillon et al. (2016) for WHT and TRAPPIST, in Ducrot et al. (2018) for SSO and LT, and in Gillon et al. (2017) and Burdanov et al. (2019) for AAT, UKIRT, and VLT.

For all ground-based observations, a standard calibration (bias, dark, and flat-field correction) was applied to each image, and fluxes were measured for the stars in the field with the DAOPHOT aperture photometry software (Stetson 1987). Differential photometry was then performed using an algorithm developed by Murray et al. (2020) to automatically choose and combine multiple comparison stars, optimized to use as many stars as possible, weighted appropriately (accounting for variability, color and distance to target star), to reduce the noise levels in the final differential light curves. This reduction and photometry was followed by an Markov Chain Monte Carlo (MCMC) analysis to retrieve transit parameters.

### 2.3. K2 and HST Observations

The K2 mission (Howell et al. 2014) observed the TRAPPIST-1 system over campaigns 12 and 19 (Luger et al. 2017b), in both long- and short-cadence imaging modes. We

<footnote n="21">²¹ http://sha.ipac.caltech.edu</footnote>

<footnote n="22">²² https://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/iracinstrumenthandbook/</footnote>
<table><caption>**Table 1**
Number of Transits from Ground-based and Space-based Observations</caption><thead><tr><th>Planet</th><th>HCT</th><th>SSO/TS/TN</th><th>LT</th><th>WHT</th><th>VLT/AAT/UKIRT</th><th>HST</th><th>Spitzer</th><th>K2</th><th>Duplicates</th><th>Total (N<sub>i</sub>)</th></tr></thead><tbody><tr><td>b</td><td>1</td><td>45</td><td>7</td><td>1</td><td>10</td><td>1</td><td>64</td><td>48</td><td>17</td><td>160</td></tr><tr><td>c</td><td>0</td><td>28</td><td>8</td><td>0</td><td>7</td><td>1</td><td>47</td><td>30</td><td>14</td><td>107</td></tr><tr><td>d</td><td>0</td><td>11</td><td>1</td><td>1</td><td>5</td><td>2</td><td>23</td><td>17</td><td>7</td><td>53</td></tr><tr><td>e</td><td>0</td><td>18</td><td>4</td><td>0</td><td>3</td><td>2</td><td>18</td><td>11</td><td>7</td><td>49</td></tr><tr><td>f</td><td>0</td><td>9</td><td>2</td><td>0</td><td>4</td><td>2</td><td>16</td><td>7</td><td>6</td><td>34</td></tr><tr><td>g</td><td>0</td><td>11</td><td>0</td><td>0</td><td>3</td><td>2</td><td>13</td><td>5</td><td>4</td><td>30</td></tr><tr><td>h</td><td>0</td><td>3</td><td>2</td><td>0</td><td>0</td><td>0</td><td>7</td><td>4</td><td>2</td><td>14</td></tr><tr><td>Total</td><td>1</td><td>125</td><td>24</td><td>2</td><td>32</td><td>10</td><td>188</td><td>122</td><td>57</td><td>447</td></tr></tbody></table>

**Note.** Duplicates indicate the excess planet transits observed simultaneously with two or three distinct observatories (as indicated in Table 14). Details on the corresponding observations can be found in Gillon et al. (2016, 2017), Grimm et al. (2018), de Wit et al. (2016, 2018), Delrez et al. (2018b), Ducrot et al. (2018, 2020), and Burdanov et al. (2019).

only use the short-cadence data from campaign 12 for this analysis, with ~1 minute sampling. We use our own photometric pipeline to track the star and produce a light curve from the Target Pixel Files (TPF). To model and correct TRAPPIST-1's stellar variability and K2's pointing-drift-correlated systematic noise, we use a Gaussian process with a quasi-periodic kernel, following the procedure described in Grimm et al. (2018). The campaign 12 data contain 48, 30, 17, 11, 7, 5, and 4 transits of planets b, c, d, e, f, g, and h, respectively.

Transit times for Hubble Space Telescope observations were utilized, as described in Grimm et al. (2018), de Wit et al. (2016, 2018), and Wakeford et al. (2019).

## 2.4. Transit-time Measurements and Analysis

Gathering together the heterogeneous sample of transits obtained from a variety of ground- and space-based telescopes, we transformed the time stamps to the BJD<sub>TDB</sub> time standard prior to photometric analysis. We analyzed the data sets together with a global photometric analysis of all single-planet transits, as described in Ducrot et al. (2020), with a separate analysis of the overlapping transits once the single-transit analysis was completed.

For each planet, a fixed time of transit for epoch zero ($T_0$) and fixed period ($P$) were used, but with timing offset ("TTV") as a fitted parameter for each transit as described by Ducrot et al. (2020). To derive $T_0$ and $P$, a linear regression of the timings as a function of their epochs was performed for each planet to derive an updated mean transit ephemeris; their exact values can be found in Table 4 of Ducrot et al. (2020). The timing offsets are then added back to the ephemeris to obtain the measured transit times and uncertainties.

The final observed data set for the transit-timing analysis is given by $y = \{t_{obs,ij}, \sigma_{ij}, j = 1, ..., N_i\}$; $i = 1, ..., 7)$, where $i$ labels each of the seven planets, $N_i$ is the number of transits for the $i$th planet (Table 1), and $j$ labels each transit for the $j$th planet, so that $t_{obs,ij}$ is the $j$th observation of the $i$th planet, and $\sigma_{ij}$ is the corresponding measurement error. The total number of transits is $N_{trans} = \sum_{i=1}^{N_p} N_i = 447$, where $N_p$ is the number of transiting planets.

Table 14 lists the complete set of transit times and uncertainties, which were utilized in the present analysis.

With this sample of transit times collected, we proceed to describe our dynamical analysis, starting with the likelihood function and evidence for outliers.

## 3. Excess of Outliers and Robust Likelihood Model

We first carried out a preliminary seven-planet, plane-parallel N-body model fit to the transit times using a $\chi^2$ log likelihood function, i.e., assuming a Gaussian uncertainty for each transit time given by the derived timing uncertainty, which we optimized using the Levenberg-Marquardt algorithm. We found that the residuals of the fit contain many more outliers than is probable assuming a Gaussian distribution for the timing uncertainties.

Figure 1 shows the cumulative distribution function (CDF) and a histogram of the normalized residuals versus a single Gaussian probability distribution function (PDF) with unit variance (orange line). This CDF disagrees with the Gaussian CDF in the wings for $P(>z) \lesssim 0.1$ and $P(>z) \gtrsim 0.9$, where $z = (t_{obs,ij} - t_{ij}(x_{dyn}))/\sigma_{ij}$ are the normalized residuals, with the model time, $t_{ij}(x_{dyn})$, as a function of the dynamical model parameters, $x_{dyn}$, described below. This indicates that there is a significant excess of outliers with large values of $|z|$ relative to a Gaussian distribution. The histogram in Figure 1 also demonstrates this clearly: there are eight data points with $z < -3$ and 7 with $z > 3$. With 447 transit-time measurements, we would only expect $\approx 1.2$ data points with $|z| > 3$ if the distribution were Gaussian with accurately estimated uncertainties. This excess is even more apparent at $|z| > 4$.

We have examined the individual transits that show these discrepancies, and there is nothing unusual about their light curves, such as flares, overlapping transits, or other anomalies. The outliers appear for each of the planets (save h), in both ground- and space-based data, and for measurements with different sizes of uncertainties. We do not think that our N-body model is in error (and we have tried to fit with an extra planet, without a significant improvement in the number of outliers; see Section 8). Consequently, we believe that these outliers are due to variations in the measured times of transits which are not associated with the dynamics of the system.

We suspect instead that these outliers are a result of some systematic error(s) present in the data. There are a variety of possibilities: uncorrected instrumental/observational systematics, time-correlated noise due to stellar variability, stellar flares (which may be too weak to be visible by eye, but might still affect the times of transit), or stellar spots (Oshagh et al. 2013; Ioannidis et al. 2015). Again, our examination of the light curves did not point to a single culprit, so we are unable to model and/or correct for any of these effects. Our data are not unique in this respect: similar outliers have been seen in other transit-timing analyses, as described in Jontof-Hutter et al. (2016).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Probability distribution of normalized residuals. Left: cumulative distribution function of the normalized residuals, *z*. The blue and brown lines are a sequence of normalized residuals. The orange line is the CDF of a Gaussian distribution. The dotted green line is the CDF of a Student's t-distribution. Right: histogram of the normalized residuals. The blue and brown data points are a histogram of the normalized residuals with Poisson uncertainties. The other lines have the same meaning as the left panel for the probability distribution function (PDF), scaled to match the histogram. In both panels, the >3σ outliers are indicated in brown. </></figcaption></figure>

Our transit-timing model will be affected by these timing outliers, which make an excessive contribution to the $\chi^2$ of the model, and thus can affect the inference of the model parameters. This can cause both the parameters and the uncertainties to be misestimated. To make progress, we have modified the likelihood model to account for outliers.

We use a heavy-tailed likelihood function, which better describes the residual distribution: a Student's t-distribution (Jontof-Hutter et al. 2016). We fit the normalized residuals to a model in which the width of the distribution was allowed to vary, which we parameterize with an additional factor multiplying the variance, which we refer to below as $V_1$. For the Student's t-distribution, there is only one additional free parameter: the number of degrees of freedom, $\nu$, which we treat as a continuous parameter.

Figure 1 shows a histogram of the outliers of the best-fit transit-timing model (described below) and shows that the Student's t-distribution gives a much higher probability for outliers.

With the description of the data set complete, we next describe our efforts to model the data.

# 4. Transit-timing Analysis

In this section, we describe our transit-timing analysis in detail, starting first with a description of our dynamical model.

## 4.1. N-body Integration

We integrate the N-body dynamics in Cartesian coordinates with a novel symplectic integrator, NbodyGradient, which is based on the algorithm originally described in Hernandez & Bertschinger (2015), derived from the nonsymplectic operator

of Gonçalves Ferrari et al. (2014).<fnref n="23" /> The time-evolution operator of the integrator is a succession of Kepler two-body problems and simple “kick” and “drift” operators. The advantage over traditional symplectic methods (Wisdom & Holman 1991) is that the dominant error is due to three-body interactions, while in the standard methods, the dominant error is due to two-body interactions, meaning close encounters between nonstellar bodies are treated poorly (Hernandez & Dehnen 2017). The Kepler problem for each pair is solved with an efficient universal Kepler solver (Wisdom & Hernandez 2015). The symplectic integrator is made to be time-symmetric to yield second-order accuracy (Hernandez & Bertschinger 2015). Then, a simple operator is introduced to double the order of the method (Dehnen & Hernandez 2017). We have found that numerical cancellations occur between Kepler steps and negative drift operators, and so we have introduced an analytic cancellation of these terms to yield an algorithm that is numerically stable, which converges for small time steps (E. Agol & D. M. Hernandez 2021, in preparation).

The initial conditions are specified with Jacobi coordinates (Hamers & Zwart 2016), and we use a set of orbital elements for each planet given by $x_{dyn} = (\{m_i, P_i, t_{0,i}, e_i \cos \omega_i, e_i \sin \omega_i\}; i = 1, ..., N_p)$, where $N_p$ is the number of planets for a total of $5N_p$ dynamical parameters. In addition, we take the star to have a mass, $m_0 = M_*/M_\odot$, which we fix to 1. The units of time for the code are days, while the length scale of the code is taken to be $m_0^{1/3}$ au.<fnref n="24" /> The initial orbital ephemeris, ($P_i$,

<footnote n="23">²³ The code may be found at https://github.com/ericagol/NbodyGradient.</footnote>

<footnote n="24">²⁴ Note that as we take $m_0 = 1$ in our simulations, we need to multiply the output of positions and velocities from the code by $(M_*/M_\odot)^{1/3}$ to scale to a stellar mass $M_*$.</footnote>
$t_{0,i}$), consists of the period and initial time of transit which each planet would have if it orbited a single body with a mass of the sum of the star and the interior planets, unperturbed by the exterior planets. We use these variables (in lieu of the initial semimajor axis and mean longitude) as they are well constrained by the observed times of transits. We convert these analytically to the time of periastron passage, once the Kepler equation is solved, to yield the initial eccentricity anomaly for each initial Keplerian. Finally, the eccentricity, $e_i$, and longitude of periastron, $\omega_i$, for each Keplerian, we parameterize in terms of $e_i \cos \omega_i$ and $e_i \sin \omega_i$ to avoid the wrapping of the angle $\omega_i$. We transform from Jacobi coordinates to Cartesian coordinates to complete the initial conditions.

For our transit-timing analysis, we assume that the planets are plane parallel and edge on in their orbits, allowing us to ignore the inclination and longitude of nodes for each planet.

A symplectic integration time step, $h$, is selected to be small, <5%, compared with the orbital period of the innermost planet (Wisdom & Holman 1991). For most of our integrations, we use a time step of $h = 0.06$ days, or about 4% of the orbital period of planet b.

The model transit times are found by tracking the positions of each planet relative to the star across a time step. Then, when the dot product of the relative velocity of the planet and star with their relative position goes from negative to positive, and the planet is between the star and observer, we flag a routine that iterates with Newton's method to find the model transit time, which is taken to be when this dot product equals zero (Fabrycky 2010), corresponding to the midpoint of the transit if acceleration is negligible over the duration of the transit. The resulting model we obtain is for the $j$th transit of the $i$th planet, giving each model transit time as a function of the initial conditions, $t_{ij}(\mathbf{x}_{dyn})$, which can then be compared to the observed times, $t_{obs,ij}$.

Once the model transit times have been found for every planet over the duration of the time integration, these are then matched with the observed transit times to compute the likelihood using Student's t-probability distribution. The log likelihood function for each data point is given by

$$ \log \mathcal{L}_{ij}(\mathbf{x}_{dyn}, \nu, V_i) = -\frac{\nu+1}{2} \log \left( 1 + \frac{(t_{obs,ij} - t_{ij}(\mathbf{x}_{dyn}))^2}{\nu \sigma_{ij}^2 V_i} \right) - \frac{1}{2} \log(\pi \nu V_i) + \log \Gamma \left( \frac{\nu+1}{2} \right) - \log \Gamma \left( \frac{\nu}{2} \right), \quad (2) $$

where $\Gamma(x)$ is the Gamma function (Fisher 1925).

The total log likelihood function which we optimize is given by

$$ \log \mathcal{L}(\mathbf{x}_{dyn}, \nu, V_i) = \sum_{i=1}^{N_p} \sum_{j=1}^{N_i} \log \mathcal{L}_{ij}(\mathbf{x}_{dyn}, \nu, V_i), \quad (3) $$

where $N_p$ is the number of planets; we use $N_p = 7$ for most of our analysis.

Note that we assume that the timing errors are uncorrelated. Most transits are well separated in time, and thus this is an accurate assumption as the noise should be uncorrelated on these timescales. There are a small number of transits (about 6%) that overlap in time and thus may have correlated uncertainties; we do not account for this in the likelihood function.

## 4.2. Uncertainty Analysis

We carried out the uncertainty analysis on the model parameters with three different approaches:

1. Laplace approximation

2. Likelihood profile

3. MCMC

First, in our Laplace approximation analysis, we assume a uniform prior on the model parameters and expand the likelihood as a multidimensional normal distribution. We maximize the likelihood model using the Levenberg-Marquardt algorithm, which requires the gradient and Hessian of the negative log likelihood. Once the maximum likelihood is found, we compute an approximate Hessian at the maximum likelihood (see Appendix A). The inverse of the Hessian matrix yields an estimate of the covariance among the parameters at the maximum likelihood, whose diagonal components provide an initial estimate for the parameter uncertainties; we will also use the Hessian for more efficient sampling of the Markov chain.

The second approach we use is to compute the likelihood profile for each model parameter. In this case, each parameter is varied over a grid of values over a range given by $\pm 3\sigma_{xi}$, where $\sigma_{xi}$ equals the square root of the diagonal component for the $i$th model parameter from the covariance matrix. At each value along the grid for each parameter, we optimize the likelihood with a constraint that keeps the parameter pinned at the grid point. This results in a profile of the maximum likelihood of each parameter, optimized with respect to all other parameters, which yields a second estimate for the uncertainties on the parameters. The likelihood-profile approach does not assume a normal distribution and is useful for checking for a multimodal probability distribution which can trip up Markov chain analysis.

However, both of these error estimates are incomplete as they do not account for nonlinear correlations between parameters, for the non-Gaussian shape of the posterior probability, nor for the prior probability distribution.<fnref n="25" /> Nevertheless, the agreement between the two estimates gives a starting point for evaluating our Markov chain analysis and for gauging the convergence of the chains, which we describe below.

In our initial Markov chain sampling, we found that the parameters of Student's t-distribution, $\nu$ and $V_i$, were strongly nonlinearly correlated and displayed a likelihood profile which was non-Gaussian. After experimenting with reparameterization, we found that $\log \nu$ and $V_i e^{1/(2\nu)}$ gave a parameterization that showed a nearly Gaussian likelihood profile in each parameter and also showed more linear correlations between these two parameters. Accordingly, we chose to sample in these transformed parameters so that our set of model parameters is $\mathbf{x} = (\mathbf{x}_{dyn}, \log \nu, V_i e^{1/(2\nu)})$.

In Appendix B, we define the prior function $\Pi(\mathbf{x})$, which multiplies the likelihood to give the posterior probability distribution,

$$ P(\mathbf{x}) \propto \Pi(\mathbf{x}) \times \mathcal{L}(\mathbf{x}), \quad (4) $$

<footnote n="25">²⁵ In principle, we could include a prior in the Laplace and likelihood-profile analyses.</footnote>
so that we can proceed to discussing the Markov chain sampling of the posterior probability of the model parameters given the data.

### 4.3. Markov Chain Sampler

We sample our posterior probability, $P(x)$, with a Markov chain sampler. There are 37 free parameters—4 orbital elements and 1 mass ratio for each planet, and 2 parameters for the Student's t-distribution. Given the high dimensionality of our model, we chose to use a Markov chain sampler that efficiently samples in high dimensions: Hamiltonian Monte Carlo (HMC; Duane et al. 1987; Neal 2011; Monnahan et al. 2016; Betancourt 2017).<fnref n="26" /> This sampler requires the gradient of the likelihood function with respect to the model parameters. The gradient of the likelihood requires the gradient of each model transit time with respect to the initial conditions of the N-body integrator.

We have written a module for our N-body integrator that computes the gradient of each model transit time by propagating a Jacobian for the positions and velocities of all bodies across every time step throughout the N-body integration (E. Agol & D. M. Hernandez 2021, in preparation). This is multiplied by the Jacobian of the coordinates at the initial time step computed with respect to the initial Keplerian elements and masses, which specify the initial conditions and comprise the N-body model parameters.

When a transit time is found during the N-body integration with NbodyGradient, we compute the derivative of each transit time with respect to the coordinates at the preceding time step, which we multiply with the Jacobian at that step to obtain the gradient of each transit time with respect to the initial conditions. The gradient of the prior with respect to the model parameters and the gradient of the likelihood with respect to the model times and the Student's t-distribution parameters are each computed with automatic differentiation, using forward-mode derivatives (Revels et al. 2016). The gradient of the likelihood with respect to the dynamical model parameters is found by applying the chain rule to the automatic derivatives of the likelihood with respect to the model times with the derivatives computed in the N-body model (from NbodyGradient).

For our HMC analysis, we augment the simulation parameters with a set of conjugate momenta, $p$, with the same dimension. We sample from the probability distribution, $e^{-H(x,p)}$, where $H$ is a Hamiltonian given by the negative log posterior,

$$H(p, x) = \frac{1}{2} p^T M^{-1} p - \log \mathcal{L}(x) - \log \Pi(x), \quad (5)$$

where $p$ is defined from Hamilton's equations,

$$\dot{p} = -\frac{\partial H}{\partial x}. \quad (6)$$

We take the mass matrix, $M$, to be the approximate Hessian matrix evaluated at the maximum likelihood, $M = \mathcal{H}(x_0)$

(Equation (A5)). Similarly, the Hamiltonian can be used to compute the evolution of the parameter “coordinates,”

$$\dot{x} = + \frac{\partial H}{\partial p}. \quad (7)$$

The dot represents the derivative with respect to an artificial “time” coordinate, which can be used to find a trajectory through the $(x, p)$ phase space that conserves the “energy” defined by this Hamiltonian.

We carry out a Markov chain using the standard approach for HMC. First, we draw the initial momentum from the multivariate Gaussian distribution defined by the kinetic energy term in the Hamiltonian,

$$p = M^{1/2} Z, \quad (8)$$

where $Z_n \sim N(0,1)$ is an element of a vector of random normal deviates for $n = 1, \dots, N_{\text{param}}$. We then carry out a leapfrog integration of Hamilton's equations for $N_{\text{leap}}$ steps from the starting point with a “time” step $\epsilon$ to obtain a proposal set of parameters $(x_{\text{prop}}, p_{\text{prop}})$. Because energy is not conserved precisely due to the finite differencing of the leapfrog integration, we then apply a Metropolis rejection step to accept the proposal step with probability

$$p_{\text{accept}} = \min(\exp(-(H(x_{\text{prop}}, p_{\text{prop}}) - H(x, p))), 1) \quad (9)$$

to determine whether to accept the proposed step and add it to the Markov chain, or to reject it and copy the prior step to the chain.

We carried out some trial integrations to tune two free parameters: $\epsilon_0$ and $N_{\text{leap},0}$. We draw the “time” step, $\epsilon$, for each integration from the absolute value of a normal distribution with width $\epsilon_0$, i.e., $\epsilon \sim |N(0, \epsilon_0)|$. We draw the number of leapfrog steps for each integration from a uniform probability, $N_{\text{leap}} \sim \text{round}(N_{\text{leap},0}\mathcal{U}(0.8, 1.0))$. We found that a choice of $\epsilon_0 = 0.1$ and $N_{\text{leap},0} = 20$ results in a proposal for which the Metropolis rejection gives a high average acceptance rate of 70%.

We ran 112 HMC chains for 2000 steps each (i.e., 2000 leapfrog integrations). Each leapfrog integration averaged about 7 minutes and so the chains took 9 days and 4 hours to complete.<fnref n="27" /> We found a minimum mean effective sample size of 57 over all chains for a total number of independent samples of 6409.

### 4.4. Results

The TTVs are shown in Figure 2, along with our best-fit model. The model is a very good description of the data, although a few outliers are clearly visible by eye. As advertised, the outer five planets show large-amplitude oscillations with the timescale $P_{\text{TTV}}$. We have created a second figure in which a polynomial with an order between 5 and 30 is fit and removed from the data, and the resulting differences are shown in Figure 3. The result shows high-frequency variations that are associated with the synodic periods of pairs of adjacent planets, typically referred to as “chopping.” The chopping TTVs encode the mass ratios of the companion planets to the star without the influence of the eccentricities and thus provide

<footnote n="26">²⁶ a.k.a. “Hybrid Monte Carlo.” Note that the “Hamiltonian” referred to in HMC is not a physical Hamiltonian but an artificial one used for treating the negative log probability as a potential energy function and adding a kinetic energy term, with an artificial momentum conjugate to each model parameter (“coordinate”). For a description of HMC and a discussion of applications to cosmology, including N body, see Leclercq et al. (2014) and Jasche & Kitaura (2010) and references therein.</footnote>

<footnote n="27">²⁷ These were run on 4 Broadwell Xeon Processors with 28 cores and 128 GB of memory, where each processor is a node in the Hyak Mox cluster at the University of Washington.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Transit-time variation measurements (orange/red error bars) and best-fit transit-time model (blue/green lines) for a subset of our Spitzer/K2/ground-based data set. The TTVs are the transit times for each planet with a best-fit linear ephemeris removed. Brown error bars indicate >3σ outliers. </></figcaption></figure>

a constraint on the planet-star mass ratios which are less influenced by degeneracies with the orbital elements (Deck & Agol 2015). The chopping variations are clearly detected for

each planet (except planet d), which contributes to the higher precision of the measurements of the planet masses in this paper.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Observed transit times with a polynomial subtracted (orange error bars) compared with the short-timescale chopping variations of the best-fit model (blue model; same polynomial removed). Green dots show the analytic chopping relation from Deck & Agol (2015) due to adjacent planets, also with a low-order polynomial removed. For the inner four planets, we have only plotted data with uncertainties smaller than the chopping semiamplitude (many observations have large uncertainties that would obscure the plot).</figcaption></figure>

The results of the posterior distribution of our transit-timing analysis are summarized in Table 2 with the mean and ±34.1% confidence intervals (1σ) computed from the standard deviation of the Markov chains. The correlations between parameters are depicted in Figure 29. There are 35 parameters that describe the planets, in addition to two parameters for the Student's t-distribution, log ν = 1.3609 ± 0.2337 and V<sub>1</sub>e<sup>1/(2ν)</sup> = 0.9688 ± 0.1166 (Figure 4). The posterior mass ratios and ephemerides are consistent with nearly Gaussian distributions. The eccentricity vectors show deviations from a Gaussian distribution for the inner two planets b and c, as shown in Figure 5. The Laplace approximation covariance uncertainty

estimates that are overplotted as Gaussian distributions very closely match the likelihood profile for each parameter. This agreement is reassuring: it indicates that the likelihood distribution is closely approximated by a multidimensional normal distribution near the maximum likelihood. In the eccentricity vector coordinates, the prior probability distribution is peaked at zero to ensure that the volume of phase space at larger eccentricities does not dominate the probability distribution, as shown in the lower-right panel of Figure 5. For the planets that have a likelihood distribution that overlaps strongly with zero, the prior distribution causes the Markov chain posterior to have a significantly different distribution
<table><caption>**Table 2**

Parameters of the TRAPPIST-1 System from Transit-timing Analysis and Their 1σ Uncertainties</caption><thead><tr><th></th><th>μ<br>=<sup>M<sub>⊕</sub></sup>/<sub>M<sub>⊕</sub></sub><br>(0.09M<sub>☉</sub> / M<sub>☉</sub>)</th><th>M<sub>p</sub><br>[10<sup>-5</sup>M<sub>☉</sub>]</th><th>σ<sub>μ</sub><br>μ</th><th>P<br>(day)</th><th>t<sub>0</sub><br>(BJD<sub>TDB</sub>-2,450,000)</th><th>e cos ω</th><th>e sin ω</th></tr></thead><tbody><tr><td>b</td><td>1.377 1 ± 0.0593</td><td>4.596 ± 0.198</td><td>4.3</td><td>1.510 826 ± 0.000006</td><td>7 257.550 44 ± 0.00015</td><td>-0.00215 ± 0.00332</td><td>0.002 17 ± 0.00244</td></tr><tr><td>c</td><td>1.310 5 ± 0.0453</td><td>4.374 ± 0.151</td><td>3.5</td><td>2.421 937 ± 0.000018</td><td>7 258.587 28 ± 0.00027</td><td>0.000 55 ± 0.00232</td><td>0.000 01 ± 0.00171</td></tr><tr><td>d</td><td>0.388 5 ± 0.0074</td><td>1.297 ± 0.025</td><td>1.9</td><td>4.049 219 ± 0.000026</td><td>7 257.067 68 ± 0.00067</td><td>-0.00496 ± 0.00186</td><td>0.002 67 ± 0.00112</td></tr><tr><td>e</td><td>0.693 2 ± 0.0128</td><td>2.313 ± 0.043</td><td>1.8</td><td>6.101 013 ± 0.000035</td><td>7 257.827 71 ± 0.00041</td><td>0.004 33 ± 0.00149</td><td>-0.00461 ± 0.00087</td></tr><tr><td>f</td><td>1.041 1 ± 0.0155</td><td>3.475 ± 0.052</td><td>1.5</td><td>9.207 540 ± 0.000032</td><td>7 257.074 26 ± 0.00085</td><td>-0.00840 ± 0.00130</td><td>-0.00051 ± 0.00087</td></tr><tr><td>g</td><td>1.323 8 ± 0.0171</td><td>4.418 ± 0.057</td><td>1.3</td><td>12.352 446 ± 0.000054</td><td>7 257.714 62 ± 0.00103</td><td>0.003 80 ± 0.00112</td><td>0.001 28 ± 0.00070</td></tr><tr><td>h</td><td>0.326 1 ± 0.0186</td><td>1.088 ± 0.062</td><td>5.7</td><td>18.772 866 ± 0.000214</td><td>7 249.606 76 ± 0.00272</td><td>-0.00365 ± 0.00077</td><td>-0.00002 ± 0.00044</td></tr></tbody></table>

Note. Note that the mass ratios, $\mu = M_p/M_*$, of the planets are computed relative to the star, which is assumed to have a mass of 0.09 $M_\odot$ (this is later combined with the estimate of stellar mass to give our estimates of the planet masses). We also report $\mu$ in units of $10^{-5}$, and the fractional precision on the measurement of $\mu$, $\sigma_\mu/\mu$. The parameters $P$, $t_0$, $e \cos \omega$, and $e \sin \omega$ describe the osculating Jacobi elements at the start of the simulation, on date BJD<sub>TDB</sub>-2,450,000 = 7257.93115525 days.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Likelihood profile (dark line) and Gaussian distribution with Laplace approximation uncertainty (light line) for $\log \nu$ (left) and $V_1 e^{1/(2\nu)}$ (right). The posterior probability distributions are shown with blue histograms. </></figcaption></figure>

from the likelihood profile. This is not due to the prior favoring small eccentricities; rather, it is simply a correction for the bias that results from using $e_i \cos \omega_i$ and $e_i \sin \omega_i$ as Markov chain parameters that favor higher eccentricities (Ford 2006).

The marginalized posterior distributions of the ratio of the planet masses to the star masses, scaled to a stellar mass of 0.09 $M_\odot$, are given in Table 2 and shown in Figure 6. The likelihood profile of the planet-to-star mass ratios is also plotted in Figure 6 and appears to be well behaved. These likelihood profiles are also approximately Gaussian in shape and track the inverse Hessian evaluated at the maximum likelihood to estimate the covariance (also plotted). Compared with the mass estimates from Grimm et al. (2018), the masses of each planet have increased with the exception of planet e, which has decreased, and planet h, which remains the same (Table 3). The mass ratios of the posterior distribution from the Markov chain are slightly shifted to smaller values than the likelihood profile and Laplace approximation probabilities for all planets save b and g.

The Student's t-distribution parameters show a posterior distribution that is shifted from the likelihood profile/Laplace probability distribution (Figure 4). This bias is due to the fact that the likelihood distribution of these parameters shifts

upwards whenever the transit-timing model parameters deviate from their maximum-likelihood values. The peak of the posterior distribution of these parameters corresponds to $\nu = 3.9$ and $V_1^{1/2} = 0.87$, which indicates that the core of the distribution is narrower than the transit-timing uncertainties indicate, while the wings of the distribution are close to $\nu = 4$, which was the value used by Jontof-Hutter et al. (2016).

### 4.5. Independent N-body TTV Analysis

In addition to the N-body code described above, we use the GPU hybrid symplectic N-body code GENGA (Grimm & Stadel 2014) with a Differential Evolution MCMC Method (DEMCMC; ter Braak 2006) as described in Grimm et al. (2018) to perform an independent TTV analysis. The parameters for the MCMC analysis are $\mathbf{x} = \{m_i, P_i, t_{0,i}, e_i, \omega_i\}; i = 1, \dots, N_p$. The mass of the star is taken to be $M_* = 0.09 M_\odot$, and the time step of the N-body integration is set to $h = 0.05$ days. The likelihood is assumed to be a normal distribution with the timing errors derived from the timing analyses. For comparison, we have rerun the likelihood-profile computation described above using a normal distribution in place of a Student's t-distribution. The derived masses from the two different analyses agree well with a maximal deviation of the median masses of better than 0.4%, while the mass-ratio uncertainties agree to better than 13%. The eccentricities and longitudes of periastron at the initial time agree as well. We interpret this as a validation of the numerical techniques being employed in this paper.

With the transit-timing analysis completed, we now use the N-body model to improve the estimate of the stellar density and the planet-to-star radius ratios. To do so we create a photodynamic model, described next.

### 5. Photodynamical Analysis

With the mass ratios and orbital parameters derived from the transit-timing analysis, we wish to improve our derivation of the planet and stellar parameters from Spitzer photometry. The transit depth, transit duration, and ingress/egress duration combined with the orbital period constrain the impact parameters and density of the star (Seager & Maller-Ornelas 2003). Combining these constraints for each of the planets enables a more precise constraint on the density of the star (Kipping et al. 2012). The transit durations are affected by the (small) eccentricities but to a lesser extent. We account for the dynamical constraints on the transit-timing model to
<figure><img src="image_8.png" /><figcaption>Figure 5. Eccentricity vector probability distribution for each planet (y-axes are the relative probability). Thick histograms are the marginalized posterior distributions from the Markov chain analysis. Thin light lines are the Laplace approximation. Thin dark lines are the likelihood profiles. The lower-right panel shows the distribution of e cos ω or e sin ω for a uniform prior on e ∈ U(0, 0.1) and ω ∈ U(0, 2π). </></figcaption></figure>

improve the photometric constraints on these parameters, albeit with the dynamical parameters fixed at the maximum likelihood.

We fit a “photodynamical” model (Carter et al. 2012) to the data with the following procedure. From the best-fit, plane-parallel, edge-on transit-time model, we compute the sky velocity at each of the midtransit times, t₀, from the model (in N-body code units). We then convert the code units to physical units using the density of the star, obtaining the sky velocity, v<sub>sky</sub>, in units of R<sub>*</sub>/day. We account for quadratic limb darkening of the star with parameters (q<sub>1,Ch1</sub>, q<sub>2,Ch1</sub>, q<sub>1,Ch2</sub>, q<sub>2,Ch2</sub>) in the two Spitzer channels, and for each planet, we specify a planet-to-star radius ratio (R<sub>p</sub>/R<sub>*</sub>) and we assume a midtransit impact parameter (b<sub>0</sub>), which is constant for all transits of a given planet. We assume that the limb-darkening parameters are a function of wavelength for the two Spitzer channels, while we treat the planet radius ratios as identical in both wave bands based on their consistency across all planets in Ducrot et al. (2020), giving a total of 19 free parameters for the photodynamical model.

We ignore acceleration during the transits, treating the impact parameters as a function of time as

$$b(t) = \sqrt{(v_{\text{sky}}(t - t_0))^2 + b_0^2} \quad (10)$$

in units of the stellar radius, R<sub>*</sub>. Although this expression ignores the curvature and inclination of the orbits, as well as the acceleration of the planet, the star is so small compared with the orbital radius that this approximation is extremely accurate. The transit model is integrated with an adaptive Simpson rule over each Spitzer exposure time (which has a uniform duration binned to 2.15 minutes), as described in Agol et al. (2020), yielding a light curve computed with a precision of better than 10<sup>-7</sup> for all cadences.

We compute a photometric model for all seven planets for all of the Spitzer data in selected windows around each of the observed transits. Starting with Spitzer photometric data, which were already corrected for systematic variations based on the analysis by Ducrot et al. (2020), we fit each transit window with the transit model multiplied by a cubic polynomial, whose coefficients are solved for
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Probability distribution of the planet-to-star mass ratios, scaled to a stellar mass of $M_* = 0.09M_\odot$; panels range from small masses to large. Thick histograms show the posterior probability distribution of the Markov chain analysis. Horizontal error bars show the mean and $1\sigma$ intervals for the mass ratios from Grimm et al. (2018). Dark solid bell curves are the likelihood profiles; light, dotted bell curves are the Laplace approximation. <>/</figcaption></figure>

via regression at each step in the Markov chain. We transform the $q_1, q_2$ limb-darkening parameters to $u_1, u_2$ in each band using the formalism of Kipping (2013) to compute the transit model from Agol et al. (2020). After carrying out an initial optimization of the model, we take the photometric error to be the scatter in each observation window to yield a reduced $\chi^2$ of unity in that window. With this photometric scatter, we compute a $\chi^2$ of the model with respect to the Spitzer photometric data, and we optimize the model using a Nelder-Mead algorithm.

## 5.1. Photodynamic Results

To compute the uncertainties on the photodynamical model parameters, we use an affine-invariant MCMC algorithm

(Goodman & Weare 2010).<fnref n="28" /> We used a uniform prior with bounds on each parameter given in Table 4. The posterior distributions of the results of the fit are given in Table 5, while the correlations between parameters are shown in Figure 30. We utilized 100 walkers run for 50,000 generations, discarding the first 1500 generations for burn in. We computed the effective sample size using the integrated autocorrelation length, finding a minimum effective sample size of 6000 over all 19 parameters.<fnref n="29" />

<footnote n="28">²⁸ As implemented in the package https://github.com/madsjulia/AffineInvariantMCMC.jl.</footnote>

<footnote n="29">²⁹ Using https://github.com/tpapp/MCMCDiagnostics.jl.</footnote>
<table><caption>**Table 3**

Planet-to-star Mass Ratios in Units of $M_{\odot}/(0.09M_{\odot})$ from Grimm et al. (2018) and Planet-to-star Radius Ratios $R_p/R_*$ from Delrez et al. (2018b) Compared with the Results from This Paper</caption><thead><tr><th>Source</th><th>Quantity</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th><th>g</th><th>h</th></tr></thead><tbody><tr><td>Grimm</td><td><sup>M<sub>p</sub> 0.09M<sub>☉</sub></sup>/<sub>M<sub>☉</sub> M<sub>*</sub></sub></td><td>1.017<sup>+0.154</sup><sub>-0.143</sub></td><td>1.156<sup>+0.142</sup><sub>-0.131</sub></td><td>0.297<sup>+0.039</sup><sub>-0.035</sub></td><td>0.772<sup>+0.079</sup><sub>-0.075</sub></td><td>0.934<sup>+0.080</sup><sub>-0.078</sub></td><td>1.148<sup>+0.098</sup><sub>-0.095</sub></td><td>0.331<sup>+0.056</sup><sub>-0.049</sub></td></tr><tr><td>This paper</td><td><sup>M<sub>p</sub> 0.09M<sub>☉</sub></sup>/<sub>M<sub>☉</sub> M<sub>*</sub></sub></td><td>1.3771 ± 0.0593</td><td>1.3105 ± 0.0453</td><td>0.3885 ± 0.0074</td><td>0.6932 ± 0.0128</td><td>1.0411 ± 0.0155</td><td>1.3238 ± 0.0171</td><td>0.3261 ± 0.0186</td></tr><tr><td>Delrez</td><td>R<sub>p</sub>/R<sub>*</sub></td><td>0.0853 ± 0.0004</td><td>0.0833 ± 0.0004</td><td>0.0597 ± 0.0006</td><td>0.0693 ± 0.0007</td><td>0.0796 ± 0.0006</td><td>0.0874 ± 0.0006</td><td>0.0588 ± 0.0012</td></tr><tr><td>This paper</td><td>R<sub>p</sub>/R<sub>*</sub></td><td>0.0859 ± 0.0004</td><td>0.0844 ± 0.0004</td><td>0.0606 ± 0.0005</td><td>0.0708 ± 0.0006</td><td>0.0804 ± 0.0005</td><td>0.0869 ± 0.0005</td><td>0.0581 ± 0.0009</td></tr></tbody></table>
<table><caption>**Table 4**
Prior Bounds on Photodynamic Parameters</caption><thead><tr><th>Parameter</th><th>Units</th><th>Prior</th></tr></thead><tbody><tr><td>b<sub>0</sub></td><td>R<sub>*</sub></td><td>U(0, 1)</td></tr><tr><td>R<sub>p</sub>/R<sub>*</sub></td><td>...</td><td>U(0, 0.2)</td></tr><tr><td>ρ<sub>*</sub></td><td>ρ<sub>☉</sub></td><td>U(0, 100)</td></tr><tr><td>(q<sub>1,Ch 1</sub>, q<sub>2,Ch 1</sub>)</td><td>...</td><td>U(0, 1)</td></tr><tr><td>(q<sub>1,Ch 2</sub>, q<sub>2,Ch 2</sub>)</td><td>...</td><td>U(0, 1)</td></tr></tbody></table>

Note. Note that the same bounds on the impact parameter, $b_0$, and radius ratio, $R_p/R_*$, are placed on all seven planets.

To help visualize the model, a photodynamical model with the best-fit parameters is shown in Figure 7 computed over 1600 days. Planets b and c have short periods and are far from a $j:j+1$ period ratio. Hence, both of these planets show weak TTVs and straighter, but still slightly meandering, river plots. The outer five planets are pairwise close to a series of $j:j+1$ resonances, showing strong TTVs on the timescale of the TTV period of ≈490 days. The other prominent feature for the outer four planets is the slight zigzag of the transits due to chopping (shown in Figure 3).

Table 3 shows the radius ratios from Delrez et al. (2018b) alongside those from the present analysis. The precision of the measurements did not improve significantly, while the radius ratios shifted by $1\sigma-2\sigma$. Figure 8 shows the posterior probability distribution of impact parameters in units of the stellar radius, $b_0$, derived from the photodynamical model. Figure 9 shows the probability distribution of stellar density. The density correlates with the impact parameters of each planet, reaching a tail of lower values for the higher impact parameters of each planet. The tail of the density probability distribution has an approximately exponential scaling with the density below the peak and cuts off as a normal distribution above. In Table 5, we report the median and 68.3% confidence interval of the stellar density. The inferred density is both slightly larger and more precise than prior analyses (Delrez et al. 2018b), which we discuss below.

Combining the measured density with the measured orbital periods of the planets, we derive the semimajor axis of each planet in units of the stellar radius,

$$ \frac{a}{R_{*}} = \left( \frac{P^{2}GM_{\odot} \rho_{*}}{4\pi^{2}R_{\odot}^{3} \rho_{\odot}} \right)^{1/3} \quad (11) $$

With the measured impact parameters, we compute the inclinations of the planets from Winn (2010),

$$ I = \cos^{-1} \left[ b_0 \left( \frac{a}{R_*} \right)^{-1} \right], \quad (12) $$

where we have ignored the eccentricity in this formula due to the extremely small values of the eccentricities of the planets from the transit-timing analysis (see Table 2). The resulting inclination posterior distribution is displayed in Figure 10. Although the inclination is derived from the impact parameters, which we constrain to be positive, in practice the photodynamical model cannot distinguish between inclinations of $I$ and $180 - I$ (Figure 10), and so we created a histogram of these two options with equal probability.

## 5.2. Mutual Inclinations and Stellar Density

The outer four planets, e through h, have inclinations that are more precisely determined, and, remarkably, their peak probabilities are aligned very closely, to less than $0.1^\circ$, save for the degeneracy of $I$ versus $180 - I$. The inner three planets have poorer constraints on their inclinations due to the larger uncertainty of their impact parameters (as seen in Figure 8). Yet, their inclination posteriors have significant overlap with the outer four planets.

As just mentioned, because each inclination may only be inferred relative to the center of the star, the derived distribution is reflected through $180 - I$. However, if some of the planets orbited above and some below the plane of the disk of the star, it would be very improbable for the outer four planets to show such a precise alignment. We conclude that it may be likely that all of the planets transit the same hemisphere of the star as shown in Luger et al. (2017a): the planets' 3D orbital inclinations are likely precisely aligned. This also implies that their longitudes of ascending node are likely aligned as well, and so in principle we can place a prior on the scatter of the mutual inclinations of the planets. We have rerun a photodynamic Markov chain with an inclination prior such that the planets' inclinations are drawn from a Gaussian about their mean value, with a standard deviation of $\sigma_\theta$, which is allowed to freely vary in the chain. We find a very tightly aligned distribution of inclinations under this assumption, shown in Figure 11. We also find that very small values of $\sigma_\theta$ are preferred, with $\sigma_\theta = 0.041^{+0.031}_{-0.016}$. If the outer and inner planets are in fact derived from a common inclination distribution, this implies that the TRAPPIST-1 planetary orbits are extremely flat, even flatter than the Galilean moons that have a dispersion in inclination of $0.25^\circ$.

The inclination prior also enables a more precise and symmetric estimate of the density of the star, $\rho_*/\rho_{\odot} = 53.22 \pm 0.53$. Why is this? Well, the inclination prior tightens the distribution of the impact parameters of planets b and c (as can be seen by comparing Figures 10 and 11). These inner two planets have deep and frequent transits and the sharpest ingress and egress, and hence they provide the tightest constraint upon the density of the star of all seven planets (Ducrot et al. 2020). Thus, given that the inclination prior tightens the distributions of inclinations of these two planets, the stellar density posterior is correspondingly tighter, and the low stellar density tail of the posterior is eliminated (see Figure 9). Despite this tighter constraint upon the stellar density, we decide to forego its use in computing the densities of the planets given the assumptions inherent in the inclination prior.

The coplanarity of the planets may be used to constrain the presence of a more distant, inclined planet given the scatter in their mutual inclinations induced by gravitational perturbations (Jontof-Hutter et al. 2018). Such an analysis should be carried out, but we leave this to future work.

## 6. Planet Densities and Mass–Radius Relation

With the completion of the transit-timing analysis and photodynamic analysis, we are now ready to revisit the mass–radius relation of the TRAPPIST-1 planets.

The only component missing is a constraint upon the mass of the host star. We use the recent analysis by Mann et al. (2019), who have constructed a sample of nearby M-dwarf binaries to
<table><caption>**Table 5**

Parameters Derived from the Photodynamic Model</caption><thead><tr><th>Parameter:</th><th>ρ<sub>*</sub>/ρ<sub>☉</sub></th><th>q<sub>1,Ch1</sub></th><th>q<sub>2,Ch1</sub></th><th>q<sub>1,Ch2</sub></th><th>q<sub>2,Ch2</sub></th></tr></thead><tbody><tr><td>Value:</td><td>53.17<sup>+0.72</sup><sub>-1.18</sub></td><td>0.133 ± 0.052</td><td>0.26 ± 0.19</td><td>0.059 ± 0.024</td><td>0.49 ± 0.20</td></tr><tr><td>Parameter:</td><td>ρ<sub>*</sub> (g cm<sup>-3</sup>)</td><td>u<sub>1,Ch1</sub></td><td>u<sub>2,Ch1</sub></td><td>u<sub>1,Ch2</sub></td><td>u<sub>2,Ch2</sub></td></tr><tr><td>Value:</td><td>75.05<sup>+1.02</sup><sub>-1.66</sub></td><td>0.161 ± 0.093</td><td>0.20 ± 0.15</td><td>0.218 ± 0.056</td><td>0.021 ± 0.098</td></tr></tbody></table><table><thead><tr><th>Planet:</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th><th>g</th><th>h</th></tr></thead><tbody><tr><td>R<sub>p</sub>/R<sub>*</sub></td><td>0.08590 ± 0.00037</td><td>0.08440 ± 0.00038</td><td>0.06063 ± 0.00052</td><td>0.07079 ± 0.00055</td><td>0.08040 ± 0.00047</td><td>0.08692 ± 0.00053</td><td>0.05809 ± 0.00087</td></tr><tr><td>Depth (%)</td><td>0.7378 ± 0.0064</td><td>0.7123 ± 0.0064</td><td>0.3676 ± 0.0063</td><td>0.5012 ± 0.0078</td><td>0.6465 ± 0.0076</td><td>0.7555 ± 0.0092</td><td>0.3375 ± 0.0101</td></tr><tr><td>T (minutes)</td><td>36.06 ± 0.11</td><td>42.03 ± 0.13</td><td>48.87 ± 0.24</td><td>55.76 ± 0.26</td><td>62.85 ± 0.25</td><td>68.24 ± 0.28</td><td>76.16 ± 0.56</td></tr><tr><td>τ (minutes)</td><td>2.889 ± 0.046</td><td>3.320 ± 0.054</td><td>2.816 ± 0.044</td><td>3.825 ± 0.071</td><td>5.158 ± 0.089</td><td>6.310 ± 0.109</td><td>4.846 ± 0.113</td></tr><tr><td>b/R<sub>*</sub></td><td>0.095<sup>+0.065</sup><sub>-0.061</sub></td><td>0.109<sup>+0.059</sup><sub>-0.061</sub></td><td>0.063<sup>+0.063</sup><sub>-0.043</sub></td><td>0.191<sup>+0.041</sup><sub>-0.041</sub></td><td>0.312<sup>+0.023</sup><sub>-0.018</sub></td><td>0.379<sup>+0.018</sup><sub>-0.014</sub></td><td>0.378<sup>+0.024</sup><sub>-0.023</sub></td></tr><tr><td>a/R<sub>*</sub></td><td>20.843<sup>+0.094</sup><sub>-0.155</sub></td><td>28.549<sup>+0.129</sup><sub>-0.212</sub></td><td>40.216<sup>+0.182</sup><sub>-0.299</sub></td><td>52.855<sup>+0.239</sup><sub>-0.392</sub></td><td>69.543<sup>+0.314</sup><sub>-0.516</sub></td><td>84.591<sup>+0.382</sup><sub>-0.628</sub></td><td>111.817<sup>+0.505</sup><sub>-0.830</sub></td></tr><tr><td>I(°)</td><td>89.728 ± 0.165</td><td>89.778 ± 0.118</td><td>89.896 ± 0.077</td><td>89.793 ± 0.048</td><td>89.740 ± 0.019</td><td>89.742 ± 0.012</td><td>89.805 ± 0.013</td></tr></tbody></table>

**Note.** Top: stellar density (in units of solar density), limb-darkening parameters ($q_1, q_2$) in Spitzer channels 1 and 2, and stellar density in cgs units and limb-darkening parameters $u_1$ and $u_2$. Bottom: planet-to-star radius ratio, $R_p/R_*$; transit depth, $(R_p/R_*)^2$; transit duration, $T$ (from first to fourth contact); ingress/egress duration, $\tau$ (from first to second contact or third to fourth contact); impact parameter in units of stellar radius, $b_0$ (assumed to be positive); ratio of semimajor axis to stellar radius, $a/R_*$; and inclination $I$ in degrees (for $b_0 > 0$).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** River plots showing every transit over 1600 days for one planet per panel (left to right are b-h, as labeled; the transits of companion planets are omitted from each panel). The x-axis ranges over 200/400 × 30 s exposures centered on the mean ephemeris for the nth transit for b-d/e-h, respectively (note the 30 s exposures have a higher resolution than the binned Spitzer time resolution). Each row contains a transit model, with green being out of transit, and blue in transit. There are (1059, 661, 395, 262, 173, 129, 85) transits of planets b-h, respectively. Planets d and h have the smallest sizes, and hence the shallowest depths, causing a lighter color during transit. </></figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Probability of planet impact parameters using the photodynamic model described in the text. </></figcaption></figure>

calibrate the mass-luminosity ($M_* - M_{K_S}$) relation of M dwarfs down to a mass of 0.075 $M_⊙$.30 Given the precise parallax measurement available for TRAPPIST-1 thanks to Gaia (Lindegren et al. 2018), the relation yields an estimated mass of $M_* = 0.0898 \pm 0.0023 M_⊙$.

To derive the masses of the planets, we draw planet-to-star mass ratios from the posterior distribution of the transit-timing analysis (Section 4), which we multiply by the mass of the star drawn from a normal distribution with $M_* = 0.0898 \pm 0.0023 M_⊙$. We then draw the planet-to-star radius ratios and stellar density from the posterior distribution from the photodynamic analysis (Section 5). With the same mass draw,

we compute the stellar radius as

$$R_* = \left( \frac{M_* \rho_\odot}{M_\odot \rho_*} \right)^{-1/3} R_\odot, \qquad (13)$$

which we multiply by each of the radius ratios drawn from the same sample to obtain the planet radii. We carry this out for a large number of samples to derive the probability distribution of the masses and radii of the entire posterior probability sample of the planets.

The probability distribution for the masses and radii of the seven planets are shown in Figure 12. The maximum-likelihood values and the posterior distributions (for $1\sigma$ and $2\sigma$ confidence) are both plotted in this figure. We postpone to Section 9 a detailed analysis of the densities and resulting constraints on the bulk compositions of the planets.

<footnote n="30">30 Note that “M” is being used in three ways here: spectral category (M dwarf), stellar mass ($M_*$), and absolute magnitude in the $K_S$ band, $M_{K_S}$.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Stellar density derived from the photodynamic model relative to the solar density, with no prior (blue solid line) and with the relative inclination prior (orange dashed line). </></figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Posterior distribution of inclination angles of the planets given the photodynamical model. </></figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Posterior distribution of inclination angles of the planets from the photodynamical model assuming a prior on the mutual inclinations of $\Pi_i(2\pi\sigma_0^2)^{-1/2}e^{-(h-l)^2/(2\sigma_0^2)}$. </></figcaption></figure>

In addition to masses and radii, we also derive other planetary properties, given in Table 6. Each of the planets has a density intermediate between Mars ($\rho_{\oplus} = 3.9335 \text{ g cm}^{-3} = 0.713 \rho_{\oplus}$) and Earth ($\rho_{\oplus} = 5.514 \text{ g cm}^{-3}$). The surface gravities span a range from 57% of Earth (planet h) to 110% of Earth (planet b).

## 7. Stellar Parameters

A by-product of our analysis is a revision of the properties of the host star. The empirically based mass estimate for the star based on Mann et al. (2019) is consistent with the mass derived by Van Grootel et al. (2018), who first proposed that the mass of the TRAPPIST-1 star is $\approx 0.09 M_{\oplus}$ based upon stellar evolution models and a ground-based parallax measurement. Ducrot et al. (2020) find a luminosity for the star of $L = (5.53 \pm 0.19) \times 10^{-4} L_{\odot}$, which, when compared with stellar evolution models, yields a mass of $M = 0.09016 \pm 0.0010 M_{\oplus}$, which is also consistent with the Mann et al. (2019) value. Burgasser & Mamajek (2017) found an older age for the host star, $7.6 \pm 2.2$ Gyr, which implies an inflated radius for the star compared with evolutionary models.

Our analysis differs slightly from our prior Spitzer analyses (Delrez et al. 2018b; Ducrot et al. 2020) in that we do not place a prior upon the quadratic limb-darkening coefficients of the TRAPPIST-1 host star. This is motivated by the fact that late-M-dwarf atmospheres are very complex to model and have yet to match observed spectra precisely (Allard et al. 2011, 2012; Juncher et al. 2017), and thus, it is possible that limb-darkening predictions may not be reliable. We investigated using a higher-order quartic limb-darkening law and found that this was disfavored by the Bayesian information criterion (BIC) and that the best-fit model differed negligibly in the model parameters. We also simulated more realistic limb-darkening models based on 3D stellar atmospheres (Claret 2018) and found that a quadratic law was sufficient to recover the correct model parameters with negligible systematic errors.

The TRAPPIST-1 system has the advantage that the planets sample different chords of the stellar disk (Figure 8; also see Delrez et al. 2018b), and given the large number of transiting planets, we are afforded multiple constraints on the stellar limb-darkening parameters. Figure 13 shows our posterior constraints upon the limb-darkening parameters of the star based on our photodynamical model, which are reported in Table 5.

Based on the updated stellar density, we have updated the physical parameters of the star. We adopt the luminosity from Ducrot et al. (2020) and the mass from Mann et al. (2019) given the complete and careful analysis from both of those papers. With our updated constraint on the density of the star, we rederive the other parameters of the star, which are summarized in Table 7. In this table, the stellar effective temperature was computed from the stellar luminosity and radius, with errors computed via Monte Carlo.

## 8. Search for an Eighth Planet

With the detection of multiple transits of the six inner planets in TRAPPIST-1 and a single transit of planet h, a clue to the orbital period of planet h was the series of GLRs found between adjacent triplets of planets (Papaloizou 2014). This relation was then used to predict candidate periods of planet h, based on different integer pairs for its commensurability with planets f and g, and a search through the prior data eliminated all but one possibility at 18.766 days. A subsequent observation of the TRAPPIST-1 system with the K2 spacecraft
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Mass–radius relation for the seven TRAPPIST-1 planets based on our transit-timing and photodynamic analysis. Each planet’s posterior probability is colored by the equilibrium temperature (see color bar), with the intensity proportional to probability, while the $1\sigma$ and $2\sigma$ confidence levels from the Markov chain posterior are plotted with solid lines. Theoretical mass–radius relations are overplotted using the model in Dorn et al. (2016) for an Earth-like molar Fe/Mg = 0.83 ratio with a core (black dashed) and core free (red), and a range of cored models with molar Fe/Mg = 0.75 ± 0.2 (gray). U18 refers to Unterborn et al. (2018; see text). The solid black line was calculated for a 5% water composition, for irradiation low enough (i.e., for planets e, f, g, and h) that water is condensed on the surface (assuming a surface pressure of 1 bar and a surface temperature of 300 K). The umber dashed and solid lines were calculated for a 0.01% and a 5% water composition, respectively, for irradiation high enough (i.e., for planets b, c, and d) that water has fully evaporated in the atmosphere, with the U18 interior model with Fe/Mg = 0.83 and Mg/Si = 1.02 (Turbet et al. 2020). Earth, Venus, and Mars are plotted as single points, also colored by their equilibrium temperatures. <//></figcaption></figure>

<table><caption>**Table 6**
Planetary Parameters from Combining the Transit-timing and Photodynamic Analysis</caption><thead><tr><th>Planet:</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th><th>g</th><th>h</th></tr></thead><tbody><tr><td>R (R<sub>⊕</sub>)</td><td>1.116<sup>+0.014</sup><sub>-0.012</sub></td><td>1.097<sup>+0.014</sup><sub>-0.012</sub></td><td>0.788<sup>+0.011</sup><sub>-0.010</sub></td><td>0.920<sup>+0.013</sup><sub>-0.012</sub></td><td>1.045<sup>+0.013</sup><sub>-0.012</sub></td><td>1.129<sup>+0.015</sup><sub>-0.013</sub></td><td>0.755<sup>+0.014</sup><sub>-0.014</sub></td></tr><tr><td>M (M<sub>⊕</sub>)</td><td>1.374 ± 0.069</td><td>1.308 ± 0.056</td><td>0.388 ± 0.012</td><td>0.692 ± 0.022</td><td>1.039 ± 0.031</td><td>1.321 ± 0.038</td><td>0.326 ± 0.020</td></tr><tr><td>ρ (ρ<sub>⊕</sub>)</td><td>0.987<sup>+0.048</sup><sub>-0.050</sub></td><td>0.991<sup>+0.040</sup><sub>-0.043</sub></td><td>0.792<sup>+0.028</sup><sub>-0.030</sub></td><td>0.889<sup>+0.030</sup><sub>-0.033</sub></td><td>0.911<sup>+0.025</sup><sub>-0.029</sub></td><td>0.917<sup>+0.025</sup><sub>-0.029</sub></td><td>0.755<sup>+0.059</sup><sub>-0.055</sub></td></tr><tr><td>g (g<sub>⊕</sub>)</td><td>1.102 ± 0.052</td><td>1.086 ± 0.043</td><td>0.624 ± 0.019</td><td>0.817 ± 0.024</td><td>0.951 ± 0.024</td><td>1.035 ± 0.026</td><td>0.570 ± 0.038</td></tr><tr><td>v<sub>esc</sub> (v<sub>esc,⊕</sub>)</td><td>1.109 ± 0.026</td><td>1.092 ± 0.022</td><td>0.701 ± 0.010</td><td>0.867 ± 0.012</td><td>0.997 ± 0.012</td><td>1.081 ± 0.013</td><td>0.656 ± 0.020</td></tr><tr><td>S (S<sub>⊕</sub>)</td><td>4.153<sup>+0.161</sup><sub>-0.159</sub></td><td>2.214<sup>+0.086</sup><sub>-0.085</sub></td><td>1.115<sup>+0.043</sup><sub>-0.043</sub></td><td>0.646<sup>+0.025</sup><sub>-0.025</sub></td><td>0.373<sup>+0.015</sup><sub>-0.014</sub></td><td>0.252<sup>+0.010</sup><sub>-0.010</sub></td><td>0.144<sup>+0.006</sup><sub>-0.006</sub></td></tr><tr><td>a (10<sup>-2</sup> au)</td><td>1.154 ± 0.010</td><td>1.580 ± 0.013</td><td>2.227 ± 0.019</td><td>2.925 ± 0.025</td><td>3.849 ± 0.033</td><td>4.683 ± 0.040</td><td>6.189 ± 0.053</td></tr><tr><td>R (10<sup>8</sup> cm)</td><td>7.119<sup>+0.087</sup><sub>-0.077</sub></td><td>6.995<sup>+0.086</sup><sub>-0.077</sub></td><td>5.026<sup>+0.071</sup><sub>-0.066</sub></td><td>5.868<sup>+0.082</sup><sub>-0.075</sub></td><td>6.664<sup>+0.085</sup><sub>-0.077</sub></td><td>7.204<sup>+0.094</sup><sub>-0.085</sub></td><td>4.817<sup>+0.091</sup><sub>-0.088</sub></td></tr><tr><td>M (10<sup>27</sup> g)</td><td>8.211 ± 0.412</td><td>7.814 ± 0.335</td><td>2.316 ± 0.074</td><td>4.132 ± 0.130</td><td>6.205 ± 0.184</td><td>7.890 ± 0.226</td><td>1.945 ± 0.122</td></tr><tr><td>ρ (g cm<sup>-3</sup>)</td><td>5.425<sup>+0.265</sup><sub>-0.272</sub></td><td>5.447<sup>+0.222</sup><sub>-0.235</sub></td><td>4.354<sup>+0.156</sup><sub>-0.163</sub></td><td>4.885<sup>+0.168</sup><sub>-0.182</sub></td><td>5.009<sup>+0.138</sup><sub>-0.158</sub></td><td>5.042<sup>+0.136</sup><sub>-0.158</sub></td><td>4.147<sup>+0.322</sup><sub>-0.302</sub></td></tr><tr><td>g (10 m s<sup>-2</sup>)</td><td>1.080 ± 0.051</td><td>1.065 ± 0.042</td><td>0.611 ± 0.019</td><td>0.801 ± 0.024</td><td>0.932 ± 0.024</td><td>1.015 ± 0.025</td><td>0.558 ± 0.037</td></tr><tr><td>v<sub>esc</sub> (km s<sup>-1</sup>)</td><td>12.400 ± 0.292</td><td>12.205 ± 0.241</td><td>7.839 ± 0.110</td><td>9.694 ± 0.133</td><td>11.145 ± 0.137</td><td>12.087 ± 0.142</td><td>7.335 ± 0.227</td></tr><tr><td>S (10<sup>6</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</td><td>5.652<sup>+0.220</sup><sub>-0.216</sub></td><td>3.013<sup>+0.117</sup><sub>-0.115</sub></td><td>1.518<sup>+0.059</sup><sub>-0.058</sub></td><td>0.879<sup>+0.034</sup><sub>-0.034</sub></td><td>0.508<sup>+0.020</sup><sub>-0.019</sub></td><td>0.343<sup>+0.013</sup><sub>-0.013</sub></td><td>0.196<sup>+0.008</sup><sub>-0.008</sub></td></tr><tr><td>a (10<sup>11</sup> cm)</td><td>1.726 ± 0.015</td><td>2.364 ± 0.020</td><td>3.331 ± 0.028</td><td>4.376 ± 0.037</td><td>5.758 ± 0.049</td><td>7.006 ± 0.060</td><td>9.259 ± 0.079</td></tr></tbody></table>

**Note.** The units are given with respect to Earth first and cgs second.

revealed four more transits of planet h occurring at precisely the period that was predicted (Luger et al. 2017b). The existence of the GLRs among the seven known planets has been used to forecast the possible existence of an eighth planet

interior (Pletser & Basano 2017) and exterior (Kipping 2018) to the seven known transiting planets. There is yet to be a definitive detection of an eighth transiting planet based upon the currently available data (Ducrot et al. 2020).
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Limb-darkening constraints, $1\sigma$ and $2\sigma$ confidence contours. Red is Spitzer IRAC Channel 1 (3.6 µm), while green is Channel 2 (4.5 µm). Error bars indicate the limb-darkening parameters and uncertainties used as priors in Ducrot et al. (2020). </></figcaption></figure>

<table><caption>**Table 7**
Updated Stellar Parameters Based on the Combined Analysis</caption><thead><tr><th>Parameter</th><th>Value</th><th>References</th></tr></thead><tbody><tr><td>M(M<sub>☉</sub>)</td><td>0.0898 ± 0.0023</td><td>Mann et al. (2019)</td></tr><tr><td>R(R<sub>☉</sub>)</td><td>0.1192 ± 0.0013</td><td>This paper</td></tr><tr><td>L(L<sub>☉</sub>)</td><td>0.000553 ± 0.000019</td><td>Ducrot et al. (2020)</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>2566 ± 26</td><td>This paper</td></tr><tr><td>log<sub>10</sub>(g(cm s<sup>-2</sup>))</td><td>5.2396<sup>+0.0056</sup><sub>-0.0073</sub></td><td>This paper</td></tr></tbody></table>

<table><caption>**Table 8**
Predictions for a GLR of Planets g and h with an Eighth Planet, Planet i, with Period $P_i$</caption><thead><tr><th>p</th><th>q</th><th>P<sub>i</sub> (day)</th><th>P<sub>i</sub>/P<sub>h</sub></th><th>j</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>39.029</td><td>2.08</td><td>1</td></tr><tr><td>1</td><td>2</td><td>25.347</td><td>1.35</td><td>3</td></tr><tr><td>1</td><td>3</td><td>22.695</td><td>1.21</td><td>4</td></tr><tr><td>2</td><td>3</td><td>28.701</td><td>1.53</td><td>2</td></tr></tbody></table>

**Note.** The ratio with the period of planet h is given, as well as the value of j for which $P_i/P_h \approx (j+1)/j$.

It may be possible to detect an exterior eighth planet via TTVs induced on the inner seven planets. Planet h should experience the strongest perturbations by an exterior eighth planet due to the fact that TTVs are a very strong function of the proximity of planets to one another and also to resonance. Table 8 shows the predictions for the period of planet “i,” $P_i$, assuming a GLR configuration with planets g and h given by

$$P_i = q / (-pP_g^{-1} + (p+q)P_h^{-1}) \quad (14)$$

for a range of $1 \le p, q \le 3$, which is the same range of integers for the GLRs among the inner seven planets. Interestingly, these cases are all close to a $j:j+1$ period ratio with planet h and thus should strongly perturb planet h due to forcing at this frequency.

We carried out a transit-timing search for an eighth planet by placing planets with mass ratios between $2 \times 10^{-6}$ and $5 \times 10^{-5}$ at these four trial orbital periods in a coplanar

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Limits on an eighth planet, “i,” for a search near the periods in Table 8. The eight-planet models are only plotted if they led to an improvement in log likelihood. Only two of the optimized likelihoods reach the difference in BIC >0 indicated on the plots; however, these two cases have an eighth planet with a relatively large eccentricity and are distant from a GLR with g and h. Orange points have eccentricities smaller than 0.01; light-blue points have larger eccentricities. </></figcaption></figure>

configuration with the other seven planets drawn from a random orbital phase at the initial time and with eccentricity vector elements drawn from a random normal of width 0.005. We placed a Gaussian prior on the eccentricity vector elements of the eighth planet with a standard deviation of 0.14 to avoid unstable configurations. We then optimized the likelihood with the eight-planet model, carrying out 11,200 optimizations on 112 CPUs with 100 optimizations per CPU, lasting seven days each for about 20,000 CPU hours.

We then carried out a search for evidence of perturbations by planet i by determining if the optimized likelihood of the transiting planets was improved by adding an eighth planet to the transit-timing model, using the BIC to penalize the additional degrees of freedom of the eight-planet model (Wit et al. 2012). We searched for a change to BIC for the eight-planet model over the seven-planet model with a difference of better than $5 \log N_{\text{trans}} = 30.5$. Given that the inner seven planets show orbital eccentricities with values $\lesssim 0.01$, we only considered an eighth planet candidate plausible if it shows an eccentricity less than this cutoff.

In all 11,200 trial optimization cases, we found that only two of the eight-planet models did exceed the BIC criterion, but both significantly exceed an eccentricity of 0.01. Figure 14 shows the change in BIC versus orbital period and mass for planet “i,” assuming a mass of the star of $M_* = 0.09 M_☉$. These two cases with $\Delta \text{BIC} > 0$ do not appear to be plausible planet candidates: they only just exceed the BIC criterion, they both have large eccentricities, and they are not in close proximity to a GLR with planets g and h (even though the initial parameters of the optimization were started near a GLR).

We also carried out a search for an eighth planet interior to planet b and found even smaller improvements in the log likelihood than in the exterior case.

We have not carried out an exhaustive search for eight-planet models at other orbital periods due to the significant volume of parameter space to search. However, it is still possible that an exterior eighth planet is perturbing planet h and may modify its transit times to a point that affects the posterior masses we infer from our seven-planet model. In principle, one could include the effect of an eighth planet on the mass inference by adding it
to the Markov chain modeling; in practice, this would be a challenging model to sample due to the multimodal nature of the parameter space. We defer such analysis to future work.

## 9. Interior Compositions

In this section, we present a theoretical interpretation of the planets' interior properties based upon the mass-radius relation we inferred in Section 6. As there is significant degeneracy in the possible interior compositions, we present a menu of different possibilities in Section 9.2. However, we start with an approach that is less dependent upon the assumption of interior composition, which we term the “normalized density.”

### 9.1. Initial Analysis of Planet Densities across the System

The probability distribution for the masses and radii of the seven planets are shown in Figure 12 alongside several theoretical mass-radius relationships added for comparison. We have added three rocky mass-radius relationships with different bulk Fe/Mg compositions: (1) molar Fe/Mg = 0.75 ± 0.2 as suggested by Unterborn et al. (2018) to represent the rocky interior of all TRAPPIST-1 planets with a $1\sigma$ range of Fe/Mg ratios consistent with local stellar abundances, (2) the Sun-like value of molar Fe/Mg = 0.83 (Lodders et al. 2009), and (3) a core-free model with Earth-like refractory ratios, but in which all of the iron is oxidized in the mantle (Elkins-Tanton & Seager 2008). Rocky interiors are calculated similar to the models of Dorn et al. (2016) with two adaptations: we are using the equation of state of Hakim et al. (2018) for pure iron and Sotin et al. (2007) for silicates. We have also added the theoretical mass-radius relationships for planets endowed with a water layer, both for planets that are irradiated less (black line; water) and more (umber lines; steam) than the runaway greenhouse irradiation threshold (Turbet et al. 2020).

The comparison of measured masses and radii with theoretical mass-radius relationships reveals several striking results. First, all seven TRAPPIST-1 planets appear to be consistent with a line of interior isocomposition at the $1\sigma$ level. There are multiple theoretical mass-radius curves that overlap with all seven planets' mass-radius probability distributions (Figure 12), which may be a good indication that the composition varies little from planet to planet. Second, all of the TRAPPIST-1 planets have lower uncompressed densities than solar system terrestrial planets. This likely means that the TRAPPIST-1 planets either have a lighter interior (e.g., lower iron content) or are enriched with volatiles (e.g., water).

We next searched for variations of density across the planets. For this, we took each planetary density calculated from $10^4$ samples and divided by the density of the closest pair of mass and radius of a fully differentiated 20 wt% iron, 80 wt% silicate ($MgSiO_3$) interior planet with no surface layers, which is less iron rich than Earth. A planet with a normalized density of 1 has exactly the same density as the reference model, while a normalized density >1 (<1) is denser (lighter), than the reference model, respectively. Figure 15 shows the resulting histograms of the posterior probability of the normalized TRAPPIST-1 planet densities. We then plot in Figure 16 the normalized densities (along with their $1\sigma$ uncertainty) as a function of the orbital periods of the planets. The normalized planet density appears very uniform across the seven planets, with perhaps a slight decrease with the increase of the orbital period (or the distance to the host star). We fit a line to the normalized density, y, versus orbital

<figure><img src="image_18.png" /><figcaption>Figure 15. Probability density function of the normalized density of all seven planets in the system. </></figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 16. Normalized planet densities (with $1\sigma$ error bars) vs. planet orbital periods. The light-blue band is the 68% confidence interval of the weighted mean normalized density of all seven planets. The orange lines show the 68.3% confidence intervals of linear fits to the normalized densities computed from $10^4$ draws from the posterior. The mean fit to the normalized density vs. period is $y = aP + b$, where $a = 1.042 \pm 0.034$ and $b = -0.0043 \pm 0.0036$. </></figcaption></figure>

period, *P*, for $10^4$ posterior samples and found a relation of $y = (1.042 \pm 0.034) - (0.0043 \pm 0.0036)P$, where the coefficients are the 68.3% confidence interval. There is only weak evidence for a declining trend of normalized density with orbital period: 88% of the fits to the $10^4$ posterior samples have slopes with a negative value, while 12% of the slopes fit have a positive value. If in the future more precise data strengthen this trend, then this may indicate that either (i) the outer planets are depleted in heavy elements (e.g., iron) compared to the inner ones, or (ii) the outer planets are enriched in volatiles (e.g., water) compared to the inner ones. However, based on the current data, we suggest that the planets' compositions could be rather uniform in nature.

The interpretation of these observations in terms of internal compositions is discussed in more detail next.
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Iron core mass fraction vs. the planetary orbital periods for a fully differentiated model with molar Mg/Si = 1.02 and no surface layer. The approximate values for Earth (McDonough 2014), Mars (Khan et al. 2018), the Moon (Barr 2016), and common chondrites (Palme et al. 2014) are indicated, as well as the 1σ confidence intervals of the TRAPPIST-1 planets. The light-blue box is the 68.3% confidence region of the weighted mean of all seven planets. The orange lines show the median and 68.3% confidence interval for linear fits to the 10⁴ posterior values for all seven planets. </></figcaption></figure>

## 9.2. Range of Possible Interior Compositions and Volatile Contents

In this subsection, we discuss a range of possible compositions of the planets based on their measured densities, starting with a volatile-poor model in which the densities are fit by varying the CMF (Section 9.2.1) and followed by an analysis in which the solid planets are taken to have an Earth-like composition, to which is added a water fraction needed to create the observed densities (Section 9.2.2). Alternatively, the planets might be explained with an enhanced oxygen content by which all of the iron is oxidized, making the planets core free (Section 9.2.3).

### 9.2.1. Core Mass Fraction

If we assume that the planets' atmospheres contribute a negligible amount to their total radius and that the planets are fully differentiated, composed of rocky mantles ($MgSiO_3$) and iron cores only, then the densities may be used to constrain the portion of the planets' mass that is contained within their cores.

We evaluated the CMFs of the TRAPPIST-1 planets as follows. For each mass/radius pair in our posterior distribution, we have estimated the CMF by linearly interpolating between precalculated mass-radius relationships with our employed interior model. We arbitrarily set each mass/radius pair lighter than a pure silicate ($MgSiO_3$) planet to a CMF of 0. Alternatively, we repeated the same procedure but discarded all CMF values lower or equal to 0. However, we found that the estimate of the CMF is only marginally changed (and only for planets g and h).

Our CMF estimates are provided in Figure 17 and Table 9. Estimates range from $16.1_{-4.2}^{+3.5}$ wt% for planet g up to $26.6_{-5.1}^{+4.6}$ wt% for planet c, which, despite the different central values, have considerable overlapping probability distributions. Figure 17 shows that within the uncertainties, the CMF/iron fractions of the planets are very consistent with one another, with the mean of all planets of $21 \pm 4$ wt% (taking into account the correlations between the planets' CMFs).

There may be a slight trend of the inferred CMF, which decreases with increasing orbital period. The trend is qualitatively similar to that reported on the normalized density (see Figure 16), with similarly weak support: only 88% of the linear fits to the 10⁴ posterior CMF values have a slope with orbital period that is negative, while 12% are positive.

### 9.2.2. Surface Water Content

The observed (weak) variation in the planet densities among all seven planets may instead be due to their differing volatile (e.g., water) inventories.

If we assume a rocky Earth-like interior (CMF = 32.5%, fully differentiated) and only allow an additional condensed<fnref n="31" /> water layer to contribute to the total radius, we can estimate the water mass fractions of the seven planets (b: $2.8_{-1.9}^{+2.1}$ wt%, c: $2.3_{-1.7}^{+1.8}$ wt%, d: $4.4_{-1.5}^{+2.0}$ wt%, e: $2.9_{-1.5}^{+1.7}$ wt%, f: $4.5_{-1.2}^{+1.8}$ wt%, g: $6.4_{-1.6}^{+2.0}$ wt%, h: $5.5_{-3.1}^{+4.5}$ wt%). The lower densities of planets d, f, g, and h can allow for two to three times as much water than for planets b, c, and e. For this simple estimate, we assumed a water layer with a surface temperature of 300 K at 1 bar.

Actual surface conditions and assumed iron content can, however, lead to much larger differences in the estimated water budgets between the inner three and outer four planets. This stems from the fact that the inner three planets are more irradiated than the runaway greenhouse irradiation limit (Kopparapu et al. 2013; Wolf 2017; Turbet et al. 2018) for which all water is vaporized, forming a thick $H_2O$-dominated steam atmosphere. Taking into account the expectation that water should be vaporized for the three inner TRAPPIST-1 planets (Turbet et al. 2019, 2020), their water mass fractions drop drastically to less than 0.01 wt%, i.e., more than several times lower than the water ocean mass fraction of the Earth.

Figure 18 shows the expected water mass fractions for each of the TRAPPIST-1 planets and for four distinct interior compositions (18, 25, 32.5, and 50 wt% iron content). It shows that the same qualitative trend of water versus orbital period is relatively robust across a large range of assumptions on the interior composition thanks to the transition from runaway greenhouse for planets b–d to surface liquid water for planets e–h.

Higher estimated water budgets for the outer three or four planets could be a clue that they formed beyond the water condensation line at $\approx 0.025$ au (Unterborn et al. 2018). This could also be due to the significant differences in water loss (through atmospheric escape) arising from variations of irradiation and gravity among the TRAPPIST-1 planets (Lissauer 2007; Bolmont et al. 2017; Bourrier et al. 2017). However, again, we caution that trends in the planetary volatile content are only weakly supported by the current data.

### 9.2.3. Core-free Planets

Given that the data may be consistent with an isocomposition mass-radius relation, we next consider another intriguing possibility: that the interiors of the planets are fully oxidized. If, instead of forming a core, all of the iron is oxidized and remains in the mantle, the size of a planet may increase by a few percent (Elkins-Tanton & Seager 2008). This turns out to be about the amount of radius inflation necessary to match the

<footnote n="31">³¹ Note that it is likely unwarranted to assume condensed surface water for the inner three planets given their location within the runaway greenhouse zone (Turbet et al. 2020).</footnote>
<table><caption>**Table 9**
Core Mass Fractions, Molar Fe/Mg Ratio (for a Fully Differentiated Model), and Water Mass Fractions Inferred for Each TRAPPIST-1 Planet, as well as the Weighted Means</caption><thead><tr><th>Planet:</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th><th>g</th><th>h</th><th>Avg b-h</th></tr></thead><tbody><tr><td>CMF (wt%)</td><td>25.2<sup>+5.3</sup><sub>-6.0</sub></td><td>26.6<sup>+4.6</sup><sub>-5.1</sub></td><td>19.7<sup>+4.7</sup><sub>-5.1</sub></td><td>24.6<sup>+4.3</sup><sub>-4.9</sub></td><td>20.1<sup>+3.5</sup><sub>-4.2</sub></td><td>16.1<sup>+3.5</sup><sub>-4.2</sub></td><td>16.5<sup>+9.3</sup><sub>-10.0</sub></td><td>20.9 ± 3.6</td></tr><tr><td>Fe/Mg molar ratio</td><td>0.60<sup>+0.18</sup><sub>-0.18</sub></td><td>0.64<sup>+0.16</sup><sub>-0.16</sub></td><td>0.44<sup>+0.14</sup><sub>-0.13</sub></td><td>0.58<sup>+0.14</sup><sub>-0.14</sub></td><td>0.45<sup>+0.10</sup><sub>-0.11</sub></td><td>0.34<sup>+0.09</sup><sub>-0.10</sub></td><td>0.35<sup>+0.27</sup><sub>-0.23</sub></td><td>0.47 ± 0.07</td></tr><tr><td>H<sub>2</sub>O (wt%) for:</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>CMF = 18%</td><td>&lt;10<sup>-3</sup></td><td>&lt;10<sup>-3</sup></td><td>&lt;10<sup>-3</sup></td><td>0.0<sup>+0.0</sup><sub>-0.0</sub></td><td>0.0<sup>+0.0</sup><sub>-0.0</sub></td><td>0.72<sup>+1.3</sup><sub>-0.72</sub></td><td>0.6<sup>+3.4</sup><sub>-0.6</sub></td><td></td></tr><tr><td>CMF = 25%</td><td>&lt;10<sup>-3</sup></td><td>&lt;10<sup>-3</sup></td><td>&lt;10<sup>-3</sup></td><td>0.3<sup>+1.8</sup><sub>-0.3</sub></td><td>1.9<sup>+1.5</sup><sub>-1.3</sub></td><td>3.5<sup>+1.6</sup><sub>-1.3</sub></td><td>3.0<sup>+3.8</sup><sub>-3.0</sub></td><td></td></tr><tr><td>CMF = 32.5%</td><td>&lt;10<sup>-3</sup></td><td>&lt;10<sup>-3</sup></td><td>&lt;10<sup>-3</sup></td><td>2.9<sup>+1.7</sup><sub>-1.5</sub></td><td>4.5<sup>+1.8</sup><sub>-1.2</sub></td><td>6.4<sup>+2.0</sup><sub>-1.6</sub></td><td>5.5<sup>+4.5</sup><sub>-3.1</sub></td><td></td></tr><tr><td>CMF = 50%</td><td>0.05<sup>+0.08</sup><sub>-0.03</sub></td><td>0.03<sup>+0.05</sup><sub>-0.02</sub></td><td>0.002<sup>+0.002</sup><sub>-0.0009</sub></td><td>9.4<sup>+2.2</sup><sub>-1.8</sub></td><td>12<sup>+2.0</sup><sub>-1.7</sub></td><td>14<sup>+2.0</sup><sub>-1.7</sub></td><td>12<sup>+4.4</sup><sub>-3.9</sub></td><td></td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Theoretical water content estimates (along with $1\sigma$ error bar) vs. planetary orbital periods. Colors depict different compositions for the rocky interior (18, 25, 32.5, and 50 wt% CMF). For high CMF, estimated water contents are larger in order to fit the total mass and radius. </></figcaption></figure>

TRAPPIST-1 planets when compared with our solar system planets.

If we assume that the refractory ratios match a solar composition and that all seven planets lack an atmosphere, then it turns out that all seven planets are consistent with a core-free, oxidized composition (Figure 12; red line). For this model, the bulk mass abundance ratios for Fe/Si/Mg/O are 29.2/17.3/15.3/38.2 wt% with a magnesium number of 0.55 (Mg/(Mg +Fe)) mol fraction; this model has a significant increase in oxygen compared to the bulk Earth with 29.7 wt% (McDonough 2014). Such a scenario would likely require the formation of the planets at large distances from the star in a highly oxidizing environment (Elkins-Tanton & Seager 2008) and a lower devolatization temperature intermediate between that of Earth and chondrites (Wang et al. 2019). Hence, although this hypothesis efficiently explains the TRAPPIST-1 data, it remains to be seen whether a geochemical model that results in the high oxidation of iron throughout the processes of planet formation and evolution (Kite et al. 2020) can be constructed.

## 10. Discussion

Here we discuss some of the implications of the results in the foregoing sections.

### 10.1. Timing Uncertainties

As reported in Section 3, the transit-timing measurements we have made show an excess of outliers with respect to the

measurement uncertainties of each transit. We were unable to identify a culprit (or culprits) for these discrepancies but wish to speculate on what may be the origin of these outliers. The cumulative distribution of these outliers (Figure 1) indicates that about 10% of transits are affected at some level. It is also interesting to note that the core of the distribution has a slightly smaller width of about 87% of the measurement errors, indicating that for about 90% of the transits, the uncertainties may be overestimated. This may be a consequence of inflating the uncertainties to account for correlated noise rather than modeling the data with, for example, a Gaussian process; further reanalysis of the data will be needed to check this hypothesis.

Could the timing outliers be due to stellar flares? In Vida et al. (2017) and Ducrot et al. (2020), the frequency distribution of stellar flares is shown to be rising toward smaller flare energies. This could mean that the more frequent, but lower energy, flares occur at a level that is swamped by the photon noise and thus not visible to an observer. We used the spectrum and energy calibration of Spitzer flares measured by Ducrot et al. (2020) to extrapolate the frequency of lower energy flares (which are not detected in Spitzer due to photon noise). As an example, for planet h, the transit time can be affected by a flare that occurs at ingress or egress (duration $2\tau \approx 10$ minutes). We estimate that a flare of energy $10^{31}$ erg could cause a $1.5\sigma$ timing outlier if it occurs during ingress or egress. This has a probability of only $\approx 0.3\%$ of occurring during the 10 minutes of ingress or egress, and thus cannot be responsible for 10% of outliers for planet h. We carried out a similar estimate for the other planets, and we conclude that low-level flaring activity cannot be the cause of the timing outliers.

Other possible causes of the timing outliers are correlated stellar variability, starspot crossings, or instrumental systematics. We do not yet have an estimate of the magnitudes of these effects and so cannot reach a conclusion about where the origin of the timing outliers lies.

### 10.2. Possible Systematic Errors

In this section, we consider possible factors that might affect our inference of the densities of the planets. Simulated planetary densities predict CMFs that are similar to Earth, with a very small scatter (Scora et al. 2020). Hence, the fact that the TRAPPIST-1 planets have inferred planetary densities that are less than this could be due to systematic uncertainties that are not captured by our modeling.

The transit depths determine the planet-to-star radius ratios, but these measurements are affected by the nonuniform surface brightness of the star. Fortunately, the multiple impact
parameters of the planets yield a constraint on the infrared limb darkening, which is fairly weak compared with optical bands. However, starspots can also affect the inferred transit depths (Czesla et al. 2009; Kipping 2012; Oshagh et al. 2013, 2014; McCullough et al. 2014; Rackham et al. 2018). If spots are present on an active latitude that is not on the same hemisphere as the planetary transit chords, this can cause all of the planet radii to be misinferred by a similar factor.

TRAPPIST-1 may have complex surface inhomogeneities, including regions brighter or darker than the mean photosphere (Morris et al. 2018a; Zhang et al. 2018; Wakeford et al. 2019). It is possible that bright or dark regions could bias the apparent transit depths toward larger or smaller measurements, depending on which type of inhomogeneity dominates. Time-variable contamination should average out with many observations, while time-steady inhomogeneity will not, such as active latitudes, polar spots, or even hemispheric asymmetry (Yadav et al. 2015; Brown et al. 2020). We modeled the transit transmission in the K2, SPECULOOS, LT, near-infrared, and Spitzer bands from Ducrot et al. (2020) for all seven planets using the contamination formula from Rackham et al. (2018) with a time-steady, three-temperature model with the temperatures of the three components ranging from 2000 to 2980 K and the covering fraction varying from 0 to 1. The mean effective temperature is constrained by our stellar model parameters (Table 7). We assumed that all seven planets transit the region with the larger covering fraction and that their transit depths are achromatic. We ran a Markov chain fit to the transmission spectra, interpolating the fluxes in the bands between the effective temperature grid points that were spaced by 20 K; we find that the posterior parameters with maximum likelihood are temperatures of (2980, 2331, 2071) K with covering fractions of (0.8, 82.1, 17.1)% . We then computed the expected impact on the transit depths in the two IRAC channels. The constraints are tight: we find that the observed radii should only change by a factor of 1.0072 ± 0.0097 in Channel 1 and 1.0071 ± 0.0108 in Channel 2 (these are the ratios of the observed radii to the actual radii). These factors are consistent with unity at better than 1σ and have uncertainties that are comparable to or smaller than the uncertainties on the absolute planetary radii. We conclude that this form of self-contamination does not greatly influence our results but should lead to caution in the interpretation. This constraint is much stronger than the analysis of Morris et al. (2018b).

Our mass precisions are predicated on a complete model of
the dynamics of the system. We ignore tides and general
relativity, which are too small in amplitude to affect our results
at the current survey duration and timing precision (Bolmont
et al. 2020). Should an eighth planet be lurking at longer orbital
periods, which has yet to reveal itself via significant TTVs or
transits, this may modify our timing solution and shift the
masses slightly. In our timing search for an additional planet,
however, we found that such a planet might only cause shifts at
the ≈1σ level. This possibility begs for caution in interpreting
the potential variation of the iron fraction with orbital period:
should an eighth planet be present beyond planet h, its timing
impact would likely affect the masses of the exterior planets
more significantly than the interior planets. Drawing stronger
conclusions about the variation of the planet iron/CMFs will
likely require longer-term monitoring, especially of planet h,
and/or higher precision timing measurements such as are

expected with JWST, to place tighter constraints on an eighth
planet.

10.3. Planet Masses and Radii in Context

In our current analysis of the transit-timing data for TRAPPIST-1, we have found larger mass ratios for all planets save planet e compared with our most recent analysis in Grimm et al. (2018). Even though most of the planets have shifted by 1σ or more, this does not indicate that the prior analysis was in error. In fact, the masses of all of the planets are strongly correlated, and thus when one planet shifts in the transit-timing solution, they all shift. With the more extensive data set analyzed here, we provide a better constraint over the transit-timing timescale and can also better account for outliers thanks to some redundancy in our measurements. Given the high precision of the Spitzer timing measurements, we expect that our current analysis may remain the most reliable constraint on the masses of the planets until the transit times can be measured with JWST.

In Figure 19, we compare our measurements for the seven TRAPPIST-1 planets with our solar system planets and with exoplanets with radii <1.7R⊕ and masses measured to >5σ retrieved from the NexSci database on 2020 February 26 (Akeson et al. 2013; Christiansen 2018), as well as planet parameters reported in Dai et al. (2019) and Kepler-93b from Dressing et al. (2015).<fnref n="32" /> The uncertainties on the other planets’ masses are the best available to date from radial-velocity (RV) measurements and yet they are much larger than the uncertainties for the TRAPPIST-1 planets, whether considered in a relative or absolute sense. The larger uncertainties of the RV planets make the CMFs difficult to constrain for these more massive planets—core-free and cored models are consistent with most of these planets’ parameters at the 1σ level (Figure 19). Nevertheless, it is notable that the rocky planets for which we currently have data seem to be similar in composition to Earth (Dressing et al. 2015); however, the actual range of bulk rock compositions of rocky exoplanets relative to their host stars is currently debated. This also appears consistent with the observation that the evaporation valley requires rocky planets and their gaseous brethren to have a composition that is a mix of silicates and iron (Owen & Wu 2017).

10.4. Comparison with Radial Velocities

Given the measurements of the masses we have made with
transit timing, this brings up the question: what RV
uncertainties would be required to make mass measurements
of similar precision?

The precision of the mass measurements may be placed in
context by comparing with current RV capabilities. The
predicted semiamplitudes for the seven planets are given in
Table 10. The predicted RV variation of the star induced by the
TRAPPIST-1 planets is plotted in Figure 20, also based upon
our mass measurements from transit timing. The sums of the
semiamplitudes of the planets equal ≈12.7 m s⁻¹, which is
close to the peak amplitude when the planets are all orbiting on
the same side of the star (near 218 days in the plotted figure).
How does this compare with current RV measurements?

Recently Hirano et al. (2020) were able to make high-
precision measurements of the RV of the TRAPPIST-1 host
star, achieving a constraint on the linear variation of the star to

<footnote n="32">³² Note: we corrected Kepler-105b with the Gaia DR2 revised radius of the host star (Berger et al. 2018; Fulton & Petigura 2018).</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Radius vs. mass for solar system terrestrial planets (green dots), TRAPPIST-1 (orange error bars), and other potentially rocky exoplanets from the NExSci database, Dressing et al. (2015), and Dai et al. (2019; red error bars). Planets with a smaller mass uncertainty are shown in a darker red color. Also plotted is a mass-radius relation with a core-mass fraction compatible with Earth (blue) and a core-free model in which the refractory elements retain the solar abundance ratios (purple). </></figcaption></figure>

<table><caption>**Table 10**
RV Semiamplitudes, $K_p$, for the TRAPPIST-1 Planets Predicted from Our Measured Masses</caption><thead><tr><th>Planet</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th><th>g</th><th>h</th></tr></thead><tbody><tr><td>$K_p$ (cm s<sup>-1</sup>)</td><td>382.0</td><td>310.7</td><td>77.6</td><td>120.7</td><td>158.1</td><td>182.2</td><td>39.1</td></tr><tr><td>RV equivalent precision for TRAPPIST-1 host (cm s<sup>-1</sup>)</td><td>19</td><td>13</td><td>2.5</td><td>3.8</td><td>4.7</td><td>5.2</td><td>2.4</td></tr><tr><td>RV equivalent precision for 1 $M_\odot$ host at 1 au (cm s<sup>-1</sup>)</td><td>0.62</td><td>0.50</td><td>0.11</td><td>0.20</td><td>0.28</td><td>0.34</td><td>0.18</td></tr></tbody></table>

**Note.** Equivalent RV precision required to measure the masses to the same precision as measured with TTVs around TRAPPIST-1. Also, equivalent RV precision required if each planet were placed around a solar twin at one astronomical unit.

a precision of 2.5 m s⁻¹, which they ascribe to stellar variability. To compare this with our transit-timing results, the semiamplitude precision that would be needed to achieve the same mass error bars that we have achieved with transit timing ranges from 2.4 to 19 cm s⁻¹, up to 100 times more precise than the RV measurements. Future observations may be able to achieve higher precision RV measurements of TRAPPIST-1 but will continue to contend with stellar variability (Klein & Donati 2019).

Were these planets orbiting a Sun-like star, the semiamplitude RV error would need to be even smaller to achieve the same mass precision we have achieved with transit timing. Table 10 lists what semiamplitude precisions would be required if each one of these planets was placed around a solar twin at one astronomical unit.

The required precision ranges from 1 to 6 mm s⁻¹. This is nearly two orders of magnitude more precise than the highest precision RV measurements for short-period exoplanets reported to date, such as Tau Ceti g, which has a reported RV semiamplitude precision of 11 cm s⁻¹ (Feng et al. 2017). We conclude that the mass precisions of Earth-sized, Earth-insolation planets based on RV must be improved by two orders of magnitude to match our TTV precision for the TRAPPIST-1 system.

## 10.5. Planetary Dynamics

In this section, we discuss some of the dynamical aspects of the planetary system: the eccentricities, the longitudes of periastron, and the GLR angles.
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Predicted radial-velocity variation of the TRAPPIST-1 host star induced by its seven known transiting planets, as well as the current measurement error bar reported by Hirano et al. (2020), which they interpret as an upper limit, thanks to stellar variability. Also plotted are the equivalent semiamplitudes for the seven planets that would be required to achieve the same mass precision as measured with TTVs. <>/</figcaption></figure>

### 10.5.1. Eccentricities

The posterior distribution of the initial eccentricities of the planets is shown in Figure 21. In prior analyses of the TTVs of the TRAPPIST-1 system, we found that the inner two planets, b and c, had significant eccentricities (Grimm et al. 2018). In contrast, with the current analysis, we find that the eccentricity probability distributions of these two planets are significant near zero eccentricity. This is consistent with *N*-body models that include tidal damping of the orbits, which predict that planets b and c should have low eccentricities, $\lesssim 10^{-3}$ (Luger et al. 2017b; Turbet et al. 2018). The other planets are all consistent with the predictions of the tidal evolution model (Luger et al. 2017b).

Figure 22 shows the posterior probability distribution for the eccentricity vectors of each planet. The only two planets consistent with zero eccentricity at $1\sigma$ confidence are planets b and c (blue and orange contours). The other five planets have nonzero eccentricities.

Now, the eccentricity vectors plotted in Figure 22 show the values at the initial time. However, over time, the eccentricity vector of each planet can be decomposed into two components: the mean eccentricity vector (over some timescale) and the variable component (which is time variable, with multiple oscillation timescales driven by the mutual planetary perturbations). Figure 23 shows the eccentricity over a single oscillation for all seven planets. The outer five planets are close to first-order resonances with adjacent planets, and the super-period for each of these planets is close to $P_{TTV} \approx 490$ days thanks to the

near-GLR commensurability for all triplets of planets. This leads to a nearly circular oscillation over this timescale due to circulation of the first-order resonances driving oscillations in the eccentricity vectors of each of these planets. The inner two planets are close to second- and third-order resonances with adjacent planets (b and c are close to 8:5, which is third order, while c and d are close to 5:3, which is second order). Because the strength of these interactions scales as a higher power of eccentricity, these planets show much smaller variation in the time-variable components of their eccentricity vectors. Because planets b and c are close to a third-order resonance, their eccentricity vectors show a three-fold symmetry. On longer timescales, these patterns precess, filling a circular pattern over time. The time-variable eccentricity vector patterns are very similar over the range of posterior values, indicating that it is primarily this component that is constrained by the TTVs of the planets.

The total eccentricity vectors show a wider range of behavior, thanks to a wider variation of the mean eccentricity, as shown in Figure 24. It is clear from this figure that each planet executes an eccentricity vector oscillation about a mean value (which was subtracted off for Figure 23). Unfortunately, the mean eccentricity is less constrained by the TTVs (Linial et al. 2018), and so there is a much wider range of eccentricity vectors that is allowed that manifests as strong correlations among the eccentricity vectors of pairs of planets (Figure 29).
<figure><img src="image_24.png" /><figcaption>Figure 21. Probability distribution of the eccentricities of the planets at the initial time based upon the transit-timing model. </></figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 23. Variable component of the osculating eccentricity vectors plotted from a simulation over 12 days for planets b and c, and over ≈490 days for planets d–h, with initial parameters drawn from the posterior distribution. </></figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 22. Posterior probability distribution for the eccentricity vectors at the initial time for each of the planets. Contours are the 1σ and 2σ confidence limits. The maximum-likelihood parameters are shown as solid points. </></figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 24. Osculating eccentricity vectors computed from a simulation for all seven planets shown for three different draws from the posterior: the first with eccentricities nearest the median of the posterior distribution, the second with eccentricities farthest from the median, and a third drawn randomly from the posterior. As with Figure 23, planets b and c are plotted over 12 days, while planets d–h are plotted over ≈490 days. </></figcaption></figure>

### 10.5.2. Laplace Angles

A remarkable property of the TRAPPIST-1 system is the near-commensurability of adjacent triplets of planets (Luger et al. 2017b), akin to Laplace resonances, with GLR angles given by

$$ \phi_{i,i+1,i+2} = p\lambda_i - (p+q)\lambda_{i+1} + q\lambda_{i+2}, \quad (15) $$

where $\lambda_i$ is the mean longitude of the ith planet, and $p$ and $q$ are small integers. In the case of an isolated triplet of planets, a stable configuration takes on $\phi = 180^{\circ}$, but when planets are captured into a series of GLR commensurabilities, their mutual torques displace the stable configuration (Delisle 2017).

Long-term dynamical simulations show that these GLR angles can take on stable values for extended durations and sometimes can quickly jump in value, flipping symmetrically about 180° (Mah 2018; Brasser et al. 2019), resulting in two possible angles for each triplet of stars, $\phi$ and $360^{\circ} - \phi$. Based on the prior measured planet-to-star mass ratios, Mah (2018) predicted the value of the three-body resonance angles resulting from the values at the end of the simulation.
<figure><img src="image_28.png" /><figcaption>**Figure 25.** GLR angles plotted over 100 yr for three draws from the posterior distribution: one with low eccentricities, one with high, and one randomly chosen. These are compared with the predictions from Mah (2018), shown by dashed horizontal lines, with the values for $\phi_{bcd}$ and $\phi_{cde}$ flipped about 180° (i.e., changed from $\phi$ to 360° - $\phi$). </></figcaption></figure>

In Figure 25, we show the GLR angles for the following triples:

$$ 
\begin{aligned}
\phi_{bcd} &= 2\lambda_b - 5\lambda_c + 3\lambda_d, \\
\phi_{cde} &= \lambda_c - 3\lambda_d + 2\lambda_e, \\
\phi_{def} &= 2\lambda_d - 5\lambda_e + 3\lambda_f, \\
\phi_{efg} &= \lambda_e - 3\lambda_f + 2\lambda_g, \\
\phi_{fgh} &= \lambda_f - 2\lambda_g + \lambda_h.
\end{aligned} 
\quad (16) 
$$

Differences between the predicted and observed angles agree within 0°5–10°, where the predicted values for $\phi$ are taken from Mah (2018), but allowing $\phi_{bcd}$ and $\phi_{cde}$ to be flipped about 180°. It is possible with the updated mass ratios from our analysis that the predictions will be more accurate, which awaits further simulation.

### 10.5.3. Long-term Stability

Prior studies of the TRAPPIST-1 system by Tamayo et al. (2017) found long-lived configurations for systems that had formed via migration. Quarles et al. (2017) examined the stability of the TRAPPIST-1 system, refining the large uncertainties from prior measurement (Gillon et al. 2017) to further constrain the masses of the system. Given the much tighter constraints we have placed on the masses of the planets and the orbital eccentricities, here we reexamine the long-term stability of our posterior distribution.

We have used the GPU N-body integrator GENGA (Grimm & Stadel 2014) to carry out long-term simulations of a set of $10^4$ posterior samples from the timing analysis. These simulations were carried out for $10^7$ yr, which corresponds to 2.4 billion orbital periods of planet b and 195 million orbital periods of planet h. We used a time step of 0.06 days, which gives a total number of $6.1 \times 10^{10}$ integration steps. We find that 100% of these posterior samples are stable over this entire timescale. To check the stability of the samples, we analyzed the evolution of the semimajor axis, $a$, and eccentricity, $e$, of all samples and planets. We compared the

<table><caption>**Table 11**
Evolution of the Semimajor Axes, $a$, and Eccentricities, $e$, from $10^4$ Samples over 10 Myr</caption><thead><tr><th>Planet</th><th>Δā</th><th>Δē</th><th>max(Δa)</th><th>max(Δe)</th></tr></thead><tbody><tr><td>b</td><td>-6.52e-09</td><td>1.73e-04</td><td>7.12e-07</td><td>0.0020</td></tr><tr><td>c</td><td>-1.44e-08</td><td>1.62e-04</td><td>2.51e-06</td><td>0.0018</td></tr><tr><td>d</td><td>-1.06e-08</td><td>1.45e-05</td><td>4.07e-06</td><td>0.0009</td></tr><tr><td>e</td><td>2.05e-08</td><td>4.44e-05</td><td>8.53e-06</td><td>0.0008</td></tr><tr><td>f</td><td>2.45e-07</td><td>5.13e-05</td><td>3.00e-05</td><td>0.0012</td></tr><tr><td>g</td><td>8.24e-08</td><td>5.01e-05</td><td>2.19e-05</td><td>0.0011</td></tr><tr><td>h</td><td>-1.23e-07</td><td>2.11e-04</td><td>3.00e-05</td><td>0.0035</td></tr></tbody></table>

**Note.** For each sample and planet, the difference of the average of $a$ and $e$ over the first and last Myr are computed as $\Delta\bar{a}$ and $\Delta\bar{e}$; we report the maximum over all samples. These numbers show that all samples remain stable over 10 Myr.

average values over the first Myr and the last Myr. Table 11 gives the average over all samples, and the maximum differences between the first and the last Myr. In all cases, the variations are small, $\le 0.002$. These results suggest that the simulations could be stable even on a much longer timescale. In addition, we have carried out long-term (50 Myr) integrations with tidal damping for two posterior samples, one with low and one with high values of the eccentricity of planet b. Using a range of values of tidal damping (from 1/10 to 100 times Earth's), we find in all cases that the system remained stable (using Posidonius; Bolmont et al. 2020).

More interesting is the evolution of the five GLR resonant angles, shown in Figure 26. In order to describe the evolution of the GLR angles, we define three categories:

1. Category I: remaining in GLR for 10 Myr, with a maximum difference from the initial value of less than 45°

2. Category II: remaining in GLR for 10 Myr, with a maximum difference from the initial value of more than 45°. In this category, the GLR angles can jump between different states.

3. Category III: not remaining in GLR for 10 Myr.

The threshold of 45° is chosen arbitrarily but is found to be practical to distinguish simulations where the GLR angles jump between different states (Category II), or remain in the same state (Category I). Figure 26 shows the three different categories in different colors, as well as a histogram of all 10,000 simulations over 10 Myr for all five GLR angles. The exact number of simulations in the three categories is given in Table 12. The GLR angles from planets b, c, and d as well as planets d, e, and f show a unique resonant state. Planets c, d, and e and planets e, f, and g have a dominant state and a subdominant state, while planets f, g, and h have a dominant state and two symmetric subdominant states. Our new samples show a better conservation of the GLR angles than was found in Grimm et al. (2018), where the longest resonance time was found to be 2 Myr.

## 10.6. Forecasts for JWST

### 10.6.1. Forecast Transit Times

With our transit-timing model, we can forecast the probabilities of future transit times and hence better help to plan transit observations with JWST. This is important for both optimizing the efficient use of the telescope and determining concurrent transits (i.e., two or more planets crossing the face of the star at the same time). This is especially important for transit transmission spectroscopy as the signal will be small, and hence many transits may need to be observed. With the
<figure><img src="image_29.png" /><figcaption>**Figure 26.** Evolution of the GLR angles $\phi$ for 10,000 samples over 10 Myr. The simulations can be split into three categories: x=category I remaining in resonance (black), II remaining in resonance but jump between states (dark blue), and III not remaining in resonance (light blue). </></figcaption></figure>

<table><caption>**Table 12**
Number of Posterior Samples Falling into the Three Resonant Categories for the Five GLR Angles</caption><thead><tr><th></th><th>bcd</th><th>cde</th><th>def</th><th>efg</th><th>gh</th></tr></thead><tbody><tr><td>Category I</td><td>130</td><td>178</td><td>755</td><td>7943</td><td>6462</td></tr><tr><td>Category II</td><td>1571</td><td>874</td><td>7653</td><td>1855</td><td>1578</td></tr><tr><td>Category III</td><td>8299</td><td>8948</td><td>1592</td><td>202</td><td>1960</td></tr></tbody></table>

Note. The total number of posterior samples is 10,000.

observation of initial transits with JWST, the ephemerides can be refined/updated; however, our current forecasts provide the starting point for planning JWST observations.

Table 15 gives our forecast for the upcoming times of transit through 2023 October to cover the first 2 yr of the JWST mission (six months after the end of Cycle 1, given the present launch date of 2021 October).

### 10.6.2. Simulated JWST TTV Analysis

Based on the measured properties of TRAPPIST-1, we have carried out a preliminary analysis forecasting future transit observations with the JWST. Already there are several JWST Guaranteed Time Observation (GTO) programs that plan to observe the TRAPPIST-1 planetary system, primarily for the purposes of spectroscopic characterization (GTO programs 1177, 1201, 1279 and 1331).<fnref n="33" /> It is very likely that additional observations will be scheduled during guest observing time throughout the duration of the JWST mission as the detection of spectroscopic features requires observations of multiple transits for each of the planets (Barstow & Irwin 2016; Morley et al. 2017; Lustig-Yaeger et al. 2019; Fauchez et al. 2020). An effort to coordinate these observations among the exoplanet and planetary science communities is underway via the TRAPPIST-1 JWST Community Initiative (Gillon et al. 2020). All to say, long-term studies of TRAPPIST-1 for spectroscopy will also

yield transit times for each transit observed, enabling a transit-timing analysis of the results.

To estimate the maximum possible precision of observations with JWST, we have simulated a five-year program in which every transit of every planet in TRAPPIST-1 is observed with NIRSPEC (Birkmann et al. 2016). The NIRSPEC instrument was chosen as its prism mode covers 0.5–5 µm, covering the peak of the SED of the star and thus maximizing the number of photons detected, which is about two orders of magnitude per transit greater than that collected by Spitzer. Although such a complete set of transits will be impossible to collect (thanks to limits due to scheduling and time-allocation), this analysis yields an estimate of the most optimistic results we might expect from JWST.

We have carried out simulations of transits of each of the planets as observed by NIRSPEC. We include realistic estimates of photon noise and correlated stellar variability based on the pattern of variations detected with the Spitzer Space Telescope, using a Gaussian Process model created with *celerite* (Foreman-Mackey et al. 2017). We do not include instrumental systematics under the assumption that over the timescales of ingress/egress, which are what limit the timing precision, the noise contribution will be dominated by photon noise and stellar variations. From these simulations, we found that the posterior timing precision ranges from 0.6 to 1.7 s per transit, much more precise than the measurements reported in the present paper.

Next, we created a simulated set of transit-timing observations at the two windows each year when the TRAPPIST-1 system is observable with JWST (Figure 27). For each transit time, we drew the time from the distribution of uncertainties from the posteriors of the simulated transit data.

Finally, we utilized our code for transit-timing analysis to optimize a plane-parallel model with seven planets. At the maximum likelihood of the fit, we computed the Hessian to estimate the uncertainties on the model parameters. Figure 27 shows the simulated transit-timing observations with JWST. This includes about 600 transits observed with the telescope (again, the maximum possible over the nominal 5 yr JWST mission). Figure 28 shows the results of the mass measurements in the simulations. We find that the masses can be recovered to better than 0.02% for planets d–h and to 0.1% for planets b and c.

Of course, it will be impossible to arrange such a large number of transit observations of this system. But, even if the number of observations is an order of magnitude smaller, we expect that the signal-to-noise should scale with the square root of the number of measurements made, and thus the outer planets will still have mass measurements precise to the order of a part per thousand.

### 10.7. Stellar Parameters

The stellar density we derive using the photodynamic model, $\rho_{*} = 53.17^{+0.72}_{-1.18}\rho_{\odot}$, is in $1\sigma$ agreement with prior analyses. Most recently, Delrez et al. (2018b) found a density of $\rho_{*} = (52.3 \pm 2.2)\rho_{\odot}$, twice as uncertain as our analysis. Our approach yields a density of superior precision due to several factors. The transit times in the Spitzer data are constrained by all of the measured transits in the photodynamic model so that fewer degrees of freedom are needed to fit the times (37 free parameters in the N-body model versus 447 transit times fit to each transit).

The stellar mass we take from the analysis by Mann et al. (2019), $M_{*} = 0.0898 \pm 0.0023 M_{\odot}$.<fnref n="34" /> This mass has a

<footnote n="33">³³ For specifications of these programs, see https://www.stsci.edu/jwst/observing-programs/approved-gto-programs.</footnote>

<footnote n="34">³⁴ https://github.com/awmann/M_-M_K-</footnote>
<figure><img src="image_30.png" /><figcaption>**Figure 27.** Simulated observations of all of the transits of TRAPPIST-1 detectable with JWST. Each transit has an uncertainty of ≈0.6–1.7 s, assumed to be observed with NIRSPEC (which maximizes the number of photons collected by any JWST instrument). From retrieval, we obtain ≲0.1% mass precision. <//></figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure 28.** Simulated planet masses based on 5 yr of JWST observations of every TRAPPIST-1 transit with NIRSPEC. The recovered mass ($M_{\text{out}}$) minus the input mass vs. the input mass ($M_{\text{in}}$). The masses relative to the star can be recovered to better than 0.1% precision. <//></figcaption></figure>

precision of 2.6%, which limits the mass precision for several of the planets. We are at the point that to improve the mass measurements of the planets we will need to improve the measurements of the star.

We used the luminosity estimate from Ducrot et al. (2020), which is slightly lower than that estimated by Gonzales et al. (2019) due to a difference in the measured bolometric flux. We are consistent with Gonzales et al. (2019) for the reported value of $R^2T_{\text{eff}}^4$ at $1\sigma$, while our $T_{\text{eff}}$ is more precise (28 K versus 42 K), R is 2.5 times more precise, and our log g is more precise by an order of magnitude.

## 11. Conclusions

The Spitzer discovery of seven transiting planets orbiting the TRAPPIST-1 star by Gillon et al. (2017) promised the determination of the interior compositions of these planets via dynamical analysis. We have now analyzed the complete set of transit-time measurements of the TRAPPIST-1 planets from Spitzer, augmented by additional transits from the ground, K2, and HST. Our primary conclusions are:

1. We have measured the masses, radii and densities to high fractional precision, 1%–8%, based on an *N*-body model and a photodynamical model with seven planets. This improves upon RV current precision by up to two orders of magnitude.

2. The pattern of masses and radii may be consistent with a uniform planetary composition for all seven planets that have lower uncompressed densities than Earth, Mars, or Venus, with weaker evidence for a declining normalized density with orbital period (88% confidence). The planet
properties may either be consistent with a CMF of
$21 \pm 4$ wt%, or an Earth-like core and mantle with a
surface water content that varies from <0.001% for the
inner three planets to ≈5% for the outer four, or core-free
planets with highly oxidized iron in the mantle that
elevates the interior light-element content. These are not
unique explanations.

3. The planets appear to be dynamically cold, with eccentricities less than ≈1%, and inclinations that may be coplanar to a few hundredths of a degree.

4. The system is stable on long timescales and shows a pattern of generalized Laplace resonances with angles that match predictions from migration simulations of Mah (2018).

5. We provide a forecast of the future times of transit for the planets (Table 15) to help in planning observations with JWST, which may yield more precise constraints upon the planets' masses.

6. We have yet to find strong evidence for an eighth planet.

Based upon these properties, we next speculate on some
possible scenarios for the formation and evolution of the system.

11.1. *Expectations for the Compositions of the TRAPPIST-1
Planets from Formation Scenarios*

As mentioned, our analysis suggests that the TRAPPIST-1
planets have somewhat lower uncompressed bulk densities than
Earth (see Table 6 and Figure 12). It is possible that these lower
densities result from a deficit of high-density material (e.g., less
iron) relative to Earth, or an excess of low-density material
(e.g., having more water), or both; in this section, we speculate
about formation scenarios that may be consistent with these
planets' bulk densities.

In general, planets that formed within the same protoplanet-
tary disk are expected to have similar budgets in relative
refractory elements (Bond et al. 2010; Elser et al. 2012) but can
have very different volatile element budgets (Öberg &
Bergin 2016). Similar relative refractory elements (Fe, Mg,
Si) implies similar CMFs for all seven planets, assuming full
differentiation. As suggested by Dorn et al. (2018), the
refractory composition may best be described by studying the
densest planet of the system, planet c with 22%-31% CMF.
Thus, with this assumption, all of the planets may likely have a
22%-31% CMF but different light-element mass fractions
(which may increase slightly with orbital period, Figure 19).

Is an overall CMF of 22%-31% realistic for terrestrial planet
interiors? This range of CMF implies lower Fe/Mg and Fe/Si
values compared to Earth (and the Sun). Elemental abundances of
rocky interiors are expected to be reflected in the photospheric
abundance of the host star as argued by Unterborn et al. (2018)
and Dorn et al. (2018). Unfortunately, measuring the photospheric
abundances of this cool and active host star remains very
challenging. However, Unterborn et al. (2018) estimated the
stellar molar Fe/Mg number ratio to be 0.75 ± 0.2 by analyzing
Sun-like stars of similar metallicity to TRAPPIST-1, which may
be slightly lower than the solar value. This corresponds to a CMF
of 24%-35% for a fully differentiated model. The corresponding
mass-radius curve for a rocky interior of this range of Fe/Mg
values is plotted in Figure 12 (gray curve and shaded region). It
overlaps well with the densest planets, c and b. This means that
the expected range of stellar abundances supports a possible
overall CMF value of 22%-31%, assuming full differentiation.

Could there be a variation of Fe/Mg ratios among the planets? Rocky planet accretion should preserve the integrated iron/rock ratio. Consider a population of planetary embryos and planetesimals that accrete into a system of rocky planets. Giant collisions between growing planetary embryos can change the iron/rock ratios of individual objects by preferentially stripping the outer, rock-dominated layers from differentiated embryos (e.g., Benz et al. 1988; Marcus et al. 2010; Asphaug & Reufer 2014). But from a system-wide perspective, it is a zero-sum game unless rock or iron is preferentially lost from all of the planets. Rock is the major component of loosely bound impact debris and more likely to be lost either by differential aerodynamic drag (Weidenschilling 1977) or solar wind drag (Spalding & Adams 2020), and so the integrated iron/rock ratio should only increase. Hypothetical variations in Fe/Mg can otherwise be caused if large portions of planetary building blocks condense at different high temperatures (>1200 K). During planet formation, such temperatures are only reached in a tiny region very close to the ultracool dwarf star. Consequently, both Unterborn et al. (2018) and Dorn et al. (2018) have assumed that all seven planets have similar refractory element ratios (i.e., Fe/Si, Fe/Mg). Whether rocky planets can have a wider compositional distribution than stars remains to be seen (Plotnykov & Valencia 2020).

Alternatively, the lower measured bulk densities of the TRAPPIST-1 planets relative to Earth-like composition might be explained by core-free interiors (Elkins-Tanton & Sea- ger 2008) in which the oxygen content is high enough such that all iron is oxidized. If the refractory elements (Mg, Fe, Si) follow solar abundances, a fully oxidized interior would contain about 38.2 wt% of oxygen, which lies between the value for Earth (29.7 wt%) and CI chondrites (45.9 wt%). Such an interior scenario can easily describe the observed bulk densities (red line in Figure 12), and this may bolster the long- range migration scenario in which the planets formed in a highly oxidizing environment that enabled the iron to remain in the mantle even after migration. Based on the elemental composition, these models have an oxygen fugacity of ΔIW = -0.91,<fnref n="35" /> which is more oxidized than Earth or even Mars but is comparable to the oxidation state of small bodies, both in our solar system and accreted by white dwarfs (Doyle et al. 2019).

However, the evidence for a core-free planet may rest on
knowing the refractory abundances of the TRAPPIST-1 host
star, which have yet to be constrained. Alas, our interpretation
of the planets' compositions may be limited by our imprecise
knowledge of the host star: its radius, its mass, its photospheric
inhomogeneity, and its refractory abundances all affect our
measurement and interpretation of the masses, radii, and
compositions of the TRAPPIST-1 planets. In this paper, our
measurements of the relative planetary radii and masses have
reached such a precision that the fault may now lie in the star.

11.2. Future Work

We conclude by pointing out directions for building upon
the work described in this paper:

1. We have yet to identify the origin of timing outliers that show an excess relative to a normal distribution. This

<footnote n="35">³⁵ Oxygen fugacity is stated relative to the iron–wüstite equilibrium reaction Fe + 0.5O₂ = FeO (wüstite) such that ΔIW = log(fO₂/rock – log(fO₂/Iw).</footnote>
may be addressed with higher precision measurements that may be able to identify a source of noise responsible for these outliers.

2. Our analysis assumes a plane-parallel system with seven planets and does not yet couple the dynamical and photometric analysis (our photodynamics held the dynamical model fixed). Future analysis with a fully coupled photodynamical model with 3D orbits and more than seven planets may be valuable.

3. We need more transits measured for planets d and h in order to better measure the amplitude and phase on the TTV timescale, as well as to better constrain the presence of planets beyond h.

4. The interpretation of the compositions of the planets is limited by the unknown composition of the host star. A measurement of the Mg/Fe and Fe/Si ratios would help to interpret the core and mantle compositions. Both sets of constraints would help to limit the range and break degeneracies of possible interior compositions of the planets (Dorn et al. 2015; Bitsch & Battistini 2019).

5. Without a constraint on the detailed abundance ratios of the host star, a Bayesian interpretation of the bulk densities of the planets should be warranted (Dorn et al. 2016) to better quantify the range of possible compositions.

6. More detailed spectral analysis of the stellar photosphere to ascertain the impact of an inhomogeneous stellar atmosphere on the radius ratios would be warranted.

We anticipate that once JWST launches, we will obtain higher precision constraints upon the dynamics of the system, yielding much improved constraints upon the planets' bulk densities, which will further improve the interpretation of their interior compositions.

This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Support for this work was provided by NASA through an award issued by JPL/Caltech. E.A. was supported by a Guggenheim Fellowship and NSF grant AST-1615315. This research was partially conducted during the Exostar19 program at the Kavli Institute for Theoretical Physics at UC Santa Barbara, which was supported in part by the National Science Foundation under grant No. NSF PHY-1748958. We also acknowledge support from NASA's NExSS Virtual Planetary Laboratory, funded under NASA Astrobiology Institute Cooperative Agreement Number NNA13AA93A, and the NASA Astrobiology Program grant 80NSSC18K0829. This work was facilitated though the use of the advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington. TRAPPIST is a project funded by the Belgian Fonds (National) de la Recherche Scientifique (F. R.S.-FNRS) under grant FRFC 2.5.594.09.F. TRAPPIST-North is a project funded by the University of Liège, in collaboration with Cadi Ayyad University of Marrakech (Morocco). B.-O.D., C.D., and J.H. acknowledge support from the Swiss National Science Foundation (grants PP00P2-163967, PZ00P2_174028, and 200020_19203, respectively). Calculations were performed on UBELIX (http://www.id.unibe.ch/hpc), the HPC cluster at the University of Bern. This work has been carried out in the framework of the PlanetS

National Centre of Competence in Research (NCCR) supported by the Swiss National Science Foundation (SNSF). This research has made use of NASA's Astrophysics Data System and of services produced by the NASA Exoplanet Science Institute at the California Institute of Technology. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 832738/ESCAPE and European Research Council (ERC; grant agreement Nos. 679030/WHIPLASH and 803193/BEBOP). M.T. thanks the Gruber Foundation for its generous support to this research. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) ERC Grant Agreement No. 336480 and from the ARC grant for Concerted Research Actions, financed by the Wallonia-Brussels Federation. V.V.G. is an F.R.S.-FNRS Research Associate. M.G. and E.J. are F.R.S.-FNRS Senior Research Associates. R.M. acknowledges support for this research from NASA (grant 80NSSC18K0397). Z.L. acknowledges support from the Washington NASA Space Grant Consortium Summer Undergraduate Research Program.

We thank Trevor Branch for discussions about HMC and automatic differentiation. We thank Pramod Gupta, Diana Windemuth, and Tyler Gordon for help in using Hyak. We thank Vardan Adibekyan, Rory Barnes, Jim Davenport, Natalie Hinkel, Dave Joswiak, and Sarah Millholland for useful discussions. Finally, we thank the referees for valuable feedback which improved this paper.

*Software:* Matplotlib (Hunter 2007; Caswell et al. 2020), Julia (Bezanson et al. 2017), Limbdark.jl (Agol et al. 2020).

## Appendix A
### Approximate Hessian Matrix

Here we approximate the posterior probability distribution as a multidimensional Gaussian, assuming a uniform prior. The log likelihood for each data point with indices i and j may be written as a function of the observed transit times and uncertainties, the modeled transit times, and the Student's t-distribution model parameters, such that

$$ \mathcal{L}_{ij}(\mathbf{x}) = \mathcal{L}_{ij}(t_{ij}(\mathbf{x}_{\text{dyn}}); \log \nu, V_1 e^{1/2\nu}; t_{ij,\text{obs}}, \sigma_{ij}), \quad (A1) $$

where all of the dependence on the dynamical model parameters enters through $t_{ij}(\mathbf{x}_{\text{dyn}})$. The maximum posterior probability also corresponds to the maximum likelihood in this limit, in which case we expand the log likelihood for the ith planet and jth transit as a Taylor series:

$$ \log \mathcal{L}_{ij}(\mathbf{x}) \approx \log \mathcal{L}_{ij}(\mathbf{x}_0) + \frac{1}{2} \sum_{k,l} \frac{\partial^2 \log \mathcal{L}_{ij}(\mathbf{x})}{\partial x_k \partial x_l} \Big|_{\mathbf{x}_0} (x_k - x_{0,k})(x_l - x_{0,l}), \quad (A2) $$

where we have used the fact that the gradient of the log likelihood vanishes at the maximum-likelihood value of the model parameters, $\mathbf{x}_0$, and the indices $k, l = 1, \dots, 5N_p + 2$ for $x_k$ and $x_l$, where the first $5N_p$ parameters are the dynamical parameters, $\mathbf{x}_{\text{dyn}}$, and the last two parameters are the Student's t-distribution likelihood parameters, $\log \nu$ and $V_1 e^{1/2\nu}$. Now, the width of the Gaussian distribution at the maximum likelihood is governed by the Hessian matrix, with elements
given by

$$ \mathcal{H}_{kl}(\mathbf{x}_0) = -\sum_{i,j} \left. \frac{\partial^2 \log \mathcal{L}_{ij}(\mathbf{x})}{\partial x_k \partial x_l} \right|_{\mathbf{x}_0}, \quad (A3) $$

which involves second derivatives of the negative log likelihood with respect to the model parameters. The derivatives of $t_{ij}$ with respect to $\mathbf{x}_{dyn}$ we compute with the NbodyGradient code; however, the second derivatives of the transit times with respect to the dynamical model parameters are not computed with our N-body code. We drop these transit-time second derivative terms, which we justify as follows.

For the Hessian matrix elements that involve second derivatives with respect to both dynamical model parameters, $1 \le k, l \le 5N_p$, we can write

$$ \frac{\partial^2 \log \mathcal{L}_{ij}(\mathbf{x})}{\partial x_k \partial x_l} = \frac{\partial^2 \log \mathcal{L}_{ij}(\mathbf{x})}{\partial^2 t_{ij}} \times \frac{\partial t_{ij}}{\partial x_k} \frac{\partial t_{ij}}{\partial x_l} + \frac{\partial \log \mathcal{L}_{ij}(\mathbf{x})}{\partial t_{ij}} \frac{\partial^2 t_{ij}}{\partial x_k \partial x_l}, \quad (A4) $$

where $t_{ij} = t_{ij}(\mathbf{x}_{dyn})$ is implied in this and subsequent equations.

Now, at the maximum likelihood, there is a balance of residuals that are both positive and negative, such that the second component of this equation has terms with positive and negative signs for different values of i and j. This causes the second term in this equation to average to a small value compared with the first term when the sum is carried out over i and j (the planet and transit indices). So, we drop the second term in this equation.

Adding in the cases of the Hessian matrix elements that involve the likelihood parameters, ($x_{5N_p+1}$, $x_{5N_p+2}$) = ($\log \nu$, $V_1 e^{1/2\nu}$), we compute the Hessian as

$$ \mathcal{H}_{k,l}(\mathbf{x}) = -\sum_{i,j} \left\{ \begin{array}{ll} \frac{\partial^2 \log \mathcal{L}_{ij}(\mathbf{x})}{\partial t_{ij} \partial t_{ij}} & 1 \le k, l \le 5N_p \\ \frac{\partial^2 \log \mathcal{L}_{ij}(\mathbf{x})}{\partial t_{ij} \partial x_l} \frac{\partial t_{ij}}{\partial x_k} & 1 \le k \le 5N_p; l > 5N_p \\ \frac{\partial^2 \log \mathcal{L}_{ij}(\mathbf{x})}{\partial t_{ij} \partial x_k} \frac{\partial t_{ij}}{\partial x_l} & k > 5N_p; 1 \le l \le 5N_p \\ \frac{\partial^2 \log \mathcal{L}_{ij}(\mathbf{x})}{\partial x_k \partial x_k} & k, l > 5N_p \end{array} \right. \quad (A5) $$

where the partial derivatives with respect to $t_{ij}(\mathbf{x}_{dyn})$, $x_{5N_p+1} = \log \nu$, and $x_{5N_p+2} = V_1 e^{1/2\nu}$ are computed with automatic differentiation.

The inverse of the Hessian matrix is used in the Levenberg-Marquardt optimization, and when evaluated at the maximum likelihood, is used to estimate the covariance matrix from which the square root of the diagonal components are used to estimate the widths of the posterior distribution for each model parameter, $\mathbf{x} = (\mathbf{x}_{dyn}, \log \nu, V_1 e^{1/2\nu})$, which are plotted in Figures 6, 5, and 4. This approximated Hessian is also used to define the mass matrix for the HMC simulations.

## Appendix B: Transit-timing Prior

We use a uniform prior for each mass and orbital element, with smooth bounds on each, with the exception of the initial eccentricity vectors. Because we sample in the eccentricity vector of each planet, $e_i = (e_i \cos \omega_i, e_i \sin \omega_i)$, the volume of parameter space scales $\propto e_i$, and so a $1/e_i$ prior is needed to yield a posterior that has a uniform probability with eccentricity, $e_i$, for the ith planet (Eastman et al. 2013).

In addition to the eccentricity prior, we place smooth bounds on the parameters. For each bound, we choose upper and lower limits for which the prior starts to transition from 1 to 0 with a cubic dependence. For the bound on a function of our parameters of value $\xi$, we specify

$$ \Pi_{\text{bound}}(\xi) = \begin{cases} 0 & \xi \le \xi_1 \\ 3y^2 - 2y^3; \quad y = \frac{\xi - \xi_1}{\xi_2 - \xi_1} & \xi_1 < \xi < \xi_2 \\ 1 & \xi_2 \le \xi \le \xi_3, \\ 3y^2 - 2y^3; \quad y = \frac{\xi_4 - \xi}{\xi_4 - \xi_3} & \xi_3 < \xi < \xi_4 \\ 0 & \xi \ge \xi_4 \end{cases} \quad (B6) $$

so that the total prior is given by

$$ \Pi(\mathbf{x}) = \prod_{i=1}^{N_p} e_i^{-1} \times \prod_{j=1}^{N_{\text{bound}}} \Pi_{\text{bound}}(f_j(\mathbf{x})), \quad (B7) $$

where the values of $\xi_1 - \xi_4$ and each transformation of parameters, $f = \{f_j(\mathbf{x})\; j=1, \dots, N_{\text{bound}}\}$ are given in Table 13, where $N_{\text{bound}} = 4N_p + 2$. The prior probability, then, is given by $\Pi(\mathbf{x})$, which we multiply by the likelihood function before sampling.
<table><caption>**Table 13**
Prior Probability Boundary Limits for the TRAPPIST-1 Planet Parameters</caption><thead><tr><th>Parameter</th><th>Bounds Function f<sub>j</sub></th><th>Lower Bound</th><th>Lower Transition</th><th>Upper Transition</th><th>Upper Bound</th></tr></thead><tbody><tr><td>Mass ratio</td><td>log<sub>10</sub> μ</td><td>-8</td><td>-7</td><td>-3</td><td>-2</td></tr><tr><td>Eccentricity</td><td>e</td><td>...</td><td>...</td><td>0.2</td><td>0.3</td></tr><tr><td>Period of b</td><td>P<sub>b</sub> (days)</td><td>1.49</td><td>1.50</td><td>1.52</td><td>1.53</td></tr><tr><td>Period of c</td><td>P<sub>c</sub> (days)</td><td>2.40</td><td>2.41</td><td>2.43</td><td>2.44</td></tr><tr><td>Period of d</td><td>P<sub>d</sub> (days)</td><td>4.03</td><td>4.04</td><td>4.06</td><td>4.07</td></tr><tr><td>Period of e</td><td>P<sub>e</sub> (days)</td><td>6.08</td><td>6.09</td><td>6.11</td><td>6.12</td></tr><tr><td>Period of f</td><td>P<sub>f</sub> (days)</td><td>9.18</td><td>9.19</td><td>9.22</td><td>9.23</td></tr><tr><td>Period of g</td><td>P<sub>g</sub> (days)</td><td>12.33</td><td>12.34</td><td>12.36</td><td>12.37</td></tr><tr><td>Period of h</td><td>P<sub>h</sub> (days)</td><td>18.75</td><td>18.76</td><td>18.78</td><td>18.79</td></tr><tr><td>Initial transit time of b</td><td>t<sub>0,b</sub> - 2,457,257 (days)</td><td>0.53</td><td>0.54</td><td>0.56</td><td>0.57</td></tr><tr><td>Initial transit time c</td><td>t<sub>0,c</sub> - 2,457,257 (days)</td><td>0.57</td><td>0.58</td><td>0.60</td><td>0.61</td></tr><tr><td>Initial transit time d</td><td>t<sub>0,d</sub> - 2,457,257 (days)</td><td>0.05</td><td>0.06</td><td>0.08</td><td>0.09</td></tr><tr><td>Initial transit time e</td><td>t<sub>0,e</sub> - 2,457,257 (days)</td><td>0.81</td><td>0.82</td><td>0.84</td><td>0.85</td></tr><tr><td>Initial transit time f</td><td>t<sub>0,f</sub> - 2,457,257 (days)</td><td>0.05</td><td>0.06</td><td>0.08</td><td>0.09</td></tr><tr><td>Initial transit time g</td><td>t<sub>0,g</sub> - 2,457,257 (days)</td><td>0.70</td><td>0.71</td><td>0.73</td><td>0.74</td></tr><tr><td>Initial transit time h</td><td>t<sub>0,h</sub> - 2,457,249 (days)</td><td>0.58</td><td>0.59</td><td>0.61</td><td>0.62</td></tr><tr><td>Degrees of freedom</td><td>ν</td><td>0.5</td><td>1.0</td><td>50</td><td>100</td></tr><tr><td>Log variance factor</td><td>log V<sub>1</sub></td><td>-2</td><td>-1</td><td>5</td><td>10</td></tr></tbody></table>

Note. The bounds are chosen so as to not affect the parameters as much as possible.
## Appendix C
### Corner Plots

Figures 29 and 30 show corner plots of the variables from the transit-timing and photodynamical analyses, respectively.

<figure><img src="image_32.png" /><figcaption>**Figure 29.** Corner plot of variables in the transit-timing analysis with 1σ and 2σ confidence contours. Lagrangian orbital elements are defined as $k_b = e_b \cos \omega_b$ and $h_b = e_b \sin \omega_b$, and similarly for the other planets. Planet masses are defined relative to the star. </></figcaption></figure>
<figure><img src="image_33.png" /><figcaption>**Figure 30.** Corner plot of variables in the photodynamical analysis with 1σ and 2σ confidence contours. Planet radii are relative to star. </></figcaption></figure>
## Appendix D
Tables

Tables of the best-fit transit times (Table 14) and the forecast times (Table 15) are given in this appendix.

<table><caption>**Table 14**
Observed Transit Times with Uncertainties, along with the Mean, $t_{\text{post}}$, and Standard Deviation, $\sigma_{\text{post}}$, of the Times from the Posterior Sample</caption><thead><tr><th>Planet</th><th>Epoch</th><th>t<sub>obs</sub></th><th>σ<sub>obs</sub></th><th>t<sub>post</sub></th><th>σ<sub>post</sub></th><th>Source</th></tr></thead><tbody><tr><td>b</td><td>0</td><td>7322.515 310</td><td>0.000 710</td><td>7322.517 902</td><td>0.000 127</td><td>TS</td></tr><tr><td>c</td><td>0</td><td>7282.805 700</td><td>0.001 400</td><td>7282.805 840</td><td>0.000 248</td><td>TS</td></tr><tr><td>d</td><td>0</td><td>7560.797 300</td><td>0.002 300</td><td>7560.801 847</td><td>0.000 384</td><td>TS</td></tr><tr><td>e</td><td>0</td><td>7312.713 000</td><td>0.002 700</td><td>7312.713 920</td><td>0.000 360</td><td>TS</td></tr><tr><td>f</td><td>0</td><td>7321.525 200</td><td>0.002 000</td><td>7321.522 338</td><td>0.000 781</td><td>TS</td></tr><tr><td>g</td><td>0</td><td>7294.786 000</td><td>0.003 900</td><td>7294.772 215</td><td>0.000 905</td><td>TS</td></tr><tr><td>h</td><td>0</td><td>7662.554 360</td><td>0.002 000</td><td>7662.550 913</td><td>0.001 498</td><td>Spitzer</td></tr></tbody></table>

**Note.** Times are in BJD<sub>TDB</sub>-2,450,000 while uncertainties are in days.

(This table is available in its entirety in machine-readable form.)

<table><caption>**Table 15**
Mean, $t_{\text{post}}$, and Standard Deviation, $\sigma_{\text{post}}$, of Forecast Times from the Posterior Sample</caption><thead><tr><th>Planet</th><th>Epoch</th><th>t<sub>post</sub></th><th>σ<sub>post</sub></th></tr></thead><tbody><tr><td>b</td><td>1</td><td>7259.061 221</td><td>0.000 152</td></tr><tr><td>b</td><td>2</td><td>7260.572 662</td><td>0.000 163</td></tr><tr><td>b</td><td>3</td><td>7262.083 452</td><td>0.000 159</td></tr><tr><td>b</td><td>4</td><td>7263.594 099</td><td>0.000 150</td></tr><tr><td>b</td><td>5</td><td>7265.104 973</td><td>0.000 151</td></tr><tr><td>b</td><td>6</td><td>7266.615 845</td><td>0.000 153</td></tr><tr><td>b</td><td>7</td><td>7268.126 697</td><td>0.000 155</td></tr><tr><td>b</td><td>8</td><td>7269.637 314</td><td>0.000 148</td></tr><tr><td>b</td><td>9</td><td>7271.148 093</td><td>0.000 146</td></tr><tr><td>b</td><td>10</td><td>7272.659 543</td><td>0.000 157</td></tr><tr><td>b</td><td>11</td><td>7274.170 344</td><td>0.000 152</td></tr><tr><td>b</td><td>12</td><td>7275.680 980</td><td>0.000 144</td></tr><tr><td>b</td><td>13</td><td>7277.191 799</td><td>0.000 143</td></tr><tr><td>b</td><td>14</td><td>7278.702 670</td><td>0.000 145</td></tr><tr><td>b</td><td>15</td><td>7280.213 581</td><td>0.000 150</td></tr><tr><td>b</td><td>16</td><td>7281.724 209</td><td>0.000 143</td></tr><tr><td>b</td><td>17</td><td>7283.234 973</td><td>0.000 139</td></tr><tr><td>b</td><td>18</td><td>7284.746 420</td><td>0.000 150</td></tr><tr><td>b</td><td>19</td><td>7286.257 235</td><td>0.000 146</td></tr><tr><td>b</td><td>20</td><td>7287.767 862</td><td>0.000 137</td></tr><tr><td>b</td><td>21</td><td>7289.278 632</td><td>0.000 136</td></tr><tr><td>b</td><td>22</td><td>7290.789 504</td><td>0.000 138</td></tr><tr><td>b</td><td>23</td><td>7292.300 472</td><td>0.000 144</td></tr><tr><td>b</td><td>24</td><td>7293.811 108</td><td>0.000 137</td></tr><tr><td>b</td><td>25</td><td>7295.321 861</td><td>0.000 133</td></tr><tr><td>b</td><td>26</td><td>7296.833 297</td><td>0.000 143</td></tr><tr><td>b</td><td>27</td><td>7298.344 125</td><td>0.000 140</td></tr><tr><td>b</td><td>28</td><td>7299.854 748</td><td>0.000 131</td></tr><tr><td>b</td><td>29</td><td>7301.365 471</td><td>0.000 129</td></tr><tr><td>b</td><td>30</td><td>7302.876 339</td><td>0.000 131</td></tr></tbody></table>

**Note.** Times are in BJD<sub>TDB</sub>-2,450,000 while uncertainties are in days. Thirty lines are previewed; full table is available in machine-readable format.

(This table is available in its entirety in machine-readable form.)
ORCID iDs

Eric Agol https://orcid.org/0000-0002-0802-9145

Caroline Dorn https://orcid.org/0000-0001-6110-4610

Simon L. Grimm https://orcid.org/0000-0002-0632-4407

Martin Turbet https://orcid.org/0000-0003-2260-9856

Elsa Ducrot https://orcid.org/0000-0002-7008-6888

Laetitia Delrez https://orcid.org/0000-0001-6108-4808

Michaël Gillon https://orcid.org/0000-0003-1462-7739

Brice-Olivier Demory https://orcid.org/0000-0002-9355-5165

Artem Burdanov https://orcid.org/0000-0001-9892-2406

Khalid Barkaoui https://orcid.org/0000-0003-1464-9276

Zouhair Benkhaldoun https://orcid.org/0000-0001-6285-9847

Emeline Bolmont https://orcid.org/0000-0001-5657-4503

Adam Burgasser https://orcid.org/0000-0002-6523-9536

Sean Carey https://orcid.org/0000-0002-0221-6871

Julien de Wit https://orcid.org/0000-0003-2415-2191

Daniel Fabrycky https://orcid.org/0000-0003-3750-0183

Daniel Foreman-Mackey https://orcid.org/0000-0002-9328-5652

Jonas Haldemann https://orcid.org/0000-0003-1231-2389

David M. Hernandez https://orcid.org/0000-0001-7648-0926

James Ingalls https://orcid.org/0000-0003-4714-1364

Emmanuel Jehin https://orcid.org/0000-0001-8923-488X

Zachary Langford https://orcid.org/0000-0001-7574-4440

Jérémy Leconte https://orcid.org/0000-0002-3555-480X

Susan M. Lederer https://orcid.org/0000-0003-2805-8653

Rodrigo Luger https://orcid.org/0000-0002-0296-3826

Renu Malhotra https://orcid.org/0000-0002-1226-3305

Victoria S. Meadows https://orcid.org/0000-0002-1386-1710

Brett M. Morris https://orcid.org/0000-0003-2528-3409

Francisco J. Pozuelos https://orcid.org/0000-0003-1572-7707

Didier Queloz https://orcid.org/0000-0002-3012-0316

Sean N. Raymond https://orcid.org/0000-0001-8974-0758

Franck Selsis https://orcid.org/0000-0001-9619-5356

Marko Sestovic https://orcid.org/0000-0002-8124-8360

Amaury H. M. J. Triaud https://orcid.org/0000-0002-5510-8751

Valerie Van Grootel https://orcid.org/0000-0003-2144-4316

References

Agol, E., & Fabrycky, D. C. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte (Cham: Springer), 7

Agol, E., Luger, R., & Foreman-Mackey, D. 2020, AJ, 159, 123

Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989

Allard, F., Homeier, D., & Freytag, B. 2011, in ASP Conf. Ser. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. C. M. Johns-Krull, M. K. Browning, & A. A. West (San Francisco, CA: ASP), 91

Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

Asphaug, E., & Reufer, A. 2014, NatGe, 7, 564

Barkaoui, K., Burdanov, A., Hellier, C., et al. 2019, AJ, 157, 43

Barr, A. C. 2016, IGRE, 121, 1573

Barstow, J. K., & Irwin, P. G. J. 2016, MNRAS, 461, L92

Batalha, N. E., Lewis, N. K., Line, M. R., Valenti, J., & Stevenson, K. 2018, ApJL, 856, L34

Benz, W., Slattery, W. L., & Cameron, A. G. W. 1988, Icar, 74, 516

Berger, T. A., Huber, D., Gaidos, E., & van Saders, J. L. 2018, ApJ, 866, 99

Betancourt, M. 2017, arXiv:1701.02434

Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2017, SIAMR, 59, 65

Birkmann, S. M., Ferruit, P., Rawle, T., et al. 2016, Proc. SPIE, 9904, 99040B

Bitsch, B., & Battistini, C. 2019, A&A, 633, A10

Bolmont, E., Demory, B.-O., Blanco-Cuaresma, S., et al. 2020, A&A, 635, A117

Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728

Bond, J. C., Lauretta, D. S., & Orien, D. P. 2010, Icar, 205, 321

Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121

Brasser, R., Barr, A. C., & Dobos, V. 2019, MNRAS, 487, 34

Brown, B. P., Oishi, J. S., Vasil, G. M., Lecoanet, D., & Burns, K. J. 2020, ApJL, 902, L3

Burdanov, A., Delrez, L., Gillon, M., & Jehin, E. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte (Cham: Springer), 130

Burdanov, A. Y., Lederer, S. M., Gillon, M., et al. 2019, MNRAS, 487, 1634

Burgasser, A. J., & Mamajek, E. E. 2017, ApJ, 845, 110

Carey, S. J., Lacy, M. D., Laine, S. J., et al. 2004, Proc. SPIE, 5487, 211

Carter, J. A., Agol, E., Chaplin, W. J., et al. 2012, Sci, 337, 556

Casali, M., Adamson, A., de Oliveira, C. A., et al. 2007, A&A, 467, 777

Caswell, T. A., Droettboom, M., Lee, A., et al. 2020, matplotlib/matplotlib: REL: v3.3.1, Zenodo, doi:10.5281/ZENODO.592536

Christiansen, J. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte (Cham: Springer), 150

Claret, A. 2018, A&A, 618, A20

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

Dai, F., Masuda, K., Winn, J. N., & Zeng, L. 2019, ApJ, 883, 79

de Wit, J., Wakeford, H. R., Gillon, M., et al. 2016, Natur, 537, 69

de Wit, J., Wakeford, H. R., Lewis, N. K., et al. 2018, NatAs, 2, 214

Deck, K. M., & Agol, E. 2015, ApJ, 802, 116

Dehnen, W., & Hernandez, D. M. 2017, MNRAS, 465, 1201

Delisle, J.-B. 2017, A&A, 605, A96

Delrez, L., Gillon, M., Queloz, D., et al. 2018a, Proc. SPIE, 10700, 107001I

Delrez, L., Gillon, M., Triaud, A. H. M. J., et al. 2018b, MNRAS, 475, 3577

Dorn, C., Khan, A., Heng, K., et al. 2015, A&A, 577, A83

Dorn, C., Mosegaard, K., Grimm, S. L., & Alibert, Y. 2018, ApJ, 865, 20

Dorn, C., Venturini, J., Khan, A., et al. 2016, A&A, 597, A37

Doyle, A. E., Young, E. D., Klein, B., Zuckerman, B., & Schlichting, H. E. 2019, Sci, 366, 356

Dressing, C. D., Charbonneau, D., Dumusque, X., et al. 2015, ApJ, 800, 135

Duane, S., Kennedy, A., Pendleton, B. J., & Roweth, D. 1987, PhLB, 195, 216

Ducrot, E., Gillon, M., Delrez, L., et al. 2020, A&A, 640, A112

Ducrot, E., Sestovic, M., Morris, B. M., et al. 2018, AJ, 156, 218

Eastman, J., Gaudí, B. S., & Agol, E. 2013, PASP, 125, 83

Eastman, J., Siverd, R., & Gaudí, B. S. 2010, PASP, 122, 935

Elkins-Tanton, L. T., & Seager, S. 2008, ApJ, 688, 628

Elser, S., Meyer, M. R., & Moore, B. 2012, Icar, 221, 859

Fabrycky, D. C. 2010, in Exoplanets, ed. S. Seager (Tucson, AZ: Univ. Arizona Press), 217

Fauchez, T. J., Turbet, M., Villanueva, G. L., et al. 2019, ApJ, 887, 194

Fauchez, T. J., Turbet, M., Wolf, E. T., et al. 2020, GMD, 13, 707

Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10

Feng, F., Tuomi, M., Jones, H. R. A., et al. 2017, AJ, 154, 135

Fisher, R. A. 1925, Metron, 5, 90

Ford, E. B. 2006, ApJ, 642, 505

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220

Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264

Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, SSRv, 123, 485

Gillon, M. 2018, NatAs, 2, 344

Gillon, M., Jehin, E., Fumel, A., Magain, P., & Queloz, D. 2013, in EPJ Web of Conf. 47, Hot Planets and Cool Stars (Garching: EDP Sciences), 03001

Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Natur, 533, 221

Gillon, M., Jehin, E., Magain, P., et al. 2011, in EPJ Web of Conf. 11, Detection and Dynamics of Transiting Exoplanets (Garching: EDP Sciences), 06002

Gillon, M., Meadows, V., Agol, E., et al. 2020, arXiv:2002.04798

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Natur, 542, 456

Gonçalves Ferrari, G., Boekholt, T., & Portegies Zwart, S. F. 2014, MNRAS, 440, 719

Gonzales, E. C., Faherty, J. K., Gagné, J., et al. 2019, ApJ, 886, 131

Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65

Grimm, S. L., Demory, B.-O., Gillon, M., et al. 2018, A&A, 613, A68

Grimm, S. L., & Stadel, J. G. 2014, ApJ, 796, 23

Hakim, K., Rivoldini, A., Hoolst, T. V., et al. 2018, Icar, 313, 61

Hamers, A. S., & Zwart, S. F. P. 2016, MNRAS, 459, 2827
He, M. Y., Triaud, A. H. M. J., & Gillon, M. 2016, MNRAS, 464, 2687

Hernandez, D. M., & Bertschinger, E. 2015, MNRAS, 452, 1934

Hernandez, D. M., & Dehnen, W. 2017, MNRAS, 468, 2614

Hirano, T., Gaidos, E., Winn, J. N., et al. 2020, ApJL, 890, L27

Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398

Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90

Ingalls, J. G., Krick, J. E., Carey, S. J., et al. 2016, AJ, 152, 44

Ioannidis, P., Huber, K. F., & Schmitt, J. H. M. M. 2015, A&A, 585, A72

Jasche, J., & Kitaura, F. S. 2010, MNRAS, 407, 29

Jehin, E., Gillon, M., Queloz, D., et al. 2011, Msngr, 145, 2

Jehin, E., Gillon, M., Queloz, D., et al. 2018, Msngr, 174, 2

Jontof-Hutter, D., Ford, E. B., Rowe, J. F., et al. 2016, ApJ, 820, 39

Jontof-Hutter, D., Truong, V. H., Ford, E. B., Robertson, P., & Terrien, R. C. 2018, AJ, 155, 239

Juncher, D., Jørgensen, U. G., & Helling, C. 2017, A&A, 608, A70

Khan, A., Liebske, C., Rozel, A., et al. 2018, JGRE, 123, 575

Kipping, D. 2018, RNAAS, 2, 136

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

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

Kipping, D. M., Dunn, W. R., Jasinski, J. M., & Manthri, V. P. 2012, MNRAS, 421, 1166

Kite, E. S., Jr., B., F., Schaefer, L., & Ford, E. B. 2020, ApJ, 891, 111

Klein, B., & Donati, J.-F. 2019, MNRAS, 488, 5114

Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., & Megeath, S. T. 2008, ApJ, 673, 526

Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131

Krissansen-Totton, J., Garland, R., Irwin, P., & Catling, D. C. 2018, AJ, 156, 114

Leclercq, F., Pisani, A., & Wandelt, B. 2014, in Proc. Int. School of Physics, Enrico Fermi, Vol. 186, New Horizons for Observational Cosmology (Amsterdam: IOS Press), 189

Lichtenberg, T., Golabek, G. J., Burn, R., et al. 2019, NatAs, 3, 307

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2

Linial, I., Gilbaum, S., & Sari, R. 2018, ApJ, 860, 16

Lissauer, J. J. 2007, ApJL, 660, L149

Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122

Lodders, K., Palme, H., & Gail, H.-P. 2009, in Solar System, Landolt-Börnstein - Group VI Astronomy and Astrophysics, Vol. 4B, ed. J. Trümper (Berlin: Springer), 712

Luger, R., Lustig-Yaeger, J., & Agol, E. 2017a, ApJ, 851, 94

Luger, R., Sestovic, M., Kruse, E., et al. 2017b, NatAs, 1, 0129

Lustig-Yaeger, J., Meadows, V. S., & Lincowski, A. P. 2019, AJ, 158, 27

Mah, J. 2018, Master's thesis, Univ. Hong Kong, https://hub.hku.hk/handle/10722/265344

Mann, A. W., Dupuy, T., Kraus, A. L., et al. 2019, ApJ, 871, 63

Marcus, R. A., Sasselov, D., Hernquist, L., & Stewart, S. T. 2010, ApJL, 712, L73

McCullough, P. R., Crouzet, N., Deming, D., & Madhusudhan, N. 2014, ApJ, 791, 55

McDonough, W. 2014, in Treatise on Geochemistry, ed. H. D. Holland & K. K. Turekian (2nd ed.; Oxford: Elsevier), 559

Monnahan, C. C., Thorson, J. T., & Branch, T. A. 2016, Methods in Ecology and Evolution, 8, 339

Morley, C. V., Kreidberg, L., Rustamkulov, Z., Robinson, T., & Fortney, J. J. 2017, ApJ, 850, 121

Morris, B. M., Agol, E., Davenport, J. R. A., & Hawley, S. L. 2018a, ApJ, 857, 39

Morris, B. M., Agol, E., Hebb, L., et al. 2018b, ApJL, 863, L32

Murray, C. A., Delrez, L., Pedersen, P. P., et al. 2020, MNRAS, 495, 2446

Neal, R. M. 2011, in Handbook of Markov Chain Monte Carlo, ed. S. Brooks et al. (Boca Raton, FL: CRC Press), 113, https://www.mcmchandbook.net/HandbookChapter5.pdf

Öberg, K. I., & Bergin, E. A. 2016, ApJL, 831, L19

Oshagh, M., Santos, N. C., Boisse, I., et al. 2013, A&A, 556, A19

Oshagh, M., Santos, N. C., Ehrenreich, D., et al. 2014, A&A, 568, A99

Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29

Palme, H., Lodders, K., & Jones, A. 2014, in Treatise on Geochemistry, ed. H. D. Holland & K. K. Turekian (2nd ed.; Oxford: Elsevier), 15

Papaloizou, J. C. B. 2014, IJAsB, 14, 291

Pletser, V., & Basano, L. 2017, arXiv:1703.04545

Plotnykov, M., & Valencia, D. 2020, MNRAS, 499, 932

Quarles, B., Quintana, E. V., Lopez, E., Schlieder, J. E., & Barclay, T. 2017, ApJ, 842, L5

Rackham, B. V., Apai, D., & Giampapa, M. S. 2018, ApJ, 853, 122

Revels, J., Lubin, M., & Papamarkou, T. 2016, arXiv:1607.07892

Sagear, S. A., Skinner, J. N., & Muirhead, P. S. 2020, AJ, 160, 19

Scora, J., Valencia, D., Morbidelli, A., & Jacobson, S. 2020, MNRAS, 493, 4910

Seager, S., & Mallen-Ornelas, G. 2003, ApJ, 585, 1038

Sestovic, M., & Demory, B.-O. 2020, A&A, 641, A170

Siebenmorgen, R., Carraro, G., Valenti, E., et al. 2011, Msngr, 144, 9

Sotin, C., Grasset, O., & Mocquet, A. 2007, Icar, 191, 337

Spalding, C., & Adams, F. C. 2020, PSJ, 1, 7

Steele, I. A., Smith, R. J., Rees, P. C., et al. 2004, Proc. SPIE, 5489, 679

Stetson, P. B. 1987, PASP, 99, 191

Storrie-Lombardi, L. J., & Dodd, S. R. 2010, Proc. SPIE, 7737, 77370L

Tamayo, D., Rein, H., Petrovich, C., & Murray, N. 2017, ApJL, 840, L19

ter Braak, C. 2006, Statistics and Computing, 16, 239

Tinney, C. G., Ryder, S. D., Ellis, S. C., et al. 2004, Proc. SPIE, 5492, 998

Turbet, M., Bolmont, E., Ehrenreich, D., et al. 2020, A&A, 638, A41

Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86

Turbet, M., Ehrenreich, D., Lovis, C., Bolmont, E., & Fauchez, T. 2019, A&A, 628, A12

Unterborn, C. T., Desch, S. J., Hinkel, N. R., & Lorenzo, A. 2018, NatAs, 2, 297

Valencia, D., Sasselov, D. D., & O'Connell, R. J. 2007, ApJ, 665, 1413

Van Grootel, V., Fernandes, C. S., Gillon, M., et al. 2018, ApJ, 853, 30

Vida, K., Kövári, Z., Pál, A., Oláh, K., & Kriskovics, L. 2017, ApJ, 841, 124

Wakeford, H. R., Lewis, N. K., Fowler, J., et al. 2019, AJ, 157, 11

Wang, H. S., Lineweaver, C. H., & Ireland, T. R. 2019, Icar, 328, 287

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

Winn, J. N. 2010, in Exoplanets, ed. S. Seager (Tucson, AZ: Univ. Arizona Press), 55

Wisdom, J., & Hernandez, D. M. 2015, MNRAS, 453, 3015

Wisdom, J., & Holman, M. 1991, AJ, 102, 1528

Wit, E., van den Heuvel, E., & Romeijn, J.-W. 2012, Statistica Neerlandica, 66, 217

Wolf, E. T. 2017, ApJL, 839, L1

Wunderlich, F., Godolt, M., Grenfell, J. L., et al. 2019, A&A, 624, A49

Yadav, R. K., Gastine, T., Christensen, U. R., & Reiners, A. 2015, A&A, 573, A68

Zhang, Z., Zhou, Y., Rackham, B. V., & Apai, D. 2018, AJ, 156, 178