It Takes Two Planets in Resonance to Tango around K2-146

K2-146 was first observed in K2 Campaign 5. The $\sim 2.2{R}_{\oplus }$ mini-Neptune, K2-146 b, was validated independently by Hirano et al. (2018) and Livingston et al. (2018). K2-146 b orbits an M3.0V dwarf and was reported to have an orbital period of 2.645 days. The system was also independently flagged as a planetary candidate by Pope et al. (2016), Libralato et al. (2016), and Dressing et al. (2017). Hirano et al. reported strong TTVs with amplitudes of over 30 minutes. The observed orbital perturbation of the planet is likely caused by either a massive object in its vicinity or an additional object orbiting in or close to an MMR. The stellar parameters of K2-146 are summarized in Table 1.

Table 1.  Stellar Parameters and Photometric Magnitudes of K2-146

Parameter	Value and Uncertainty	Source
EPIC	211924657	a
2MASS	2MASS J08400641+1905346	b
Gaia	661192902209491456	c
R.A.	08 40 06.42	c
Decl.	+19 05 34.42	c
μR.A. $(\mathrm{mas}\,{\mathrm{yr}}^{-1})$	−15.92 ± 0.12	c
μdecl. $(\mathrm{mas}\,{\mathrm{yr}}^{-1})$	−129.02 ± 0.07	c
Parallax (mas)	12.582 ± 0.075	c
RV (km s−1)	40.6 ± 0.2	This work
Spectral type	M3.0V	d
Teff (K)	3385 ± 70	d
$[\mathrm{Fe}/{\rm{H}}]$ (dex)	−0.02 ± 0.12	d
log g	4.906 ± 0.041	d
Ms (M⊙)	0.358 ± 0.042	d
Rs (R⊙)	0.350 ± 0.035	d
Ls (L⊙)	0.015 ± 0.003	d
v  sin i $(\mathrm{km}\,{{\rm{s}}}^{-1})$	<2.00	This work
Photometric Magnitudes
Kep	15.03	a
Gaia G	14.98	c
Johnson B	17.69	e
Johnson V	16.18	e
J	12.18	b
H	11.60	b
K	11.37	b

Note. References for the sources: (a) EXOFOP-K2: https://exofop.ipac.caltech.edu/k2/; (b) The Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006); (c) Gaia DR2 (Gaia Collaboration et al. 2016, 2018); (d) Hirano et al. (2018); (e) The AAVSO Photometric All-sky Survey (APASS; Henden et al. 2009).

Download table as:  ASCIITypeset image

This paper is organized as follows. Section 2 describes the K2 photometric observations, data reduction, and planet detection. Section 3 describes the spectroscopic observations, which are used to rule out the presence of an outer planet companion. Section 4 describes the extraction of transit times for K2-146 b and c, and transit analyses. The TTV model and analysis are described in Section 5. The transit parameters of the mini-Neptune pair are refined following the TTV analysis in Section 6. In Section 7 we evaluate the stability, orbital resonance, and interior composition of the planet paird, and interpret their possible implications for the evolution history of the system. We present our conclusions in Section 8.

We searched the K2 light curves for transit signals using the DST algorithm (Cabrera et al. 2012), which optimizes the fit to the transit shapes with a parabolic function. Figure 2 shows the periodograms of the DST statistics measured in all light curves. The top left panel of Figure 2 shows that the ∼2.6 day signal of K2-146 b was detected in C05, and subsequently in C16 and C18.

Zoom In Zoom Out Reset image size

The 2.6 day signal was then filtered and the light curves were analyzed with the DST algorithm again. Strong peaks at ∼4 days are found in the periodograms of the C16 and C18 data, as shown in the bottom left panel of Figure 2. However, no significant detection is found in C05, and transits of the outer planet were not observed upon visual inspection due to the noise level of the C05 light curve. We also ran our transit search algorithm on the K2SFF C05 data since it has a slightly lower scatter. Although we detected hints of a transit signal at ∼4 days, the detection was not significant. We attribute this to a precessing orbital plane of this outer planet, which we discuss in later sections.

The characterization of the multi-planet system K2-146 follows the approach outlined here: the transit parameters are derived iteratively. We first performed a global analysis to extract the transit times and transit parameters of planet b and planet c (Section 4). The transit parameters of the two planets were analyzed independently using a stacked transit light curve. We then model the transit times of the planets to derive their respective orbital elements (Section 5). Finally, we use information from the TTV-deduced orbital elements to model the stacked transits of planet b and planet c, and improve the precision of the system parameters (Section 6).

To measure the absolute and and projected rotational velocities of K2-146, we use the empirical spectral matching routines described in Stefansson et al. (2020), which closely follow the SpecMatch-emp algorithm described in Yee et al. (2017), adapted for use on a high signal-to-noise ratio (S/N) library of as-observed HPF spectra. Briefly, the spectral matching algorithm compares an as-observed spectrum of the target star to a library of high S/N slowly rotating as-observed spectra with known stellar parameters using a ${\chi }^{2}$ figure of merit. The algorithm first aligns the target and library spectra in wavelength by accounting for the barycentric shift of the target star. The absolute RV of the target star is calculated and accounted for by cross-correlating the target spectrum by a binary mask as described in Stefansson et al. (2020). In doing so, we obtain an absolute RV value of γ = 40.6 ± 0.2 km s⁻¹, where the scatter is obtained by calculating the scatter of the absolute RV independently determined for 8 different HPF orders clean of tellurics.

After aligning the target and library spectra in wavelength, the spectra are resampled to a common wavelength grid. Residuals and a corresponding χ² value are then calculated by subtracting from the target spectrum an artificially broadened version of the library spectrum that is also scaled by a fifth-order Chebychev polynomial to account for any residual systematics. A best-fit v sin i value for each library spectrum is then obtained by minimizing the χ². This optimization is repeated for all of the library stars discussed in Stefansson et al. (2020). After this loop, the top five best-fit spectra are then further optimized to create a final composite spectrum that results in an even better fit to the observed target spectrum.

To arrive at a final v sin i value, and its associated scatter, we independently estimated the v sin i for five HPF orders clean of tellurics, resulting in five independent point estimates of the v sin i. In doing so, we formally obtained a median value of v sin i =1.1 km s⁻¹ with an rms scatter of 200 m s⁻¹. However, as discussed in Stefansson et al. (2020), this value is at the lower limit we estimate to be able to reliable measure the v sin i of the star, given the resolution of HPF of R = 55,000, corresponding to FWHM = 6 km s⁻¹. Even though the derived values are all consistent at the 200 m s⁻¹ level, at this resolution, we estimate that we should be able to detect rotational velocities upwards of ∼2 km s⁻¹. We thus conclude that K2-146 has a v sin i <2 km s⁻¹ from the HPF spectra.

Precise RVs are extracted from the HPF spectra via the method described in Metcalf et al. (2019) and Stefansson et al. (2020). Briefly, our RV extraction code is a modified version of SERVAL (Zechmeister et al. 2018), which determines RV by template-matching (using a χ² metric) the stellar spectrum at a single epoch with a "master" template made by coadding all spectra of the target. This method has proven especially useful for M dwarfs (Anglada-Escudé & Butler 2012), whose complex spectra are not as amenable to cross-correlation with a binary mask (e.g., Baranne et al. 1996) as hotter stars. The RV measurements of K2-146 are shown in Figure 3 and the values are listed in Appendix A.

Zoom In Zoom Out Reset image size

We used the Lomb–Scargle (LS) periodogram (Lomb 1976; Scargle 1982) by implementing the fast algorithm of Press & Rybicki (1989) to search for possible periodic signals in the RV data. The resulting power spectrum is shown in the upper right panel of Figure 3. We computed the corresponding false-alarm-probabilities (FAP) using the approximation method proposed by Baluev (2008), and found no significant peaks with a FAP of less than 10%. The rms of the RV data is $\sigma =25.95$ m s⁻¹. This value was used as the 1σ upper limit in the RV semi-amplitude, K. We determined the upper mass limit (M_p sin i) of a possible outer planet by computing ${M}_{p}\sin {iK}/P$ over a range of periods (P). The bottom panel of Figure 3 shows the 1σ and 3σ mass upper limits for an outer planet in a circular orbit. We are able to exclude the presence of an outer planet with M_p sin i > 86 M_⊕ in a circular orbit for periods up to ∼76 days (or a < 0.25 au).

We extracted the transit times using the Python Tool for Transit Variations (PyTTV; J. Korth 2020, in preparation). This tool uses PyTransit (Parviainen 2015) for transit modeling, PyDE³⁷ for optimization and emcee (Foreman-Mackey et al. 2013) for posterior sampling.

Before modeling, the light curves around each transit were detrended by subtracting a second order polynomial fit to the out-of-transit light curve. These cutout segments were then the input for the detailed modeling with PyTTV. Some transits were excluded from the analysis (see Figure 1): transit numbers 3, 8, 384, 424, and 427 for planet b and transit numbers 2, 10, 14, and 17 for planet c.³⁸

The transit time extraction and transit fits were performed for each planet independently. The individual transits are fitted collectively with the Mandel & Agol (2002) model, each with their own transit center but sharing the rest of the transit parameters (individual transit fitting models are found in Figures B.3 and B.4). The optimization was done by computing the log-likelihood, log L, in our code to estimate transit parameters:

$$\begin{eqnarray}&&\mathrm{log}L=-\displaystyle \frac{1}{2}\sum _{i=0}^{N}\left[\displaystyle \frac{{({x}_{i}-{\mu }_{i})}^{2}}{{\sigma }_{i}^{2}}+\mathrm{ln}(2\pi {\sigma }_{i}^{2})\right],\end{eqnarray} \tag{ 1 }$$

where x_i is the model, μ_i is the data, and σ_i is the error in the data for the ith point, respectively. The fitted parameters are radius ratio k, impact parameter b, and stellar density ρ_s, all with uniform priors. We used a quadratic limb-darkening model with the "triangular sampling" parameterization presented by Kipping (2013). Gaussian priors were imposed on the limb-darkening coefficients, u₁ and u₂, where the central values of the coefficients were calculated with PyLDTK (Parviainen & Aigrain 2015), which utilizes the spectrum library of Husser et al. (2013). The transits of planet c is grazing, the a/R_s and hence the mean stellar density is less well constraint. Thus the stellar density posterior from the planet b analysis was used as a stellar density prior in the planet c analysis. To account for the long exposure times of the K2 observation we applied a supersampling (n = 10) as suggested in Kipping (2010).

For posterior sampling, we ran five MCMC chains with 5000 steps whereby the previous run was used as a burn-in for the current MCMC run. The chains were checked for convergence visually. The posteriors and the derived planetary parameters are shown in Appendix B. The resulting corner plots are shown in Figures B.1 and B.2. The transit times of both planets are used for a detailed TTV analysis in the following section.

The fitted transit times of planet b (red) and planet c (blue) are shown in the O − C diagram in Figure 4. An anti-correlation between the O − C values is clearly visible in campaigns 16 and 18. This shows that the two planets are orbiting in the same system and that they are gravitationally interacting with each other, confirming the planetary nature of the signals. For a better visualization, the inset in Figure 4 shows a magnification of the different campaigns. The colored error bars mark the 1σ and 3σ uncertainties of the fitted transit times.

Zoom In Zoom Out Reset image size

An independent transit analysis of K2-146 b and c was performed using stacked transit light curves. The light curves were cut such that only data within eight transit durations, centered on each transit, were used in our analysis. The selected transits of K2-146 b and K2-146 c used in our analysis are indicated by red and blue lines, respectively, in Figure 1.

We employed a Markov-chain Monte Carlo (MCMC) approach to derive the system parameters of K2-146. The stacked transits of K2-146 b and K2-146 c were modeled using the analytical functions by Mandel & Agol (2002). The transit model was implemented using the package PyTransit (Parviainen 2015), and a quadratic stellar limb-darkening law was applied. The fitted transit parameters are the planet-to-star radius ratios, k_b and k_c, the orbital inclinations, i_b and i_c, stellar density, ρ_s, and the triangle sampling (Kipping 2010) of the quadratic stellar limb-darkening coefficients u₁ and u₂, where uniform priors were used. The orbital periods of the planets are kept fixed.

The MCMC method was implemented using the Python package emcee (Foreman-Mackey et al. 2013) for Bayesian parameter estimation. χ² statistics was used for likelihood estimation in our model. We computed the log-likelihood, log L, following Equation (1). An initial burn-in phase of 20 MCMC chains ×10,000 steps was implemented to optimize the convergence of the fit. To obtain reasonable uncertainties in the transit parameters, we rescaled the error bars such that the value of the reduced χ² equals 1. We then initiated 100 MCMC chains of 5 × 10⁴ steps to sample the posterior space. We checked for convergence and discarded the first 2000 steps, then adopted the median, 16th, and 84th percentiles of the samples in the marginalized posterior distributions as the fitted values and their 1σ uncertainties. The results of the fitted transit parameters of K2-146 b and K2-146 c are presented in Appendix B. The best-fitted transit parameters of K2-146 b and K2-146 c obtained here are generally consistent within ∼1σ with those derived by PyTTV.

The stacked transit light curves of K2-146 b and c were reanalyzed incorporating the eccentricity information from the TTV analysis. The transits were modeled with Transit and Light Curve Modeller (TLCM; Csizmadia 2020), a software tool for joint RV and transit light-curve fits, or transit light-curve fits only. It utilizes the Mandel & Agol (2002) and Eastman et al. (2013) subroutines to calculate the transit light-curve shapes for every moment when we have an observation. A wavelet-filter (Carter & Winn 2009) can be applied to model the red-noise effects. Contamination is also taken into account. The light-curve part uses a quadratic limb-darkening law. The stellar radius, based on the value measured by Hirano et al. (2018), as well as the spectroscopic log g values, can be used as priors for the fit. Parameter estimation is done via the Genetic Algorithm, refined by Simulated Annealing (Geem 2001). The final parameter estimation is done using several chains of MCMC with at least 10⁵ steps. The median and the width of the chains will define the final adopted solutions and their uncertainty ranges. The Gelman–Rubin statistic (e.g., Croll 2006) is used to check the convergence of chains. For detailed descriptions of TLCM, we refer the reader to the following works: Csizmadia et al. (2011, 2015), Smith et al. (2017).

The fitted parameters are the epoch of mid-transit, T₀, the scaled semimajor axis (a/R_s), planet-to-stellar radius ratio (R_p/R_s), the impact parameter, b, and the quadratic stellar limb-darkening coefficients u₊ = u₁ + u₂, and u₋ = u₁ − u₂. The orbital periods, P, of the planets were kept fixed. In addition, we used the eccentricity and argument of periastron derived from Section 5 as priors to perform our analysis. Table 3 presents the best-fit transit parameters of K2-146 b and K2-146 c. The resulting best-fit transit models of the inner and outer planets are shown in Figure 6.

Zoom In Zoom Out Reset image size

Table 3.  Best-fit Planet Parameters of K2-146 b,c from a Stacked Transit Analysis, and Their Corresponding 1σ Uncertainties

Parameter	Description (unit)	Values and Uncertainties
		K2-146 b	K2-146 c
P	Period (day)	2.6698	3.9663
T0	Epoch (day from transit center)	0.00009 ± 0.00038	−0.00015 ± 0.00076
Rp	Radius (${R}_{\oplus }$)	2.25 ± 0.10	${2.59}_{-0.39}^{+1.81}$
a	Semimajor axis (au)	0.0248 ± 0.0002	0.0327 ± 0.0006
b	Impact parameter	0.391 ± 0.069	0.930 ± 0.097
i	Inclination (°)	88.5 ± 0.3	87.3 ± 0.3
Rp/Rs	Scaled planet radius	0.0589 ± 0.0014	${0.0680}_{-0.0254}^{+0.1226}$
a/Rs	Scaled semimajor axis	15.250 ± 0.126	20.064 ± 0.412
${u}_{+}={u}_{1}+{u}_{2}$	Combined limb-darkening coefficient	0.575 ± 0.171	0.678 ± 0.169
${u}_{-}={u}_{1}-{u}_{2}$	Combined limb-darkening coefficient	0.022 ± 0197	0.072 ± 0.210
ρp	Density (${\rm{g}}\,{\mathrm{cm}}^{-3}$)	2.702 ± 0.494	${2.246}_{-1.846}^{+1.883}$

Note. The orbital periods P of the planets were kept fixed at the values derived from the TTV analysis.

Download table as:  ASCIITypeset image

The outer planet K2-146 c was only found to transit in the latter campaigns C16 and C18. This suggests orbital plane precession is at play. We searched for seasonal changes in the impact parameters of both planets. The transits of the inner and outer planets were stacked separately in each K2 campaign. The mean orbital elements of each campaign were obtained from the posterior samples of the TTV analysis. We fixed $\sqrt{e}\sin \omega$ and $\sqrt{e}\cos \omega$ values on these mean values for the analysis of separate campaigns. We then ran TLCM on the campaign-stacked data; the free parameters in each campaign were the scaled semimajor axis, impact parameter, radius ratio, and epoch. The radius ratio was the same from campaign to campaign.

For planet b, we obtained an impact parameter of b = 0.42 ± 0.11, b = 0.25 ± 0.10, and b = 0.41 ± 0.15, for C05, C16, and C18 respectively. While the impact parameter is practically the same in C05 and C18, it differs by approximately 1σ in C16 in comparison to the two other campaigns. This indicates that we do not see a significant change in the impact parameter of planet b in these data. In the case of planet c, we obtained an impact parameter of b = 0.94 ± 0.10 and b = 0.92 ± 0.10 for C16 and C18, respectively. Again, we do not see significant changes in the impact parameter of the planet, although the null detection of transits in C05 does suggest that it has been drifting. We attribute this to a shorter time baseline between C16 and C18 ($\sim 200\,\mathrm{days}$) compared to the $\sim 1000\,{\rm{d}}{\rm{a}}{\rm{y}}$ separation between C05 and C16.

We investigated the dynamical stability for the K2-146 multi-planet system, independent of the TTV analysis, to obtain the mass limits on the two Neptune-size planets and to establish the stability of the system.

7.1.1. Hill Radius of K2-146

7.1.2. N-body Simulations of K2-146

The N-body simulation code Mercury6 (Chambers 1999) was used to study the orbital evolution and determine the dynamic stability of the system. We chose the "hybrid symplectic and Bulirsch–Stoer integrator" mode of Mercury6 to compute close encounters in the system. We adopted the stellar mass and radius reported in Hirano et al. (2018) for the central star.

We performed some initial simulations to place upper mass limits on the two planets based on their dynamic stability. We employed 5105 numerical integrations, each with an integration period of 2 Myr. An initial step size of 2.7 days was selected; subsequent step sizes were adjusted by the variable time-step algorithms in the program to maintain integration accuracy. The orbital parameters of the system were recorded every 2 yr. For each integration, the orbital periods (P_b and P_c) were chosen from a Gaussian distribution using the center and $1\sigma$ reported by PyTTV in Section 4. The eccentricities of the planets (e_b and e_c) were drawn from a uniform distribution between e = 0 and e = 0.3. The two planets were assumed to be coplanar with fixed inclinations and longitudes of the ascending node where i_b, i_c = 90° and Ω_b, Ω_c = 0°. The arguments of periastron (${\omega }_{{\rm{b}}}$ and ω_c) and mean anomalies of the planets were drawn from a uniform sample where $0^\circ \lt {\omega }_{{\rm{b}}},{\omega }_{{\rm{c}}}\lt 360^\circ$ and $0^\circ \lt {{ \mathcal M }}_{{\rm{b}}},{{ \mathcal M }}_{{\rm{c}}}\lt 360^\circ$. The planetary masses of K2-146 b and K2-146 c (M_b and M_c respectively) were randomly drawn from a half-normal distribution with a width of 1000 ${M}_{\oplus }$. In each simulation, the system becomes unstable when (1) the planets collide with one another or with the central star; or (2) the planets are ejected from the system, i.e., ${e}_{{\rm{b}}},{e}_{{\rm{c}}}\gt 1$ and/or a > 10 au. Our results showed that a significant fraction of systems are dynamically stable for 2 Myr for planets with masses up to ∼30 ${M}_{\oplus }$. We found that over 90% of these systems remained stable for 2 Myr if the mass ratio of K2-146 b and K2-146 c is close to unity.

We further used a subset of the posterior samples from our TTV analysis in Section 5 to check the stability of the system. We performed 153 simulations where the starting orbital parameters of the system were defined by randomly selecting posterior samples from the TTV analysis. Again, each simulation was integrated for a period of 2 Myr, the step sizes were automatically adjusted by the algorithm, and the orbital parameters were recorded for every 2 yr. All 153 of our simulations were shown to be stable over at least 2 Myr. Figure 7 shows the planet masses drawn from our TTV posterior sample.

Zoom In Zoom Out Reset image size

The orbital periods of K2-146 b and K2-146 c show a 3:2 commensurability. To assess whether the planet pair is trapped in a 3:2 MMR, we monitored the orbital evolution simulation of the so-called resonance angles over a 2000 yr long interval.

The resonance arguments, Θ₁ and Θ₂, for a pair of planets in a 3:2 MMR are defined as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Theta }}}_{1} & = & 3\cdot {\lambda }_{{\rm{c}}}-2\cdot {\lambda }_{{\rm{b}}}-{\varpi }_{{\rm{b}}},\\ {{\rm{\Theta }}}_{2} & = & 3\cdot {\lambda }_{{\rm{c}}}-2\cdot {\lambda }_{{\rm{b}}}-{\varpi }_{{\rm{c}}}.\end{array}\end{eqnarray} \tag{ 2 }$$

The quantity $\lambda ={ \mathcal M }+\varpi$ is often called the mean longitude, and ${ \mathcal M }$ is the mean anomaly. The longitude of pericenter is $\varpi =\omega +{\rm{\Omega }}$ (e.g., Murray & Dermott 1999). The angle between two planets at conjunction is 0°. One can think of each resonant angle as tracking the angle at which conjunctions happen, relative to each planet's pericenter (e.g., when λ_c = λ_b, the conjunction occurs at λ_b − ϖ_b = Θ₁).

We performed 20,000 simulations to evaluate the resonance angles of the planet pair. In each simulation, the planet-to-star mass ratios and orbital elements of the planets were drawn from the posterior samples of the TTV analysis. The stellar mass in Table 1 was used to convert the mass ratios into planetary masses. Then the orbits were numerically integrated for 2000 yr. The values of the relevant parameters were saved for every 36 days of the integration, then the resonant angles were calculated. For each integration, we recorded the maximum amplitude of the resonant angles. We took the modulo 360° values of the resonant variables.

Figure 8 shows the histograms of the maximum half-amplitude of the resonant angles. The wide peak at around Θ₁/2 = 135° is because the location of conjunctions librates with a half-amplitude of ≈135° around a location roughly anti-aligned with the pericenter of planet b's orbit. We found approximately 74% of simulations show libration with a half-amplitude of around 135°. Figure 9 shows the simulation that gave the smallest libration angle. For the remaining cases, the resonance angles either always circulate, or chaotically switch from libration to circulation. The special cases can be tested in more detailed dynamics studies in the future.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The widths of the peaks in these histograms are caused by the uncertainties of the masses and orbital elements derived from TTVs. To increase their precisions, more transit observations are needed from this system. Despite the faintness of the host star, which leads to lower accuracy in TTVs, it would be worthwhile to try to observe it with high cadence with Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) and PLAnetary Transits and Oscillations (PLATO of stars; Rauer et al. 2014). Alternatively, observations may also be obtained with the less precise (∼10 minutes) timing values from CHaracterising ExOPlanet Satellite (CHEOPS; Broeg et al. 2013), because the precision could be compensated by many hours of TTV-values in this case.

Our TTV analysis in Section 5 and transit parameters refined in Section 6 revealed that K2-146 b has a mass and radius of 2.25 ± 0.10 ${R}_{\oplus }$ and 5.6 ± 0.7 ${M}_{\oplus }$, respectively, corresponding to a bulk density of 2.702 ± 0.494 ${\rm{g}}\,{\mathrm{cm}}^{-3}$. The mass and radius of the outer planet K2-146 c are ${2.59}_{-0.39}^{+1.81}$ ${R}_{\oplus }$ and 7.1 ± 0.9 ${M}_{\oplus }$, respectively, which correspond to a bulk density of ${2.246}_{-1.846}^{+1.883}$ ${\rm{g}}\,{\mathrm{cm}}^{-3}$. K2-146 b and K2-146 c are orbiting at distances of 0.0248 au and 0.0327 au, respectively. We assumed the planets have an albedo of 0, and a reradiation factor of one-fourth, where atmospheric circulation redistributes the energy around the planetary atmosphere, then re-radiates the energy back into space. Under these assumptions, the equilibrium temperatures of planet b and planet c are approximately 590 K and 520 K, respectively.

Figure 10 shows a mass–radius plot of known planets with M_p < 30 ${M}_{\oplus }$ where the mass of the planets are determined with a precision better than approximately 30%. The solid, dashed, and dashed dot lines represents the mass–radius relations for different planetary compositions as derived in Zeng et al. (2016, 2019). The color of each data point indicates a planet's equilibrium temperature corresponding to the color bar. The masses and radii of K2-146 b and c are consistent with cases of a 100% H₂O interior, a water-rich core with the addition of an H₂O-gaseous atmosphere, or an Earth-like rocky interior with a small fraction of H₂ envelope. The large radius uncertainty of K2-146 c means that it could have a more massive H₂ envelope.

Zoom In Zoom Out Reset image size

Observational evidence of planet orbital architectures allows one to place constraints on the formation and dynamical evolution of the system. In the case of K2-146, we observed a number of interesting traits:

1. Masses and radii of the planets are consistent with a water-dominated core and the presence of water or an H₂ envelope—it was previously proposed that in situ formation of mini Neptunes and super Earths is possible if 50–100${M}_{\oplus }$ of rocky material is delivered to the inner disk for planet assembly (Hansen & Murray 2012). However, the large fraction of solids would drift toward the host star on a relatively short timescale, preventing in situ formation. Furthermore, in situ formation is unlikely to produce a significant fraction of atmospheric masses for close-in planets (Schlichting 2014; Inamdar & Schlichting 2015). Therefore, close-in planets with an atmosphere were likely formed at larger separation from the star in the presence of a gas disk. Subsequently, planets accrete gaseous envelopes as they migrate inward toward their current locations.

2. The evolution of resonance arguments for the planet pair suggests that K2-146 b and c are likely trapped in a 3:2 mean motion resonance—during formation, planets interact with protoplanetary disks. This drives the migration of planets inward through the disk due to exchange of angular momentum (Goldreich & Tremaine 1979, 1980; Lin & Papaloizou 1979). Convergent migration can occur in one of two ways: (1) planets formed at wide separations can move toward one another with different migration speeds; (2) planets formed in close proximity are massive enough to form a gap in the disk where inner and outer disks push the planets toward each other. When the planet orbital periods approach a commensurability, dynamical interactions that follow cause planets to migrate collectively inward while preserving period commensurability (Snellgrove et al. 2001). Under favorable disk parameters, planet masses, and migration speeds, the planet pair can enter a 3:2 MMR after breaking the 2:1 MMR barrier (e.g., Kley 2000; Nelson & Papaloizou 2002), such as the case of HD 45364 (Correia et al. 2009) and KOI-1599 (Panichi et al. 2019).

3. Our TTV model revealed that both the inner and outer planets have moderate eccentricities of ${e}_{b}=0.14\pm 0.07\,{\rm{a}}{\rm{n}}{\rm{d}}{{\rm{e}}}_{{\rm{b}}}=0.16\pm 0.07,$ and are apsidally anti-aligned, i.e., ${\rm{\Delta }}\omega \,={\omega }_{{\rm{b}}}-{\omega }_{{\rm{c}}}\approx 180$°—convergent migration could have played a role in the observed eccentricity in the K2-146 planet pair. After the planets are captured in resonance, the planet pair migrates inward while maintaining resonance, which leads to orbital eccentricity excitation (Lee & Peale 2002; Batygin & Morbidelli 2013). While the protoplanetary disk is present, the planets could experience eccentricity damping as they migrate. The moderate eccentricities observed in K2-146 b and c imply that the migration process must be fast enough in order to minimize the damping efficiency. The anti-aligned apsides are a natural result of such a migration process. The conjunctions occur when K2-146 c is near periapse and K2-146 b is near apoapse. The longitudes of periapse of both planets are required to precess at the same rate for a stable configuration, such that the lines of apsides are locked in an anti-aligned state (Lee & Peale 2002). Due to the close proximity of the two planets, such a mechanism is present to avoid close encounters.

4. The outer planet K2-146 c showed a change in impact parameter—gravitational interaction between planets can give rise to apsidal and nodal precession around the host star. The change in impact parameter of the outer planet observed from K2 Campaign 5 to Campaign 16 is likely an indication of orbital precession. A misalignment between the planet orbits can be suspected, resulting in a precession of the line of nodes of planet c (Miralda-Escudé 2002). The precession would then lead to a change in the length of the transit chord. In our TTV analysis, we assumed the planets have coplanar orbits because, if the mutual inclination is large, the two planets are unlikely to transit simultaneously even if their orbits are precessing (see Mills & Fabrycky 2017). This assumption could be tested directly via a joint modeling of TTVs and TDVs, or a photodynamical model of the light curves, which will enable a measurement of the mutual orbital inclination through the constraint on the nodal precession rate (e.g., Kepler-117 Almenara et al. 2015, Kepler-108 Mills & Fabrycky 2017).