It takes two planets in resonance to tango around K2-146

TL;DR

This study discovers and characterizes a second planet in the K2-146 system, revealing a resonant, gravitationally interacting pair with specific orbital and physical properties, based on transit timing variations and numerical simulations.

Contribution

The paper reports the detection of a second planet in the K2-146 system and analyzes their resonant interaction using TTVs and orbital modeling, which is a novel insight into this system.

Findings

Two planets are likely in a 3:2 mean motion resonance.

The planets have moderate eccentricities around 0.14 and 0.16.

The system is dynamically stable for at least 2 million years.

Abstract

K2-146 is a cool, 0.358 M_sun dwarf that was found to host a mini-Neptune with a 2.67-days period. The planet exhibited strong transit timing variations (TTVs) of greater than 30 minutes, indicative of the presence of a further object in the system. Here we report the discovery of the previously undetected outer planet, K2-146 c, in the system using additional photometric data. K2-146 c was found to have a grazing transit geometry and a 3.97-day period. The outer planet was only significantly detected in the latter K2 campaigns presumably because of precession of its orbital plane. The TTVs of K2-146 b and c were measured using observations spanning a baseline of almost 1200 days. We found strong anti-correlation in the TTVs, suggesting the two planets are gravitationally interacting. Our TTV and transit model analyses revealed that K2-146 b has a radius of 2.25 $\pm$ 0.10 \R_earth and a mass of 5.6 $\pm$ 0.7 M_earth, whereas K2-146 c has a radius of $2.59_{-0.39}^{+1.81}$ R_earth and a mass of 7.1 $\pm$ 0.9 M_earth. The inner and outer planets likely have moderate eccentricities of $e = 0.14 \pm 0.07$ and $0.16 \pm 0.07$, respectively. Long-term numerical integrations of the two-planet orbital solution show that it can be dynamically stable for at least 2 Myr. The evaluation of the resonance angles of the planet pair indicates that K2-146 b and c are likely trapped in a 3:2 mean motion resonance. The orbital architecture of the system points to a possible convergent migration origin.

Figures

Tables

Table 1. Table 1: Stellar parameters and photometric magnitudes of K2-146.

\topruleParameter	Value and uncertainty	Source
EPIC	211924657	a
2MASS	2MASS J08400641+1905346	b
Gaia	661192902209491456	c
RA	08 40 06.42	c
DEC	+19 05 34.42	c
${\mu }_{RA}$ $\left[\mathrm{mas}/\mathrm{yr}\right]$	$-15.92\pm 0.12$	c
${\mu }_{Dec}$ $\left[\mathrm{mas}/\mathrm{yr}\right]$	$-129.02\pm 0.07$	c
Parallax $\left[\mathrm{mas}\right]$	$12.582\pm 0.075$	c
Spectral type	M3.0V	d
$T_{\mathrm{eff}}$ $\left[K\right]$	$3385\pm 70$	d
$\left[\mathrm{Fe}/H\right]$ $\left[dex\right]$	$-0.02\pm 0.12$	d
$\log g$	$4.906\pm 0.041$	d
Ms $\left[M_{\odot }\right]$	$0.358\pm 0.042$	d
Rs $\left[R_{\odot }\right]$	$0.350\pm 0.035$	d
Ls $\left[L_{\odot }\right]$	$0.015\pm 0.003$	d
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

Table 2. Table 2: Masses and orbital elements for K2-146 b and c determined from TTV modeling.

	K2-146 b	K2-146 c
Fitted parameters
$M_p/M_s$ [$\times {10}^{-5}$]	$4.7\pm 0.2$	$6.0\pm 0.2$
$P$ [days]	$2.6698\pm 0.0001$	$3.9663\pm 0.0002$
$\sqrt{e}\cos \omega$	$-{0.36}_{-0.08}^{+0.11}$	${0.40}_{-0.10}^{+0.08}$
$\sqrt{e}\sin \omega$	$-{0.07}_{-0.08}^{+0.09}$	$0.01\pm 0.07$
$t_0$ [$\mathrm{BJD}-2454833$]	$3467.4345\pm 0.0007$	$3466.6019\pm 0.0009$
Derived parameters
$M_p$ [$M_{\oplus }$]aaDerived using $M_s=0.358\pm 0.042M_{\odot }$.	$5.6\pm 0.7$	$7.1\pm 0.9$
$e$	$0.14\pm 0.07$	$0.16\pm 0.07$
$\omega$ [deg]	${191}_{-15}^{+11}$	$2_{-10}^{+12}$

Table 3. Table 3: Best-fit planet parameters of K2-146 b,c from a stacked transit analysis, and their corresponding 1 − σ 1 𝜎 1-\sigma uncertainties. The orbital periods, P 𝑃 P of the planets were kept fixed at the values derived from the TTV analysis.

Parameter	Description [unit]	Values and uncertainties
		K2-146 b	K2-146 c
$P$	Period [day]	$2.6698$	$3.9663$
$T_0$	Epoch [day from transit center]	$0.00009\pm 0.00038$	$-0.00015\pm 0.00076$
$R_p$	Radius [$R_{\oplus }$]	$2.25\pm 0.10$	${2.59}_{-0.39}^{+1.81}$
$a$	Semi-major axis [AU]	$0.0248\pm 0.0002$	$0.0327\pm 0.0006$
$b$	Impact parameter	$0.391\pm 0.069$	$0.930\pm 0.097$
$i$	Inclination [∘]	$88.5\pm 0.3$	$87.3\pm 0.3$
$R_p/R_s$	Scaled planet radius	$0.0589\pm 0.0014$	${0.0680}_{-0.0254}^{+0.1226}$
$a/R_s$	Scaled semi-major axis	$15.250\pm 0.126$	$20.064\pm 0.412$
$u_{+}=u_1+u_2$	Combined limb-darkening coefficient	$0.575\pm 0.171$	$0.678\pm 0.169$
$u_{-}=u_1-u_2$	Combined limb-darkening coefficient	$0.022\pm 0197$	$0.072\pm 0.210$
${\rho }_p$	Density [$g{\mathrm{cm}}^{-3}$]	$2.702\pm 0.494$	${2.246}_{-1.846}^{+1.883}$

Table 4. Table 4: Transit parameters of K2-146 b and K2-146 c derived from PyTV analysis and stacked transit analysis. The reported posterior values reported are the median and 1- σ 𝜎 \sigma uncertainties.

Parameter	PyTV		Stacked transit analysis
	K2-146 b	K2-146 c	K2-146 b	K2-146 c
P [day]	$2.65044\pm 0.00007$	$3.98974\pm 0.00069$	-	-
$R_p/R_s$	$0.057\pm 0.001$	$0.084\pm 0.013$	$0.0527\pm 0.0004$	$0.0604\pm 0.0046$
b	$0.27\pm 0.19$	$0.96\pm 0.02$	$0.03\pm 0.03$	$0.81\pm 0.03$
$a/R_s$	$15.06\pm 1.16$	$20.44\pm 1.19$	$15.51\pm 0.03$	$20.38\pm 0.07$
$i$ [∘]	$88.9\pm 0.9$	$87.3\pm 0.2$	$89.9\pm 0.11$	$87.7\pm 0.09$
$R_p$ [$R_{\oplus }$]	$2.18\pm 0.22$	$3.16\pm 0.59$	$2.01\pm 0.02$	$2.31\pm 0.18$
a [AU]	$0.0246\pm 0.0031$	$0.0332\pm 0.0039$	$0.0252\pm 0.0166$	$0.0332\pm 0.0001$

Equations

$$\log L=-\frac{1}{2}\sum _{i=0}^N\left[\frac{(x_i-{\mu }_i{)}^2}{{\sigma }_i^2}+\ln (2\pi {\sigma }_i^2)\right],$$

$${\Theta }_1$$

$${\Theta }_2$$

$$\frac{n_c}{n_b}=\frac{P_b}{P_c}=\frac{p}{p+q}$$

Full text

methods: observational — techniques: photometric — planets and satellites: detection — stars: individual (K2-146)

††journal: ApJ

1 Introduction

The Kepler (Borucki et al., 2010) and K2 (Howell et al., 2014) missions have brought about many exciting discoveries since the spacecraft was launched in 2009. Statistical studies using the Kepler planet sample revealed that sub-Neptune size planets with $R_{p}<4$ $\rm R_{\oplus}$  are by far the most common type of planets in the galaxy (e.g. Borucki et al. (2011); Howard et al. (2012); Batalha et al. (2013); Dressing & Charbonneau (2013); Petigura et al. (2013); Fressin et al. (2013)). Previous works have also shown that short-period planets with radii between 1.5-6 $\rm R_{\oplus}$ are common in near co-planar multi-planet systems (Lissauer et al., 2011; Fabrycky et al., 2014). More recently, Weiss et al. (2018) used precisely determined stellar and planetary parameters to show that, multi-planet systems are dynamically packed and that adjacent planets in the same system are likely to have similar sizes.

Multi-planet systems are of particular interest because these systems can provide insights on the formation and evolution of our own Solar system. Gaining more knowledge about these systems is important for understanding the planetary system dynamics. To address these topics we need to know the planetary parameters, in particular their radii and masses. Unfortunately, only a small number of planets with precisely measured masses are known because there are inadequate telescope resources for sufficient spectroscopic measurements, or the stars are simply too faint. The mass determination for Earth- or Neptune-sized planets around Sun-like stars or faint stars by radial velocity (RV) is particularly difficult with currently available telescopes and instrumental technique because of the small Doppler reflex motion of the host star.

In multi-planet systems, planets can experience mutual gravitational interactions that perturb their orbits. One of the consequences of these is that individual transits vary periodically around a mean orbital period. This effect is referred to Transit Timing Variations (TTVs, e.g.  Holman & Murray, 2005; Agol et al., 2005). This effect is most prominent when the orbital periods of the planets are close to a Mean-Motion Resonance (MMR), and it can be measured even for low mass planets. Thus TTVs are sometimes the only chance to characterize the planetary system. For example, a TTV analysis of KOI-142 revealed a pair of planets orbiting in a near 2:1 resonance (Nesvorný et al., 2013). Extensive RV observations were obtained for K2-19 b and c, a two planet system in a near 3:2 MMR (Armstrong et al., 2015; Narita et al., 2015; Barros et al., 2015; Nespral et al., 2017). The measured RV masses of the Neptune-sized planets were found to be consistent with the TTV-derived masses. Precise mass determination via TTVs showed a pair of planets, Kepler-36 b and c, that have distinctly different bulk densities, hinting at different formation origin of the planets (Carter et al., 2012). Dynamical modelling of both the TTV and transit duration variations (TDVs) can also reveal mutual inclination of a pair of planets (e.g. Kepler-108; Mills & Fabrycky, 2017) and uncover the presence of an additional non-transiting companion in some cases (e.g. KOI-872 system; Nesvorný et al. 2012, Kepler-448 b and Kepler-693 b; Masuda 2017).

1.1 The cool star K2-146

K2-146 was first observed in the K2 Campaign 5. The $\sim 2.2R_{\oplus}$ mini-Neptune, K2-146 b, was validated independently by Hirano et al. (2018) and Livingston et al. (2018). K2-146 b orbits around 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 amplitude 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 a MMR. The stellar parameters of K2-146 are summarised in Table 1.

This paper is organised as follows. Section 2 describes the K2 photometric observations, data reduction and planet detection. Section 3 describes the extraction of transit times for K2-146 b and c. The TTV model and analysis are described in Section 4. Subsequent derivation of transit parameters of the mini-Neptune pair are presented in Section 5. In Section 6 we evaluate the stability, orbital resonance, and interior composition of the planets pair, and interpret their possible implications on the evolution history of the system. Finally, we draw our conclusions in Section. 7.

2 K2 photometry

K2-146 was observed in Campaign 5, 16 and 18 (hereafter, C05, C16, and C18, respectively) in the long cadence mode. The photometric observations were obtained between 2015 April 27 and 2018 July 02, spanning a baseline of almost 1200 days. The K2 target pixel data was downloaded from the Mikulski Archive for Space Telescope111http://archive.stsci.edu/kepler/data_search/search.php (MAST). A custom pipeline was implemented for light curve reduction and is described below.

The photometric analysis was conducted for each campaign separately. For each campaign, the timestamps were combined and the quality of the light curve was tested using different thresholds to the number of counts per pixel. An optimal aperture was selected using the 100 counts per pixel threshold. Using this aperture we calculated the flux for each frame. To correct for possible correlation with the movement of pointing, we cut the light curves into segments with a length of 4 days. Between adjacent segments, there is an overlapping region of 0.8 days. The overlapping regions help us to avoid edge effects when fitting the data in the time domain to remove stellar variability.

For each light curve segment, outliers (including transits) were identified and masked before fitting a multidimensional polynomial, over the POS_CORR columns (which measures the relative motion of the star), to the data. To avoid influence from time-dependant variability, we simultaneously obtained a third order polynomial fit over the time. This fit was then applied to the whole segment, including the transits, to correct for possible correlation with the telescope pointing. After the correlation to the POS_CORR columns have been removed, we masked outliers and transit events again before fitting a seventh order polynomial to the light curve segments to correct for stellar variability. Finally, all light curve segments are normalized and stacked together.

We compared light curves generated from our custom pipeline with ones that are publicly available from the Vanderburg & Johnson (2014) pipeline (K2SFF) and from Luger et al. (2018) (EVEREST; kindly provided by Luger). We found that the noise level of the light curves from the three pipelines are comparable. Light curves from Campaign 16 and Campaign 18 generated from our pipeline have a slightly lower overall scatter, whereas the scatter in the Campaign 5 lightcurve is slightly higher than the K2SFF pipeline. For a consistent analysis, we opted to use the light curves obtained from our pipeline (as shown in Figure 1) for light curve modelling and TTV analysis.

2.1 Planet detection

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 $\sim 2.6$ days signal of K2-146 b was detected in the C05, and subsequently in C16 and C18.

The 2.6-day signal was then filtered and the light curves were analyzed with the DST algorithm again. Strong peaks at $\sim 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 transit signal at $\sim 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 3). The transit parameters of the two planets were analysed independently using a stacked transit light curve. We then model the transit times of the planets to derive their respective orbital elements (Section 4). 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 5).

3 Transit time extraction

3.1 PyTV

We extracted the transit times using the Python Tool for Transit Variations (PyTV, Korth 2019, in prep.). This tool uses PyTransit (Parviainen, 2015) for transit modelling, PyDE222 https://github.com/hpparvi/PyDE for optimisation and emcee (Foreman-Mackey et al., 2013) for posterior sampling.

Before modelling, the light curve around each transit were detrended by subtracting a second order polynomial fit to the out-of-transit light curve. These cut-out segments were then the input for the detailed modelling with PyTV. Some transits were excluded from the analysis (see Figure 1): transit number 3, 8, 384, 424 and 427 for planet b and transit numbers 2, 10, 14 and 17 for planet c.333We checked the light curve flags and found that transit number 3 and 424 both have flags 2048 (impulsive outlier) and 1048576 (thruster firing), transit number 8 has flag 8192 (cosmic ray), transit number 384 has flag 1024 (sudden sensitivity dropout), and transit number 427 has flag 32768 (no fine point). For planet c the transits with numbers 2 and 10 have flag 1048576. For transit number 14 and 17 we found no non-zero flags but they might not be visible because of the long cadence observations (30 minutes) and are therefore missed in the light curve.

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 appendix Figures 12 and 13). The optimization was done by computing the log-likelihood, $\log L$, in our code to estimate transit parameters:

[TABLE]

where $x_{i}$ is the model, $\mu_{i}$ is the data, and $\sigma_{i}$ is the error in the data for the $i^{\rm th}$ point, respectively. The fitted parameters are radius ratio $k$, impact parameter $b$, stellar density $\rho_{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_{1}$ and $u_{2}$, where the central values of the coefficients were calculated with PyLDTK (Parviainen & Aigrain, 2015) which utilises 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 5 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 Table 4. The resulting corner plots are shown in the appendix (Figures 10 and 11). 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 3. 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 visualisation, the inset in figure 3 shows a magnification of the different campaigns. The coloured error bars mark the 1-$\sigma$ and 3-$\sigma$ uncertainties of the fitted transit times.

3.2 Stacked transit analysis

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 8 transit durations, centred on each transits 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 modelled 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, $\rho_{s}$, and the triangle sampling (Kipping, 2010) of the quadratic stellar limb-darkening coefficients $u_{1}$ and $u_{2}$, where uninformative priors was used. The orbital periods or the planets are kept fixed.

The MCMC method was implemented using the Python package emcee (Foreman-Mackey et al., 2013) for Bayesian parameter estimation. A $\chi^{2}$ 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 $\times$ 10000 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 $\chi^{2}$ equals to 1. We then initiated 100 MCMC chains of $5\times 10^{4}$ steps to sample the posterior space. We checked for convergence and discarded the first 2000 steps, then adopted the median, $16^{\rm th}$, and $84^{\rm th}$ percentiles of the samples in the marginalized posterior distributions as the fitted values and their 1-$\sigma$ uncertainties. The results of the fitted transit parameters of K2-146 b and K2-146 c are presented in appendix Table 4. The best-fitted transit parameters of K2-146 b and K2-146 c obtained here are generally consistent within $\sim$1-$\sigma$ with those derived by PyTV .

4 TTV Analysis and results

We modeled the TTVs for the two planets using the TTVFast code (Deck et al., 2014), considering Newtonian gravitational interactions between the host star and the two planets alone. For each planet, we fitted planet-to-star mass ratio, orbital period $P$, eccentricity and argument of periastron parameterization ($\sqrt{e}\cos\omega$ and $\sqrt{e}\sin\omega$, so that the uniform priors on these parameters correspond to the priors flat in $e$ and $\omega$), and time $t_{0}$ of inferior conjunction closest to the dynamical epoch $t_{\rm epoch}(\mathrm{BJD}-2454833)=3467.8$. The elements are osculating Jacobi elements defined at $t_{\rm epoch}$, and the time $t_{0}$ is related to the time of periastron passage $\tau$ via $2\pi(t_{0}-\tau)/P=E_{0}-e\sin E_{0}$, where $E_{0}=2\arctan\left[\sqrt{{1-e}\over{1+e}}\tan\left({\pi\over 4}-{\omega\over 2}\right)\right]$. Considering that both planets are transiting in the C18 data, which are close to the dynamical epoch, the inclination and longitude of ascending node at the epoch were fixed to be $\pi/2$ and [math], respectively. The likelihood was defined using the usual $\chi^{2}$ as $\exp(-\chi^{2}/2)$, where we directly adopted the errors from PyTV because the scatter in the data around the best model was found to be consistent with the assigned values. We adopted uniform priors for all these parameters and used emcee (Foreman-Mackey et al., 2013) to sample from their posterior distribution.

Figure 4 shows the TTV models generated with 20 sets of parameters randomly drawn from the posterior distribution. Table 2 summarizes the median and 68% credible interval of the marginal posterior distribution for each parameter: the upper parts show the fitted parameters, and the lower part shows the derived parameters. The planet-to-star mass ratios combined with the host star mass yield planetary masses of $5.6\pm 0.7\,M_{\oplus}$ for K2-146 b and $7.1\pm 0.9\,M_{\oplus}$ for K2-146 c. Moderate eccentricities are favored for both planets, but they show a strong negative correlation and one of the planets may have a nearly circular orbit. The $99.7\%$ upper limit for the eccentricity is $0.3$ for both planets. The observed difference in the arguments of periastron is consistent with anti-alignment of the apses. The implications of these features will be discussed in Section 6.

5 Transit modelling

5.1 Updating transit parameters

The stacked transit light curves of K2-146 b and c were re-analyzed incorporating the eccentricity information from the TTV analysis. The transits were modelled with TLCM (Transit and Light Curve Modeller; Csizmadia 2019, MNRAS, under review), a software tool for joint radial velocity and transit light curve fit, or transit light curve fit only. It utilizes the Mandel & Agol (2002) and Eastman et al. (2013) subroutines to calculate the transit light curve shapes for every time 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 Genetic Algorithm, refined by Simulated Annealing (Geem, 2001). The final parameter estimation is done by using several chains of MCMC with at least $10^{5}$ steps. The median and the width of the chains will define the finally adopted solutions and its 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_{0}$, the scaled semi-major 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_{1}+u_{2}$, and $u_{-}=u_{1}-u_{2}$. The orbital periods, $P$, of the planets were kept fixed. In addition, we used the eccentricity and argument of periastron derived from Section 4 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 model of the inner and outer planets are shown in Figure 5.

5.2 Planet impact parameters

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 semi-major 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\pm 0.11$, $b=0.25\pm 0.10$, and $b=0.41\pm 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-$\sigma$ 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\pm 0.10$ and $b=0.92\pm 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 $\sim 1000\,\mathrm{days}$ separation between C05 and C16.

6 Discussion

6.1 Dynamic stability of the K2-146 system

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.

A Hill-sphere is the region where the planet’s gravity is dominating over the central star. If the Hill-spheres overlap then there is a chance to collide or to remove one of the planets from the system. Using the orbital parameters reported in Table 2, we find that planets $b$ and $c$ have Hill-sphere radii of 2.1% and 2.2% of their semi-major axes which suggests that the planets are likely stable.

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 employed 5105 integrations, each with an integration period of 2 Myr. An initial step size of 2.7 d 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 years. For each integration, the orbital periods ($P_{b}$ and $P_{c}$) were chosen from a Gaussian-distribution with the center and 1$-\sigma$ reported by PyTV in Section 3. 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 co-planar with fixed inclinations and longitudes of the ascending node where $i_{b},i_{c}=90^{\circ}$ and $\Omega_{b},\Omega_{c}=0^{\circ}$. The planetary masses of K2-146 b and K2-146 c ($M_{b}$ and $M_{c}$ respectively) were at first randomly drawn from a half-normal distribution with a width of $1000$ $\rm M_{\oplus}$. We further increased the number of data points in the lower mass regions by uniformly sampling the 0-30 $\rm M_{\oplus}$ range (the 0-3 $\rm M_{\oplus}$ was even more densely sampled). All samples were merged in the presented final result. The arguments of periastron ($\omega_{b}$ and $\omega_{c}$) and mean anomalies of the planets were drawn from a uniform sample where $0^{\circ}<\omega_{b},\omega_{c}<360^{\circ}$ and $0^{\circ}<\mathcal{M}_{b},\mathcal{M}_{c}<360^{\circ}$. In each simulation, the system becomes unstable when (1) the planets collides with one another or with the central star; or (2) the planets are ejected from the system, i.e. $e_{b},e_{c}>1$ and/or $a>10$ au.

Our results presented in Figure 6 shows that a significant fraction of systems are dynamically stable for 2 Myr for planets with masses up to $\sim$30 $\rm M_{\oplus}$. The stable and unstable configurations are denoted by cyan and red circles, respectively. We calculated the fraction of stable orbits within each grid with widths of 5 $\rm M_{\oplus}$. We found that over 90% of the systems remained stable for 2 Myr if the mass ratio of K2-146 b and K2-146 c is close to unity. The upper mass limits of the planets obtained from dynamical constraint is consistent with the planet masses derived from TTVs in Section 4. The mass ratio of the planets determined from the TTV analysis is $0.783\pm 0.006$. The system is thus very likely to be stable on a timescale of at least 2 Myr.

6.2 Orbital resonance of the sub-Neptunes

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

The resonance arguments, $\Theta_{1}$ and $\Theta_{2}$, for a pair of planets in a 3:2 MMR are defined as

[TABLE]

so the difference is only in the last term. The mean motion resonance is defined as

[TABLE]

therefore in our case $p=2$ and $q=1$. 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+\Omega$ (e.g. Murray & Dermott 1999). The resonant angles measure the angle between the two planets at the conjunction point. If any resonant angle librates rather than circulates, then the planets are in mean motion resonance.

We performed twenty thousand 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. Stellar mass in Table 1 was used to convert the mass ratios into planetary masses. Then the orbits were numerically integrated for 2000 years. The values of the relevant parameters were saved for every 36 days of the integration, and 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 7 shows the histograms of the maximum half-amplitude of the resonant angles. The largest peak at around $\Theta_{1}/2=150$∘ is due to libration, so planet c is librating with a half-amplitude of $\sim 150$∘ around the conjunction points at the time of predicted conjunction times (every third orbital cycle). The smaller, narrower peak at $\Theta_{1}/2=180$∘ corresponds to either horseshoe-shaped orbits or $\Theta_{1}$ circulates and then $\Theta_{2}$ librates. We found approximately 77% of simulations show a libration amplitude of around 150*∘*  which implies there is at least a 77% chance that the system is trapped in 3:2 MMR instead of a random commensurability. Figure 8 shows the simulation which gave the smallest libration angle.

The width of the peak 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 smaller accuracy in TTVs, it would be worthy to try to observe it with high cadence with TESS (Transiting Exoplanet Survey Satellite; Ricker et al., 2014) and PLATO (PLAnetary Transits and Oscillations of stars; Rauer et al., 2014). Alternatively, observations may also be obtained with the less precise ($\sim$10 minutes) timing values from CHEOPS (CHaracterising ExOPlanet Satellite; Broeg et al., 2013), because the precision could be compensated with many hours of TTV-values in this case.

6.3 Planet composition

Our light curve analysis in section 5 and TTV analysis in section 4 reveal that K2-146 b has a mass and radius of $2.25\pm 0.10$ $\rm R_{\oplus}$ and $5.6\pm 0.7$ $\rm M_{\oplus}$, respectively, corresponding to a bulk density of $2.702\pm 0.494$ $\rm g\ cm^{-3}$. The mass and radius of the outer planet K2-146 c are $2.59_{-0.39}^{+1.81}$ $\rm R_{\oplus}$ and $7.1\pm 0.9$ $\rm M_{\oplus}$, respectively, which correspond to a bulk density of $2.246_{-1.846}^{+1.883}$ $\rm g\ cm^{-3}$. K2-146 b and K2 -146 c are orbiting at a distance of 0.0248 AU and 0.0327 AU, respectively. We assumed the planets have an albedo of 0, and a re-radiation factor of 1/4, where atmospheric circulation redistributes the energy around the planetary atmosphere, then re-radiate 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 9 shows a mass-radius plot of known planets with $M_{p}<30$ $\rm 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 colour bar. The masses and radii of K2-146 b and c are consistent with cases of a 100% H2O interior, a water-rich core with the addition of H2O-gaseous atmosphere, or an Earth-like rocky interior with a small fraction of H2-envelope. The large radius uncertainty of K2-146 c means that it could have a more massive H2-envelope.

6.4 Formation and evolution

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 H2 envelope – it was previously proposed that in situ formation of mini Neptunes and super Earths is possible if 50-100$\rm 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 towards the host star in a relatively short timescale, preventing in situ formation. Furthermore, in situ formation is unlikely to produce 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 migrated inwards towards their current locations.

2. Evolution of resonance arguments for the planet pair suggests that K2-146 b and c are likely trapped in 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 in 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 towards 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 would push the planets towards each other. When the planet orbital periods approaches a commensurability, dynamical interactions that follow cause planets to migrate collectively inwards while preserving period commensurability (Snellgrove et al., 2001). Under favourable 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 $0.14\lesssim e\lesssim 0.16$, and are apsidally anti-aligned, i.e. $\Delta\omega=\omega_{b}-\omega_{c}\approx 180$∘ – convergent migration could have played a role in the observed eccentricity in K2-146 planet pair. After the planets are captured in resonance, the planet pair migrates inwards 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 implies that the migration process must be fast enough in order to minimize the damping efficiency. The misaligned apsides could have been the result of such 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 is locked in an anti-aligned state (Lee & Peale, 2002). Due to close proximity of the two planets, such 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 gives 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 (cf. 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)).

7 Summary and conclusion

The strong TTV detected in the mini-Neptune K2-146 b suggested the presence of an additional body in the system. Further photometric observations from K2 revealed an additional mini-Neptune K2-146 c orbiting at a 4-day period, forming a 3:2 mean-motion resonance with the inner planet. The long observation baseline allowed precise determination of the planet masses via TTV. This demonstrates the importance of follow-up transit observations for parameter and dynamical constraints of a TTV system.

N-body simulations of K2-146 performed in this work provided a glimpse into the possible stability and resonance configuration of the planet pair. We found that the planets are probably captured into a 3:2 MMR during migration, and that their current orbital configuration can be dynamically stable for at least 2 Myr. Furthermore, the change in the impact parameter of the outer planet suggests some orbital plane precession, resulting in the displacement of the chord of transit and hence the change in transit depth and duration of K2-146 c. This effect can be further investigated using both TTVs and TDVs to constrain the orbit precession rate, and mutual inclination in the system. Further observations with TESS, and in the future PLATO can also provide a better precision in the transit times measurement. A detailed migration model would be valuable to study different precession rates leading to a stable orbital configuration.

Small planets around M dwarfs are frequently found in multi-planet systems. In fact, occurrence studies suggests that there are typically around 2.5 small planets (R_{p}<4$$\rm R_{\oplus}) per M dwarfs with periods shorter than 200 days (Dressing & Charbonneau, 2015). Only a small handful of planets, are known around M dwarfs which have planetary radii and masses below the radius and mass of Neptune (e.g. TRAPPIST-1 system; Gillon et al. 2017, LHS-1140 system; Dittmann et al. 2017; Ment et al. 2019, L 98-59 System; Kostov et al. 2019, Gl 357 system; Luque et al. 2019). The These are laboratories to test planet formation theories and dynamical evolutions, providing clues to the processes involved in building multi-planet systems containing the smallest possible planets.

The expected RV semi-amplitudes of K2-146 b and K2-146 c are $K_{b}=5.1\rm~{}m~{}s^{-1}$ and $K_{c}=5.7\rm~{}m~{}s^{-1}$. The faintness of the host star (J$=12.18$ mag) means that mass measurement by means of RV follow up is challenging for many currently available instruments. Recent precision RV instruments (e.g. Infrared Doppler instrument (IRD); Kotani et al. (2014)) on 8m class telescopes and future generations of high-resolution infra-red (IR) spectrographs such as CRIRES+ (Dorn et al., 2014) on the Very Large Telescope (VLT) are sensitive to cool, low mass M dwarfs which radiate mostly in the IR. RV measurements of the K2-146 system may be within reach.

Acknowledgement

We thank Mareike Godolt for helpful discussion on interpretation of planet compositions. KWFL, JK, SzCs, ME, SG, APH, MP and HR acknowledge support by DFG grants PA525/18-1, PA525/19-1, PA525/20-1, HA3279/12-1 and RA714/14-1 within the DFG Schwerpunkt SPP 1992, “Exploring the Diversity of Extrasolar Planets”. M.F., I.G., and C.M.P. gratefully acknowledge the support of the Swedish National Space Agency (DNR 163/16 and 174/18). This work is partly supported by JSPS KAKENHI Grant Numbers JP18H01265 and JP18H05439, and JST PRESTO Grant Number JPMJPR1775. O.B. acknowledges support from the UK Science and Technology Facilities Council (STFC) under grants ST/S000488/1 and ST/R004846/1. SM acknowledges support from the Spanish Ministry under the Ramon y Cajal fellowship number RYC-2015-17697. PK acknowledges support by GACR 17-01752J. SA and EK acknowledge support by the Danish Council for Independent Research, through a DFF Sapere Aude Starting Grantnr. 4181-00487B. SA, EK, and ML acknowledge support by the Stellar Astrophysics Centre which funding is provided by The Danish National Research Foundation (Grant agreement no.: DNRF106). This work was supported by the KESPRINT collaboration, an international consortium devoted to the characterization and research of exoplanets discovered with space-based missions. www.kesprint.science This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of NASA’s Astrophysics Data System Bibliographic Services, the ArXiv preprint service hosted by Cornell University.

Appendix A Supplementary figures