© The Authors 2025

1 Introduction

Super-Earths and sub-Neptunes represent a class of planets that do not exist in the Solar System. Theoretical works have sought to understand why such planets do not orbit the Sun. Notably, it has been proposed that Jupiter played a role in inhibiting the formation of super-Earths and sub-Neptunes (Izidoro et al. 2015; Lambrechts et al. 2019). With estimated radii in the range 1 R_⊕ < R_p < 4 R_⊕ , the discovery of these planets around other stars was the result of large surveys, coupled with important instrumental developments. Radial velocity (RV) blind searches with high-resolution spectrographs (e.g. HARPS, Pepe et al. 2000), together with the space-based Kepler transit survey (Borucki et al. 2010) have revealed this new category of planets that is intermediate in mass and radius to Earth and Neptune (Mayor & Udry 2008). Since then, the number of characterised super-Earths and sub-Neptunes has increased dramatically thanks to the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015). The efficient coordination with ground-based RV follow-up of transiting exoplanet candidates has allowed the community to fill in the mass-radius space and has provided important insights into the composition of these planets. While super-Earths and sub-Neptunes overlap in their range of masses (e.g. Otegi et al. 2020; Parc et al. 2024), they show a more distinct separation in radius (Fulton et al. 2017; Mayo et al. 2018; Cloutier & Menou 2020). This separation, known as the radius valley, suggests also a separation in bulk density, as proposed by Luque & Pallé (2022) for planets around M-dwarf stars. This separation in density has been more recently supported by pebble-based formation models for all stellar types (Venturini et al. 2024); however, the question of a structural difference between various populations of sub-Neptunes is still heavily debated (e.g. Rogers et al. 2023; Parc et al. 2024). Understanding the composition of super-Earths and sub-Neptunes, as well as the origin for such differences, will bring important insights into the general theory of planet formation and evolution.

Several models attempt to explain the emergence of the observed differences between super-Earths and sub-Neptunes. Photo-evaporation (Owen & Wu 2013) and core-powered mass loss (Ginzburg et al. 2018) are two popular hypotheses positing that some planets would have lost most of their atmospheres to become the super-Earths we observe today. Other hypotheses explore in situ formation of super-Earths, where the planets would form at a later stage in a gas-poor environment (Lopez & Rice 2018). These various models make different predictions regarding the location of the valley in the orbital period-radius space (e.g. Cloutier & Menou 2020). Therefore, a continued effort to detect and characterise super-Earths and sub-Neptunes is of utmost interest. Among the newly discovered planets, multi-planet systems are particularly valuable. They provide an opportunity to compare planet properties within the same formation environment and around the same star and to compare these characteristics to dynamical processes. Multi-planet systems are key laboratories to interpret population-level studies and constrain formation and evolution models.

Within this observational effort, the high-resolution echelle spectrograph HARPS-N (Cosentino et al. 2012) plays a key role. With a proven stability at the sub-m s⁻¹ level (Cosentino et al. 2014; John et al. 2023), HARPS-N has played a significant role in measuring the masses of small planet candidates identified with data from the Kepler (Borucki et al. 2010), K2 (Howell et al. 2014), and TESS (Ricker et al. 2015) missions. As of the 14th of May 2024, the HARPS-N Collaboration has measured the mass of more than a third of the transiting small planets (R_p < 4 R_⊕) with a precision better than 30%, playing a major role in populating the mass-radius diagram (e.g. Pepe et al. 2013; Malavolta et al. 2018; Mortier et al. 2020; Cloutier et al. 2020; Lacedelli et al. 2021; Rajpaul et al. 2021; Bonomo et al. 2023; Palethorpe et al. 2024).

In this paper, we report on the discovery and characterisation of a super-Earth and a sub-Neptune orbiting the K-dwarf TOI-1453 (TIC 198390247). The system is located in the TESS northern continuous viewing zone (CVZ) and benefits from a wealth of photometric data. We carried out follow-up RV observations with HARPS-N during two seasons. The combination of these two data sets allowed us to constrain the planetary masses and radii. The planets span the radius valley and have a ratio of orbital periods of nearly 3 to 2. Additionally, the system is orbited by a faint stellar companion. This system therefore contains the ingredients to set important constraints on the formation and evolution of super-Earths and sub-Neptunes.

In Sect. 2, we describe the ensemble of observations that we used to both validate the planet candidates and characterise them. Section 3 presents the astrometric and atmospheric characterisation of the star, as well as an analysis of its activity. We then tackle the data analysis in Sect. 4 so to constrain the planet physical parameters. In Sect. 5, we study the dynamics of this two-planet system and investigate the transit timing variations (TTVs). We explore the possible structures of the sub-Neptune TOI-1453 c in Sect. 6. Finally, we discuss the place of TOI-1453 in the context of an exoplanet population and present the main takeaways of this peculiar planetary system in Sect. 7.

2 Observations

This section summarises the ensemble of TOI-1453 observations used in the present work.

2.1 TESS photometry

The TESS observations of this star started on 2019 July 18 (sector 14), and have since then covered a large number of sectors. As of April 2024, TESS has observed TOI-1453 in sectors 14–17, 19–26, 40–41, 47, 49–54, 56–57, 59–60, and 73–76, for a total of 29 sectors and a baseline of 1711 days. All the data are available with a 2 minute cadence. We note also that 20 second cadence data are available for this star on the sectors comprised between 40 and 60. These data were processed by the Science Processing Operations Center at NASA Ames Research Center (SPOC; Jenkins et al. 2016).

An early report from the SPOC based on sectors 14, 15, and 16 identified a candidate transiting planet with an orbital period of 13.18 days. This finding was obtained using the Transiting Planet Search (TPS) pipeline module (Jenkins 2002; Jenkins et al. 2010, 2020). However, the phase-folded light curve presented a significant flux decrease at exactly half the period, suggesting that the true orbital period could be 6.59 days. The orbital period of this candidate was updated to 6.59 days in the SPOC report for sector 21, with a corresponding transit depth of 894 ± 228 ppm. MIT’s Quick Look Pipeline (QLP, Huang et al. 2020a,b) also performed a transiting planet search of sectors 14– 16 using FFI data and detected the signature of TOI-1453 c at the correct period of 6.59 days. The SPOC fitted an initial limb- darkened transit model to the signature (Li et al. 2019) which passed a suite of diagnostic tests to help make or break the planetary nature of the signal (Twicken et al. 2018). The TESS Science Office (TSO) reviewed the results and issued an alert for TOI-1453 c on 14 November 2019 (Guerrero et al. 2021). In the meantime, the SPOC Data Validation report for sectors 14– 19 detected the signature of a small second planetary candidate with an orbital period of 4.31 days and a transit depth of 243±243 ppm, and the TSO issued an alert for TOI-1453 b on 21 February 2020. Beyond detecting the signatures of TOI-1453 b and c, the SPOC difference image centroiding tests conducted for the search of sectors 14–55 constrained the location of the host star to within 5.0±5.5 arcsec of the transit source of TOI-1453 b, and to within 1.4±3.4 arcsec of the transit source of TOI-1453 c, respectively.

Figure 1 shows a sample TESS target pixel file from sector 24, illustrating the low spatial resolution of TESS (the pixel scale is 21″). As a result, blending with other sources is expected. The pipeline photometric aperture is selected to maximise the signal- to-noise ratio (S/N) of the target star. This is highlighted in Fig. 1 with light gray thick contours to the pixels. In this work, we employed the pipeline photometric aperture to retrieve the light curves. Some blended sources may be hidden by the target on the image. We therefore identified all the known nearby sources from the Gaia DR3 catalogue and marked these close sources in red on the figure. In particular, a very close source appears inside the pipeline photometric aperture. It has a relative Gaia magnitude of ∆G=5.24 with respect to the target and a similar proper motion to TOI-1453, indicative of a gravitationally bound object.

To retrieve the light curves, we employed a custom extraction from the pipeline aperture of the target pixel files. Using lightkurve (Lightkurve Collaboration 2018), we removed the background scattered light by applying a principal components analysis on the out-of-aperture pixels matrix and decorrelating the target flux from the five main components via a linear regression. Additionally, we removed the most common systematic trends observed in each TESS CCD and for each sector via the co-trending basis vectors (CBVs), which are generated in the pre-search data conditioning (PDC) component of the TESS pipeline. In particular, we employed a combination of multi-scale and spike CBVs, in order to correct for the systematics on different wavelet-based band passes as well as the short impulsive systematics. For each band pass of the multi-scale CBVs, we used eight CBVs. We extracted and corrected the light curve from systematics with the joint background correction and the CBV detrending, using the function CBVCorrector.correct_gaussian_prior in lightkurve. The resulting light curve did not need further detrending. Indeed, we performed a model comparison on each sector, between models with and without Gaussian process (GP) detrending. These analyses were carried out with the Juliet Python package (Espinoza et al. 2019), which employs the code batman (Kreidberg 2015) to fit for the two-planet transit model. We tested different GP kernels in our model comparison, each of them implemented in the code celerite (Foreman-Mackey et al. 2017) which Juliet depends on. The model parameters of the fit were explored using the nested sampling algorithm dynesty (Speagle 2020), which also delivers the Bayesian evidence of the model. Unambiguously, our analyses have provided more evidence to support use of the simpler model without GP detrending.

We also compared the resulting light curve with the PDC simple aperture photometry (PDCSAP; Stumpe et al. 2012; Smith et al. 2012; Stumpe et al. 2014). We noticed an overall improvement with our custom extraction, with a smaller, or at worst equivalent, scatter in the individual sectors. Furthermore, because of the quality flag applied in the PDCSAP light curve, a significant number of transits are missed. With our custom extraction, we could retrieve 21 additional transits of the inner planet candidate and 10 additional transits of the outer candidate over all the sectors with a minor impact on the light curve quality, as noted with the scatter comparison. In the following, we use the light curve obtained from our custom extraction. We provide a synthesis of each sector of our dataset in Table A.1.

Fig. 1 Example of TESS target pixel file centred on TOI-1453 (TIC 198390247). The target is highlighted with a white cross at the centre of the grid. The zone inside of which we perform aperture photometry is contoured in bold gray. All the sources in the neighbourhood of TOI- 1453 identified with the Gaia DR3 catalogue are shown in red, and the arrows illustrate their proper motions on the sky. This plot was designed with tpfplotter (Aller et al. 2020).

Fig. 2 AO image of TOI-1453 and corresponding sensitivity curve obtained with NIRC2 installed on the Keck2 telescope (cf. https://exofop.ipac.caltech.edu/tess/view_tag.php?tag=19377, PI: Gonzales). The filter used for this image is Brγ. The faint companion at about 2″ away from TOI-1453 is revealed by this image.

2.2 High-resolution imaging

High-resolution images were taken from the ground to search for nearby stars that might be blended with TOI-1453 in the TESS images. Their identification is important not only to identify potential false positive transits, but also to update the light dilution correction factor when needed. High-resolution images of TOI-1453 and its surrounding field were obtained with three different instruments and telescopes.

Two adaptive optics (AO) images were obtained on 2020 May 28 with the NIRC2 instrument coupled with the Keck2 10m telescope located at Mauna Kea observatory, Hawaii. One of the images was taken in the J_cont photometric band (1.21 µm), the other in the B_rγ filter (2.17 µm). These images have a resolution of 9.942 milliarcsec/pixel, and have a contrast detection threshold at 0.5″ angular separation of ∆mag = 6.43 and 7.41, respectively, for the J_cont and B_rγ filters. Figure 2 presents the contrast curve together with the resulting image, obtained in the B_rγ filter. TOI-1453 appears at the center of the image, and its faint stellar companion is seen about 2″ away.

Two speckle images were obtained on 2020 June 6 with the Alopeke instrument installed on the Gemini 8m telescope at Mauna Kea observatory, Hawaii. The achieved resolution of this data is 10 mas/pixel. Each image was obtained in a distinct narrow photometric band, around 562 nm and 832 nm, for which the contrast threshold at 0.5″ is measured to ∆mag = 4.83 and 6.62, respectively.

Finally, lucky images of TOI-1453 were taken on 2021 March 22 and 2021 September 29, using the AstraLux instrument on the 2.2m telescope at Calar Alto observatory. The observations were carried out in the SDSSz bandpass, and present a resolution of 23.27 milliarcsec/pixel. The first observation taken in March 2021 has a contrast detection threshold at 1″ of ∆mag = 5.6, while the second observation obtained in September 2021 has a threshold contrast of 4.3.

All these high-resolution observations jointly confirm that besides the known companion 1.9″ away from TOI-1453, it is most likely that the star has no other close stellar companion. A mismatch was identified by Lillo-Box et al. (2024) between the position of the stellar companion in the AstraLux lucky images and the Gaia data. The latter indicates a separation of 3.02″, while they report a separation of 1.85″ with AstraLux (consistent with the other high-resolution images). The origin of this discrepancy is yet unclear, and might indicate that the companion star moved between the Gaia epoch and 2020–2021.

2.3 TRES reconnaissance spectroscopy

Early on after the planet candidate alert released by the TSO, two spectra of TOI-1453 were obtained with the Tillinghast Reflector Echelle Spectrograph (TRES) installed on the 1.5 m telescope at Fred L. Whipple Observatory (FLWO), Arizona. The two observations were performed on 2020 January 30 and 2020 March 16. The aim of those observations was to confirm the single-line nature of the spectrum, assess the activity level of the star and its convenience for precise RV follow-up, and to reject false positive hypotheses. TRES presents a spectral resolution of 44 000 in the broad visible range (from 385 to 909.6 nm). The first observation revealed a single-line spectrum, appropriate for precise RV follow-up. Because of the initial confusion regarding the true orbital period of TOI-1453.01, the second TRES spectrum was taken at the same phase as the first one. A negligible velocity shift was noted.

2.4 Ground-based photometry

The TESS pixel scale is ~21″ pixel⁻¹ and photometric apertures typically extend out to roughly 1 arcminute, generally causing multiple stars to blend in the TESS photometric aperture. In our attempt to determine the true source of the TESS detection, we acquired follow-up ground-based photometric observations of the fields around TOI-1453 as part of the TESS Follow-up Observing Program (TFOP; Collins 2019)¹.

2.4.1 LCOGT

We used the TESS Transit Finder, which is a customised version of the Tapir software package (Jensen 2013), to schedule our transit observations. We observed two full transit windows of TOI-1453.01 (candidate planet c) in Pan-STARRS ɀ-short band on UTC 2020 May 17 and 2020 August 24 from the Las Cumbres Observatory Global Telescope (LCOGT) (Brown et al. 2013) 1 m network node at McDonald Observatory near Fort Davis, Texas, United States (McD). The observation on 2020 August 24 was simultaneously observed with two 1 m telescopes at McD. The 1 m telescopes are equipped with a 4096 × 4096 SINISTRO camera with an image scale of 0″.389 per pixel, resulting in a 26′ × 26′ field of view.

The images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018), and differential photometric data were extracted using AstroImageJ (Collins et al. 2017). We used circular photometric apertures with radii 5″.1 to 8″.2. The target star aperture included all of the flux from the nearest known neighbour in the Gaia DR3 catalogue (Gaia DR3 1433062331333157632), which is 1″.9 northeast of TOI-1453 and has a delta T-band magnitude of 5.4, making it ~100 times fainter than TOI-1453.

We visually compared the light curves of the 16 nearby stars within 2′.5 of the target star that are bright enough to potentially cause the TESS detection. No obvious nearby eclipsing binaries (NEBs) were detected in 15 of 16 of the stars, but the 1″.9 neighbour result is inconclusive due to heavy blending with the target star in the photometric aperture. We conclude that the transit is occurring either in TOI-1453, or the 1″.9 neighbour star. We investigate this further in Sect. 4.2. The light curve data are available on the EXOFOP-TESS website².

2.4.2 Maury Lewin Observatory

A full transit window of TOI-1453.02 was observed from Maury Lewin Observatory in Ic band on UTC 2020 May 09. The Maury Lewin Astronomical Observatory 0.36 m telescope is located near Glendora, CA. It is equipped with a 3326 × 2504 SBIG STF8300M camera having an image scale of 0″.84 per pixel, resulting in a 23′ × 17′ field of view. The images were calibrated and differential photometric data were extracted using AstroImageJ. The shallow ~250 ppm transit signal would be difficult to detect from ground-based observatories, so we instead checked for NEBs in the 29 nearby stars within 2′.5 from the target star that are bright enough to potentially cause the TESS detection. We found no obvious NEB signals, but the results do not rule out NEBs at high confidence in most of the stars, including the 1″.9 neighbour Gaia DR3 1433062331333157632. The NEB check light curves are available on the EXOFOP-TESS website.

2.4.3 Adams Observatory

Another full transit window of TOI-1453.02 was observed from Adams Observatory in the Rc band on UTC 2020 June 04. The Adams Observatory is located at Austin College in Sherman, TX. The 0.6 m telescope is equipped with a FLI ProLine PL16803 detector that has an image scale of 0″.38 pixel⁻¹, resulting in a 26′ × 26′ field of view. The observations from this facility were also aimed at ruling out NEBs in stars in close proximity to TOI-1453 at the predicted ephemerides of the candidate planet b. We did not detect any clear NEB signals, although again we cannot rule out this possibility with a high level of confidence. The results are available on the EXOFOP-TESS website.

2.4.4 MuSCAT2

MuSCAT2 (Narita et al. 2019) observed TOI-1453 on the night of March 13, 2024. The instrument is a multi-band optical imager mounted on the 1.5 m Telescopio Carlos Sánchez (TCS) at the Teide Observatory, Spain. MuSCAT2 is equipped with four cameras capable of taking simultaneous images in the g′, r′, i′, and ɀ_s bands. Each CCD has 1024 × 1024 pixels with a field of view of 7 4 × 7 4 arcmin² .

The telescope was defocused to avoid saturation of the target. The exposure times were set to 10 s, 10 s, 5 s and 5 s in g′, r′, i′ and ɀ_s, respectively. The data were reduced by the MuS- CAT2 pipeline (Parviainen et al. 2019); the pipeline performs dark and flat field calibrations, aperture photometry, and transit model fitting including instrumental systematics. Due to the faint reference stars, the light curve scatter was too large to detect the transit with any significance. Our analysis eliminated 3 stars as possible EBs during the observed window.

2.5 HARPS-N high precision spectroscopy

All the above follow-up observations increased our confidence in the likely planetary nature of the candidates TOI-1453.01 and TOI-1453.02. Early, we began a precise RV follow-up campaign with HARPS-N (Cosentino et al. 2012, 2014). HARPS-N is a high-resolution (R = 115 000) echelle spectrograph mounted on the 3.6 m Telescopio Nazionale Galileo (TNG) at Roque de los Muchachos observatory, La Palma, Spain. The HARPS-N instrument is stabilised in temperature and pressure so to achieve the submeter-per-second accuracy and its spectra cover a wavelength range between 383 and 693 nm.

2.5.1 DRS

We began the observations of TOI-1453 with HARPS-N on 2020 March 12 and gathered spectra for two seasons until 2021 September 19. In total, we obtained 100 RV measurements as part of the HARPS-N Collaboration Time (HCT). Each observation had an exposure time of 1800 s, which led to a median S/N of 58 in order 50. The RVs were extracted with version 3.0.1 of the data reduction software (DRS, Dumusque et al. 2021), adapted from the ESPRESSO pipeline with additional HARPS- N systematics corrections. We employed a G9 spectral mask to derive the cross-correlation function (CCF) and extract the RVs. Our DRS RV time series presents a median RV uncertainty per data point of 1.56 m/s and a root mean square (rms) deviation of 5.04 m/s.

2.5.2 s-BART

We also extracted the RVs with s-BART (Silva et al. 2022). This is a template-matching pipeline to extract achromatic RV variations. After removing from the spectra wavelength domains sensitive to stellar activity and telluric features, the pipeline creates a spectral template of the star. The latter is built from the series of observations obtained on the given star and assembling the template order-by-order. Finally, it employs a semi-Bayesian approach to undertake the template-matching, assuming each spectrum differs from the template only through a certain Doppler shift. The posterior RV is derived from each observation using a Markov chain Monte Carlo (MCMC) algorithm. This technique makes it possible to obtain reliable RV uncertainties. The s-BART time series presents a median RV uncertainty of 1.28 m/s and a rms deviation of 4.03 m/s.

2.5.3 YARARA

Finally, the RV extraction was performed with YARARA (Cretignier et al. 2021). YARARA is a post-process applied to the spectra derived by the DRS meant to further clean them from instrumental and Earth atmosphere signatures. The spectra are thoroughly corrected for cosmic rays, tellurics, ghosts, fringing, and instrumental focus variations. Additionally, YARARA can also correct for some stellar activity features via a multi-linear regression approach with the S-index, for each wavelength of the spectral time series.

However, this reduction step did not lead to a smaller scatter in the RV time series of TOI-1453, and no clear activity signal was detected. Therefore, we omitted this correction. The merged-order 1d spectra were continuum-normalised with RASSINE (Cretignier et al. 2020b), and the RVs were extracted from the CCF computed on a careful line selection (Cretignier et al. 2020a). Together with the RVs, we also extracted various stellar activity indicators such as the full width at half maximum (FWHM) of the CCF, the H_α, and the log $\left({R^{\prime }}_{HK}\right)$. The YARARA RV time series has a median RV uncertainty of 1.04 m/s and a rms deviation of 2.73 m/s, which are significantly reduced compared to the DRS due to the powerful correction of additional systematics and improved selection of stellar lines to compute the CCF.

Because both the RV uncertainties and rms are the smallest with the YARARA time series, we use the latter in the analyses presented in this paper. Among the 100 observations, 6 were rejected due to anomalous CCFs, leaving us with 94 measurements. The data have been made available via the Centre de Données astronomiques de Strasbourg (CDS) and Table A.2 illustrates the format of this dataset with the five first entries.

3 Stellar characterisation

In this section we describe our analyses to characterise TOI- 1453. We refer to Table 1 for the adopted stellar parameters.

3.1 Astrometry

According to the third Gaia data release (Gaia Collaboration 2016, 2023), the star draws a yearly parallax motion of 12.676±0.010 milliarcsec. This corresponds to a distance of 78.89±0.06 pc from the Sun. Coupling this distance measurement and the systemic RV with the measured proper motions (µ_α, µ_δ) = (–70.974, –0.085) mas yr⁻¹, we calculated the U_LSR, V_LSR, W_LSR velocity components of the star in the Galaxy with respect to the Local Standard of Rest (LSR). This estimation includes the correction of the solar motion in the LSR, which we took from Dehnen & Binney (1998). We find reasonably small velocities: (U_LSR, V_LSR, W_LSR) = (–58.99, –5.33, –10.61) km/s. Based on those velocities, we computed the probabilities that TOI-1453 belongs to the thin disc, thick disc and halo following the work of Reddy et al. (2006). We obtained P_thin =98.5%, P_thick =1.5%, P_halo < 0.1%. Our analysis of TOI-1453 kinematics hence indicates that the star belongs to the thin Galactic disc.

As was shown in Sect. 2.1, TOI-1453 has a close stellar companion with a magnitude difference of ΔG = 5.24. They share similar proper motions, according to the Gaia data, which suggests that they form a stellar binary system. As a result, TOI-1453.01 and TOI-1453.02 are classified as S-type planet candidates, meaning that they presumably orbit a single component of a multi-star system. The angular separation between the two stellar components is estimated to 1.9″, which translates into a projected separation of ∼150 AU given the distance estimated from the Gaia parallax. In their recent work, Christian et al. (2024) estimated the mass of the stellar companion to TOI-1453, based on an empirical relation between the mass and the magnitude in the K band for low-mass stars (Mann et al. 2019). They found a mass of 0.238±0.029 M_⊙. They also performed a fit of the Gaia astrometric data, from which they constrained the semi-major axis of the binary to be ${171}_{-58}^{+127}$ AU and a (poorly constrained) orbital eccentricity of ${0.63}_{-0.39}^{+0.32}$. Additionally, they measured an orbital inclination of ${91.98}_{-1.27}^{+3.86}$ deg, and concluded that the orbit of the binary is aligned with the planetary orbits. These results have important implications for the disc-binary interactions that took place in this system. In Sect. 4.2, we discuss our analyses aimed at identifying the true host to the transits.

Based on the mass and orbital elements of the stellar companion to TOI-1453, we could estimate its expected impact on the RV time series to be ∼2.3 ms⁻¹ over the timespan of the HARPS-N observations, under the hypothesis of a circular orbit. We fit the RV data with a linear trend and obtained a slope of –0.37 ms⁻¹/year, i.e. less than 1 ms⁻¹ drift over the timespan of the HARPS-N observations. This contrast between the observed and expected drifts is likely due to the large orbital eccentricity which induces asymmetric RV variations. Another possibility is the existence of an underlying additional RV drift that counterbalances the effect of the binary. This additional drift could be instrumental or stellar in origin. As HARPS-N is known to be very stable on such timescales, we rejected the instrumental hypothesis. A drift due to the stellar magnetic cycle would have an opposite direction and similar amplitude to the drift due to the binary. The combination of these conditions makes the stellar origin less likely than that of an eccentric binary orbit. The high eccentricity is supported by the work of Christian et al. (2024) and would create a significantly smaller RV variation along most of the orbit.

3.2 Stellar atmospheric parameters

The high-resolution HARPS-N spectra were used to constrain various stellar parameters. To proceed, we employed three independent characterisation techniques: SPC, ARES+MOOG, and CCFPams. Ultimately, we combined the results from the three methods in order to obtain final stellar parameters that account for model systematics. This process of stellar characterization was described in detail in Mortier et al. (2020). Here we briefly synthesise the three methods.

The Stellar Parameter Classification tool (SPC, Buchhave et al. 2012, 2014) is a spectrum synthesis method that works on the individual spectra and merges the results into a weighted average. Each HARPS-N spectrum is crosscorrelated with a library of model template spectra. The stellar parameters for TOI-1453 derived by SPC are T_eff = 5044±50K, log g = 4.68±0.1, [m/H] = –0.21±0.08. Additionally, the SPC method constrains the projected rotational velocity of the star to υ sini < 2 km/s.

The ARES+MOOG method is a curve-of-growth technique based on a line list of neutral and ionised iron lines (Sousa 2014). The equivalent widths of those lines were computed on the HARPS-N stacked spectra using ARESv2 (Sousa et al. 2015)³. These results are then injected into the radiative transfer code MOOG to extract parameters such as T_eff, log g and iron abundance [Fe/H]. We note that the HARPS-N wavelength range contains a large number of iron lines. Therefore, we consider the measurement of [Fe/H] to be equivalent to the overall metal- licity of the star, [m/H]. Following this procedure, we obtained T_eff = 4932±86K, log g = 4.87±0.16, [Fe/H]=−0.35±0.05.

Finally, CCFPams (Malavolta et al. 2017)⁴ estimates T_eff, log g and the iron abundance from the measurement of the equivalent width of the CCFs, using an empirical calibration relation. This technique provides T_eff = 4948 ±68 K, log g = 4.693±0.228, [Fe/H] = –0.35±0.05. The surface gravity that we report includes a correction factor coming from a survey of transit-constrained log g, as presented by Mortier et al. (2014).

Noteworthy, we found interstellar absorption in the HARPS- N spectra in the wide NaD lines around 5890 and 5896 Å. The absorption features are consistently blue-shifted by ∼15.5 km s⁻¹ with respect to the centre of the lines. We provide further details on this observation in Fig. B.1.

3.3 Stellar radius, mass, and age

We ran the isochrones code (Morton 2015) with the MESA Isochrones and Stellar Tracks set of isochrones (MIST, Dotter 2016) to estimate the stellar mass, radius, and age. The input parameters for MIST were the effective temperature and metal- licity, the Gaia parallax, as well as the stellar apparent magnitude in eight photometric bands: B, V, J, H, K, W1, W2, and W3. We ran the MIST code three times, for each of the three input T_eff and [Fe/H], with the purpose of alleviating systematic differences in stellar modelling (e.g. Tayar et al. 2022). Finally, we combined the three posteriors and report the results in Table 1. For the stellar mass and radius, we add a relative uncertainty of 4% in quadrature to the uncertainties reported by our isochrones fits, in order to account for the systematic differences between various stellar evolution models (Tayar et al. 2022).

In this table, we report the value for log gfrom our isochrones fits, however, we note that T_eff and [Fe/H] are better constrained from our spectroscopic analyses. For these two parameters, we report the weighted averages from the three methods: T_eff=4975±68K, [Fe/H]=−0.31±0.06. The results indicate that TOI-1453 is slightly metal-poor, while from our isochrones fits we estimate the age to be ${12.0}_{-3.8}^{+1.1}$ Gyr. Those results have to be considered in regards of the kinematics of the star in the Galaxy, which unambiguously revealed that the star belongs to the thin disc (cf. Sect. 3.1). In a recent work, Gallart et al. (2024) reconstructed the star formation history in the solar neighbourhood within 100 pc, using Gaia data. They showed that at this location, the disc of the Milky Way started to form ∼11 Gyr ago at roughly solar metallicity, and that the metallicity of the stellar population followed a slightly decreasing trend until ∼6 Gyr ([m/H] in-between [–0.5,0]). The metallicity, age, and kinematics estimations of TOI-1453 are in agreement with this global picture.

3.4 Stellar activity

The HARPS-N spectra provide us with valuable information about the level of activity of TOI-1453 during the observations. We extracted the YARARA log $\left({R^{\prime }}_{HK}\right)$time series. To proceed, we fitted a profile of plage on the CaII H and K lines (see Cretig- nier et al. 2024b) but adapted for the K spectral type (see Fig. 10 in Cretignier et al. 2024a), which leads to a precise but inaccurate S-index time series. Then, we calibrated this time series via a linear transformation on the $\log \left({R^{\prime }}_{HK}\right)$ estimated from the DRS and our log S-index. Finally, we used the coefficients of this transformation to build a precise and accurate YARARA $\log \left({R^{\prime }}_{HK}\right)$ time series. The latter indicates $\text{log\hspace{0.17em}}\left(R_{HK}^{\prime }\right)=-4.91\pm 0.06$, pointing towards a low activity level. We applied the empirical relation of Noyes et al. (1984) to obtain a rough estimate on the stellar rotation period. We obtained P_rot=42.9 days, while we undertook the same process for the relation of Mamajek & Hillenbrand (2008) and found P_rot =43.8 days.

We report in Fig. 3 the time series, periodograms and correlation plots of the YARARA RVs together with a few spectroscopic activity indicators: $\log \left({R^{\prime }}_{HK}\right)$, the FWHM and bisector span of the CCF, and H_α . The plots reveal little to no correlation between our RV time series and the activity indicators. Additionally, we observe that the RV time series is contained within ± 10 m s⁻¹ . To compute the periodogram of $\log \left({R^{\prime }}_{HK}\right)$, FWHM and H_α , we first subtracted a quadratic trend from the time series. The concave shape of these trends suggests a phase of low stellar magnetic activity. We observe a significant signal at 40.3 days in the residual periodogram of $\log \left({R^{\prime }}_{HK}\right)$, and 19.0 days in the periodogram of FWHM. The strongest (but not significant) periodicity in the RVs is 17.7 days.

In an attempt to further constrain stellar rotation in the HARPS-N data, we used SPLEAF (Delisle et al. 2020, 2022) to fit a multi-dimensional Gaussian Process (GP) simultaneously on the RVs and $\log \left({R^{\prime }}_{HK}\right)$. To proceed, we employed the Exponential-sine periodic (ESP) kernel, which is an approximation of the squared-exponential periodic kernel. It is an efficient kernel to capture stellar rotation features (e.g. Stalport et al. 2023). We fitted the GP together with a Keplerian around P=6.59 days to account for the signal of the outer planet candidate (there is no serious hint of the inner candidate in the RV data alone). However, we could not reach a satisfying convergence despite a vast exploration of different priors. While the period of the GP consistently peaked at ∼39 days, the other hyperparameters were strongly unconstrained. We also tried to focus the analyses only on the second season of observations, which shows a slightly higher activity level. The results were similarly inconclusive. We also tested the multi-dimensional GP on FWHM, bisector span, and H_α, all leading to bad convergence. Finally, we tested the simpler case of a GP applied on the RV only, which might ease the convergence when stellar activity is weak. Again, we fitted this GP simultaneously with one Keplerian. We were unsuccessful in constraining the hyperparameters, this time including also the period of the GP. These results confirm the low activity of the star during our HARPS-N observations.

We analysed the TESS light curve in search for periodicity in the photometry. To proceed, we downloaded the SAP fluxes and removed the flux measurements strongly contaminated by the Earth or Moon scattered light. Thanks to the extensive timespan of the TESS observations, we defined three datasets roughly corresponding to the three cycles of observations: sectors 14 to 26 (12 sectors), sectors 49 to 60 (9 sectors), and sectors 73 to 76 (4 sectors). For each dataset, we undertook a period search by computing generalised Lomb-Scargle (GLS) periodograms. We present the results in Fig. 4. This analysis shows an increasing activity over the timespan of the observations, as is also hinted in the $\log \left({R^{\prime }}_{HK}\right)$ time series from the HARPS-N observations. Beyond this trend, we have not observed any significant periodicity. Noteworthy, the main periodicity in the last cycle of TESS data peaks at 39.4d, similar to the period found in the $\log \left({R^{\prime }}_{HK}\right)$. However, more than 800 days separate these two datasets, so this should not help our multi-dimensional GP converge.

Our analyses suggest a tentative stellar rotation period around 40 days. However, we remind that this result is based on observations taken during a phase of low stellar activity, which affects the significance of our conclusion.

4 Data analysis

4.1 Photometric analyses

We decided to utilise exclusively the TESS data for our photometric analyses, and to not include the ground-based LCOGT photometry or the other photometric datasets listed in Sect. 2.4. The reason is that the planetary transits in the ground-based data are very shallow (the inner planet is not detected, the outer planet is at the limit of detection). Therefore, these data do not help constrain the transits shape and depth. Adding such ground-based photometry can still be very useful when it significantly extends the observational baseline of the system, which can serve to refine the orbital periods. However in our case, the ground-based photometry does not extend the timespan of the extensive TESS dataset.

We started our planet analyses by searching for the transit signals in our extracted TESS light curves. Using the Python package Transit Least Squares (TLS, Hippke & Heller 2019), we consistently found the transit signals of the two planet candidates throughout the TESS observations. The linear ephemerides obtained with this preliminary search were then used as priors to fit a two-planet model to the light curve with Juliet. Juliet models photometric data with the batman package (Kreidberg 2015). Our model consists of the sum of two transit components and a noise component. In our case, the noise component simply consists of a jitter term σ_TESS added in quadrature to the photometric measurements (this noise model is satisfactory as we explained in Sect. 2.1).

The transit model for each planet is parametrised with the orbital period, P, the reference time at inferior conjunction, T₀, the orbital eccentricity, e, the argument of periastron, ω, the relative radius, R_p/R_⋆), and the impact parameter, b. The latter two parameters were included via the (r1, r2) parametrisation suggested by Espinoza (2018), and with uniform priors in [0,1]. We fixed e to 0, and set Gaussian priors on P and T₀ based on the results of TLS (for planet b: N(4.31, 0.1) and N(2459392.91, 0.1), for planet c: N(6.59, 0.1) and N(2459395.37, 0.1), where N indicates a normal distribution). In this transit model, the stellar density ρ is used together with P to constrain the scaled semi-major axis a/R_⋆ of each planet, while reducing the number of parameters – from the two parameters a_b /R_⋆ and a_c /R_⋆, we now only fit for the stellar density (Espinoza et al. 2019). We set priors on ρ according to our stellar characterisation analyses – cf. Sect. 3. The transit model also contains instrument-dependent parameters such as an out-of-transit flux offset, δ_flux (prior: N(0, 0.1)), a dilution parameter, D, and the limb darkening (LD) coefficients. The dilution parameter accounts for the light contamination from the neighbouring stars, which results in an underestimation of the transit depth. This parameter is defined in Juliet as: $$D=\frac{1}{1+{\displaystyle {\sum }_k{F}_k}/F_T},$$

where F_k is the flux of each star k in the photometric aperture except the target star, whose flux is F_T. Given that we used the same aperture than the SAP and PDCSAP to extract the light curve, the dilution factor is simply given by the crowdsap parameter indicated in the FITS files of the TPF. This metric provides the ratio between the target flux and the total flux inside the aperture, which is equivalent to D. Its average value over all the sectors is 0.98435. To account for the effect of LD, we used a quadratic LD law with coefficients q₁ and q₂ defined in Kipping (2013). The stellar parameters were used to constrain those coefficients for the TESS passband with the help of LDCU⁵. LDCU is a modified version of the Python routine implemented by Espinoza & Jordán (2015) that computes the LD coefficients and their corresponding uncertainties using a set of stellar intensity profiles accounting for the uncertainties on the stellar parameters (e.g. Deline et al. 2022). The stellar intensity profiles are generated based on two libraries of synthetic stellar spectra: ATLAS (Kurucz 1979) and PHOENIX (Husser et al. 2013). The code takes as input various stellar parameters T_eff, log g, [m/H], which we estimated in Sect. 3. For each intensity profile, LDCU performs three fits for different atmosphere model data selection processes (all, limb-excluded, and interpolated data). To obtain the most conservative LD estimations, we extracted the merged distributions from all the intensity profiles, and further increased the uncertainties by a factor ~3. We used these results to set Gaussian priors on the LD coefficients (N(0.3911, 0.05) and N(0.3358, 0.05) for q₁ and q₂, respectively).

Having defined the model, we proceeded to fit the data from 29 TESS sectors. The model parameters were explored with the nested sampling algorithm dynesty using a number of live points n_live ≥ N², where N is the number of model parameters to be explored (not accounting for the fixed parameters). We assumed that the fit had converged when the variation in logevidence from one iteration to the next went below 0.01. In Fig. 5, we present the TESS light curve together with the best-fit transit model.

Fig. 3 HARPS-N YARARA RVs, $\log \left(R_{HK}^{\prime }\right)$, CCF FWHM, bisector span, and Hαtime series. Left: the various time series. Middle: GLS periodograms, together with the 10, 1, and 0.1% FAP levels. Right: correlation plots with the RVs. In the plots, we also indicate the Spearman correlation coefficient.

Fig. 4 GLS periodograms of the TESS SAP photometry separated into three datasets: sectors 14 to 26 (light blue), sectors 49 to 60 (dark blue), and sectors 73 to 76 (gray). In the legend, the periods at which the power is maximum are indicated inside parentheses for each dataset.

Fig. 5 TESS light curve of TOI-1453 (29 sectors comprised between sector 14 and sector 76) together with the transit model in black. (a) Full TESS light curve. The zone highlighted in purple illustrates the timespan of the HARPS-N observations. The three deeper transit features in the model correspond to simultaneous transits of planets b and c. (b) TESS light curve phase-folded on the transits of TOI-1453 b (TOI-1453.02). The orange dots are binned TESS data for better visibility. (c) Same for TOI-1453 c (TOI-1453.01), where the green dots are binned TESS data.

4.2 Validation of the planet candidates

To validate the two planets as transiting TOI-1453, we explored all the possible scenarios that could create our observations. To proceed, we used TRICERATOPS⁶ (Giacalone & Dressing 2020; Giacalone et al. 2021). This publicly available code employs a Bayesian probabilistic approach to test the astrophysical false positive hypotheses such as nearby eclipsing binaries. TRICERATOPS takes as input the normalised light curve phase-folded onto the transit feature of interest. We phase-folded the entire TESS light curve onto the period and phase of TOI-1453.01, and binned the resulting data into 400 data points to increase the computation speed. Additionally, we also provided the contrast curve of the speckle image in the 562 nm narrow photometric band from Gemini, which helps constrain the false positive probabilities (FPPs). We ran the TRICERATOPS test for TOI-1453.01, and obtained FPP=0.0023±0.0017. We repeated the same process for TOI- 1453.02, and found FPP=0.0211 ±0.0020. While this test validates TOI-1453.01 as a planet orbiting TOI-1453, namely TOI- 1453 c, the result is not so clear for TOI-1453.02 since the FPP is slightly above 1%. In particular, we cannot firmly rule out the nearby stellar companion 1.9″ away as the source of the transit. We note however that this scenario seems unlikely, because that would imply a transit depth of ~ 10% on this faint star and subsequently a large planet around a small star. In contrast, the scenario of two small planets around TOI-1453 is more likely. As was already noted in the context of the Kepler space mission, the probability of FP decreases sharply in multi-planet candidate systems compared to singles (Lissauer et al. 2012, 2014). These studies revealed a correction factor larger than 10 on the FPPs between single and multi-planet candidate systems. Taking the system multiplicity into account in our FPP estimates above, we conclude that both planet candidates present FPP below 1%.

To proceed further, we performed a second validation test based on the computation of the stellar density. In Table 1, we report our estimation of ρ_⋆ based on our fits of isochrones constrained from the HARPS-N spectra (cf. Sect. 3). Additionally, the transit shape independently constrains the stellar density as well. The comparison between the transit-based and the model-based stellar densities represents a good test of the stellar host. The stellar companion to TOI-1453, DR3 1433062331333157632, is much fainter than TOI-1453. Provided they both are at the same distance from the Sun, this implies that the spectral types of the two stars are also very distinct. From the mass estimate of Christian et al. (2024), M_B = 0.24 M_⊙, we expect an effective temperature T_eff<3500 K. At such small T_eff, the stellar density is expected to be larger than TOI-1453 by a factor of at least 5 (Baraffe et al. 2015). The densities of the two stars being significantly different from each other, their comparison to the transit-based densities could bring valuable information.

The stellar density can be derived from the transit parameters following the formula (Seager & Mallén-Ornelas 2003): $${\rho }_{\star }=\left(\frac{4{\pi }^2}{P^{2}G}\right){\left[\frac{{\left(1+\frac{R_p}{R_{\star }}\right)}^2-b^2\left(1-{\sin }^2\left(T_{14}\pi /P\right)\right)}{{\sin }^2\left(T_{14}\pi /P\right)}\right]}^{3/2},$$

where P is the orbital period of the considered planet, R_p and R_⋆ are the planetary and stellar radii, b is the transit impact parameter, and T₁₄ is the total transit duration. Therefore, we performed a fit unconstrained by the stellar priors so to measure the stellar density from the transit parameters. We set wide uninformative priors on the stellar density and LD coefficients (log U[500, 10 000] kg m⁻³ and U[0,1] for ρ_⋆ and LD coefficients q₁, q₂, respectively), and fixed the dilution factor D to 1. From the posteriors of that analysis, we computed the stellar densities based on the transit parameters of TOI-1453.01 and TOI-1453.02, respectively. We found ${\rho }_{\star }(.01)={1.04}_{-0.84}^{+1.13}{\rho }_{\odot }$ and ${\rho }_{\star }(.02)={1.94}_{-1.61}^{+1.66}{\rho }_{\odot }$. These two stellar densities are compatible with the model-based stellar density of TOI-1453 within 1σ. Furthermore, they are not compatible with the expected high density of the stellar companion DR3 1433062331333157632, which validates both planets. In the following, we will refer to TOI-1453.02 as TOI-1453b, and TOI-1453.01 asTOI-1453c.

4.3 Joint TESS-HARPS-N analysis

The results of the transit fit described in Sect. 4.1 were used to constrain a global fit to the transits of each planet and the RVs. Notably, we narrowed down the priors on P and T₀, and restricted the relative radius exploration of both planets to values within the range R_p/R_⋆=[0.01, 0.05] so as to reduce the computation cost⁷. Regarding the RVs, we included the YARARA data presented in Fig. 3. The addition of the RV data leads to the inclusion of four more parameters in the global model: the RV semi-amplitude K_b and K_c of planets b and c, an RV offset µ_HARPSN and a jitter term µ_HARPSN for the HARPS-N data. We placed wide priors on these parameters. We also tested the inclusion of a RV linear drift in the model, but both slope and intercept were compatible with zero so we decided to not include this trend in the final fit. The set of parameters and priors used in our joint fit is described in Table C.1.

While GP detrending of YARARA RVs can efficiently account for stellar red noise (e.g. Stalport et al. 2023), our observations do not reveal significant stellar activity, as already discussed in Sect. 3.4. To further test whether GP detrending was needed, we compared joint transit+RV fits with different input RV data using Juliet. In one case, we input the original YARARA RVs, in the other case we use the RVs detrended from the multi-dimensional GP $\left(\text{RV},\log \left(R_{HK}^{\prime }\right)\right)$. We employed the same two-planet model in both fits with uniform priors on the RV semi-amplitudes between 0 and 100 m/s, and we compared the posterior planet parameters after a Bayesian exploration of the parameter space with dynesty (to decrease the computation time, we performed those analyses on a reduced TESS dataset composed of the sectors 24, 25, 26, 40, 41). The RV semiamplitudes have negligible differences between the two fits, with ∆K ~ 0.02 m/s for both planets b and c. Therefore, there is no measurable impact of our GP correction on the planet masses. Additionally, we obtain slightly smaller uncertainties on K_b and K_c , and hence higher-significance detections, for the fits that use the original YARARA datasets without GP detrending. In conclusion, for the global joint transit+RV fit we opt for the original YARARA RVs, without further detrending.

We undertook a joint fit of the HARPS-N RVs and 29 TESS sectors. We explored the model parameters with dynesty using the same constraint on the number of live points (n_live ≥ N²) and the same criterion on convergence than with the fit on the photometric data alone. The results of our global joint fit are reported in Table C.1. We obtain nearly identical orbital inclinations between both planets. This is an expected outcome of tidal dissipation, which aligns the planetary orbits with the stellar spin given the stellar rotation period is not small compared to the orbital periods (e.g. Bolmont et al. 2014; Ogilvie 2020). Therefore, the observed orbital alignment further supports both planets orbiting the same star. Regarding LD coefficients, the posteriors are essentially identical to the priors, despite the conservative approach from LDCU. This highlights here that the LD coefficients are essentially constrained from the stellar parameters rather than the transits.

Thanks to the extensive TESS dataset, we constrain well the transit depths of the inner and outer planets to 224± 12.5 ppm (2.8% uncertainty on R_b/R_⋆) and 802±20.7 ppm (1.3% uncertainty on R_c/R_⋆), respectively. With our set of HARPS-N observations, we measure RV semi-amplitudes of K_b=0.61 m/s (with 1.9σ significance), K_c=1.26 m/s (non-zero with 3.5σ significance), respectively. We present the best fit RV curves of TOI-1453 b and c in Fig. 6. Therefore, the RV detection of the inner planet remains insignificant with our current dataset. However, we set reliable constraints on the mass of the outer planet: $M_c={2.95}_{-0.84}^{+0.83}$ M_⊕. With an estimated radius of 2.224±0.095 R_⊕, we conclude that this sub-Neptune presents a very low bulk density of ${0.267}_{-0.080}^{+0.087}$ρ_⊕.

Fig. 6 HARPS-N RV phase-folded on the results of the joint fit for both TOI-1453 b (left) and TOI-1453 c (right). The black line represents the best-fit model RV curve of each planet.

5 Dynamical state and TTVs

5.1 Proximity of the planets to the 3:2 MMR

The planets, with their orbital periods of P_b =4.314 days and P_c=6.589 days, present a period ratio of P_c/P_b=1.527, close to the exact 3:2 commensurability. Therefore, we carried out a dynamical analysis to assert if the planets are dynamically bound into a state of mean motion resonance (MMR). Furthermore, with these analyses we also investigated the stability of the planetary system. We did not include the stellar companion to TOI-1453 in our dynamical analyses. Indeed, our investigations focus on short-term dynamics (and short-term instability), enough to reveal MMR. Within these short timescales, the stellar companion only has a marginal effect on the MMR structure (e.g. Stalport et al. 2022b).

Chaos maps are powerful tools to visualise MMR in a planetary system (e.g. Goździewski et al. 2001; Correia et al. 2009; Stalport et al. 2022b; Daquin & Charalambous 2023). These are 2-dimensional sections of the parameter space, inside which the dynamical chaos is evaluated. We focused our investigations on the short-term chaos which is a direct outcome of the interactions between the MMR (Chirikov 1979). To estimate the amount of chaos in a planetary system, we coupled the results of numerical evolution with the Numerical Analysis of Fundamental Frequencies (NAFF) chaos indicator (Laskar 1990, 1993). NAFF is a process used to precisely compute the orbital mean motion of a planet based on a time series of its mean-longitudes over time. From this time series, a refined Fourier frequency analysis is performed in order to identify the main frequency of the orbital evolution, which is associated with the mean motion. Separating the mean-longitudes time series in two halves, we can apply this technique on each half of the integration and compute the quantity: $$\text{NAFF}=\underset{j}{\max }\left[{\log }_{10}\frac{\left|n_{j,2}-n_{j,1}\right|}{n_{j,0}}\right],$$

where n_j,1 and n_j,2 are the mean motions of planet j over the first and second half of the integration, respectively, and n₀ is the initial mean motion of that same planet and serves as a normalisation factor. The temporal variation of the mean motion offers direct information on the amount of chaos in the planetary orbit. Indeed, the classical (non-chaotic) theory of dynamics of Laplace and Lagrange predicts constant orbital mean motion over secular timescales. This is no longer true in the presence of chaos. Therefore, the drift in the orbital mean motion n_j,2 − n_j,1 is a measure of chaos in a planetary orbit. We estimated the chaos in each of the planetary orbits, j, and defined the level of chaos in the whole system as the level of chaos of the planetary orbit that displays the highest degree of variation.

In practice, we generated a grid of 201×201 systems exploring the initial period ratio P_c / P_b and the orbital eccentricity e_c of the outer planet, while all the other parameters had their initial conditions fixed to the best-fit solution (cf. Table C.1). We simulated the evolution of each system in the grid with the adaptive time-step high-order integrator IAS15 (Rein & Spiegel 2015) included in the Python package REBOUND⁸ (Rein & Liu 2012). The timespan of each simulation was set to 10kyr (i.e. about 554 000 revolutions of planet c). We included an orbital correction due to General Relativity (GR) effects following the framework of Anderson et al. (1975) and implemented in REBOUNDx⁹ (Tamayo et al. 2020b). During these simulations, the planet mean longitudes were evenly recorded to build a time series of twenty thousand values. These time series were split in two halves, from which we computed the mean motion of TOI-1453 b and c on each half and derived the NAFF indicator. This process was repeated for each of the 40401 system configurations.

In Fig. 7, we present the resulting grid of system configurations. The colour is indicative of NAFF, with redder regions denoting more chaotic systems. White areas indicate strongly unstable configurations, that did not reach the end of the numerical simulation (escape or close encounter as defined by the mutual Hill radius between planets b and c). Higher eccentricity orbits naturally lead to stronger dynamical chaos and potentially short-term instability. The larger the period ratio, the weaker the planet-planet gravitational interactions. As such, towards the right end of the map, the systems harbour a weakly chaotic behaviour for a wider range of initial orbital eccentricities. This is illustrated in Fig. 7 by the positive slope of the limit between blue and red, and red and white systems. The vertical line informs about the position of the system on the x-axis, which is well-constrained by the transit observations. At that position, we observe a rough upper limit of e_c<0.12 to keep a weak chaos. We stress however that this limit is approximate. While calibrating NAFF on the actual orbital stability is possible (Couetdic et al. 2010; Stalport et al. 2022a), the purpose of chaos maps is to highlight the dynamical topology of the parameter space. The blue areas are feasible solutions, while the other more chaotic configurations will ultimately lead to orbital instability within short timescales. There exists rigorous dynamical techniques to refine the planet parameters based on reliable parameter explorations (e.g. Tamayo et al. 2020a; Stalport et al. 2022a), but this is beyond the scope of the present investigation. The 3:2 MMR appears in Fig. 7 as the vertical V-shape feature of stable systems extending to higher eccentricities. As this chaos map shows, the planets in TOI-1453 lie outside of the 3:2 MMR.

Fig. 7 Chaos map of the TOI-1453 system, with a resolution of 201×201. This map estimates the chaos in a 2D section of the parameter space, along Pc/Pb (x-axis) and ec (y-axis). The colour code depicts the NAFF chaos estimation, where redder areas are more chaotic. The white regions were unstable before the end of the short-term numerical simulations (escape or close encounter). The actual period ratio of the planet pair constrained from our photometric fit is represented by the vertical line.

5.2 TTV investigations

In an attempt to further constrain the mass of TOI-1453 b, we used the large TESS dataset to explore potential TTVs of the outer planet TOI-1453 c. Due to the gravitational interactions between the two planets, their transit times should deviate from the pure Keplerian ephemerides with an amplitude that depends chiefly on the eccentricities, the mass ratios, and the orbital periods. While the transits of TOI-1453 b are too shallow for measuring their individual timings, we investigated potential TTVs on the transits of TOI-1453 c as a result of the planet-planet gravitational interactions.

The planets are close to, but outside of the 3:2 MMR. Therefore, their TTV amplitudes should be modest. We carried out numerical simulations to estimate the expected TTV amplitudes in the TESS dataset. To proceed, we explored the influence of the mass and the orbital eccentricity of the inner planet on the TTVs of the outer one, via the design of a 101 × 101 grid of simulations. For each combination of parameters, we used TTVFast (Deck et al. 2014) to integrate the system on the timespan of our TESS observations (i.e., 1711 days) with a time step of 0.05 day. We verified the validity of this time step by performing the same integration with smaller time steps, and obtained similar results. During the integration, we recorded the transit times of TOI-1453 c and fitted this time series with a line, the slope of which gives the orbital period of the linear ephemerides. We retained the maximum deviation among the transit times from this linear ephemerides, and repeated this process for each point in the grid. The resulting map of simulated TTV amplitudes is presented in Fig. 8. While a large range of TTV amplitudes is observed throughout the grid, most of the parameter space is consistent with TTVs smaller than 20 minutes in amplitude. Notably, within the 68.27% range of masses M_b derived from our joint fit, we observe small amplitude (≲20 min) TTV up to an orbital eccentricity of e_b ≲ 0.15.

With these prospects in mind, we investigated TTVs in the TESS data via a per-sector fit of the transit times. Our transit model includes TOI-1453 b with linear ephemerides, while the outer planet has varying transit times. Indeed, instead of providing the parameters P_c and T_0c, we added each transit time as a free parameter. Therefore, we included as many parameters as there are transits of planet c in the considered TESS sector. For both planets, we used the results of our joint global fit to impose narrow Gaussian priors on the r1 and r2 parameters, as well as P_b and T_{0, b}. Our analyses were carried out with Juliet and the dynesty nested sampling algorithm. We undertook a two-step process to precisely constrain the transit times of TOI-1453 c. For each sector, we performed a first fit with wide Gaussian priors on the individual transit times (0.1d of standard deviation), with the mean of the Gaussian derived from the linear ephemerides prediction. These results were used to narrow down the priors on the transit times for a second fit. We retained the results of that second fit to estimate the individual transit times. At the end of this process, we obtained the individual transit times of TOI-1453 c in all the sectors. In Fig. 9, we present the deviation of each transit time with respect to the linear ephemeris, well-known as the O-C diagram (Observed timings minus the Calculated ones). It does not show large variations in the transit times: if any, their amplitudes are smaller than 20 minutes. This is in line with the results of our dynamical simulations.

Even if the amplitude of the TTV signal is expected to be small, we checked if the time series of O–C transit times could harbour a significant periodicity. We computed the GLS periodogram (GLS; Zechmeister & Kürster 2009) of these transit times, and find the absence of significant periodicity as shown in the middle panel of Fig. 9. The highest peak slightly exceeds a False Alarm Probability (FAP) of 10% only. However, due to the proximity of the system to the 3:2 MMR, we note that one of the peaks lies close to the expected TTV super-period associated with this resonance. The TTV super-period is a longterm periodic variation of the TTV pattern often dominant in amplitude, and caused by the proximity of a pair of planets to a MMR. Close to a j:k MMR, the angles jλ_c − kλ_b − (j − k)ϖ_b and jλ_c − kλ_b - (j − k)ϖ_c vary slowly, where λ and ϖ depict the mean longitude and the longitude of periastron, respectively (in the co-planar case, the longitude of periastron is equivalent to the argument of periastron ω). Consequently, the two planets regularly end up at the same relative position with respect to each other, which creates a constructive effect. Due to these strong mutual gravitational interactions, the planets will deviate significantly from their Keplerian orbits on a timescale given by the TTV super-period (e.g. Agol & Fabrycky 2018): $$P_{\text{TTV}}=1/\left|j/P_c-k/P_b\right|.$$

The further away the planet pair from the resonance, the shorter in time and smaller in amplitude is this effect. In the case of TOI- 1453, the planet pair has an expected super-period of 120 days due to its proximity to the 3:2 MMR. This is represented by the vertical dotted line in Fig. 9.

To further investigate the tentative TTV periodicity close to the expected super-period, we also analysed the 20-second cadence TESS photometry when available. Our motivation was to attempt to reduce the uncertainties on the transit times and maybe increase the significance of the observed signals. However, the 20 s TESS photometry is not available for all the sectors. Particularly, the data products are available in 20s cadence for the following sectors: 40, 41, 49, 50, 51, 52, 53, 54, 56, 57, 59, and 60 (12 sectors out of 29). For each of them, we repeated the process outlined above to retrieve precise transit times, and added these timings to the sectors that only contain 120s data. The new periodogram is presented in Fig. 9 as a dashed line. It unambiguously shows that the significance of the signal decreases close to the TTV super-period. Additionally, we note that all the peaks are located below the FAP 10% line, making them particularly insignificant. In conclusion, with our current dataset we do not detect significant TTVs on TOI-1453 c. The use of advanced TTV detection and characterisation techniques, such as RIVERS (Leleu et al. 2021), could potentially unveil a TTV signal and constrain the mass of TOI-1453 b.

The mass of TOI-1453 b is likely smaller than TOI-1453 c, given its smaller radius and non-detection in the RV data. Therefore, the TTV amplitude on TOI-1453 b due to planet-planet interactions should be larger than on TOI-1453 c. As we already indicated, the individual transits of TOI-1453 b are too shallow to constrain their TTVs. As a last investigation, we explored an alternative route by stacking the transit times of TOI-1453 b in all the pairs of sectors, and searched for potential orbital period variations as a proxy for TTVs. To proceed, we performed linear ephemerides fits on all the pairs of consecutive sectors (sectors 14–15, 16–17, 19–20, 21–22, 23–24, 25–26, 40–41, 49–50, 51–52, 53–54, 56–57, 59–60, 73–74, 75–76) to derive a total of 14 period estimates. With datasets of two sectors, we can stack enough transits of the inner planet to obtain a robust period estimation. For each pair of sectors, we undertook a two-step fitting procedure again. First, we carried out a fit with a wide prior on P_b (normal prior with standard deviation of 0.1 day) and centred on the linear ephemerides expected value obtained from our global joint fit. Secondly, we used these results to set narrow priors on P_b for a new fit, and saved the results of this second fit. We present the time series of the period variations ∆ P_b compared to the linear ephemerides in the bottom panel of Fig. 9. Throughout the TESS observations, we observe a maximum point-to-point period variation of 7 min. Compared to the linear ephemerides of our global fit, we do not observe period deviations larger than 4 min. Furthermore, there is no visible correlation between these small variations and the TTV diagram of the outer planet. Hence again, we do not find significant TTV signal. The absence of significant orbital period variations also highlights the robustness of our planet radius estimations. The TTV, if any, are not large enough to blur the signal of stacked transits.

Fig. 8 Grid of simulated TTVs of TOI-1453 c as a function of the mass of the inner planet Mb and its orbital eccentricity eb. The vertical black line depicts the best-fit mass estimate from our joint fit, while the surrounding white band represents the 68.27% confidence interval on this estimation.

Fig. 9 TTV searches in the TESS LCs. The top plot shows the O- C diagram of the individual transits of planet c, and the middle panel presents the GLS periodogram of those observations. The solid line is the periodogram derived from all the 2 min cadence data. The dashed line was obtained after replacing the 2 min by 20s cadence data when available, and proceeding to new transit fits (cf. the orange points in the top plot). We do not find significant periodicity. The bottom plot presents the variation of the orbital period of planet b as measured from pairs of consecutive sectors (with 2 min cadence).

Fig. 10 Posterior distributions of the interior structure of TOI-1453 c, showing mass fractions of the inner core, mantle, and envelope, as well as the water fraction in the envelope. The top row corresponds to a water prior compatible with a formation scenario beyond the ice line, while the bottom row assumes a water-poor composition, in agreement with formation inside the ice line. Colours represent three priors for the planetary Si/Mg/Fe ratios: stellar (purple, A1 and B1), iron-enriched (pink, A2 and B2), and uniformly sampled from a simplex (blue, A3 and B3). Dashed lines mark the median of each distribution, with priors indicated by dotted lines.

6 Interior modelling

We applied the publicly available Bayesian modelling framework plaNETic¹⁰ (Egger et al. 2024) to model the interior structure of TOI-1453 c. In this framework, the planet is modelled as a combination of an inner iron core with up to 19% of sulphur, a silicate mantle made up of SiO₂, MgO and FeO, and a volatile layer made up of a uniform mixture of H/He and water. The plaNETic code is split into an inverse model that samples possible internal structure parameters using the chosen priors, and a forward model, that calculates the transit depth of each sampled structure, which is then compared to the observed value. This process is sped up by replacing the forward model (the planetary structure model of the BICEPS code, Haldemann et al. 2024) with a neural network, which allows for a fast and reliable characterisation of the planetary interior.

We ran six different models assuming different priors. First, we vary the priors for the planetary Si/Mg/Fe values, assuming they are stellar (e.g. Thiabaud et al. 2015), iron-enriched compared to the host star (e.g. Adibekyan et al. 2021), or leaving them unconstrained by sampling uniformly from a simplex with an upper limit of 0.75 for the iron fraction. On the other hand, we also vary the priors for the water content in the planet, once assuming that water was readily available for accretion during the planets formation, and once that water could only be accreted through the accreted gas, compatible with a formation inside the water iceline.

Figure 10 shows the resulting posterior distributions of the most important interior structure parameters, while a summary of the results for the full set of parameters can be found in Table D.1. We find that, for both chosen water priors, we do not constrain the mass fractions of the inner core and mantle layers, for which the posteriors largely agree with the priors. For the water-rich prior (top row), we find that a wide range of envelope mass fractions is compatible with the observed planetary properties, with a preference for higher envelope mass fractions. Most of these inferred envelopes are water-rich, with a median envelope water mass fraction of around 90%. Conversely, for the water-poor prior (bottom row), we inferred tightly constrained metal-poor envelopes with mass fractions of ${0.3}_{-0.1}^{+0.2}\%,{0.6}_{-0.2}^{+0.3}\%$, and ${0.2}_{-0.1}^{+0.3}\%$, for the stellar, iron-enriched and free Si/Mg/Fe priors respectively.

Finally, we also applied plaNETic to model TOI-1453b. However, we could not set any constraints on its structure due to the large uncertainties onto the planet mass.

7 Discussion and conclusion

In this paper, we validated two planets orbiting TOI-1453 based on TESS observations, ground-based photometric, spectroscopic, and high-resolution imaging follow-up. We then confirmed the outer planet with 100 HARPS-N spectra. We further analysed the 29 TESS sectors and the HARPS-N data to characterise this system. We summarise our findings below.

• TOI-1453 is a metal-poor K-dwarf star with a mass and radius of 0.715±0.035 M_⊙ and 0.720±0.029 R_⊙, respectively. It belongs to the Galactic thin disc at a distance of ~78.9 pc. This star is orbited by a stellar companion with a current angular separation of 1.9″. This corresponds to a projected distance of ~ 150 AU. Christian et al. (2024) found that the binary orbit is eccentric and co-planar with the planetary orbits and constrained the semi-major axis of this binary to ${171}_{-58}^{+127}$ AU. TOI-1453 appears very inactive in the time span of our observations.

• TOI-1453 b orbits the host star with an orbital period of P_b =4.314 days. It is a super-Earth with estimated mass and radius of $M_b={1.24}_{-0.65}^{+0.66}$ M_⊕ and $R_b={1.173}_{-0.056}^{+0.059}$ R_⊕, respectively. Its bulk composition is poorly constrained given the currently low precision on the planet mass (1.9σ detection). However, it is likely rocky with a bulk density of ${\rho }_b={0.76}_{-0.40}^{+0.44}$ ρ_⊕·

• TOI-1453 c orbits the host star with an orbital period of P_c=6.589 days. It is a sub-Neptune with estimated mass and radius of $M_c={2.95}_{-0.84}^{+0.83}$M_⊕ and R_c=2.224±0.095 R_⊕, respectively. It harbours a very low bulk density $\left({\rho }_c={0.267}_{-0.080}^{+0.087}{\rho }_{\oplus }\right)$, compatible with either a water-rich composition or the presence of an extended H atmosphere on top of a rocky planet. TOI-1453 c is one of the rare sub-Neptunes below 3M_⊕ known to date.

TOI-1453 b and c span the radius valley, as illustrated in Fig. 11. These plots present the current sample of precisely characterised planets orbiting FGK stars (T_eff >4000 K), obtained with the DACE platform¹¹. Projected along the period (left) or incident flux (right) versus radius spaces, the planet population reveals a characteristic occurrence gap between 1.7 and 2 R_⊕ (Fulton et al. 2017; Van Eylen et al. 2018; Fulton & Petigura 2018). TOI-1453 b, represented by the orange star symbol, clearly stands below the gap. While TOI-1453 c (the green star symbol) stands above the radius valley close to the edge. It is noteworthy to observe that there are very few planets with radius ~1 R_⊕, and TOI- 1453 b occupies a yet-empty area. This is likely the result of the low S/N that these planets have around FGK stars; whereas significantly more detections were made around M dwarfs (e.g. Kemmer et al. 2022; Gillon et al. 2024). Constraining and understanding the composition of TOI-1453 c and the possible contrast with TOI-1453 b further are of particular importance. Currently, there are two interpretations for the population of sub-Neptunes. On one side, they could be ‘dry’, water-poor, puffy worlds. In this scenario, the super-Earths were formerly sub-Neptunes that lost their atmospheres due either to photo-evaporation (Owen & Wu 2013) or core-powered mass loss (Ginzburg et al. 2018). On the other hand, sub-Neptunes could be water-rich worlds, as a result of formation beyond the iceline followed by inward migration (e.g. Luque & Pallé 2022).

The proximity to MMR is indicative of migration in the protoplanetary disk. Chen et al. (2024) identified in the Kepler sample a larger proportion of sub-Neptunes close to a MMR compared to super-Earths. They use this observation and planetary formation models to propose that both ex situ and in situ formation mechanisms shaped the sub-Neptunes population (some formed beyond the ice line, others inside, following distinct formation routes). As was highlighted in Leleu et al. (2024), close-to-MMR sub-Neptunes are less dense. The bulk density of TOI-1453 c supports this trend, yet further detailed investigations are needed, such as atmospheric characterisation, to lift the compositional degeneracy. Additionally, stellar metallicity is expected to impact the planet compositions (e.g. Santos et al. 2017), with metal-poor stars expected to harbour planets with small core-mass fractions. While the internal structure modelling was unable to constrain the core mass fraction, the low metallic- ity estimated in this work could partially contribute to the small bulk density of TOI-1453 c.

The fact that TOI-1453 b and c are close to but still beyond a MMR is likely a consequence of tidal dissipation. At small orbital distances, tidal torques can move planet pairs away from the resonance by increasing their period ratio (e.g. Delisle et al. 2014; Louden et al. 2023). The Kepler sample supports this tidal MMR breaking. Indeed, Delisle & Laskar (2014) observe a net preference for MMR offsets in short-period planet pairs as opposed to pairs with P>15d. The planetary tides are expected to have synchronised the planet spins with their orbital frequencies to reach a state of tidal locking (e.g. Bolmont et al. 2014). However, the circularisation timescale can be significantly larger. There is no hint of orbital eccentricity in the TESS and HARPS-N observations. According to our dynamical analysis, any potential eccentricity should not exceed 0.12 to enable an orbital stability (cf. Fig. 7). Small orbital eccentricities are also expected from tidal dissipation, given the small orbital separation of the planets. This tidal eccentricity damping counteracts the excitation from planet-planet gravitational interactions, leading to a small but non-zero equilibrium value (Mardling 2007). Using the analytical estimation of Matsumura et al. (2008), we estimated the timescale to reach tidal circularisation. Assuming an initial small orbital eccentricity of 0.01, a constant tidal quality factor Q of 10⁴ (Louden et al. 2023) and a Love number identical to the Earth (k₂=0.3), we obtain τ_circ,_c ~ 5−6 × 10⁹ yr for TOI-1453 c. Regarding TOI-1453 b, the same exercise leads to τ_circ,b ~ 9−10 × 10⁹ yr. While these timescales are approximative, they suggest that a small orbital eccentricity could be maintained.

Furthermore, the nearby stellar companion is likely to have played a role in shaping the architecture of the system we observe today. To some extent, it could have had a long-term impact on the planetary dynamics, such as inducing a spinorbit misalignment or low-frequency variations of the planetary orbital elements. Recent population studies revealed an alignment of binary orbits with the orbits of small S-type planets (e.g. Behmard et al. 2022; Christian et al. 2024) for binary separations comprised between 100 and 700 AU. This result supports theoretical investigations predicting that the binary companions can align the protoplanetary disc inside which the planets form (Batygin 2012). But also, the gravitational influence of this companion could have truncated the protoplanetary disk around TOI-1453 and subsequently influenced the formation of the planets (e.g. Ronco et al. 2021). This exact influence is yet poorly understood, as systematic studies on S-type planets (which orbit one component of a multi-star system) are just beginning. Sullivan et al. (2023) analysed the population of S-type planets from the Kepler mission and found no radius valley. However, the number of precisely characterised S-type planets is yet very small. New detections, such as that of TOI-1453, are crucial to increase the sample size and firmly answer questions related to the radius valley in S-type planets.

In Fig. 12, we present the current mass-radius diagram of small planets with a precision on the mass and radius better than 25% and 8%, respectively (cf. Otegi et al. 2020; Parc et al. 2024), and obtained via the DACE platform¹². We overplot in this diagram TOI-1453 b (in orange) and c (in green). Model compositional curves from Zeng et al. (2016) are also highlighted in the diagram with thick solid lines. TOI-1453 b has likely a rocky composition as shown in the figure, even though this observation must be taken with caution given the large uncertainties on the planet mass. The dashed lines represent the mass-radius relations for Earth-like objects (33% Iron core plus 67% silicates mantle) with H/He atmospheres of different mass fractions. These relations, taken from Lopez & Fortney (2014), take into account the thermal evolution of the planetary atmosphere, which is an important factor to interpret the mass and radius of sub-Neptunes (Rogers et al. 2023). In Fig. 12, they are given for an irradiation of 10 F_⊕, a Solar metallicity and an age of 10 Gyr. Compared to these curves, the mass and radius of TOI-1453 c are compatible with an extended H₂ atmosphere. Naturally, these dashed lines are dependent on the stellar age, metallicity, incoming irradiation, and they do not take into account the photo-evaporation effect. Nevertheless, they illustrate the degeneracy in the mass-radius diagram between rocky cores with extended H/He atmospheres and water worlds. The distinction between these two compositional possibilities cannot be resolved from the mass and radius alone. Our modelling of the interior structure of TOI-1453 c stressed this degeneracy (cf. Sect. 6). However, it allowed us to constrain the water mass fraction in the envelope both for a water-rich and a water-poor formation scenario. If planet c formed beyond the ice-line, we would expect a water content of ~90%. On the contrary, in a water-poor scenario we would expect an envelope mass fraction of ~0.5%, among which water would make up ~0.5% only.

An investigation of the atmosphere of the planet would provide important information in favour of one scenario or the other. To quantify the detectability of the atmospheres of these two planets, we computed the transmission spectroscopy metric (TSM, Kempton et al. 2018), which is correlated with the expected S/N of the transmission spectroscopy signals. While the TSM is negligible for TOI-1453 b given its small transit depth (TSM=5.1), we find a TSM of 78.8 for TOI-1453 c. This is below the threshold of 90 recommened by Kempton et al. (2018) for planets of this radius range. However, we stress that this threshold serves as an indication only. Observations with JWST could provide high quality spectra of TOI-1453 and shed some light on the atmospheric composition of planet c. In Table 2, we report additional planet properties useful for atmospheric follow-up.

All these different characteristics (planet bulk densities, stellar host parameters, dynamics, stellar companion, atmospheric prospects, etc.) constitute as many constraints that could be added to the various formation and evolution scenarios. They make TOI-1453 a particularly valuable system which, with its long photometric baseline, likely represents a PLATO archetype. PLAnetary Transits and Oscillations of stars (PLATO) is the forthcoming ESA space-based mission that will photometrically survey the same portion of sky for at least two consecutive years, with exquisite precision (Rauer et al. 2014). Thanks to this continuous monitoring, it will accumulate many transits of the short-period small-planet systems, similarly to TOI-1453. The extreme precision of PLATO will constrain further those systems, notably with a good estimate of the stellar age and rotation cycles, and very precise transit timings for TTV observations. The PLATO telescope will reveal more small-planet systems similar to TOI-1453 and amenable to RV follow-ups. This will populate the still sparse mass-radius region of small-mass subNeptunes, offering insights into the origin of these intriguing worlds.

Fig. 11 Population of precisely characterised planets (uncertainty on density smaller than 25%) orbiting FGK stars. The left panel presents the period-radius diagram, while the right panel shows the incident flux versus radius diagram. The shaded areas indicate the planet occurrence rate. The orange (green) star symbol represents TOI-1453 b (c) in these 2D spaces.

Fig. 12 Mass-radius diagram of all exoplanets smaller than 5 R⊕ and with a precision on the mass and radius better than 25% and 8%, respectively.

Data availability

Table A.2 is only available in electronic form at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/696/A86

Acknowledgements

The authors thank the referee for a constructive report. This work is based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundación Galileo Galilei (FGG) of the Istituto Nazionale di Astrofisica (INAF) at the Observatorio del Roque de los Muchachos (La Palma, Canary Islands, Spain). The HARPS-N project was funded by the Prodex Program of the Swiss Space Office (SSO), the Harvard University Origin of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory (SAO), the Italian National Astrophysical Institute (INAF), University of St. Andrews, Queen’s University Belfast, and University of Edinburgh. This paper made use of data collected by the TESS mission and are publicly available from the Mikulski Archive for Space Telescopes (MAST) operated by the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by NASA’s Science Mission Directorate. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch. M.S. thanks the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number C4000140754. A.M. acknowledges funding from a UKRI Future Leader Fellowship, grant number MR/X033244/1 and a UK Science and Technology Facilities Council (STFC) small grant ST/Y002334/1. MC acknowledges the SNSF support under grant P500PT_211024. KAC and CNW acknowledge support from the TESS mission via subaward s3449 from MIT. This research has made use of the Exoplanet Follow-up Observation Program (ExoFOP; DOI: 10.26134/ExoFOP5) website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. X.D acknowledges the support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement SCORE No 851555) and from the Swiss National Science Foundation under the grant SPECTRE (No. 200021_215200). This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. A.L. acknowledges support of the Swiss National Science Foundation under grant number TMSGI2_211697. This work makes use of observations from the LCOGT network. Part of the LCOGT telescope time was granted by NOIRLab through the Mid-Scale Innovations Program (MSIP). MSIP is funded by NSF. This article is based on observations made with the MuSCAT2 instrument, developed by ABC, at Telescopio Carlos Sánchez operated on the island of Tenerife by the IAC in the Spanish Observato- rio del Teide. This work is partly financed by the Spanish Ministry of Economics and Competitiveness through grants PGC2018-098153-B-C31. As well, this work is supported by JSPS KAKENHI Grant Number JP24H00017, JP24K00689 and JSPS Bilateral Program Number JPJSBP120249910. This work is partly supported by JSPS KAKENHI Grant Numbers JP17H04574, JP18H01265, and JP18H05439, Grant-in-Aid for JSPS Fellows Grant Number JP20J21872, JST PRESTO Grant Number JPMJPR1775, and a University Research Support Grant from the National Astronomical Observatory of Japan (NAOJ). This paper is based on observations made with the MuSCAT2 instrument, developed by ABC, at Telescopio Carlos Sánchez operated on the island of Tenerife by the IAC in the Spanish Observatorio del Teide. SNQ acknowledges support from the TESS mission via subaward s3449 from MIT. This work has been partially supported by the National Aeronautics and Space Administration under grant No. NNX17AB59G, issued through the Exoplanets Research Program. 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. M.P. acknowledge support from the European Union - NextGen- erationEU (PRIN MUR 2022 20229R43BH) and the “Programma di Ricerca Fondamentale INAF 2023”.

Appendix A TESS and HARPS-N data

In this appendix, we share a detailed summary of the TESS (Table A.1) and HARPS-N (Table A.2) data used in this paper. The light curve of each TESS sector was obtained via a custom extraction from the TPF. Regarding the HARPS-N spectroscopic data, we report the YARARA dataset which we used in this work. The full table is available on the CDS.

Appendix B Interstellar absorption

During our analysis of the HARPS-N spectra, we found narrow absorption features in the wide NaD lines around 5890 and 5896 Å. In Fig. B.1, we present the section of a spectrum where this observation was made. This feature is likely the absorption signature of the interstellar medium itself on the Sodium line since the same feature is observed on both Sodium Doublet lines and is for instance not visible in the Ni1 line between the Sodium lines.

Notably, this feature appears when the observed star is sufficiently far away (here ∼80 pc) to exhibit a large extinction (see Fig.3 in Strassmeier & Weber 2020). We note that the intensity of such absorption is sky-direction-type dependent since the gas density is not homogeneously distributed (see Fig.2 in Murga et al. 2015). By computing the RV shift of the feature compared to the wavelength in the stellar rest-frame, we estimated that the velocity of the interstellar medium along the line-of-sight was about −15.5 km/s which is compatible with the expected kinematics in the Galactic plane. The line profile itself is very different from the stellar line profiles, with a thinner profile in agreement with the expected absorption by a cool material containing small micro-turbulence.

Fig. B.1 HARPS-N spectrum zoomed onto the NaI doublet. The main absorption lines contain an additional narrow feature blueshifted with respect to their centre. These secondary lines are not visible in other spectral lines. Therefore, we exclude contamination from the nearby stellar companion as the cause of this additional absorption. These features have a constant relative position with respect to the main absorption lines of ∼15.5 km/s. We attribute this observation to interstellar absorption.

Appendix C Final system parameters

We report in Table C.1 the results of our final joint fit of the HARPS-N and TESS data.

Appendix D Internal structure of TOI-1453 c

Table D.1 details the posterior resulting from the internal structure modelling of TOI-1453 c. This table contains the mass fractions of the core, mantle and envelope together with the fraction of individual elements in the core and mantle. We report the results for both water-rich and water-poor formation scenarios, and for all three sets of Si/Mg/Fe priors.

References

All Tables

All Figures

Fig. 1 Example of TESS target pixel file centred on TOI-1453 (TIC 198390247). The target is highlighted with a white cross at the centre of the grid. The zone inside of which we perform aperture photometry is contoured in bold gray. All the sources in the neighbourhood of TOI- 1453 identified with the Gaia DR3 catalogue are shown in red, and the arrows illustrate their proper motions on the sky. This plot was designed with tpfplotter (Aller et al. 2020).

In the text

	Fig. 2 AO image of TOI-1453 and corresponding sensitivity curve obtained with NIRC2 installed on the Keck2 telescope (cf. https://exofop.ipac.caltech.edu/tess/view_tag.php?tag=19377, PI: Gonzales). The filter used for this image is Brγ. The faint companion at about 2″ away from TOI-1453 is revealed by this image.
In the text

	Fig. 3 HARPS-N YARARA RVs, $\log \left(R_{HK}^{\prime }\right)$, CCF FWHM, bisector span, and Hαtime series. Left: the various time series. Middle: GLS periodograms, together with the 10, 1, and 0.1% FAP levels. Right: correlation plots with the RVs. In the plots, we also indicate the Spearman correlation coefficient.
In the text

	Fig. 4 GLS periodograms of the TESS SAP photometry separated into three datasets: sectors 14 to 26 (light blue), sectors 49 to 60 (dark blue), and sectors 73 to 76 (gray). In the legend, the periods at which the power is maximum are indicated inside parentheses for each dataset.
In the text

Fig. 5 TESS light curve of TOI-1453 (29 sectors comprised between sector 14 and sector 76) together with the transit model in black. (a) Full TESS light curve. The zone highlighted in purple illustrates the timespan of the HARPS-N observations. The three deeper transit features in the model correspond to simultaneous transits of planets b and c. (b) TESS light curve phase-folded on the transits of TOI-1453 b (TOI-1453.02). The orange dots are binned TESS data for better visibility. (c) Same for TOI-1453 c (TOI-1453.01), where the green dots are binned TESS data.

In the text

	Fig. 6 HARPS-N RV phase-folded on the results of the joint fit for both TOI-1453 b (left) and TOI-1453 c (right). The black line represents the best-fit model RV curve of each planet.
In the text

Fig. 7 Chaos map of the TOI-1453 system, with a resolution of 201×201. This map estimates the chaos in a 2D section of the parameter space, along Pc/Pb (x-axis) and ec (y-axis). The colour code depicts the NAFF chaos estimation, where redder areas are more chaotic. The white regions were unstable before the end of the short-term numerical simulations (escape or close encounter). The actual period ratio of the planet pair constrained from our photometric fit is represented by the vertical line.

In the text

	Fig. 8 Grid of simulated TTVs of TOI-1453 c as a function of the mass of the inner planet Mb and its orbital eccentricity eb. The vertical black line depicts the best-fit mass estimate from our joint fit, while the surrounding white band represents the 68.27% confidence interval on this estimation.
In the text

Fig. 9 TTV searches in the TESS LCs. The top plot shows the O- C diagram of the individual transits of planet c, and the middle panel presents the GLS periodogram of those observations. The solid line is the periodogram derived from all the 2 min cadence data. The dashed line was obtained after replacing the 2 min by 20s cadence data when available, and proceeding to new transit fits (cf. the orange points in the top plot). We do not find significant periodicity. The bottom plot presents the variation of the orbital period of planet b as measured from pairs of consecutive sectors (with 2 min cadence).

In the text

Fig. 10 Posterior distributions of the interior structure of TOI-1453 c, showing mass fractions of the inner core, mantle, and envelope, as well as the water fraction in the envelope. The top row corresponds to a water prior compatible with a formation scenario beyond the ice line, while the bottom row assumes a water-poor composition, in agreement with formation inside the ice line. Colours represent three priors for the planetary Si/Mg/Fe ratios: stellar (purple, A1 and B1), iron-enriched (pink, A2 and B2), and uniformly sampled from a simplex (blue, A3 and B3). Dashed lines mark the median of each distribution, with priors indicated by dotted lines.

In the text

	Fig. 11 Population of precisely characterised planets (uncertainty on density smaller than 25%) orbiting FGK stars. The left panel presents the period-radius diagram, while the right panel shows the incident flux versus radius diagram. The shaded areas indicate the planet occurrence rate. The orange (green) star symbol represents TOI-1453 b (c) in these 2D spaces.
In the text

	Fig. 12 Mass-radius diagram of all exoplanets smaller than 5 R⊕ and with a precision on the mass and radius better than 25% and 8%, respectively.
In the text

Fig. B.1 HARPS-N spectrum zoomed onto the NaI doublet. The main absorption lines contain an additional narrow feature blueshifted with respect to their centre. These secondary lines are not visible in other spectral lines. Therefore, we exclude contamination from the nearby stellar companion as the cause of this additional absorption. These features have a constant relative position with respect to the main absorption lines of ∼15.5 km/s. We attribute this observation to interstellar absorption.

In the text