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

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. ²¹ 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. ²² 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 12). 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).

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 the United Kingdom Infra-Red 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.

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 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).

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_TDB 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₀) 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₀ 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 , σ_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 ith planet (Table 1), and j labels each transit for the ith planet, so that t_obs,ij is the jth observation of the ith planet, and σ_ij is the corresponding measurement error. The total number of transits is ${N}_{\mathrm{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.

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). ²³ 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 ω_i , e_i sin ω_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₀ = M_*/M_⊙, 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. ²⁴ The initial orbital ephemeris, (P_i , 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 eccentric anomaly for each initial Keplerian. Finally, the eccentricity, e_i , and longitude of periastron, ω_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 ω_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 jth transit of the ith planet, giving each model transit time as a function of the initial conditions, t_ij ( 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

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}{{ \mathcal L }}_{{ij}}({{\boldsymbol{x}}}_{\mathrm{dyn}},\nu ,{V}_{1})\\ \ =\ -\tfrac{\nu +1}{2}\mathrm{log}\left(1+\displaystyle \frac{{\left({t}_{\mathrm{obs},{ij}}-{t}_{{ij}}({{\boldsymbol{x}}}_{\mathrm{dyn}})\right)}^{2}}{\nu {\sigma }_{{ij}}^{2}{V}_{1}}\right)\\ \ \ \ \ -\ \tfrac{1}{2}\mathrm{log}\left(\pi \nu {V}_{1}\right)+\mathrm{log}{\rm{\Gamma }}\ \left(\displaystyle \frac{\nu +1}{2}\right)-\mathrm{log}{\rm{\Gamma }}\ \left(\displaystyle \frac{\nu }{2}\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where Γ(x) is the Gamma function (Fisher 1925).

The total log likelihood function which we optimize is given by

$$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }({{\boldsymbol{x}}}_{\mathrm{dyn}},\nu ,{V}_{1})=\displaystyle \sum _{i=1}^{{N}_{p}}\displaystyle \sum _{j=1}^{{N}_{i}}\mathrm{log}{{ \mathcal L }}_{{ij}}({{\boldsymbol{x}}}_{\mathrm{dyn}},\nu ,{V}_{1}),\end{eqnarray} \tag{ 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.

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 }_{{x}_{i}}$, where ${\sigma }_{{x}_{i}}$ equals the square root of the diagonal component for the ith 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. ²⁵ 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, ν and V₁, were strongly nonlinearly correlated and displayed a likelihood profile which was non-Gaussian. After experimenting with reparameterization, we found that $\mathrm{log}\nu$ and V₁ e^1/(2ν) 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 ${\boldsymbol{x}}=({{\boldsymbol{x}}}_{\mathrm{dyn}},\mathrm{log}\nu ,{V}_{1}{e}^{1/2\nu })$.

In Appendix B, we define the prior function Π( x ), which multiplies the likelihood to give the posterior probability distribution,

$$\begin{eqnarray}&&P({\boldsymbol{x}})\propto {\rm{\Pi }}({\boldsymbol{x}})\times { \mathcal L }({\boldsymbol{x}}),\end{eqnarray} \tag{ 4 }$$

so that we can proceed to discussing the Markov chain sampling of the posterior probability of the model parameters given the data.

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). ²⁶ 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,

$$\begin{eqnarray}&&H({\boldsymbol{p}},{\boldsymbol{x}})=\tfrac{1}{2}{{\boldsymbol{p}}}^{T}{{\boldsymbol{M}}}^{-1}{\boldsymbol{p}}-\mathrm{log}{ \mathcal L }({\boldsymbol{x}})-\mathrm{log}{\rm{\Pi }}({\boldsymbol{x}}),\end{eqnarray} \tag{ 5 }$$

where  p is defined from Hamilton's equations,

$$\begin{eqnarray}&&\dot{{\boldsymbol{p}}}=-\displaystyle \frac{\partial H}{\partial {\boldsymbol{x}}}.\end{eqnarray} \tag{ 6 }$$

We take the mass matrix,  M , to be the approximate Hessian matrix evaluated at the maximum likelihood, ${\boldsymbol{M}}={\boldsymbol{ \mathcal H }}({{\boldsymbol{x}}}_{0})$ (Equation (A5)). Similarly, the Hamiltonian can be used to compute the evolution of the parameter "coordinates,"

$$\begin{eqnarray}&&\dot{{\boldsymbol{x}}}=+\displaystyle \frac{\partial H}{\partial {\boldsymbol{p}}}.\end{eqnarray} \tag{ 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,

$$\begin{eqnarray}&&{\boldsymbol{p}}={{\boldsymbol{M}}}^{1/2}{\boldsymbol{Z}},\end{eqnarray} \tag{ 8 }$$

where Z_n  ∼ N(0,1) is an element of a vector of random normal deviates for n = 1, ..., N_param. We then carry out a leapfrog integration of Hamilton's equations for N_leap steps from the starting point with a "time" step to obtain a proposal set of parameters ( x _prop,  p _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

$$\begin{eqnarray}&&{p}_{\mathrm{accept}}=\min (\exp (-(H({{\boldsymbol{x}}}_{\mathrm{prop}},{{\boldsymbol{p}}}_{\mathrm{prop}})-H({\boldsymbol{x}},{\boldsymbol{p}}))),1)\end{eqnarray} \tag{ 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: ₀ and N_leap,0. We draw the "time" step, , for each integration from the absolute value of a normal distribution with width ₀, i.e.,  ∼ ∣N(0, ₀)∣. We draw the number of leapfrog steps for each integration from a uniform probability, ${N}_{\mathrm{leap}}\sim \mathrm{round}({N}_{\mathrm{leap},0}{ \mathcal U }(0.8,1.0))$. We found that a choice of ₀ = 0.1 and N_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. ²⁷ We found a minimum mean effective sample size of 57 over all chains for a total number of independent samples of 6409.

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_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 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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, $\mathrm{log}\nu =1.3609\pm 0.233\,7$ and V₁ e^1/(2ν) = 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 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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Parameters of the TRAPPIST-1 System from Transit-timing Analysis and Their 1σ Uncertainties

	$\mu \left[\tfrac{{{\boldsymbol{M}}}_{\oplus }}{0.09{{\boldsymbol{M}}}_{\odot }}\right]$	Mp	$\tfrac{{\sigma }_{\mu }}{\mu }$	P	t0	$e\cos \omega$	$e\sin \omega$
	$=\tfrac{{{\boldsymbol{M}}}_{p}}{{{\boldsymbol{M}}}_{\oplus }}\left(\tfrac{0.09{{\boldsymbol{M}}}_{\odot }}{{{\boldsymbol{M}}}_{* }}\right)$	[10−5 M*]	%	(day)	(BJDTDB-2,450,000)
b	1.377 1 ± 0.0593	4.596 ± 0.198	4.3	1.510 826 ± 0.000006	7 257.550 44 ± 0.00015	−0.00215 ± 0.00332	0.002 17 ± 0.00244
c	1.310 5 ± 0.0453	4.374 ± 0.151	3.5	2.421 937 ± 0.000018	7 258.587 28 ± 0.00027	0.000 55 ± 0.00232	0.000 01 ± 0.00171
d	0.388 5 ± 0.0074	1.297 ± 0.025	1.9	4.049 219 ± 0.000026	7 257.067 68 ± 0.00067	−0.00496 ± 0.00186	0.002 67 ± 0.00112
e	0.693 2 ± 0.0128	2.313 ± 0.043	1.8	6.101 013 ± 0.000035	7 257.827 71 ± 0.00041	0.004 33 ± 0.00149	−0.00461 ± 0.00087
f	1.041 1 ± 0.0155	3.475 ± 0.052	1.5	9.207 540 ± 0.000032	7 257.074 26 ± 0.00085	−0.00840 ± 0.00130	−0.00051 ± 0.00087
g	1.323 8 ± 0.0171	4.418 ± 0.057	1.3	12.352 446 ± 0.000054	7 257.714 62 ± 0.00103	0.003 80 ± 0.00112	0.001 28 ± 0.00070
h	0.326 1 ± 0.0186	1.088 ± 0.062	5.7	18.772 866 ± 0.000214	7 249.606 76 ± 0.00272	−0.00365 ± 0.00077	−0.00002 ± 0.00044

Note. Note that the mass ratios, $\mu ={{\boldsymbol{M}}}_{{\boldsymbol{p}}}/{{\boldsymbol{M}}}_{* }$, of the planets are computed relative to the star, which is assumed to have a mass of 0.09 M_⊙ (this is later combined with the estimate of stellar mass to give our estimates of the planet masses). We also report μ in units of 10⁻⁵, and the fractional precision on the measurement of μ, σ_μ /μ. The parameters P, t₀, $e\cos \omega$, and $e\sin \omega$ describe the osculating Jacobi elements at the start of the simulation, on date BJD_TDB-2,450,000 = 7257.93115525 days.

Download table as:  ASCIITypeset image

The marginalized posterior distributions of the ratio of the planet masses to the star masses, scaled to a stellar mass of 0.09 M_⊙, 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.

Zoom In Zoom Out Reset image size

Table 3. Planet-to-star Mass Ratios in Units of M_⊕/(0.09M_⊙) 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

Source	Quantity	b	c	d	e	f	g	h
Grimm	$\tfrac{{M}_{p}}{{M}_{\oplus }}\tfrac{0.09{M}_{\odot }}{{M}_{* }}$	${1.017}_{-0.143}^{+0.154}$	${1.156}_{-0.131}^{+0.142}$	${0.297}_{-0.035}^{+0.039}$	${0.772}_{-0.075}^{+0.079}$	${0.934}_{-0.078}^{+0.080}$	${1.148}_{-0.095}^{+0.098}$	${0.331}_{-0.049}^{+0.056}$
This paper	$\tfrac{{M}_{p}}{{M}_{\oplus }}\tfrac{0.09{M}_{\odot }}{{M}_{* }}$	1.3771 ± 0.0593	1.3105 ± 0.0453	0.3885 ± 0.0074	0.693 2 ± 0.0128	1.0411 ± 0.0155	1.3238 ± 0.0171	0.3261 ± 0.0186
Delrez	Rp /R*	0.0853 ± 0.0004	0.0833 ± 0.0004	0.0597 ± 0.0006	0.0693 ± 0.0007	0.0796 ± 0.0006	0.0874 ± 0.0006	0.0588 ± 0.0012
This paper	Rp /R*	0.0859 ± 0.0004	0.0844 ± 0.0004	0.0606 ± 0.0005	0.0708 ± 0.0006	0.0804 ± 0.0005	0.0869 ± 0.0005	0.0581 ± 0.0009

Download table as:  ASCIITypeset image

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 ν = 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 ν = 4, which was the value used by Jontof-Hutter et al. (2016).

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  x  = ({m_i , P_i , t_0,i , e_i , ω_i }; i = 1, ..., N_p ). The mass of the star is taken to be M_⋆ = 0.09M_⊙, 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.

To compute the uncertainties on the photodynamical model parameters, we use an affine-invariant MCMC algorithm (Goodman & Weare 2010). ²⁸ 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. ²⁹

Table 4. Prior Bounds on Photodynamic Parameters

Parameter	Units	Prior
b0	R*	${ \mathcal U }(0,1)$
Rp /R*	⋯	${ \mathcal U }(0,0.2)$
ρ*	ρ⊙	${ \mathcal U }(0,100)$
(q1,Ch 1, q2,Ch 1)	⋯	${ \mathcal U }(0,1)$
(q1,Ch 2, q2,Ch 2)	⋯	${ \mathcal U }(0,1)$

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

Download table as:  ASCIITypeset image

Table 5. Parameters Derived from the Photodynamic Model

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

Note. Top: stellar density (in units of solar density), limb-darkening parameters (q₁, q₂) in Spitzer channels 1 and 2, and stellar density in cgs units and limb-darkening parameters u₁ and u₂. Bottom: planet-to-star radius ratio, R_p /R_*; transit depth, (R_p /R_*)²; transit duration, T (from first to fourth contact); ingress/egress duration, τ (from first to second contact or third to fourth contact); impact parameter in units of stellar radius, b₀ (assumed to be positive); ratio of semimajor axis to stellar radius, a/R_*; and inclination I in degrees (for b₀ > 0).

Download table as:  ASCIITypeset image

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).

Zoom In Zoom Out Reset image size

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σ–2σ. Figure 8 shows the posterior probability distribution of impact parameters in units of the stellar radius, b₀, 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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,

$$\begin{eqnarray}&&\displaystyle \frac{a}{{R}_{* }}={\left(\displaystyle \frac{{P}^{2}{{GM}}_{\odot }}{4{\pi }^{2}{R}_{\odot }^{3}}\displaystyle \frac{{\rho }_{* }}{{\rho }_{\odot }}\right)}^{1/3}.\end{eqnarray} \tag{ 11 }$$

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

$$\begin{eqnarray}&&I={\cos }^{-1}\left[{b}_{0}{\left(\displaystyle \frac{a}{{R}_{* }}\right)}^{-1}\right],\end{eqnarray} \tag{ 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.

Zoom In Zoom Out Reset image size

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 01, 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 σ_θ , 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 σ_θ are preferred, with ${\sigma }_{\theta }={0\buildrel{\circ}\over{.} 041}_{-0\buildrel{\circ}\over{.} \,016}^{+0\buildrel{\circ}\over{.} \,031}$. 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 025.

Zoom In Zoom Out Reset image size

The inclination prior also enables a more precise and symmetric estimate of the density of the star, ρ_*/ρ_⊙ = 53.22 ± 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.

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σ 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σ 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⁴ 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₃) 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σ 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 period, P, for 10⁴ posterior samples and found a relation of y = (1.042 ± 0.034) − (0.0043 ± 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⁴ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The interpretation of these observations in terms of internal compositions is discussed in more detail next.

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₃) 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₃) 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 ± 4 wt% (taking into account the correlations between the planets' CMFs).

Zoom In Zoom Out Reset image size

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

Planet:	b	c	d	e	f	g	h	Avg b–h
CMF (wt%)	${25.2}_{-6.0}^{+5.3}$	${26.6}_{-5.1}^{+4.6}$	${19.7}_{-5.1}^{+4.7}$	${24.6}_{-4.9}^{+4.3}$	${20.1}_{-4.2}^{+3.5}$	${16.1}_{-4.2}^{+3.5}$	${16.5}_{-10.0}^{+9.3}$	20.9 ± 3.6
Fe/Mg molar ratio	${0.60}_{-0.18}^{+0.18}$	${0.64}_{-0.16}^{+0.16}$	${0.44}_{-0.13}^{+0.14}$	${0.58}_{-0.14}^{+0.14}$	${0.45}_{-0.11}^{+0.10}$	${0.34}_{-0.10}^{+0.09}$	${0.35}_{-0.23}^{+0.27}$	0.47 ± 0.07
H2O (wt%) for:
CMF = 18%	<10−3	<10−3	<10−3	${0.0}_{-0.0}^{+0.0}$	${0.0}_{-0.0}^{+0.0}$	${0.72}_{-0.72}^{+1.3}$	${0.6}_{-0.6}^{+3.4}$
CMF = 25%	<10−3	<10−3	<10−3	${0.3}_{-0.3}^{+1.8}$	${1.9}_{-1.3}^{+1.5}$	${3.5}_{-1.3}^{+1.6}$	${3.0}_{-3.0}^{+3.8}$
CMF = 32.5%	<10−3	<10−3	<10−3	${2.9}_{-1.5}^{+1.7}$	${4.5}_{-1.2}^{+1.8}$	${6.4}_{-1.6}^{+2.0}$	${5.5}_{-3.1}^{+4.5}$
CMF = 50%	${0.05}_{-0.03}^{+0.08}$	${0.03}_{-0.02}^{+0.05}$	${0.002}_{-0.0009}^{+0.002}$	${9.4}_{-1.8}^{+2.2}$	${12}_{-1.7}^{+2.0}$	${14}_{-1.7}^{+2.0}$	${12}_{-3.9}^{+4.4}$

Download table as:  ASCIITypeset image

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 ³¹ 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₂O-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.

Zoom In Zoom Out Reset image size

Higher estimated water budgets for the outer three or four planets could be a clue that they formed beyond the water condensation line at ≈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

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τ ≈ 10 minutes). We estimate that a flare of energy 10³¹ erg could cause a 1.5σ timing outlier if it occurs during ingress or egress. This has a probability of only ≈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.

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.007 1 ± 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.

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). ³² 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).

Zoom In Zoom Out Reset image size

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?

Zoom In Zoom Out Reset image size

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 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.

In this section, we discuss some of the dynamical aspects of the planetary system: the eccentricities, the longitudes of periastron, and the GLR angles.

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, ≲10⁻³ (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).

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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 ≈ 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.

Zoom In Zoom Out Reset image size

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).

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{\phi }_{i,i+1,i+2}=p{\lambda }_{i}-(p+q){\lambda }_{i+1}+q{\lambda }_{i+2},\end{eqnarray} \tag{ 15 }$$

where λ_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 ϕ = 180°, 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, ϕ and 360° − ϕ. 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.

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

$$\begin{eqnarray}\begin{array}{rcl}{\phi }_{\mathrm{bcd}} & = & 2{\lambda }_{{\rm{b}}}-5{\lambda }_{{\rm{c}}}+3{\lambda }_{{\rm{d}}},\\ {\phi }_{\mathrm{cde}} & = & {\lambda }_{{\rm{c}}}-3{\lambda }_{{\rm{d}}}+2{\lambda }_{{\rm{e}}},\\ {\phi }_{\mathrm{def}} & = & 2{\lambda }_{{\rm{d}}}-5{\lambda }_{{\rm{e}}}+3{\lambda }_{{\rm{f}}},\\ {\phi }_{\mathrm{efg}} & = & {\lambda }_{{\rm{e}}}-3{\lambda }_{{\rm{f}}}+2{\lambda }_{{\rm{g}}},\\ {\phi }_{\mathrm{fgh}} & = & {\lambda }_{{\rm{f}}}-2{\lambda }_{{\rm{g}}}+{\lambda }_{{\rm{h}}}.\end{array}\end{eqnarray} \tag{ 16 }$$

Differences between the predicted and observed angles agree within 05–10°, where the predicted values for ϕ are taken from Mah (2018), but allowing ϕ_bcd and ϕ_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.

Zoom In Zoom Out Reset image size

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⁴ posterior samples from the timing analysis. These simulations were carried out for 10⁷ 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 × 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 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, ≤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).

Table 11. Evolution of the Semimajor Axes, a, and Eccentricities, e, from 10⁴ Samples over 10 Myr

Planet	${\rm{\Delta }}\bar{a}$	${\rm{\Delta }}\bar{e}$	$\max ({\rm{\Delta }}a)$	$\max ({\rm{\Delta }}e)$
b	−6.52e−09	1.73e−04	7.12e−07	0.0020
c	−1.44e−08	1.62e−04	2.51e−06	0.0018
d	−1.06e−08	1.45e−05	4.07e−06	0.0009
e	2.05e−08	4.44e−05	8.53e−06	0.0008
f	2.45e−07	5.13e−05	3.00e−05	0.0012
g	8.24e−08	5.01e−05	2.19e−05	0.0011
h	−1.23e−07	2.11e−04	3.00e−05	0.0035

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

Download table as:  ASCIITypeset image

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.

Zoom In Zoom Out Reset image size

Table 12. Number of Posterior Samples Falling into the Three Resonant Categories for the Five GLR Angles

	bcd	cde	def	efg	fgh
Category I	130	178	755	7943	6462
Category II	1571	874	7653	1855	1578
Category III	8299	8948	1592	202	1960

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

Download table as:  ASCIITypeset image

10.6.1. Forecast Transit Times

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). ³³ 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

The stellar density we derive using the photodynamic model, ${\rho }_{* }={53.17}_{-1.18}^{+0.72}{\rho }_{\odot }$, is in 1σ agreement with prior analyses. Most recently, Delrez et al. (2018b) found a density of ρ_* = (52.3 ± 2.2)ρ_⊙, 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 ± 0.0023M_⊙. ³⁴ This mass has a 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}^{2}{T}_{\mathrm{eff}}^{4}$ at 1σ, while our T_eff is more precise (28 K versus 42 K), R is 2.5 times more precise, and our $\mathrm{log}g$ is more precise by an order of magnitude.

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 protoplanetary 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 & Seager 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, ³⁵ 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.

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 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.