TESS and CHEOPS discover two warm sub-Neptunes transiting the bright K-dwarf HD 15906 Open Access

ABSTRACT

1 INTRODUCTION

Exoplanet population studies have shown that small planets between the size of Earth and Neptune (the so-called super-Earths and sub-Neptunes) are the most ubiquitous in our Galaxy (Fressin et al. 2013; Kunimoto & Matthews 2020). However, there is a statistically significant drop in the occurrence rate of close-in planets (orbital period ≲ 100 d) with radii between 1.5 and 2.0 R_⊕ (Fulton et al. 2017; Fulton & Petigura 2018; Van Eylen et al. 2018). One theory is that this radius gap represents a transition between predominantly rocky planets and planets with extended H/He envelopes. There are several possible explanations for how this could arise, including gas-poor formation (Lee, Chiang & Ormel 2014; Lee & Chiang 2016; Lopez & Rice 2018; Lee, Karalis & Thorngren 2022), core-powered mass-loss (Ginzburg, Schlichting & Sari 2018; Gupta & Schlichting 2019, 2020) and photoevaporation (Owen & Wu 2013, 2017; Lopez & Rice 2018). More recently, Luque & Pallé (2022) studied small planets transiting M-dwarfs and found that the radius gap might actually be a density gap separating rocky and water-rich planets. To test these theories we need small, well-characterized planets spanning a range of equilibrium temperatures, T_eq.

Warm (defined in this paper as T_eq ≲ 700 K) sub-Neptunes transiting bright stars are particularly interesting targets for detailed characterization studies. These planets are amenable to observations to, for example, precisely measure their radii and masses and probe their atmospheres (e.g. Kreidberg et al. 2014; Benneke et al. 2019; Tsiaras et al. 2019; Delrez et al. 2021; Scarsdale et al. 2021; Orell-Miquel et al. 2022; Wilson et al. 2022). From measurements of a planet’s mass and radius, the bulk density can be calculated and its internal composition inferred. This can help distinguish between different formation mechanisms for small planets (Bean, Raymond & Owen 2021). Furthermore, since warm planets are less affected by radiation from their host star, they can retain their primordial atmospheres. Observations of these atmospheres and measurements of the carbon-to-oxygen ratio could therefore reveal their formation history (Öberg, Murray-Clay & Bergin 2011; Madhusudhan, Amin & Kennedy 2014). Multiplanet systems are especially powerful because they allow us to study planets that formed from a common protoplanetary disc, leading to additional constraints on formation and evolution models (e.g. Lissauer et al. 2011; Fang & Margot 2012; Weiss et al. 2018; Van Eylen et al. 2019; Weiss et al. 2022).

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) is an all-sky transit survey searching for exoplanets around some of the brightest and closest stars. Since its launch in 2018, it has discovered a plethora of planets orbiting bright stars, including many super-Earth and sub-Neptune planets (e.g. Gandolfi et al. 2018; Vanderburg et al. 2019; Plavchan et al. 2020; Teske et al. 2020; Leleu et al. 2021; Serrano et al. 2022a). However, due to the nature of its observing strategy, TESS is limited in its ability to discover long-period exoplanets. During its two-year primary mission, TESS observed the majority of the sky for ∼ 27 consecutive days. This means that planets with periods longer than ∼ 27 d, and some planets with periods between ∼ 13 and 27 d, would only have been observed to transit once, if at all. These single transit detections are known as ‘monotransits’ and their orbital periods are unknown, although the shape of the transit allows the period to be constrained (e.g. Wang et al. 2015; Osborn et al. 2016). In its extended mission, TESS reobserved the sky approximately two years after the first observation and, as predicted by simulations (Cooke, Pollacco & Bayliss 2019; Cooke et al. 2020, 2021), a large fraction of primary mission monotransits were observed to transit a second time. The result was a sample of ‘duotransits’ – planetary candidates with two observed transits separated by a large gap, typically two years. From the two non-consecutive transits, the period of the planet remains unknown, but there now exists a discrete set of allowed period aliases. These aliases can be calculated according to P_n = T_diff / n, where T_diff is the time between the two transit events and n ∈ {1, 2,..., n_max}. The maximum value, n_max, is dictated by the non-detection of a third transit in the TESS data.

Both monotransits and duotransits are the observational signatures of long-period planets (P ≳ 20 d). However, follow-up photometric or spectroscopic observations are required to recover their true periods. The follow-up of monotransits requires a blind survey approach (e.g. Gill et al. 2020; Villanueva et al. 2021; Ulmer-Moll et al. 2022), whereas the period aliases of a duotransit allow more targeted follow-up observations (e.g. Grieves et al. 2022; Ulmer-Moll et al. 2022). So far, the majority of these follow-up efforts have focused on giant planets, partly because their deeper transits facilitate ground-based observations. It’s vital that we also pursue follow-up of shallow duotransits to expand the sample of small, long-period planets, including warm sub-Neptunes.

The CHaracterising ExOPlanet Satellite (CHEOPS; Benz et al. 2021) is an ESA mission dedicated to the follow-up of known exoplanets. The effective aperture diameter of CHEOPS (∼ 30 cm) is about three times larger than that of TESS (∼ 10 cm), allowing it to achieve a higher per-transit signal-to-noise ratio (SNR; e.g. Bonfanti et al. 2021). Furthermore, CHEOPS performs targeted photometric observations to observe multiple transits of a planet without the need for continuous monitoring. CHEOPS is therefore very well-suited to the follow-up of small, long-period planets from TESS. We have a dedicated CHEOPS Guaranteed Time Observing (GTO) program to recover the periods of TESS duotransits, focusing on small planets that cannot be observed from the ground. We select most of our targets from the TESS Objects of Interest (TOI) Catalog (Guerrero et al. 2021) and from our specialized duotransit pipeline (Tuson & Queloz 2022). Through our CHEOPS programme, we have recovered the periods of two duotransits in the TOI 2076 system (Osborn et al. 2022), one duotransit in the HIP 9618 system (Osborn et al. 2023), one duotransit in the TOI 5678 system (Ulmer-Moll et al. 2023) and one duotransit in the HD 22946 system (Garai et al. 2023).

In this paper, we report the discovery of two warm sub-Neptunes transiting the bright (G = 9.5 mag) K-dwarf HD 15906 (TOI 461, TIC 4646810). This paper is organized as follows. In Section 2, we provide details of the photometric and spectroscopic observations used in our analyses. In Section 3, we describe our characterization of the host star and in Section 4 we describe the analyses of the system. Section 5 presents the results of our analyses and in Section 6 we validate the two planets. Finally, in Section 7, we present a discussion of our findings and outlook for future observations.

2 OBSERVATIONS

2.1 TESS photometry

HD 15906 was observed by TESS (camera 1, CCD 1) at two-minute cadence in sector 4 (2018 October 18 to November 15) and sector 31 (2020 October 21 to November 19). During both sectors, the instrument suffered from operational anomalies causing interruptions in data collection. In sector 4, no data were collected between 1418.5 and 1421.2 (BJD – 2457000) due to an instrument shutdown and sector 31 ended ∼ 2 d earlier than scheduled due to a star tracker anomaly. No more TESS observations are scheduled before the end of Cycle 6 (2024 October 1).

The TESS observations were reduced and analysed by the Science Processing Operations Center (SPOC; Jenkins et al. 2010a, 2016) at the NASA Ames Research Center. We downloaded the light curve files, created by SPOC pipeline version 5.0.20–20201120, from the Mikulski Archive for Space Telescopes (MAST) portal.¹ These files include a Simple Aperture Photometry (SAP; Twicken et al. 2010; Morris et al. 2020) light curve and a Presearch Data Conditioning Simple Aperture Photometry (PDCSAP; Smith et al. 2012; Stumpe et al. 2012, 2014) light curve that has been corrected for instrumental systematics. For our analysis, we used the PDCSAP light curves. Following the advice in the TESS Archive Manual², we rejected all data points of lesser quality using the binary digits 1, 2, 3, 4, 5, 6, 8, 10, 13, and 15. We then rejected outliers from the light curve by calculating the mean absolute deviation (MAD) of the data from the median smoothed light curve and rejecting data greater than 5 × MAD away from the smoothed data set. We repeated this process until no more outliers remained and the resulting TESS light curve is shown in Fig. 1.

Figure 1.

TESS PDCSAP light curve from sector 4 (left-hand panel) and sector 31 (right-hand panel). The 2 min cadence data (grey) has been binned to 120 min (black squares) to guide the eye. The red line is the median model from the global photometric fit, described in Section 4.2, and the red shaded region (difficult to see on this scale) is the 1σ uncertainty on the model. The blue and pink markers indicate the mid-transit times of the inner and outer planet, respectively. The lower panels show the residuals of the median model.

Open in new tabDownload slide

From the sector 4 data alone, the transiting planet search (TPS; Jenkins 2002; Jenkins et al. 2010b, 2020) performed by the SPOC pipeline identified a single planet candidate. This planet candidate was announced as TOI 461.01 in 2019 February with an epoch of 1416.3 (BJD – 2457000) and a period of 14.5 d. When the sector 31 data became available, we performed a by-eye search of the light curve and realized that TOI 461.01 was actually a combination of two planetary signals. There was one multitransiting planet candidate, with an epoch of 1416.3 (BJD – 2457000) and a period of 10.9 d, and one duotransit – a planet candidate with one transit in sector 4 and one transit in sector 31, separated by T_diff ∼ 733.8 d. When a multisector TPS was performed by SPOC in 2021 May, it correctly identified the multitransiting planet candidate and the ephemeris of TOI 461.01 was updated accordingly. This planet candidate passed all of the SPOC vetting tests (Twicken et al. 2018; Li et al. 2019), including the difference image centroid test, the odd–even depth test and the ghost diagnostic test, and the source of the transit signal was localized within 6.0 ± 4.2 arcsec of HD 15906. The duotransit did not receive a TOI designation.

The TESS data contains four transits of the inner planet candidate (TOI 461.01, hereafter called HD 15906 b) and two transits of the outer planet candidate (hereafter called HD 15906 c). From the TESS data alone, the orbital period of the outer planet candidate was ambiguous. There existed a discrete set of 36 allowed period aliases, in the range 20.4 – 733.8d (see Section 4.1), and follow-up observations were therefore required to recover the correct period.

2.2 CHEOPS photometry

To recover the period of the outer planet candidate, we observed HD 15906 through the CHEOPS GTO programme CH_PR110048 (‘Duos – Recovering long period duo-transiting planets’). Our observing strategy was informed by our analysis of the TESS data (see Section 4.1). We scheduled CHEOPS observations of the 13 highest probability period aliases (P < 31 d), giving highest priority to the four most probable period aliases (P < 22.5 d). The first and second CHEOPS visits did not reveal a transit and ruled out six period aliases in total. The third CHEOPS visit revealed a transit and uniquely confirmed a period of ∼ 21.6 d for HD 15906 c. A fourth CHEOPS visit, scheduled before the period had been confirmed, did not reveal a transit. We scheduled one additional observation of both HD 15906 b and c to improve radius precision and search for possible transit timing variations (TTVs). For all of our CHEOPS observations, we used an exposure time of 60 s with no on-board image stacking, resulting in a final light curve cadence of 60 s. A summary of our six CHEOPS observations is presented in Table 1.

Due to the fact CHEOPS is in a low-Earth orbit, with an orbital period ∼ 98.7 min, our observations suffer from interruptions caused by high levels of stray light (from the illuminated Earth limb), occultations of the target by the Earth and passage of the satellite through the South Atlantic Anomaly (SAA; Benz et al. 2021). These interruptions result in gaps in the CHEOPS light curves, reducing the observing efficiency (time spent collecting data divided by the duration of the visit). The efficiencies of our six visits are included in Table 1 and the inset of Fig. 2 shows examples of the light curve gaps.

Figure 2.

Background flux versus angle between the instrument’s line of sight (LOS) and the Earth limb for non-flagged data from all six CHEOPS observations. The red horizontal line represents our background cut of 10 000 e⁻pix⁻¹ and the red crosses correspond to the data removed from CHEOPS visit 2. Inset: CHEOPS light curve from visit 2. Only non-flagged data are plotted and the points shown as red crosses were removed by the background cut.

Open in new tabDownload slide

For each of our CHEOPS visits, sub-array images and light curves were produced by the Data Reduction Pipeline (drp 13.1.0; Hoyer et al. 2020). The sub-array images are circular, with a diameter of 200 pixels (∼ 200 arcsec), and are centred on the target star. They are calibrated and corrected for effects such as cosmic ray hits, smear trails caused by nearby stars, and variations in background flux. From these images, the drp uses aperture photometry to produce four light curves using circular apertures of different sizes. The DEFAULT, RINF, and RSUP apertures are pre-defined with radii 25, 22.5, and 30 pixels, respectively. The OPTIMAL aperture is selected per visit to minimize the effect of instrumental noise and contamination from nearby stars. We downloaded the CHEOPS sub-array images and drp light curves from the Data & Analysis Center for Exoplanets (DACE³; Buchschacher et al. 2015). Alongside the time, flux, and flux error, the drp light curves include a set of detrending vectors that can be used to model instrumental trends in the light curve. This includes the background flux, the smearing and contamination from nearby stars, the x and y centroid position of the target star, and the roll angle of the satellite. CHEOPS rolls around its pointing direction once per orbit, to maintain thermal stability, and every data point has an associated roll angle between 0 and 360 degrees.

We also extracted our own light curves from the CHEOPS sub-array images using point-spread function (PSF) photometry. This technique is complementary to the aperture photometry performed by the drp. We used the PSF Imagette Photometric Extraction (pipe) package⁴ (see description in Deline et al. 2022), which was developed specifically for CHEOPS data. pipe photometry is less sensitive to contamination from nearby stars and the effects of smear trails are removed before extracting the flux (Serrano et al. 2022b). The pipe light curves contain the time, flux, and flux error, as well as the same detrending vectors as the drp light curves, with the exception of smearing and contamination.

We performed preliminary transit fits of the drp and pipe light curves using pycheops⁵ (Maxted et al. 2021) and found that the planet parameters obtained in each case were fully compatible. We then compared the photometric precision of the drp and pipe light curves for each CHEOPS visit. First, we performed iterative outlier clipping as described in Section 2.1. Then, we calculated the MAD of each clipped light curve, see Table 2. We found that for four of the six visits, including all three transit observations, the pipe light curve had the lowest MAD. In the other two visits, the MAD of the pipe light curve was comparable to the lowest value. We therefore chose to use the pipe photometry for our analysis.

To prepare the pipe light curves for our analysis, we performed a series of cuts to the data. First, we rejected all flagged data. pipe assigns flags to data of lesser quality, for example due to outliers in centroid position or a large number of bad pixels in the frame. Next, we performed a cut to remove data with high background flux. Some of the CHEOPS light curves showed sharp spikes in the target’s flux immediately before and/or after the data gaps (see an example in the inset of Fig. 2). These spikes coincide with the target star approaching the illuminated Earth limb, causing high levels of scattered light and an increase in the background flux. This can be seen in Fig. 2, where we have plotted the background flux against the angle between the instrument’s line of sight (LOS) and the Earth limb for all six CHEOPS visits. Notice that not all of the observations with a small angle have a high background flux; it is only when the star approaches the Earth’s day side that there is a significant increase in scattered light. We removed all data with background flux > 10 000 e⁻pix⁻¹ because this adequately reduced the spikes in the light curves whilst retaining as much data as possible. After the background cut, we removed remaining outliers from the light curves using the same iterative MAD clipping described in Section 2.1. In total, these three cuts rejected 42/346 (∼ 12 per cent), 33/358 (∼ 9 per cent), 36/356 (∼ 10 per cent), 47/353 (∼ 13 per cent), 43/476 (∼ 9 per cent), and 31/375 (∼ 8 per cent) data points from each respective CHEOPS visit.

Following these steps, the pipe photometry still contained trends correlated with instrumental parameters such as background flux, centroid position, and roll angle. Rather than pre-detrending the data, we chose to fit a joint transit and detrending model, see Section 4.2.

2.3 LCOGT photometry

We conducted ground-based photometric follow-up observations of HD 15906 as part of the TESS Follow-up Observing Program⁶ (TFOP; Collins 2019) Sub Group 1.

We used the TESS Transit Finder, a customized version of the tapir software package (Jensen 2013), to schedule our transit observations. We observed full predicted transit windows of HD 15906 b in Pan-STARRS z-short band using the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network nodes at Siding Spring Observatory and McDonald Observatory on 2021 August 27 and 2021 November 1, respectively. See Table 3 for a summary of these observations. The 1.0 m telescopes are equipped with 4096 × 4096 SINISTRO cameras having an image scale of 0.389 arcsec pix⁻¹, resulting in a 26 arcmin × 26 arcmin field of view. We used an exposure time of 30 s and, with the full frame readout time of ∼ 30 s, the final image cadence was ∼ 60 s. The images were calibrated with the standard LCOGT banzai pipeline (McCully et al. 2018). The telescopes were intentionally defocused in an attempt to improve photometric precision, resulting in a typical HD 15906 full-width half-maximum (FWHM) of 6.5 arcsec. Differential photometric data were extracted using AstroImageJ (Collins et al. 2017). We used a circular photometric aperture with radius 9.3 arcsec to exclude all flux from the nearest known Gaia Data Release 3 stars (Gaia DR3; Gaia Collaboration 2016, 2022). A transit-like event was detected in both LCOGT light curves and they were included in the analysis described in Section 4.2.

2.4 WASP photometry

HD 15906 was observed 38 740 times by the Wide Angle Search for Planets at the South African Astronomical Observatory (WASP-South; Pollacco et al. 2006) between 2008 August 19 and 2014 December 19. The photometry was extracted and detrended for systematic effects following the methods described in Collier Cameron et al. (2006). Based upon a visual inspection of the light curve, we removed data with a normalized flux greater than 1.07 or less than 0.93 and we removed data with a relative flux error greater than 0.03. These cuts removed 5 231/38 740 (∼ 14 per cent) data points and the resulting light curve is shown in Fig. 3. With an average flux error of ∼ 9 ppt, we do not detect the transits of HD 15906 b or c in the WASP data. Furthermore, there were no additional transits detected in the light curve. Thanks to the long baseline, the WASP photometry is used to estimate the stellar rotation period (see Section 3.3).

Figure 3.

Archival WASP photometry. Upper: Normalized WASP light curve spanning more than 6 yr. Lower: GLS periodogram of the WASP light curve. The strongest peaks are in the range 25–30 d (red highlight), followed by 13–15 d (green highlight).

Open in new tabDownload slide

2.5 HARPS spectroscopy

The High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003) is a high-resolution (R = 115 000) fibre-fed spectrograph installed on the 3.6 m telescope at the European Southern Observatory (ESO) in La Silla, Chile. It has been operational since 2003 and the optical fibres were upgraded in 2015, leading to an offset in the measured radial velocities (RVs) (Lo Curto et al. 2015).

HARPS observed HD 15906 18 times between 2003 November 3 and 2018 February 9. There were 15 observations taken before the fibre upgrade and 3 observations taken after the upgrade. The exposure times of the observations ranged from 358 to 900 s and the average SNR at 550 nm was 53.7. The data spans ∼ 5212 d, with an average separation of ∼ 307 d between each observation. The HARPS spectra are publicly available on the ESO Science Archive Facility.

For our analysis of the HD 15906 system, we used the RVs presented in Trifonov et al. (2020). Specifically, we used the columns ‘RV_mlc_nzp’ and ‘e_RV_mlc_nzp’ for the RV and RV error, respectively. These RVs were extracted by the SpEctrum Radial Velocity AnaLyser (serval; Zechmeister et al. 2018) pipeline, where the extraction was done independently before and after the fibre upgrade and a correction was made for the nightly zero-point. The data have a root mean square (RMS) of 10.70 ms⁻¹ and the average RV uncertainty is 1.51 ms⁻¹. We present these RVs in Table 4.

2.6 FIES spectroscopy

As part of TFOP, we observed HD 15906 seven times using the FIbre-fed Échelle Spectrograph (FIES; Telting et al. 2014) at the Nordic Optical Telescope (NOT; Djupvik & Andersen 2010) between 2019 September 15 and October 12. For each observation, we used the high-resolution fibre (R ∼ 67 000) and an exposure time of 1800 s. We extracted the spectra and derived multi-order RVs following Buchhave et al. (2010). The SNR per resolution element at 550 nm ranges between 20 and 105 with a median of 97. The RMS of the RV data is 7.88 ms⁻¹ and the average uncertainty is 6.08 ms⁻¹. These FIES RVs are included in Table 4.

3 STELLAR CHARACTERIZATION

3.1 Atmospheric properties

As described in Section 2.5, HD 15906 was observed by HARPS 18 times between 2003 and 2018, with 15 observations made before the 2015 fibre upgrade. We retrieved the 15 pre-upgrade HARPS spectra from the ESO Science Archive Facility and co-added them to create a single master spectrum. This was used to perform the following spectroscopic analyses.

We performed an equivalent width (EW) analysis using ares + moog to derive the stellar atmospheric parameters (T_eff, log g, microturbulence, [Fe/H]). We followed the same methodology described in Santos et al. (2013); Sousa (2014); Sousa et al. (2021). We used the latest version of ares⁷ (Sousa et al. 2007, 2015) to measure the EWs of the iron lines in the master HARPS spectrum. We used a minimization process to find ionization and excitation equilibrium and converge to the best set of spectroscopic parameters. The iron abundances were computed using a grid of Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code moog (Sneden 1973). We also derived a more accurate trigonometric surface gravity using recent Gaia data following the same procedure as described in Sousa et al. (2021). The quoted errors for T_eff, log g, and [Fe/H] are ‘accuracy’ errors, that is they have been corrected for systematics following the discussion presented in Section 3.1 of Sousa et al. (2011). The final spectroscopic parameters and their errors are included in Table 5 and we find that HD 15906 is a K-dwarf.

We also performed an independent spectral synthesis with sme (Spectroscopy Made Easy; Valenti & Piskunov 1996; Piskunov & Valenti 2017) version 5.2.2⁸. A detailed description of the modelling can be found in Persson et al. (2018). We used the ATLAS12 stellar atmosphere grid (Kurucz 2013) and atomic and molecular line data from VALD⁹ (Vienna Atomic Line Database; Ryabchikova et al. 2015). The macro- and microturbulent velocities were held fixed to 1.5 and 0.5 kms⁻¹, respectively. The resulting T_eff, log g, and abundances were in excellent agreement with the ares + moog analysis. We additionally derived the projected rotational velocity, vsin i_⋆ = 2.7 ± 0.7 kms⁻¹.

3.2 Stellar mass and radius

We determined the stellar radius, R_⋆, of HD 15906 from the stellar angular diameter and the offset corrected Gaia DR3 parallax (Lindegren et al. 2021) using a Markov-Chain Monte Carlo Infrared Flux Method (MCMC IRFM; Blackwell & Shallis 1977; Schanche et al. 2020). We used the stellar spectral parameters as priors to construct model spectral energy distributions (SEDs) using atmospheric models from stellar catalogues. From this, we derived the stellar bolometric flux and angular diameter by comparing synthetic photometry, computed by convolving the model SEDs over broadband bandpasses of interest, to the observed data taken from the most recent data releases for the following bandpasses; Gaia G, G_BP, and G_RP, Two Micron All-Sky Survey (2MASS) J, H, and K, and Wide-field Infrared Survey Explorer (WISE) W1 and W2 (Skrutskie et al. 2006; Wright et al. 2010; Gaia Collaboration 2022). To account for systematic model uncertainties in our stellar radius error, we used stellar atmospheric models taken from a range of ATLAS catalogues (Kurucz 1993; Castelli & Kurucz 2003) and combined them in a Bayesian modelling averaging framework. Within the MCMC IRFM we attenuated the SED to correct for potential extinction and report the determined E(B–V) in Table 5. We combined the retrieved angular diameter with the offset-corrected Gaia DR3 parallax and found R_⋆ = 0.762 ± 0.005 R_⊙.

We then determined the stellar mass, M_⋆, by inputting T_eff, [Fe/H], and R_⋆ into two different stellar evolutionary models, parsec¹⁰ v1.2S (PAdova and TRieste Stellar Evolutionary Code; Marigo et al. 2017) and cles (Code Liégeois d’Évolution Stellaire; Scuflaire et al. 2008). We employed the isochrone placement algorithm (Bonfanti et al. 2015; Bonfanti, Ortolani & Nascimbeni 2016) to interpolate the input parameters within pre-computed grids of parsec isochrones and tracks and we retrieved a first estimate of the stellar mass, M_{⋆, PD} = 0.772 ± 0.037 M_⊙. A second estimate was computed through the cles code, which builds the best-fitting evolutionary track of the star by applying the Levenberg–Marquadt minimization scheme (e.g. Salmon et al. 2021) and we found M_{⋆, LG} = 0.797 ± 0.014 M_⊙. To account for model-related uncertainties, we added in quadrature an uncertainty of 4 per cent to the mass estimates obtained from each set of models (see Bonfanti et al. 2021). We note that the two outcomes are well within 1σ. We also checked their mutual consistency through the χ²-based criterion broadly presented in Bonfanti et al. (2021) and obtained a p-value = 0.49, which is greater than the normally adopted significance level of 0.05, as expected. For each mass estimate, we built the corresponding Gaussian probability density function, as described in Bonfanti et al. (2021), and we combined them to obtain a final mass value of M_⋆ = 0.790|$_{-0.036}^{+0.020} {\rm M}_{\odot }$|⁠, as presented in Table 5.

3.3 Stellar age

The isochrone fitting described in Section 3.2 also provided an estimate of the stellar age. However, the stellar mass is sufficiently low that the slow evolutionary speed of the star along its tracks led to an uninformative age of 6.8|$_{-6.3}^{+6.9}$| Gyr. To try and constrain the stellar age more precisely, we used gyrochronology, empirical |$\log R^{\prime }_{\mathrm{HK}}$| relations, and kinematics.

For the gyrochronology, we first estimate the stellar rotation period, P_rot. The TESS photometry (Fig. 1) shows flux modulation, likely caused by stellar activity, that can be used to do this. We conducted a generalized Lomb–Scargle (GLS; Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009) analysis on the TESS SAP and PDCSAP photometry and found strong peaks at 11–12 d and 25–27 d. However, this analysis is adversely affected by the short ∼ 27 d baseline of the TESS light curves. The archival WASP photometry has a much longer baseline that can be used to derive an independent estimate of the stellar rotation period. We performed a GLS analysis on the WASP light curve, the results of which are shown in Fig. 3. The strongest peaks are in the range 25–30 d, with the maximum power at 26.6 d corresponding to a best-fitting photometric amplitude of ∼ 4 ppt. The next strongest peaks are in the range 13–15 d, with a maximum power at 13.7 d and an amplitude of ∼ 3 ppt. This shorter rotation period is supported by our value of vsin i_⋆. Assuming sin i_⋆ = 1 and using the stellar radius in Table 5 leads to an upper limit of the rotation period, P_rot = 14.3 ± 3.7 d. Finally, from a GLS analysis of the HARPS and FIES RVs (see Section 4.4), we found that the peak power was at 12.27 d with a false alarm probability (FAP) of less than 1 per cent. It’s possible that this corresponds to the stellar rotation period, however, due to the very sparse sampling of the RVs, this value is unreliable. The stellar rotation period remains somewhat ambiguous, but the evidence favours a value in the range 11–15 d. Using the gyrochronological relations of Barnes (2007) and (B–V) from Table 5, these P_rot values yield a stellar age in the range 0.29–0.52 Gyr. We note that the longer P_rot values (25–30 d) would translate to an age of 1.39–1.97 Gyr. However, more recent studies have shown that the relations of Barnes (2007) might lead to an incorrect age estimate for low-mass stars because they do not account for the stalling period during spin-down (e.g. Curtis et al. 2020). Based upon a sample of benchmark stellar clusters, a rotation period of 11–15 d for a star with a similar effective temperature as HD 15906 is consistent with an age up to ∼ 1 Gyr.

Next, we computed values of |$\log R^{\prime }_{\mathrm{HK}}$| from each of the 18 HARPS spectra using actin¹¹ (Gomes da Silva et al. 2018) to extract the Ca ii index and following the method described in Gomes da Silva et al. (2021) for the |$\log R^{\prime }_{\mathrm{HK}}$| calibration. We found an average value of −4.694 ± 0.065 and, using the empirical relations of Mamajek & Hillenbrand (2008), this converts into a stellar age of 1.9 ± 0.7 Gyr.

Finally, we computed the kinematic age using the method developed in Almeida-Fernandes & Rocha-Pinto (2018) and the Galactic UVW velocities that we determined from the Gaia DR3 proper motions, offset-corrected parallax (Lindegren et al. 2021), and stellar RV, using the method outlined in Johnson & Soderblom (1987). We found a stellar age of 1.9|$^{+6.0}_{-0.7}$| Gyr, favouring an older star.

In Fig. 4, we present a comparison of the age estimates derived by our various methods. The age estimates derived from |$\log R^{\prime }_{\mathrm{HK}}$| and kinematics are consistent and they are in agreement with the gyrochronological age implied by a rotation period of 25–30 d. The favoured rotation period of 11–15 d yields a much younger age, however we reiterate that gyrochronology is not necessarily accurate for low-mass stars. We conclude that the stellar age is ambiguous based on the current data.

Figure 4.

A comparison of stellar age estimates obtained from isochrone fitting, |$\log R^{\prime }_{\mathrm{HK}}$| relations, kinematics, and gyrochronology.

Open in new tabDownload slide

4 ANALYSIS

4.1 TESS only analysis

Before pursuing CHEOPS follow-up observations of HD 15906, we used monotools¹² (Osborn 2022) to perform an analysis of the TESS data. monotools is designed for the analysis of planets with unknown periods, including duotransits. It can be used to derive the allowed period aliases and their corresponding probabilities, crucial for scheduling follow-up observations.

We built a monotools model using the stellar parameters presented in Table 5, one periodic planet and one duotransit. We defined initial guesses for transit depth, duration, and mid-transit time for the two planets using a visual inspection of the TESS light curve. Since this is a multiplanet transiting system, we selected the eccentricity distribution from Van Eylen & Albrecht (2015). We also included a Gaussian Process (GP; Rasmussen & Williams 2006; Gibson 2014) with a simple harmonic oscillator (SHO) kernel from celerite (Foreman-Mackey et al. 2017) to model the correlated noise in the light curve. We sampled the posterior probability distribution using the No-U-Turn Sampler (NUTS; Hoffman & Gelman 2014), a variant of Hamiltonian Monte Carlo, implemented via pymc3 (Salvatier, Wiecki & C. 2016).

We found that the duotransit, HD 15906 c, had 36 possible period aliases, with a minimum value, P_min, of 20.384 d. The probability of each period alias is shown in Fig. 5. These results were used to schedule our CHEOPS follow-up observations, from which we successfully determined the true period of planet c to be ∼ 21.6 d (see Section 2.2).

Figure 5.

From the TESS data alone, HD 15906 c was a duotransit with 36 possible period aliases. This plot shows these aliases and their corresponding probabilities, derived using monotools. The true period, as determined by CHEOPS follow-up observations, had the highest probability.

Open in new tabDownload slide

4.2 Global photometric analysis

Once we had confirmed the true period of HD 15906 c with CHEOPS, we performed a joint fit of the TESS, CHEOPS, and LCOGT photometric data using juliet¹³ (Espinoza, Kossakowski & Brahm 2019). This package combines transit models from batman (Kreidberg 2015) with the option to include linear models and GPs to model instrumental noise and stellar variability. We created a model consisting of two transiting planets, using the following parametrization:

• Orbital period, P, and mid-transit time, T₀, for both planets. We set broad uniform priors on P and T₀ from a visual inspection of the TESS and CHEOPS light curves.

• Planet-to-star radius ratio, p = R_P/R_⋆, and impact parameter, b, for both planets. We set uniform priors to allow exploration of all physically plausible solutions.

• Eccentricity, e, and argument of periastron, ω, for both planets. We used the eccentricity prior from Van Eylen et al. (2019) for systems with multiple transiting planets – the positive half of a Gaussian with μ = 0 and σ = 0.083. We used a uniform prior for ω, covering the full range of possible values. We decided to fit for eccentricity, rather than assuming a circular orbit, to ensure that the uncertainties on the other fitted parameters were not underestimated. We note that we repeated our final global photometric fit assuming a circular orbit, with e fixed to zero and ω fixed to 90 degrees, and all of the fitted planet parameters were consistent within 1.2σ.

• Stellar density, ρ_⋆. Using Kepler’s third law, this can be combined with P to derive a value of a/R_⋆ for each planet. This is preferred to fitting for a/R_⋆ directly; not only does it reduce the number of fitted parameters, but it also ensures a consistent value of ρ_⋆. We defined a normal prior on ρ_⋆ using the values of R_⋆ and M_⋆ presented in Table 5.

• Quadratic limb darkening parameters, q1 and q2, for each instrument. We used the Kipping (2013) parametrization of the quadratic limb darkening law and defined normal priors on q1 and q2 for each instrument. The mean was computed by interpolating tables of quadratic limb darkening coefficients (Claret & Bloemen 2011; Claret 2017, 2021), based on the stellar parameters presented in Table 5, and a standard deviation of 0.1 was used in all cases.

In addition to the transit models, we used linear models to detrend CHEOPS and LCOGT against instrumental systematics (see Sections 4.2.1 and 4.2.2). We treated each CHEOPS and LCOGT observation independently for the sake of this detrending. We also included a GP to model the variability in the TESS light curve, see Section 4.2.3. For each instrument, we included a jitter term to account for white noise and a relative flux offset term. We fixed the dilution factor to 1 due to the lack of any bright contaminating sources (see Section 6.4). We used the dynesty package to sample the posterior probability of this model with static nested sampling, using 300 live points and stopping when the difference between the evidence and the estimated remaining evidence was less than 0.01 (Speagle 2020). For a full list of the parameters and priors used in our global fit see Appendix  A and for the results of our modelling see Section 5.1.

4.2.1 CHEOPS detrending

The CHEOPS light curves contain trends that are correlated with instrumental parameters such as background flux (bg) and centroid position (x, y). There are also periodic noise features that repeat once per CHEOPS orbit due to the satellite rolling around its pointing direction. Detrending the light curve against the sine or cosine of the roll angle (ϕ) can remove these periodic instrumental effects.

The CHEOPS light curves also include stellar variability. From the TESS LC we know that HD 15906 shows stellar variability (see Fig. 1). On the shorter time-scale of a CHEOPS visit (∼ 8.3 h), this stellar variability can be modelled with a linear trend in time (t).

We included linear models in our global fit to account for these instrumental trends and stellar variability. However, for each CHEOPS observation, it was important to only select the relevant detrending parameters. To do this we used the pycheops package (Maxted et al. 2021) and the method described in Swayne et al. (2021). Briefly, we defined 10 detrending parameters: x, y, t, bg, cos(ϕ), sin(ϕ), cos(2ϕ), sin(2ϕ), cos(3ϕ), and sin(3ϕ). For each CHEOPS visit, we took the clipped light curve (see Section 2.2) and did an initial fit of a transit model with no detrending. We defined broad uniform priors on the transit parameters based on a visual inspection of the TESS and CHEOPS data. We used the RMS of the residuals from this initial fit to define normal priors on the detrending parameters, with μ = 0 and σ = RMS. We added the 10 detrending parameters to the fit one-by-one, selecting the parameter with the lowest Bayes factor at each step. When there were no remaining parameters with Bayes factor < 1, we stopped adding detrending parameters. In order to remove strongly correlated parameters, if any of the selected detrending parameters had a Bayes factor > 1, we removed the parameter with the largest Bayes factor until no more parameters with Bayes factor > 1 remained. The selected detrending parameters for each CHEOPS visit are included in Table 1.

4.2.2 LCOGT detrending

We used AstroImageJ to select the relevant detrending vectors for each LCOGT observation by jointly fitting a transit model and linear combinations of zero, one, or two detrending parameters from the available detrending vectors: airmass, time, sky background, FWHM, x-centroid, y-centroid, total comparison star counts, humidity, and exposure time. The best zero, one, or two detrending vectors were retained if they reduced the Bayesian information criterion (BIC) for a fit by at least two per detrending parameter. We found that the airmass plus FWHM detrending pair provided the best improvement to the light curve fit for both LCOGT observations. We therefore included linear models for airmass and FWHM for each LCOGT observation in our global fit.

4.2.3 TESS detrending

The TESS light curves contain correlated noise, including stellar variability and residual instrumental systematics, that we model with a GP. We initially modelled sector 4 and sector 31 jointly, using a GP with an approximate Matérn-3/2 (M32) kernel implemented via celerite (Foreman-Mackey et al. 2017). Upon a visual inspection of the results from this fit, we noticed that the TESS residuals contained a sinusoidal-like trend. We ran a GLS analysis on the TESS residuals, treating the sector 4 and sector 31 data separately, and the resulting periodograms are presented in Fig. 6. We found a significant periodic signal in the TESS sector 4 residuals, with a period of 0.22004 d and a FAP of 6.0 x 10⁻¹¹. This signal is persistent throughout the whole of sector 4 and the best-fitting sinusoidal model has an amplitude of ∼ 57 ppm. There was no corresponding detection in the TESS sector 31, CHEOPS, or LCOGT residuals. The periodic signal is present in the TESS light curve itself, it was not introduced as a result of our detrending, and we discuss its origin in Section 7.1.

Figure 6.

GLS periodogram of the TESS residuals from a global photometric fit using a GP with a Matérn-3/2 kernel to jointly model TESS sector 4 and sector 31. Upper: GLS periodogram of sector 4 residuals. Lower: GLS periodogram of sector 31 residuals. The horizontal red line is the 1 per cent false alarm probability level in each case. The significant peak at 0.22004 d in the sector 4 residuals (highlighted in green) is not present in the sector 31 residuals.

Open in new tabDownload slide

A half-cycle of the periodic signal is a similar duration to the transits and it was therefore important to check if it was affecting the fitted planet parameters. We therefore performed two additional fits, changing only the TESS detrending to account for this periodic signal. We made a custom GP kernel by adding together the M32 and SHO kernels from celerite. The M32 kernel was intended to capture the long-term variability and the SHO kernel was used to capture the short-term quasi-sinusoidal noise. We defined a normal prior on the natural frequency of the SHO kernel, ω₀, using the peak and its width from the periodogram analysis. We performed one fit where we jointly modelled the sector 4 and sector 31 data with this kernel and we also performed a fit where we decoupled the sector 4 and sector 31 data. We used the M32 plus SHO kernel for sector 4 and the M32 kernel for sector 31, motivated by the fact we only detect the periodic trend in sector 4. After performing these two fits, we checked for periodicity in the TESS sector 4 residuals. In both cases, the peak of the periodogram was still at 0.22004 d but with a FAP greater than 68 per cent. This confirms that the SHO kernel adequately removes the periodic trend from the sector 4 TESS data.

We checked the consistency of the fitted planet parameters between the three fits. The majority of the fitted planet parameters were consistent between all three of the fits within 1σ and the remaining parameters were consistent within 2σ, except for the argument of periastron of the outer planet. There was a disagreement greater than 3σ between the values from the joint M32 fit and the decoupled fit. Constraining eccentricity and the argument of periastron is challenging with photometry alone and we remind the reader that we only included them in our fit to ensure that the uncertainties on the other fitted parameters were not underestimated. We conclude that the fitted planet parameters are not significantly affected by the presence of the periodic signal in TESS sector 4.

We also compared the Bayes evidence (dlnZ) of the three fits (Table 6). We found a decisive preference for both of the fits incorporating the SHO kernel over the original fit (Kass & Raftery 1995). The joint M32 plus SHO fit had the highest evidence, preferred over the original joint M32 fit with dlnZ = 106.8, and the decoupled fit of sector 4 and 31 was preferred over the original joint M32 fit with dlnZ = 89.2.

Despite the fact the evidence favoured the model with the M32 plus SHO kernel jointly fit to sectors 4 and 31, the model where we decoupled sector 4 and sector 31 is more physically motivated. This is because we only detected the periodic signal in sector 4. We therefore chose the decoupled fit as our final global photometric fit and we present the results in Section 5.1.

We performed one last test to assess the dependence of our results on our chosen detrending model – we repeated the decoupled fit, replacing the M32 kernels with SHO kernels. For sector 4, we used an SHO kernel for the short-term quasi-periodic signal summed with a second SHO kernel for the longer term variability. For sector 31, we used a single SHO kernel. All of the fitted planet parameters were fully consistent with our final results (see Table 7) within 1σ, except for the eccentricity of the inner planet which was consistent within 2.3σ. We conclude that our results are not significantly influenced by the choice of GP kernel.

4.3 Transit timing variation analysis

From the global photometric analysis, we found that HD 15906 b and c orbit close to a 2:1 period commensurability (P_c / P_b = 1.976), an indication that the planets might be in mean motion resonance (MMR). Planets in or near a low-order period commensurability have amongst the largest amplitude TTVs (e.g. Veras, Ford & Payne 2011; Agol & Fabrycky 2018), so we therefore checked for TTVs in the HD 15906 system.

juliet can incorporate TTVs into a photometric model, however, it expects that each instrument contains at least one transit of all the planets being fit. This is not true in our case – none of the CHEOPS or LCOGT observations contain a transit of both planets. We therefore had to perform a separate TTV fit for each planet. When fitting the inner planet, we included the TESS data, CHEOPS visit 5, and both LCOGT visits. For the outer planet, we included the TESS data and CHEOPS visits 3 and 6. In total, we had seven transits of the inner planet and four transits of the outer planet.

For the fit of each planet, we used a model consisting of one transiting planet and the same detrending as described in Section 4.2. The only difference in the transit model was that we fit for the individual transit times instead of P and T₀. We set a uniform prior of width 0.1 d on each transit time based upon a visual inspection of the data. All other priors were unchanged from the global photometric analysis and we used dynesty to sample the posterior of the model with nested sampling. Our results are presented in Section 5.2.

4.4 Radial velocity analysis

From HARPS and FIES, we have 25 sparsely sampled RV data points that show a relatively large scatter (see Fig. 7). We ran a GLS periodogram on the RV data and found no significant peaks at the planetary periods. The strongest peak was at 12.27 d and the best-fitting sinusoid with this period had an amplitude of ∼ 10 ms⁻¹. It is possible that this signal is caused by stellar activity, but with such large gaps between each observation, the short-period peaks in the GLS periodogram are unreliable. We removed the best-fitting sinusoid from the RV data and re-ran the GLS periodogram – no additional peaks emerged.

Figure 7.

Upper: HD 15906 RV time-series, highlighting the sparsity of the data. HARPS data taken before/after the 2015 fibre upgrade is plotted (green squares/blue triangles) alongside the FIES data (red circles). Lower: GLS periodogram of the HARPS and FIES RV data. The photometrically derived orbital periods of the two planets, see Table 7, are indicated by the blue and pink vertical dashed lines and the red horizontal line represents the 1 per cent false alarm probability level. There are no significant peaks at the planet periods and the strongest peak is at 12.27 d.

Open in new tabDownload slide

To search for the planetary signals, we performed a series of fits to the HARPS and FIES RV data using juliet. For our first fit, we assumed that there were no planets in the system and we fit only for an offset and a white noise term for each instrument. We used a uniform prior for the offset, in the range −20 to 20 ms⁻¹, and a log-uniform prior for the white noise term, in the range 0.01 to 20 ms⁻¹. The HARPS data from before and after the fibre upgrade had to be treated as two independent instruments. However, we only had three data points from post-upgrade which was insufficient to constrain the instrumental parameters. We therefore excluded the three post-upgrade HARPS data points from our fits and we used dynesty to sample the posterior of the model.

We then added planets to our model. We performed one fit with only the inner planet, one fit with only the outer planet, and finally a fit with both planets. We used a Keplerian for each planet, generated via radvel (Fulton et al. 2018), with the following parametrization:

• Orbital period, P, and mid-transit time, T₀. We fixed these to the solution from the global photometric fit (Table 7).

• Eccentricity, e, and argument of periastron, ω. For simplicity, we fixed eccentricity to zero and ω to 90 degrees.

• Semi-amplitude, K. We used a broad uniform prior to allow exploration of the range 0 to 20 ms⁻¹.

Finally, we took the model with both planets and added a GP with a quasi-periodic kernel (Foreman-Mackey et al. 2017) to account for the stellar activity. This kernel is described by four hyperparameters: the amplitude, period, an additive factor impacting the amplitude and the scale of the exponential component. For the amplitude we used a uniform prior in the range 0 to 20 ms⁻¹ and for the period we defined a normal prior using the peak from the periodogram analysis (μ = 12.27 d, σ = 0.1 d). The other two hyperparameters were allowed to vary uniformly over a broad range. With such a small number of sparsely sampled RVs, the GP was unlikely to yield a meaningful result but we chose to include it for completeness. The results of our RV modelling are presented in Section 5.3.

We note that we also tried a joint fit of the TESS, CHEOPS, and LCOGT photometric data with the HARPS and FIES RV data using juliet. The photometric model was identical to that presented in Section 4.2 and we used the RV model with two planets but no GP. However, due to the small number of sparse RVs, the fitted planet parameters were adversely affected compared to those from the global photometric model. Therefore, we decided to present independent analyses of the photometry and RVs in this paper.

5 RESULTS

5.1 Global photometric results

In Section 4.2, we described our joint fit of the TESS, CHEOPS, and LCOGT photometry and we present the resulting fitted planetary parameters in Table 7. We also include the following derived planetary parameters: transit depth (δ = (R_P/R_⋆)²), planet radius (R_P), semimajor axis (a), orbital inclination (i), total transit duration (T_dur), insolation flux (S_P), and equilibrium temperature assuming zero bond albedo and full day-night heat redistribution (T_eq).

Figs 1, 8, and 9 show the TESS, CHEOPS, and LCOGT data alongside the global photometric model. Fig. 10 shows the detrended TESS and CHEOPS data, phase-folded on each planet with the best-fitting transit model, and Fig. 11 shows the same for the LCOGT data. For a full list of posterior values and the corner plots presenting the posterior distributions of the fitted planetary parameters, see Appendix  A.

Figure 8.

Results of the global photometric fit. This plot shows the six CHEOPS light curves, where the 60 s cadence data (grey) has been binned to 20 min (black squares) to guide the eye. The red line is the median model from the global photometric fit and the red shaded region is the 1σ uncertainty on the model. The blue and pink markers indicate the mid-transit times of the inner and outer planets, respectively. The residuals of the model are included in the panel beneath each light curve.

Open in new tabDownload slide

Figure 9.

Results of the global photometric fit. This plot shows the two LCOGT light curves, where the 60 s cadence data (grey) has been binned to 20 min (black squares) to guide the eye. The red line is the median model from the global photometric fit and the red shaded region is the 1σ uncertainty on the model. The blue markers indicate the mid-transit times of the inner planet. The residuals of the model are included in the panel beneath each light curve.

Open in new tabDownload slide

Figure 10.

Results of the global photometric fit. Upper: Phase-folded TESS (top) and CHEOPS (bottom) light curves for the inner (left-hand panel) and outer (right-hand panel) planet. The light curves have been detrended to remove the instrumental and stellar variability and the data (grey) has been binned to 20 min (black squares) to guide the eye. The median transit models for the inner (blue line) and outer (pink line) planet are included, along with 50 random samples drawn from the posterior distribution of the model. Lower: Residuals of the median transit models. Note that an arbitrary offset has been applied to the TESS data and residuals for visibility purposes.

Open in new tabDownload slide

Figure 11.

Results of the global photometric fit. Upper: Phase-folded LCOGT light curve for the inner planet. The light curve has been detrended to remove instrumental effects and the data (grey) has been binned to 20 min (black squares) to guide the eye. The median transit model for the inner planet (blue line) is included, along with 50 random samples drawn from the posterior distribution of the model. Lower: Residuals of the median transit model.

Open in new tabDownload slide

In Fig. 10, there is a small dip during the transit of the outer planet which occurs just before the mid-transit position in both the TESS and CHEOPS phase-folded light curves. Rather than being a significant feature, it is most likely a coincidence. In the CHEOPS data, there is very poor coverage of this part of the transit and the dip is exaggerated by binning. In the TESS data, the mid-transit dip is only present in the first of the two transits.

Our analysis has shown that HD 15906 b is a 2.24 R_⊕ planet orbiting its host star at a separation of 0.090 AU with a period of 10.92 d. HD 15906 c is bigger (2.93 R_⊕) and orbits the host star at a larger separation (0.141 AU) with a longer period (21.58 d). The fit favoured slightly eccentric orbits (e_b = 0.11, e_c = 0.04) with a high impact parameter (b_b = 0.86, b_c = 0.90), but the transits of both planets are non-grazing. The inner and outer planet receive 33.1 and 13.4 times the amount of flux that the Earth receives from the Sun and, assuming zero bond albedo and full day-night heat redistribution, they have equilibrium temperatures of 668 and 532 K. We remind the reader that we repeated our global photometric fit with zero eccentricity and all fitted planet parameters were consistent within 1.2σ. In this case, we derived planetary radii of 2.24 ± 0.07 R_⊕ and 2.84 ± 0.05 R_⊕ for HD 15906 b and c, respectively.

5.2 Transit timing variation results

In Section 4.3, we described our TTV analysis of the HD 15906 system. The fitted observed transit times for each planet are presented in Table 8. From these values, juliet derived the best-fitting period and mid-transit time for each planet, assuming a linear ephemeris. These values, and all other fitted planet parameters, were fully consistent with the results of the global photometric model (see Section 5.1) within 2σ.

Using the best-fitting period and mid-transit time, we computed the expected transit times for each planet. We then plotted an observed–computed (O–C) diagram, see Fig. 12, to show the TTVs. We found marginal evidence for TTVs – the maximum TTV is ∼ 10 min, but nine of the eleven transits are consistent with no TTVs within 3σ. With only eleven transits of two planets and a gap of ∼ 2 yr in the data, we did not attempt to model these TTVs. In Section 7.2, we simulate the expected TTV signals for the two planets and compare these predictions with the observations.

Figure 12.

Results of the TTV analysis. This plot shows the difference between the observed (O) transit time and the computed (C) transit time, assuming a linear ephemeris, for transits of the inner (blue) and outer (pink) planet from TESS (square), CHEOPS (diamond), and LCOGT (circle).

Open in new tabDownload slide

5.3 Radial velocity results

In an attempt to detect the two planetary signals in the HARPS and FIES data, we fit five models to the RVs (see Section 4.4). We tried a model with no planets, only the inner and outer planet, both planets and both planets plus a GP to model the stellar activity. In the fit with the GP, the posterior distributions of the GP hyperparameters were the same as the priors, which tells us the data were unable to constrain the GP model, as expected. In Table 9, we present the Bayes evidence of each fit compared to the fit with no planets. The model with no planets had the highest evidence, with a substantial or strong preference over the other models (Kass & Raftery 1995), and we therefore conclude that the two transiting planets are not detected in the current HARPS and FIES RV data. However, we can still utilize this data for validation purposes, see Section 6.1.

6 VETTING AND VALIDATION

It is important to confirm that the transits we observed with TESS, CHEOPS, and LCOGT were caused by planets orbiting HD 15906. We therefore need to rule out false positive scenarios, including:

• The target star is an eclipsing binary (EB).

• The target star has a gravitationally associated companion star that is either an EB or has transiting planets.

• There is an aligned foreground or background star, not gravitationally associated with the target star, that is either an EB or has transiting planets.

• There is a nearby star, with a small angular separation from the target star but not gravitationally associated with it, that is either an EB or has transiting planets.

Furthermore, it is important to check for nearby unresolved stars because, if not accounted for, the blended flux can lead to underestimated planetary radii and improper characterization of the host star (Ciardi et al. 2015; Furlan & Howell 2017, 2020).

As mentioned in Section 2.1, HD 15906 b passed all of the SPOC vetting tests. In addition, it has been shown that multiplanet systems are significantly less likely to be false positives than single planet systems, especially when the planets are smaller than 6 R_⊕ (Lissauer et al. 2012; Guerrero et al. 2021). In this section, we use additional observational and statistical techniques to validate the HD 15906 planetary system.

6.1 High-resolution spectroscopy

Using the HARPS and FIES data, we did not detect the RV signals induced by the two transiting objects (see Section 5.3). In this section, we use the HARPS data to rule out stellar masses for the transiting objects and place limits on the presence of a bound stellar companion.

In Fig. 13, we show the HARPS and FIES RVs folded on HD 15906 b and c using the ephemerides obtained in the global photometric analysis (Table 7).

Figure 13.

HARPS (green squares and blue triangles) and FIES (red circles) RV data folded on the inner (left-hand panel) and outer (right-hand panel) planet. The transits occur at phase zero. A Keplerian model (dotted line) has been plotted on each axis to guide the eye and the arrows illustrate that this is an upper limit. The model represents a planet on a circular orbit with a semi-amplitude equivalent to the RMS of the HARPS data (10.70 ms⁻¹), a proxy for its maximum value.

Open in new tabDownload slide

The RMS of the HARPS data (10.70 ms⁻¹) can be used as a proxy for the maximum possible semi-amplitudes of the two transiting objects. Using the stellar mass presented in Table 5, the orbital parameters presented in Table 7 and a semi-amplitude of 10.70 ms⁻¹, HD 15906 b has an upper mass limit of ∼ 32 M_⊕ and HD 15906 c has an upper mass limit of ∼ 39 M_⊕. This confirms that the two transiting objects must be of planetary mass.

Furthermore, under the assumption of a circular orbit and an orbital inclination of 90 degrees, the RMS of the HARPS data rules out a bound brown dwarf or star, with a mass greater than 13 M_Jupiter, out to ∼ 1500 AU. At the distance of HD 15906, this corresponds to an angular separation of ∼ 32 arcsec. Even down to an orbital inclination of 10 degrees, we can rule out a brown dwarf or stellar companion out to ∼ 45 AU, corresponding to an angular separation of ∼ 1 arcsec.

Finally, we checked for a linear drift in the RV data because this could be indicative of a long-period bound stellar companion. We chose the pre-upgrade HARPS data for this purpose because it has the longest baseline (> 9 yr). We used juliet to perform a fit of this data, using a model consisting of no planets, an offset, white noise, and a linear trend. The best-fitting gradient was consistent with zero within 1σ and this supports the conclusion that HD 15906 does not have a bound stellar companion.

6.2 Archival imaging

HD 15906 is a high proper motion star (μ = 195.97 mas yr⁻¹; Gaia Collaboration 2022). We therefore made use of archival imaging to check for foreground or background objects at the star’s present day position.

HD 15906 was observed on 1953 November 11 by the Oschin Schmidt Telescope, using a blue photographic emulsion (λ = 330–500 nm; Monet et al. 2003), during the first Palomar Observatory Sky Survey (POSS-I). It was observed again on 1979 September 21 by the UK Schmidt Telescope, using a blue photographic emulsion (λ = 395–540 nm; Monet et al. 2003), during the SERC-EJ survey. We downloaded these images from the Digitized Sky Survey (DSS)¹⁴ and plotted them in the first two panels of Fig. 14. HD 15906 was also observed in 2010 by the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016). We downloaded the i-filter Pan-STARRS image from the MAST and plotted it in the third panel of Fig. 14. Finally, HD 15906 was observed during TESS sector 31 in 2020. We downloaded the target pixel file (TPF) from the MAST and plotted the first good quality cadence in the final panel of Fig. 14.

Figure 14.

Images of HD 15906 spanning 67 yr, from 1953 to 2020. Left to right: POSS-I, SERC-EJ, Pan-STARRS, and TESS sector 31. All images are shown on a scale of 4 arcmin × 4 arcmin, except for the POSS-I image which is zoomed in to 1 arcmin × 1 arcmin, and centred on the 2020 position of HD 15906 (pink star). We overlaid the TESS apertures from sector 4 (blue) and sector 31 (red) on the images, as well as the 2020 positions of all known stars from Gaia DR3 (green stars; Gaia Collaboration 2022). Only one of these stars (TIC 632595010; TESS magnitude = 20.3) is within the TESS apertures.

Open in new tabDownload slide

HD 15906 moved ∼ 13 arcsec between the POSS-I observation in 1953 and TESS sector 31 in 2020. Using the POSS-I image, we rule out a foreground or background star at the TESS sector 31 position of HD 15906 down to a TESS magnitude of ∼ 18. A star this faint would be incapable of producing the transit signals we observe, even in the case of a full EB, and it would not significantly impact the derived planet parameters due to flux blending (Ciardi et al. 2015). We therefore conclude that our results are not affected by an unresolved foreground or background star.

6.3 High-resolution imaging

High-resolution imaging was used to search for nearby stars, bound or unbound, that could be contaminating the photometry. We observed HD 15906 with a combination of high-resolution resources, including near-infrared adaptive optics (NIR AO) imaging at the Keck and Lick Observatories and optical speckle imaging at Gemini-North and SOAR. While the optical observations tend to provide higher resolution, the NIR AO tend to provide better sensitivity, especially to lower mass stars. The combination of the observations in multiple filters enables better characterization of any companions that might be detected. The observations are described in detail in the following subsections and a summary is provided in Table 10. Fig. 15 shows the resulting images and contrast curves. No stellar companions were detected within the contrast and angular limits of each facility, essentially ruling out stars at least ∼ 7 magnitudes fainter than HD 15906 between 0.5 and 10 arcsec. At small angular separations, where high-resolution imaging does not achieve a high contrast, we used high-resolution spectroscopy to rule out bound companions within ∼ 1 arcsec (see Section 6.1).

Figure 15.

High-resolution imaging of HD 15906. Upper: From left to right are the speckle images from SOAR and Gemini-North and the NIR AO images from Lick and Keck. Each image is zoomed into a region of 1.6 arcsec × 1.6 arcsec centred on HD 15906. Lower: Contrast curves from each observation.

Open in new tabDownload slide

6.3.1 SOAR

We searched for stellar companions to HD 15906 with speckle imaging on the 4.1 m Southern Astrophysical Research (SOAR) telescope (Tokovinin 2018) on 2019 July 14. We observed in Cousins I-band, a similar visible bandpass to TESS. This observation was sensitive to a star 5.3 magnitudes fainter than HD 15906 at an angular distance of 1 arcsec from the target. More details of the observations within the SOAR TESS survey are available in Ziegler et al. (2020). The 5σ detection sensitivity and speckle autocorrelation functions from the observations are shown in Fig. 15. No nearby stars were detected within 3 arcsec (∼ 137 AU, if bound) of HD 15906 in the SOAR observations.

6.3.2 Lick

We observed HD 15906 on 2019 July 21 using the Shane Adaptive optics infraRed Camera-Spectrograph (ShARCS) camera on the Shane 3 m telescope at Lick Observatory (Kupke et al. 2012; Gavel et al. 2014; McGurk et al. 2014). Observations were taken with the Shane AO system in natural guide star mode in order to search for nearby, unresolved stellar companions. We collected a single sequence of observations using a Ks filter (λ₀ = 2.150 |$\mu$|m, Δλ = 0.320 |$\mu$|m). We reduced the data using the publicly available SImMER pipeline¹⁵ (Savel et al. 2020). Our reduced image and corresponding contrast curve is shown in Fig. 15. The observations rule out stellar companions ∼ 4 magnitudes fainter than HD 15906 at 0.5 arcsec (∼ 23 AU, if bound) and ∼ 9 magnitudes fainter between 2 arcsec and 10 arcsec (∼ 92–458 AU, if bound).

6.3.3 Gemini-North

HD 15906 was observed on 2019 October 15 using the ’Alopeke speckle instrument on the Gemini-North 8 m telescope (Scott et al. 2021; Howell & Furlan 2022). ’Alopeke provides simultaneous speckle imaging in two bands (562 and 832 nm) with output data products including a reconstructed image with robust contrast limits on companion detections. Three sets of 1000 × 0.06 s images were obtained and processed in our standard reduction pipeline (see Howell et al. 2011). Fig. 15 includes our final 5σ contrast curves and the 832 nm reconstructed speckle image. We find that HD 15906 has no companion stars brighter than 5–8 magnitudes below that of the target star within the angular and image contrast levels achieved. The angular region covered ranges from the 8 m telescope diffraction limit (20 mas) out to 1.2 arcsec (∼ 0.9 to 55 AU, if bound).

6.3.4 Keck

HD 15906 was observed with NIR AO high-resolution imaging at the Keck Observatory on 2020 September 9. The observations were made with the NIRC2 instrument, which was positioned behind the natural guide star AO system (Wizinowich et al. 2000), on the Keck-II telescope. We used the standard 3-point dither pattern to avoid the lower left quadrant of the detector which is typically noisier than the other three quadrants. The dither pattern step size was 3 arcsec and was repeated twice, with each dither offset from the previous dither by 0.5 arcsec. The camera was in the narrow-angle mode with a full field of view of ∼ 10 arcsec and a pixel scale of approximately 0.0099442 arcsec pix⁻¹. The observations were made in the narrow-band Br-γ filter (λ₀ = 2.1686 |$\mu$|m, Δλ = 0.0326 |$\mu$|m) with an integration time of 0.5 s with one co-add per frame for a total of 4.5 s on target. The AO data were processed and analysed with a custom set of idl tools (see description in Schlieder et al. 2021) and the resolution of the final combined image, determined from the FWHM of the PSF, was 0.048 arcsec. The sensitivity of the combined AO image was determined according to Furlan et al. (2017) and the resulting sensitivity curve for the Keck data is shown in Fig. 15. The image reaches a contrast of ∼ 7 magnitudes fainter than the host star between 0.5 and 4 arcsec (∼ 23 to 183 AU, if bound) and no stellar companions were detected.

6.4 Gaia assessment

We used Gaia DR3 (Gaia Collaboration 2022) to show that there are no nearby, resolved stars bright enough to cause the transits we observe. The images presented in Fig. 14 show that there is only one Gaia DR3 star within the TESS optimal apertures. This is TIC 632595010 with a TESS magnitude of 20.3 (> 10 mag fainter than HD 15906) and a separation of ∼ 50 arcsec from HD 15906. This star is not bright enough to be the source of the transit signals we see, even in the case of a full EB. Furthermore, as explained in Section 2.3, the LCOGT observations confirmed that the transit signals do not originate from any of the known Gaia DR3 stars.

We also searched for wide stellar companions that may be bound members of the system. Based upon similar parallaxes and proper motions (Mugrauer & Michel 2020, 2021), there are no additional widely separated companions identified by Gaia.

Finally, the Gaia DR3 astrometry provides additional information on the possibility of inner companions that may have gone undetected by either Gaia or the high-resolution imaging/spectroscopy. The Gaia Renormalized Unit Weight Error (RUWE) is a metric, similar to a reduced chi-square, where values that are ≲ 1.4 indicate that the Gaia astrometric solution is consistent with a single star whereas RUWE values ≳ 1.4 may indicate an astrometric excess noise, possibly caused by the presence of an unseen companion (e.g. Ziegler et al. 2020). HD 15906 has a Gaia DR3 RUWE value of 1.15, indicating that the astrometric fits are consistent with a single star model.

6.5 Statistical validation

We finally used triceratops (Tool for Rating Interesting Candidate Exoplanets and Reliability Analysis of Transits Originating from Proximate Stars; Giacalone et al. 2021) to statistically validate the two transiting planets in the HD 15906 system. This Bayesian tool uses the stellar and planet parameters, the transit light curve, and the high-resolution imaging to test the false positive scenarios listed at the start of Section 6 and calculate the false positive probability (FPP) and the nearby false positive probability (NFPP) of TESS planet candidates. The FPP is the probability that the observed transit is not caused by a planet on the target star and the NFPP is the probability that the observed transit originates from a resolved nearby star. To consider a planet candidate validated, it must have FPP < 0.015 and NFPP < 0.001.

We ran triceratops on both HD 15906 b and c. We used the stellar parameters presented in Table 5, the planet parameters presented in Table 7, the combined TESS, CHEOPS, and LCOGT light curve and the high-resolution imaging contrast curves from Section 6.3. triceratops only accepts one contrast curve as input, so we ran the analysis with each of the five contrast curves and compared the results. In agreement with our analysis in Section 6.4, triceratops did not identify any nearby resolved stars that were bright enough to be the source of the transits. The results confirmed that the highest probability scenario was that of two planets transiting HD 15906. The most probable form of false positive scenario for the inner planet was an unresolved background EB and for the outer planet was an unresolved bound companion that is an EB. With our archival imaging (Section 6.2) and high-resolution spectroscopy (Section 6.1), that triceratops does not consider, these scenarios become less likely. For the Gemini-North and Keck contrast curves, both planets were validated with a negligible value of NFPP and FPP < 0.015. With the SOAR and Lick contrast curves, both planets had a negligible value of NFPP, the outer planet had FPP < 0.015 and the inner planet had a FPP just greater than 0.015 (0.0159 for Lick and 0.0166 for SOAR). According to the triceratops criteria, this means that the inner signal is likely a planet. However, triceratops does not account for the fact that multiplanet systems are more likely to be real (Lissauer et al. 2012; Guerrero et al. 2021), so the fact that the outer planet was validated means the inner planet may also be considered validated. We therefore conclude that both HD 15906 b and c are validated planets according to the triceratops criteria.

7 DISCUSSION

We have presented the discovery of the HD 15906 multiplanet system. In this section, we discuss our results, compare the system to other confirmed exoplanets, and assess the feasibility of future follow-up observations.

7.1 TESS periodicity

In Section 4.2.3, we reported the detection of a sinusoidal-like signal in the TESS sector 4 light curve of HD 15906. This signal has a period of 0.22004 d (∼ 5 h) and the best-fitting sinusoidal model has an amplitude of ∼ 57 ppm, equivalent to the transit depth expected for a planet with a radius of ∼ 0.63 R_⊕. In this section, we provide a discussion of this signal and its origin.

The 0.22 d periodic signal is present in the TESS sector 4 light curve, but not the sector 31 light curve. The signal is present in the sector 4 SAP and PDCSAP flux, but not in the background flux or centroid position. We checked for a periodic signal in the nearest star of comparable magnitude (TIC 4646803; TESS magnitude = 9.51, separation = 167 arcsec). This star was observed at 30 min cadence in sector 4, so we searched the TESS-SPOC light curve (Caldwell et al. 2020) and found no periodicity at 0.22 d.

We also extracted our own HD 15906 light curves from the TESS TPFs for both sectors. This was done using a default quality bitmask and optimizing the aperture mask to reduce the combined differential photometric precision (CDPP) noise in the resulting data. The extracted target fluxes were sky-corrected using a custom background mask. Detrending was done in two steps: scattered light was corrected for using a principal component analysis and any flux modulation caused by spacecraft jitter was removed by a linear model detrending using co-trending basis vectors and the mean and average of the engineering quaternions as the basis vectors. This second method has shown promise in cleaning up TESS photometry previously (Delrez et al. 2021). Our light curves were consistent with the TESS SAP and PDCSAP flux; our sector 4 light curve contained a 0.22 d periodicity and our sector 31 light curve did not. We can therefore confirm that the periodic signal is not dependent on light curve extraction technique.

Furthermore, we performed experiments extracting light curves from apertures of different sizes and found that using an aperture of radius 1 pixel centred on HD 15906 resulted in a significantly larger amplitude variability (roughly by a factor of two) than when we used an aperture of radius 4 pixels. This is not what we would expect for a signal originating from within a pixel of HD 15906 (where we would expect the amplitude to stay roughly constant given the lack of nearby bright stars to dilute the flux) or from a blended star from larger distances (which should show larger amplitude in larger apertures).

We considered the possibility that the periodic signal is a form of stellar activity originating from HD 15906. However, a variety of arguments suggested that this was unlikely. First, the signal is strongly present in TESS sector 4, but is undetectable in any other observations. The very short period of the signal strongly disfavours it being related to the rotation period of HD 15906, given the star’s narrow spectral lines and amenability to precise RV measurements. The period (∼ 5 h) is consistent with the time-scale of granulation on the surface of a Sun-like star, but this process does not create sharp periodicities like we detected (see Fig. 6, which shows a clearly defined sharp peak in the periodogram of the sector 4 TESS residuals). Stellar pulsations can sometimes create such sharp periodicities, but main sequence stars of this type should not exhibit any pulsations on similar amplitudes or time-scales.

We finally searched for evidence that the signal originated from another star on the TESS detectors and contaminated the light curve of HD 15906 through a process other than direct overlap of the PSFs. This was a frequent occurrence during the Kepler and K2 missions (Coughlin et al. 2014) but is much less common during the TESS mission due to differences in the design of the telescopes, electronics, and optics. The bright contact binary DY Cet (TIC 441128066; TESS magnitude = 9.23) was observed on the same CCD as HD 15906 during TESS sectors 4 and 31. This EB has a period of 0.4408 d and the TESS light curves show a sinusoidal-like variability with a period of 0.2204 d (Yıldırım 2022). This is consistent with the period of the signal we detected in the TESS sector 4 light curve of HD 15906. During sector 4, DY Cet was in the same CCD columns as HD 15906, but during sector 31 it was not. We therefore conclude that the flux from DY Cet contaminated that of HD 15906 in TESS sector 4 during CCD readout, although the exact mechanism of contamination is currently unknown. We note that DY Cet cannot be the source of the transits of HD 15906 b and c, which have been independently observed by CHEOPS and LCOGT, and we reiterate that this periodic signal does not affect our fitted planet parameters (see Section 4.2.3).

7.2 Transit timing variation predictions

In Section 5.2, we reported marginal evidence for TTVs in the HD 15906 system. Here, we compute the expected TTV signals for each planet and compare them with our observations. We reiterate that modelling the TTVs is beyond the scope of this work due to the small amount of data.

We compared two methods for simulating the TTV signals of HD 15906 b and c. The first uses an approximated estimation of the TTV signal by modelling it as a linear combination of basis functions as described in Hadden et al. (2019) and implemented in the TTV2Fast2Furious (ttv2f2f)¹⁶ package. The second approach is a direct N-body simulation with the trades¹⁷ code (Borsato et al. 2014, 2019; Nascimbeni et al. 2023). For each planet, we used the orbital period and inclination presented in Table 7, the predicted mass (see Section 7.4.1) and we assumed a circular orbit for simplicity. We simulated the TTV signal for planets b and c for a time range that covers the full range of transit observations and we present the results in Fig. 16.

Figure 16.

The predicted TTV signal, as computed with TTV2Fast2Furious (ttv2f2f, crosses and solid line) and trades (open diamonds and solid line), for HD 15906 b (upper panel) and HD 15906 c (lower panel). We have also included the observed TTVs for each planet (blue circles and pink squares).

Open in new tabDownload slide

The predictions for HD 15906 b are generally in good agreement with the observations. We note that there is a slight difference between the amplitudes of the two simulated TTV signals, but both are consistent with most of the observations given their uncertainties. However, the predictions for HD 15906 c seem to be in antiphase with the observations. As expected for a two planet configuration close to a first-order period commensurability, the predicted TTVs of planets b and c are anticorrelated (Agol & Fabrycky 2018). Contrary to this expectation, the observed TTVs of HD 15906 b and c appear to be correlated. This could suggest that there is an additional, undetected planet in the system perturbing the orbits of the two observed planets. Alternatively, the TTVs might be spurious or affected by excess systematic noise from, for example, stellar activity (e.g. Oshagh et al. 2013; Ioannidis, Huber & Schmitt 2016). Given the sparse sampling of the TTV signals, future observations are required to assess the true nature of the TTVs.

7.3 Comparison with confirmed exoplanets

HD 15906 b and c have radii of 2.24 R_⊕ and 2.93 R_⊕, respectively, meaning they cannot have a purely rocky composition (Rogers 2015; Lozovsky et al. 2018). They both fall on the upper side of the radius gap and we therefore classify them as sub-Neptunes. Furthermore, with insolation fluxes of 33 S_⊕ and 13 S_⊕, and equilibrium temperatures of 668 K and 532 K, both planets are in the warm regime (T_eq ≲ 700 K).

Of more than 5300 confirmed exoplanets¹⁸, there are 66 sub-Neptune sized planets (1.75 < R_P/R_⊕ < 3.5) transiting bright stars (G ≤ 10 mag). Only 18 of these have an insolation flux less than HD 15906 b and only 5 have an insolation flux lower than HD 15906 c – GJ 143 b (Dragomir et al. 2019), ν² Lupi d (Delrez et al. 2021), HD 23472 c (Barros et al. 2022), HD 73583 c (Barragán et al. 2022), and Kepler-37 d (Marcy et al. 2014). HD 15906 c is therefore one of the most lowly irradiated sub-Neptune planets transiting such a bright star. Furthermore, there are only 5 other multiplanet systems with two warm (T_eq ≤ 700 K) sub-Neptune sized planets (1.75 < R_P/R_⊕ < 3.5) transiting a bright (G ≤ 10 mag) star – HD 108236 (Bonfanti et al. 2021), ν² Lupi (Delrez et al. 2021), HD 191939 (Orell-Miquel et al. 2023), HD 23472 (Barros et al. 2022), and TOI 2076 (Osborn et al. 2022). The HD 15906 system is therefore an interesting target for future follow-up studies, discussed further in Section 7.4.

Due to the nature of its observing strategy, TESS is biased towards the discovery of short-period planets; less than 14 per cent of planets confirmed by TESS have periods longer than 20 d, of which only half have radii smaller than 4 R_⊕. This work has demonstrated how CHEOPS can be used to follow-up TESS duotransits to expand the sample of long-period planets. Fig. 17 presents a period–radius diagram comparing the two planets in the HD 15906 system to the confirmed exoplanet population. The other TESS duotransits resolved by CHEOPS have also been included – TOI 2076 c and d (Osborn et al. 2022), HIP 9618 c (Osborn et al. 2023), TOI 5678 b (Ulmer-Moll et al. 2023), and HD 22946 d (Garai et al. 2023). Through our CHEOPS duotransit programme, we have contributed to the discovery of six planets with periods longer than 20 d, radii smaller than 5 R_⊕, and host stars brighter than G = 12 mag. There are only 18 other planets confirmed by TESS in this parameter space, illustrating the power of the TESS and CHEOPS synergy for the discovery of small, long-period planets transiting bright stars.

Figure 17.

Period–radius diagram of confirmed exoplanets with a Gaia magnitude brighter than 12, where discoveries made by TESS are highlighted in red. HD 15906 b (blue circle) and c (pink square) are included, alongside the five additional TESS duotransits resolved by CHEOPS (green squares).

Open in new tabDownload slide

7.4 Potential for future follow-up

With two warm sub-Neptunes transiting a bright (G ∼ 9.5 mag) K-dwarf, the HD 15906 system is an excellent target for future observations to measure the masses of the planets and perform atmospheric characterization. Warm sub-Neptunes are less affected by radiation from their host star than their hot counterparts, meaning their atmospheres will not have been sculpted so heavily by photoevaporation and they will more closely resemble their primordial state. Observations of these planets are therefore crucial in testing models of the formation and evolution of sub-Neptune planets. In addition, as a multiplanet system, HD 15906 will allow for comparative studies of internal structure and composition as a function of stellar irradiation.

7.4.1 Radial velocity

We were unable to detect HD 15906 b and c in the current HARPS and FIES RV data due to the small number of sparsely sampled observations (see Fig. 7 and Section 5.3). Here, we use the results of our global photometric analysis (Table 7) to predict the expected mass and semi-amplitude of the two planets. Otegi, Bouchy & Helled (2020) present a mass–radius relation that is dependent upon the density of the planet (ρ_P):

The high density case is applicable when the planet has a rocky composition and the low density case is for when the planet has a volatile-rich composition. Assuming a volatile-rich composition, the inner and outer planet would have masses of 6.21|$^{+1.56}_{-1.43}$|M_⊕ and 9.47|$^{+2.44}_{-2.21}$|M_⊕, respectively, leading to semi-amplitudes of 2.11|$^{+0.53}_{-0.49}$|ms⁻¹ and 2.54|$^{+0.67}_{-0.60}$|ms⁻¹.

The predicted semi-amplitudes of HD 15906 b and c are greater than the average HARPS RV uncertainty (∼ 1.5 ms⁻¹). This means that the planetary signals should be detectable with sufficient observations from a high-resolution spectrograph. Since HD 15906 is visible from both hemispheres, there are many instruments that would be capable of doing this. We note that there is a relatively large scatter in the current RV data (∼ 10 ms⁻¹) which could make a precise mass measurement challenging. It will require a high sampling rate and a large number of RV observations to adequately model the planetary and stellar signals.

7.4.2 Atmospheric characterization

Theory predicts a wide variety of possible chemical compositions for sub-Neptunes (e.g. Moses et al. 2013; Guzmán-Mesa et al. 2022). Atmospheric characterization can constrain their composition and, thanks to their bright host star, HD 15906 b and c are amenable to such observations. The Transmission Spectroscopy Metric (TSM; Kempton et al. 2018) can be used to rank transiting planets based on their suitability for transmission spectroscopy. It quantifies the expected SNR of the spectral features for a 10 h observation with JWST/NIRISS, assuming a cloud-free atmosphere. Using the stellar parameters from Table 5 and planet parameters from Table 7, we find that HD 15906 b and c have TSM values of 71.7 and 82.1, respectively. This puts them in the top 3 per cent of all confirmed transiting planets smaller than 4 R_⊕, where in the absence of a measured mass we computed the expected mass according to the empirical mass–radius relation used by Kempton et al. (2018). Furthermore, HD 15906 b and c have amongst the highest TSMs for such small and lowly irradiated planets, illustrated in Fig. 18. There are only six sub-Neptune sized planets (1.75 < R_P/R_⊕ < 3.5) with a higher TSM and lower irradiation than HD 15906 c, of which HD 15906 is the second brightest host star. In addition, there are only three other multiplanet systems which host two sub-Neptune sized planets with irradiation lower than 35 S_⊕ and TSM higher than 70 (HD 191939, TOI 2076 and TOI 270; Van Eylen et al. 2021; Osborn et al. 2022; Orell-Miquel et al. 2023). The HD 15906 system is therefore an interesting target for comparative studies of warm sub-Neptune composition.

Figure 18.

Transmission spectroscopy metric (TSM) as a function of stellar irradiation for all confirmed sub-Neptune sized planets (1.75 < R_P/R_⊕ < 3.5). The marker size scales with planet radius and HD 15906 b (blue) and c (pink) are plotted alongside a shaded region to show planets with a lower irradiation and higher TSM. HD 15906 is one of four multiplanet systems with two sub-Neptune sized planets with irradiation lower than 35 S_⊕ and TSM higher than 70. The other three systems are highlighted in red (HD 191939), green (TOI 2076), and purple (TOI 270).

Open in new tabDownload slide

8 CONCLUSIONS

In this paper, we have reported the discovery and validation of two warm sub-Neptune planets transiting the bright (G = 9.5 mag) K-dwarf HD 15906 (TOI 461, TIC 4646810). During TESS sectors 4 and 31, four transits of the inner planet, HD 15906 b, were observed, but there were only two transits of the outer planet, HD 15906 c, separated by ∼ 734 d. The period of the outer planet was ambiguous, with 36 possible values, and we used CHEOPS follow-up to determine the true period. Using TESS, CHEOPS, and LCOGT photometry, we precisely characterized the two planets – HD 15906 b and c have periods of 10.92 d and 21.58 d and radii of 2.24 R_⊕ and 2.93 R_⊕, respectively. We found marginal evidence for TTVs in the system and, comparing the observations to simulations, we showed that more observations are required to understand the nature of the TTV signals. Both planets are in the warm regime, with insolation fluxes of 33.1 S_⊕ and 13.4 S_⊕ and equilibrium temperatures of 668 K and 532 K. We find that HD 15906 c is one of the most lowly irradiated sub-Neptune sized planets transiting such a bright star.