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−1, 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.
Export citation and abstract BibTeX RIS
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⊙; 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 (∼0.12R⊙) 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) 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 for integers j and k for the ith and (i + 1)th planets. This proximity causes a resonant timescale for k = 1 given by
(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 PTTV ≈ 491 ± 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 tmin ≈ 2PTTV ≈ 2.7 yr. The current paper is the first with a survey time, tsurvey = 4.114 yr, such that tsurvey > tmin for the TRAPPIST-1 system.
Prior studies used the data available at the time (Delrez et al. 2018b), with tsurvey < tmin, 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σ–2σ 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 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 (BJDTDB), 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. 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. 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 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).
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 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.
Table 1. Number of Transits from Ground-based and Space-based Observations
Planet | HCT | SSO/TS/TN | LT | WHT | VLT/AAT/UKIRT | HST | Spitzer | K2 | Duplicates | Total (Ni ) |
---|---|---|---|---|---|---|---|---|---|---|
b | 1 | 45 | 7 | 1 | 10 | 1 | 64 | 48 | 17 | 160 |
c | 0 | 28 | 8 | 0 | 7 | 1 | 47 | 30 | 14 | 107 |
d | 0 | 11 | 1 | 1 | 5 | 2 | 23 | 17 | 7 | 53 |
e | 0 | 18 | 4 | 0 | 3 | 2 | 18 | 11 | 7 | 49 |
f | 0 | 9 | 2 | 0 | 4 | 2 | 16 | 7 | 6 | 34 |
g | 0 | 11 | 0 | 0 | 3 | 2 | 13 | 5 | 4 | 30 |
h | 0 | 3 | 2 | 0 | 0 | 0 | 7 | 4 | 2 | 14 |
Total | 1 | 125 | 24 | 2 | 32 | 10 | 188 | 122 | 57 | 447 |
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).
Download table as: ASCIITypeset image
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 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 BJDTDB 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 (T0) 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 T0 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 = ({tobs,ij , σij ; j = 1, ..., Ni }; i = 1, ...,7), where i labels each of the seven planets, Ni is the number of transits for the ith planet (Table 1), and j labels each transit for the ith planet, so that tobs,ij is the jth observation of the ith planet, and σij is the corresponding measurement error. The total number of transits is = 447, where Np 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 χ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) ≲ 0.1 and P(>z) ≳ 0.9, where z = (tobs,ij − tij ( x dyn))/σij are the normalized residuals, with the model time, tij ( 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 ≈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).
Our transit-timing model will be affected by these timing outliers, which make an excessive contribution to the χ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 V1. For the Student's t-distribution, there is only one additional free parameter: the number of degrees of freedom, ν, 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). 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 = ({mi , Pi , t0 ,i , ei cos ωi , ei sin ωi }; i = 1, ..., Np ), where Np is the number of planets for a total of 5Np dynamical parameters. In addition, we take the star to have a mass, m0 = 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 au. 24 The initial orbital ephemeris, (Pi , t0,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, ei , and longitude of periastron, ωi , for each Keplerian, we parameterize in terms of and 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, tij ( x dyn), which can then be compared to the observed times, tobs,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
where Γ(x) is the Gamma function (Fisher 1925).
The total log likelihood function which we optimize is given by
where Np is the number of planets; we use Np = 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 , where 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. 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, ν and V1, were strongly nonlinearly correlated and displayed a likelihood profile which was non-Gaussian. After experimenting with reparameterization, we found that and V1 e1/(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 .
In Appendix B, we define the prior function Π( x ), which multiplies the likelihood to give the posterior probability distribution,
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). 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,
where p is defined from Hamilton's equations,
We take the mass matrix, M , to be the approximate Hessian matrix evaluated at the maximum likelihood, (Equation (A5)). Similarly, the Hamiltonian can be used to compute the evolution of the parameter "coordinates,"
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,
where Zn ∼ N(0,1) is an element of a vector of random normal deviates for n = 1, ..., Nparam. We then carry out a leapfrog integration of Hamilton's equations for Nleap 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
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: 0 and Nleap,0. We draw the "time" step, , for each integration from the absolute value of a normal distribution with width 0, i.e., ∼ ∣N(0, 0)∣. We draw the number of leapfrog steps for each integration from a uniform probability, . We found that a choice of 0 = 0.1 and Nleap,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. 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 PTTV. 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.
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageThe 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, and V1 e1/(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 and as Markov chain parameters that favor higher eccentricities (Ford 2006).
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageTable 2. Parameters of the TRAPPIST-1 System from Transit-timing Analysis and Their 1σ Uncertainties
Mp | P | t0 | |||||
---|---|---|---|---|---|---|---|
[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, , 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−5, and the fractional precision on the measurement of μ, σμ /μ. The parameters P, t0, , and describe the osculating Jacobi elements at the start of the simulation, on date BJDTDB-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.
Download figure:
Standard image High-resolution imageTable 3. Planet-to-star Mass Ratios in Units of M⊕/(0.09M⊙) from Grimm et al. (2018) and Planet-to-star Radius Ratios Rp /R* from Delrez et al. (2018b) Compared with the Results from This Paper
Source | Quantity | b | c | d | e | f | g | h |
Grimm | ||||||||
This paper | 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 , 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).
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 x = ({mi , Pi , t0,i , ei , ωi }; i = 1, ..., Np ). 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.
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 & Mallen-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 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, t0, 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, vsky, in units of R*/day. We account for quadratic limb darkening of the star with parameters (q1,Ch1, q2,Ch1, q1,Ch2, q2,Ch2) in the two Spitzer channels, and for each planet, we specify a planet-to-star radius ratio (Rp /R*) and we assume a midtransit impact parameter (b0), 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
in units of the stellar radius, R*. 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−7 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 via regression at each step in the Markov chain. We transform the q1, q2 limb-darkening parameters to u1, u2 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 χ2 of unity in that window. With this photometric scatter, we compute a χ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). 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. 29
Table 4. Prior Bounds on Photodynamic Parameters
Parameter | Units | Prior |
---|---|---|
b0 | R* | |
Rp /R* | ⋯ | |
ρ* | ρ⊙ | |
(q1,Ch 1, q2,Ch 1) | ⋯ | |
(q1,Ch 2, q2,Ch 2) | ⋯ |
Note. Note that the same bounds on the impact parameter, b0, and radius ratio, Rp /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: | 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: | 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* | |||||||
a/R* | |||||||
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 (q1, q2) in Spitzer channels 1 and 2, and stellar density in cgs units and limb-darkening parameters u1 and u2. Bottom: planet-to-star radius ratio, Rp /R*; transit depth, (Rp /R*)2; 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, b0 (assumed to be positive); ratio of semimajor axis to stellar radius, a/R*; and inclination I in degrees (for b0 > 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).
Download figure:
Standard image High-resolution imageTable 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, b0, 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.
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageCombining 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,
With the measured impact parameters, we compute the inclinations of the planets from Winn (2010),
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.
Download figure:
Standard image High-resolution image5.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 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 . 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.
Download figure:
Standard image High-resolution imageThe 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.
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 calibrate the mass–luminosity (M* − ) 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 ± 0.0023M⊙.
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 ± 0.0023M⊙. 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
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σ and 2σ 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.
Download figure:
Standard image High-resolution imageIn 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 (ρ♂ = 3.9335 g cm−3 = 0.713 ρ⊕) and Earth (ρ⊕ = 5.514 g cm−3). The surface gravities span a range from 57% of Earth (planet h) to 110% of Earth (planet b).
Table 6. Planetary Parameters from Combining the Transit-timing and Photodynamic Analysis
Planet: | b | c | d | e | f | g | h |
---|---|---|---|---|---|---|---|
R (R⊕) | |||||||
M (M⊕ ) | 1.374 ± 0.069 | 1.308 ± 0.056 | 0.388 ± 0.012 | 0.692 ± 0.022 | 1.039 ± 0.031 | 1.321 ± 0.038 | 0.326 ± 0.020 |
ρ (ρ⊕ ) | |||||||
g (g⊕ ) | 1.102 ± 0.052 | 1.086 ± 0.043 | 0.624 ± 0.019 | 0.817 ± 0.024 | 0.951 ± 0.024 | 1.035 ± 0.026 | 0.570 ± 0.038 |
vesc (vesc,⊕) | 1.109 ± 0.026 | 1.092 ± 0.022 | 0.701 ± 0.010 | 0.867 ± 0.012 | 0.997 ± 0.012 | 1.081 ± 0.013 | 0.656 ± 0.020 |
S (S⊕) | |||||||
a (10−2 au) | 1.154 ± 0.010 | 1.580 ± 0.013 | 2.227 ± 0.019 | 2.925 ± 0.025 | 3.849 ± 0.033 | 4.683 ± 0.040 | 6.189 ± 0.053 |
R (108 cm) | |||||||
M (1027 g) | 8.211 ± 0.412 | 7.814 ± 0.335 | 2.316 ± 0.074 | 4.132 ± 0.130 | 6.205 ± 0.184 | 7.890 ± 0.226 | 1.945 ± 0.122 |
ρ (g cm−3) | |||||||
g (10 m s−2) | 1.080 ± 0.051 | 1.065 ± 0.042 | 0.611 ± 0.019 | 0.801 ± 0.024 | 0.932 ± 0.024 | 1.015 ± 0.025 | 0.558 ± 0.037 |
vesc (km s−1) | 12.400 ± 0.292 | 12.205 ± 0.241 | 7.839 ± 0.110 | 9.694 ± 0.133 | 11.145 ± 0.137 | 12.087 ± 0.142 | 7.335 ± 0.227 |
S (106erg cm−2 s−1) | |||||||
a (1011 cm) | 1.726 ± 0.015 | 2.364 ± 0.020 | 3.331 ± 0.028 | 4.376 ± 0.037 | 5.758 ± 0.049 | 7.006 ± 0.060 | 9.259 ± 0.079 |
Note. The units are given with respect to Earth first and cgs second.
Download table as: ASCIITypeset image
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 ≈0.09M⊙ 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 ± 0.19) × 10−4 L⊙, which, when compared with stellar evolution models, yields a mass of M = 0.09016 ± 0.0010M⊙, which is also consistent with the Mann et al. (2019) value. Burgasser & Mamajek (2017) found an older age for the host star, 7.6 ± 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.
Download figure:
Standard image High-resolution imageBased 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.
Table 7. Updated Stellar Parameters Based on the Combined Analysis
Parameter | Value | References |
---|---|---|
M(M⊙) | 0.0898 ± 0.0023 | Mann et al. (2019) |
R(R⊙) | 0.1192 ± 0.0013 | This paper |
L(L⊙) | 0.000553 ± 0.000019 | Ducrot et al. (2020) |
Teff (K) | 2566 ± 26 | This paper |
This paper |
Download table as: ASCIITypeset image
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 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).
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," Pi, assuming a GLR configuration with planets g and h given by
for a range of 1 ≤ p, q ≤ 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.
Table 8. Predictions for a GLR of Planets g and h with an Eighth Planet, Planet i, with Period Pi
p | q | Pi (day) | Pi/ Ph | j |
---|---|---|---|---|
1 | 1 | 39.029 | 2.08 | 1 |
1 | 2 | 25.347 | 1.35 | 3 |
1 | 3 | 22.695 | 1.21 | 4 |
2 | 3 | 28.701 | 1.53 | 2 |
Note. The ratio with the period of planet h is given, as well as the value of j for which Pi/Ph ≈ (j + 1)/j.
Download table as: ASCIITypeset image
We carried out a transit-timing search for an eighth planet by placing planets with mass ratios between 2 × 10−6 and 5 × 10−5 at these four trial orbital periods in a coplanar 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 Ntrans = 30.5. Given that the inner seven planets show orbital eccentricities with values ≲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.09M⊙. These two cases with Δ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).
Download figure:
Standard image High-resolution imageWe 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σ 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 104 samples and divided by the density of the closest pair of mass and radius of a fully differentiated 20 wt% iron, 80 wt% silicate (MgSiO3) 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 104 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 104 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.
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageThe interpretation of these observations in terms of internal compositions is discussed in more detail next.
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 (MgSiO3) 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 (MgSiO3) 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 wt% for planet g up to 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).
Download figure:
Standard image High-resolution imageTable 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%) | 20.9 ± 3.6 | |||||||
Fe/Mg molar ratio | 0.47 ± 0.07 | |||||||
H2O (wt%) for: | ||||||||
CMF = 18% | <10−3 | <10−3 | <10−3 | |||||
CMF = 25% | <10−3 | <10−3 | <10−3 | |||||
CMF = 32.5% | <10−3 | <10−3 | <10−3 | |||||
CMF = 50% |
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 104 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 31 water layer to contribute to the total radius, we can estimate the water mass fractions of the seven planets (b: wt%, c: wt%, d: wt%, e: wt%, f: wt%, g: wt%, h: 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 H2O-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.
Download figure:
Standard image High-resolution imageHigher 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
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 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τ ≈ 10 minutes). We estimate that a flare of energy 1031 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.
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.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.
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). 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).
Download figure:
Standard image High-resolution image10.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−1, 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?
Download figure:
Standard image High-resolution imageTable 10. RV Semiamplitudes, Kp , for the TRAPPIST-1 Planets Predicted from Our Measured Masses
Planet | b | c | d | e | f | g | h |
---|---|---|---|---|---|---|---|
Kp (cm s−1) | 382.0 | 310.7 | 77.6 | 120.7 | 158.1 | 182.2 | 39.1 |
RV equivalent precision for TRAPPIST-1 host (cm s−1) | 19 | 13 | 2.5 | 3.8 | 4.7 | 5.2 | 2.4 |
RV equivalent precision for 1 M⊙ host at 1 au (cm s−1) | 0.62 | 0.50 | 0.11 | 0.20 | 0.28 | 0.34 | 0.18 |
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.
Download table as: ASCIITypeset image
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−1, 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−1, 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−1. 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−1 (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.
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−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).
Download figure:
Standard image High-resolution imageFigure 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.
Download figure:
Standard image High-resolution imageNow, 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 PTTV ≈ 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.
Download figure:
Standard image High-resolution imageThe 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).
Download figure:
Standard image High-resolution image10.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
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:
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.
Download figure:
Standard image High-resolution image10.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 104 posterior samples from the timing analysis. These simulations were carried out for 107 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 × 1010 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 104 Samples over 10 Myr
Planet | ||||
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 and 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.
Download figure:
Standard image High-resolution imageTable 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. 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 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). 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.
Download figure:
Standard image High-resolution imageFinally, 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.
Download figure:
Standard image High-resolution imageOf 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, , 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⊙. 34 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 at 1σ, while our Teff is more precise (28 K versus 42 K), R is 2.5 times more precise, and our 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 ± 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 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, 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 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 Sklodowska-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
where all of the dependence on the dynamical model parameters enters through tij ( x 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:
where we have used the fact that the gradient of the log likelihood vanishes at the maximum-likelihood value of the model parameters, x 0, and the indices k, l = 1, ..., 5Np + 2 for xk and xl , where the first 5Np parameters are the dynamical parameters, x dyn, and the last two parameters are the Student's t-distribution likelihood parameters, and V1 e1/2ν . Now, the width of the Gaussian distribution at the maximum likelihood is governed by the Hessian matrix, with elements given by
which involves second derivatives of the negative log likelihood with respect to the model parameters. The derivatives of tij with respect to 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 ≤ k, l ≤ 5Np , we can write
where tij = tij ( 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, , we compute the Hessian as
where the partial derivatives with respect to tij ( x dyn), , and 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, , 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, , the volume of parameter space scales ∝ ei , and so a 1/ei prior is needed to yield a posterior that has a uniform probability with eccentricity, ei , 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 ξ, we specify
so that the total prior is given by
where the values of ξ1 − ξ4 and each transformation of parameters, f = {fj ( x ); j = 1, ..., Nbound} are given in Table 13, where Nbound = 4Np + 2. The prior probability, then, is given by Π( x ), which we multiply by the likelihood function before sampling.
Table 13. Prior Probability Boundary Limits for the TRAPPIST-1 Planet Parameters
Parameter | Bounds Function fj | Lower Bound | Lower Transition | Upper Transition | Upper Bound |
---|---|---|---|---|---|
ξ1 | ξ2 | ξ3 | ξ4 | ||
Mass ratio | −8 | −7 | −3 | −2 | |
Eccentricity | e | ⋯ | ⋯ | 0.2 | 0.3 |
Period of b | Pb (days) | 1.49 | 1.50 | 1.52 | 1.53 |
Period of c | Pc (days) | 2.40 | 2.41 | 2.43 | 2.44 |
Period of d | Pd (days) | 4.03 | 4.04 | 4.06 | 4.07 |
Period of e | Pe (days) | 6.08 | 6.09 | 6.11 | 6.12 |
Period of f | Pf (days) | 9.18 | 9.19 | 9.22 | 9.23 |
Period of g | Pg (days) | 12.33 | 12.34 | 12.36 | 12.37 |
Period of h | Ph (days) | 18.75 | 18.76 | 18.78 | 18.79 |
Initial transit time of b | t0,b − 2,457,257 (days) | 0.53 | 0.54 | 0.56 | 0.57 |
Initial transit time c | t0,c − 2,457,257 (days) | 0.57 | 0.58 | 0.60 | 0.61 |
Initial transit time d | t0,d − 2,457,257 (days) | 0.05 | 0.06 | 0.08 | 0.09 |
Initial transit time e | t0,e − 2,457,257 (days) | 0.81 | 0.82 | 0.84 | 0.85 |
Initial transit time f | t0,f − 2,457,257 (days) | 0.05 | 0.06 | 0.08 | 0.09 |
Initial transit time g | t0,g − 2,457,257 (days) | 0.70 | 0.71 | 0.73 | 0.74 |
Initial transit time h | t0,h − 2,457,249 (days) | 0.58 | 0.59 | 0.61 | 0.62 |
Degrees of freedom | ν | 0.5 | 1.0 | 50 | 100 |
Log variance factor | −2 | −1 | 5 | 10 |
Note. The bounds are chosen so as to not affect the parameters as much as possible.
Download table as: ASCIITypeset image
Appendix C: Corner Plots
Figures 29 and 30 show corner plots of the variables from the transit-timing and photodynamical analyses, respectively.
Download figure:
Standard image High-resolution imageDownload figure:
Standard image High-resolution imageAppendix D: Tables
Tables of the best-fit transit times (Table 14) and the forecast times (Table 15) are given in this appendix.
Table 14. Observed Transit Times with Uncertainties, along with the Mean, tpost, and Standard Deviation, σpost of the Times from the Posterior Sample
Planet | Epoch | tobs | σobs | tpost | σpost | Source |
---|---|---|---|---|---|---|
b | 0 | 7322.515 310 | 0.000 710 | 7322.517 902 | 0.000 127 | TS |
c | 0 | 7282.805 700 | 0.001 400 | 7282.805 840 | 0.000 248 | TS |
d | 0 | 7560.797 300 | 0.002 300 | 7560.801 847 | 0.000 384 | TS |
e | 0 | 7312.713 000 | 0.002 700 | 7312.713 920 | 0.000 360 | TS |
f | 0 | 7321.525 200 | 0.002 000 | 7321.522 338 | 0.000 781 | TS |
g | 0 | 7294.786 000 | 0.003 900 | 7294.772 215 | 0.000 905 | TS |
h | 0 | 7662.554 360 | 0.002 000 | 7662.550 913 | 0.001 498 | Spitzer |
Note. Times are in BJDTDB-2,450,000 while uncertainties are in days.
Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.
Download table as: DataTypeset image
Table 15. Mean, tpost, and Standard Deviation, σpost, of Forecast Times from the Posterior Sample
Planet | Epoch | tpost | σpost |
---|---|---|---|
b | 1 | 7259.061 221 | 0.000 152 |
b | 2 | 7260.572 662 | 0.000 163 |
b | 3 | 7262.083 452 | 0.000 159 |
b | 4 | 7263.594 099 | 0.000 150 |
b | 5 | 7265.104 973 | 0.000 151 |
b | 6 | 7266.615 845 | 0.000 153 |
b | 7 | 7268.126 697 | 0.000 155 |
b | 8 | 7269.637 314 | 0.000 148 |
b | 9 | 7271.148 093 | 0.000 146 |
b | 10 | 7272.659 543 | 0.000 157 |
b | 11 | 7274.170 344 | 0.000 152 |
b | 12 | 7275.680 980 | 0.000 144 |
b | 13 | 7277.191 799 | 0.000 143 |
b | 14 | 7278.702 670 | 0.000 145 |
b | 15 | 7280.213 581 | 0.000 150 |
b | 16 | 7281.724 209 | 0.000 143 |
b | 17 | 7283.234 973 | 0.000 139 |
b | 18 | 7284.746 420 | 0.000 150 |
b | 19 | 7286.257 235 | 0.000 146 |
b | 20 | 7287.767 862 | 0.000 137 |
b | 21 | 7289.278 632 | 0.000 136 |
b | 22 | 7290.789 504 | 0.000 138 |
b | 23 | 7292.300 472 | 0.000 144 |
b | 24 | 7293.811 108 | 0.000 137 |
b | 25 | 7295.321 861 | 0.000 133 |
b | 26 | 7296.833 297 | 0.000 143 |
b | 27 | 7298.344 125 | 0.000 140 |
b | 28 | 7299.854 748 | 0.000 131 |
b | 29 | 7301.365 471 | 0.000 129 |
b | 30 | 7302.876 339 | 0.000 131 |
Note. Times are in BJDTDB-2,450,000 while uncertainties are in days. Thirty lines are previewed; full table is available in machine-readable format.
Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.
Download table as: DataTypeset image
Footnotes
- 20
This refers to the condition , which is a generalization of the Laplace resonance with p = 1 and q = 2 (Papaloizou 2014).
- 21
- 22
- 23
The code may be found at https://github.com/ericagol/NbodyGradient.
- 24
Note that as we take m0 = 1 in our simulations, we need to multiply the output of positions and velocities from the code by (M*/M⊙)1/3 to scale to a stellar mass M*.
- 25
In principle, we could include a prior in the Laplace and likelihood-profile analyses.
- 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.
- 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.
- 28
As implemented in the package https://github.com/madsjulia/AffineInvariantMCMC.jl.
- 29
- 30
Note that "M" is being used in three ways here: spectral category (M dwarf), stellar mass (M*), and absolute magnitude in the KS band, .
- 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).
- 32
- 33
For specifications of these programs, see https://www.stsci.edu/jwst/observing-programs/approved-gto-programs.
- 34
- 35
Oxygen fugacity is stated relative to the iron–wüstite equilibrium reaction Fe + 0. 5O2 = FeO (wüstite) such that .