# Improved characterization of the TOI-2141 system: A dense sub-Neptune with non-transiting inner and outer companions

R. Luque, K. W. F. Lam, J. Cabrera, A. Bonfanti, Y. N. E. Eschen, G.

Olofsson, W. Benz, N. Billot, A. Brandeker, A. C. M. Correia, et al.

► To cite this version:

R. Luque, K. W. F. Lam, J. Cabrera, A. Bonfanti, Y. N. E. Eschen, et al. Improved characterization of the TOI-2141 system: A dense sub-Neptune with non-transiting inner and outer companions. Astronomy & Astrophysics - A&A, 2025, 704, pp.A174. 10.1051/0004-6361/202556515. hal-05477490

HAL Id: hal-05477490

https://hal.science/hal-05477490v1

Submitted on 26 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Improved characterization of the TOI-2141 system: A dense sub-Neptune with non-transiting inner and outer companions

R. Luque<fnref n="1" />,*, K. W. F. Lam², J. Cabrera², A. Bonfanti³, Y. Eschen⁴, G. Olofsson⁵, W. Benz⁶,⁷, N. Billot⁸, A. Brandeker⁵, A. C. M. Correia⁹, L. Fossati³, D. Gandolfi¹⁰, H. P. Osborn⁷,¹¹, C. Pezzotti¹², S. G. Sousa¹³, T. G. Wilson⁴, S. Wolf¹⁴, Y. Alibert⁷,⁶, R. Alonso¹⁵,¹⁶, J. Asquier¹⁷, T. Bárczy¹⁸, D. Barrado¹⁹, S. C. C. Barros¹³,²⁰, W. Baumjohann³, F. Biondi²¹,²², L. Borsato²², C. Broeg⁶,⁷, A. Collier Cameron²³, Sz. Csizmadia², P. E. Cubillos³,²⁴, M. B. Davies²⁵, M. Deleuil²⁶, A. Deline⁸, O. D. S. Demangeon¹³,²⁰, B.-O. Demory⁷,²⁷,⁶, A. Derekas²⁸, B. Edwards²⁹, J. A. Egger⁶, D. Ehrenreich⁸,³⁰, A. Erikson², A. Fortier⁶,⁷, M. Fridlund³¹,³², K. Gazeas³³, M. Gillon³⁴, M. Güdel³⁵, M. N. Günther¹⁷, A. Heitzmann⁸, Ch. Helling³,³⁶, K. G. Isaak¹⁷, T. Keller⁶,⁷, L. L. Kiss³⁷,³⁸, J. Korth⁸, J. Laskar³⁹, A. Lecavelier des Etangs⁴⁰, A. Leleu⁸,⁶, M. Lendl⁸, D. Magrin²², G. Mantovan⁴¹,²², P. F. L. Maxted⁴², B. Merín⁴³, C. Mordasini⁶,⁷, V. Nascimbeni²², R. Ottensamer³⁵, I. Pagano⁴⁴, E. Pallé¹⁵,¹⁶, G. Peter², D. Piazza⁴⁵, G. Piotto²²,⁴⁶, D. Pollacco⁴, D. Queloz¹¹,⁴⁷, R. Ragazzoni²²,⁴⁶, N. Rando¹⁷, H. Rauer⁴⁸,⁴⁹, I. Ribas⁵⁰,⁵¹, N. C. Santos¹³,²⁰, G. Scandariato⁴⁴, D. Ségransan⁸, A. E. Simon⁶,⁷, A. M. S. Smith², M. Stalport¹²,³⁴, S. Sulis²⁶, Gy. M. Szabó²⁸,⁵², S. Udry⁸, B. Ulmer², S. Ulmer-Moll⁵³,¹², V. Van Grootel¹², J. Venturini⁸, E. Villaver¹⁵,¹⁶, N. A. Walton⁵⁴, and D. Wolter²

*(Affiliations can be found after the references)*

Received 21 July 2025 / Accepted 29 August 2025

## ABSTRACT

**Aims.** We aim to refine the fundamental parameters of the TOI-2141 planetary system, which includes a transiting sub-Neptune orbiting a Sun-like star on a relatively long orbit of 18.26 days, by combining new photometric and spectroscopic observations.

**Methods.** We analyzed new space-based photometry from TESS and CHEOPS, as well as 61 radial velocity measurements from HARPS-N. We performed individual and joint photometric and RV analyses using several modeling tools within a Bayesian model comparison framework.

**Results.** We refined the radius and mass of the transiting planet TOI-2141 b to 3.15 ± 0.04 $R_⊕$ and 20.1 ± 1.6 $M_⊕$ to levels that are five and two times more precise, respectively, than the previously reported values. Our radial velocity analysis reveals two additional non-transiting companions with orbital periods of 5.46 and 60.45 days. Despite the innermost planet's high geometric transit probability, we found no evidence of transits in the photometric data.

**Conclusions.** The bulk properties of TOI-2141 b suggest a significant volatile envelope atop an Earth-like core. The modeling indicates a hydrogen-rich atmosphere that could have experienced mild photoevaporation over the system's history. Planets b and c are believed to exhibit a modest mutual inclination of at least ~2.4°.

**Key words.** techniques: photometric – techniques: radial velocities – planets and satellites: detection – planets and satellites: individual: TOI-2141 b – stars: solar-type

## 1. Introduction

This work revisits the TOI-2141 planetary system, discovered and characterized by Martioli et al. (2023) using space-based photometry from TESS (Ricker et al. 2014) and radial velocities (RVs) from the SOPHIE high-resolution spectrograph mounted on the 1.93-m telescope at the Observatoire de Haute-Provence (Perruchot et al. 2008; Bouchy et al. 2013). The star TOI-2141 is a relatively bright ($V = 9.46$ mag), solar analog ($T_{eff} \sim 5650$ K) hosting a transiting sub-Neptune ($R_p \sim 3 R_⊕$, $M_p \sim 24 M_⊕$) with an orbital period of 18.26 days. Previous works have argued that TOI-2141 b is likely to have a dense rocky core and, possibly, a

thick, water-rich envelope and several internal structure models could be employed to help explain its measured bulk density (as pointed out by, e.g., Rogers & Seager 2010; López & Fortney 2014; Zeng et al. 2016; Nixon & Madhusudhan 2021; Luo et al. 2024). In this work, we have added unpublished TESS observations from June 2024, four transits of TOI-2141 b with CHEOPS<fnref n="1" /> (Benz et al. 2021), and 61 RVs from HARPS-N (Cosentino et al. 2012) to refine the properties of the system, including its architecture.

<footnote n="1">¹ Observed as part of the Guaranteed Time Observation (GTO) programs CH_PR100024 and CH_PR140073. The raw and detrended photometric time-series data are available in electronic form at the CDS.</footnote>

<footnote>* Corresponding author: rluque@iaa.es</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** PDC-corrected photometry of TOI-2141 from TESS. Gray points are the original 2-min cadence TESS data, while the 15-min binned photometric data are shown in blue. Transits of the planet TOI-2141 b are marked in orange. The black line shows the best-fit model to the data described in Sect. 3.4.</figcaption></figure>

## 2. Observations

### 2.1. TESS photometry

The TESS mission observed TOI-2141 with a cadence of 2-min in Sectors 25, 26, 52, and 79. Data from Sectors 25, 26, and 52 were described and analyzed in Martioli et al. (2023), where the transiting planet TOI-2141 b was identified and confirmed. In this work, we add to our analyses the photometry from Sector 79 gathered in June 2024 and combine it with the rest of the photometry. We downloaded the TESS data products from the Mikulski Archive for Space Telescopes<fnref n="2" /> and use the Presearch

Data Conditioning (PDC) flux time series (Smith et al. 2012; Stumpe et al. 2012, 2014) computed by the TESS Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016). Figure 1 shows the complete TESS dataset used in this work.

### 2.2. CHEOPS photometry

The Characterising Exoplanets Satellite (CHEOPS) is a European Space Agency small-class mission dedicated to studying bright exoplanet host stars using high-precision photometric observations (Broeg et al. 2013; Benz et al. 2021; Fortier et al. 2024). CHEOPS observed TOI-2141 between 05 May 2023 and 16 July 2024. The observations were taken under

<footnote n="2">² https://archive.stsci.edu/</footnote>
<table><caption>**Table 1.** CHEOPS observing log.</caption><thead><tr><th>Visit</th><th>Start (UTC)</th><th>Length [h]</th><th>Archive filekey</th><th>No. Frames</th><th>RMS (ppm)</th><th>Efficiency (%)</th></tr></thead><tbody><tr><td>1</td><td>2023-05-05T03:38:19</td><td>13.29</td><td>CH_PR100024_TG016301_V0300</td><td>506</td><td>426.1</td><td>63.4</td></tr><tr><td>2</td><td>2023-06-10T15:58:19</td><td>13.37</td><td>CH_PR100024_TG016302_V0300</td><td>639</td><td>300.0</td><td>79.5</td></tr><tr><td>3</td><td>2023-06-28T22:15:19</td><td>13.34</td><td>CH_PR100024_TG016303_V0300</td><td>555</td><td>307.7</td><td>69.3</td></tr><tr><td>4</td><td>2024-06-28T05:01:58</td><td>11.04</td><td>CH_PR140073_TG002301_V0300</td><td>517</td><td>315.0</td><td>77.9</td></tr><tr><td>5</td><td>2024-07-16T10:37:59</td><td>11.84</td><td>CH_PR140073_TG002302_V0300</td><td>411</td><td>283.7</td><td>57.8</td></tr></tbody></table>

the CHEOPS Guaranteed Time Observation (GTO) programs PR100024 (MR. Improve) and PR140073 (Gas Content of Low Mass Planets). Five visits were observed with an exposure cadence of 60 s. The data were extracted with the CHEOPS DRP v14 pipeline (Hoyer et al. 2020) which performed instrument calibrations, as well as cosmic ray, background, and smearing corrections. Aperture sizes with radii between 15 and 40 pixels were extracted. We opted to use light curves extracted using aperture radii of either 22 or 23 pixels for further analyses, as these have the lowest RMS values.

Figure A.1 shows the original CHEOPS data and the detrended light curves used in our subsequent analyses. Table 1 presents an observing log summarizing the properties of each visit.

## 2.3. HARPS-N spectroscopy

We collected 61 RVs of TOI-2141 with the HARPS-N spectrograph mounted at the 3.6-m Telescopio Nazionale Galileo in La Palma between October 2022 and September 2024. The data were collected as part of the international time program titled Completing CHEOPS characterization of long-period low-mass planets" (PI: Luque; program IDs ITP22_2 and ITP23_10) granted by the International Scientific Committee of the Observatorios de Canarias of the IAC. Radial velocities and additional spectral indicators were derived using YABI (Borsa et al. 2015) and serval Zechmeister et al. (2018). YABI is an online implementation of the classical HARPS-N DRS pipeline at the INAF Trieste Observatory<fnref n="3" /> and computes cross-correlation functions (CCFs) for a given template mask following Cosentino et al. (2014). On the other hand, serval follows a template-matching technique (Anglada-Escudé & Butler 2012) to compute radial velocities from a high signal-to-noise (S/N) template built from the stellar spectra themselves. The YABI-derived radial velocities have a mean internal precision of 0.98 m s⁻¹ and show a standard deviation of 4.7 m s⁻¹. The serval radial velocities have a mean internal precision of 1.3 m s⁻¹ and show the same standard deviation. Due to their superior precision, we use the YABI RVs in our subsequent analyses.

Figure A.2 shows the HARPS-N RVs together with the SOPHIE RVs presented in Martioli et al. (2023) used in this work. A search for periodic signals in the combined dataset (subtracting for an instrumental offset) using a generalized Lomb-Scargle periodogram (Zechmeister & Kürster 2009) shows a prominent peak (false alarm probability below 10⁻⁶) at the period of the transiting planet TOI-2141 b ($P_b$ = 18.26 days, $f_b$ = 0.0547 d⁻¹). Surprisingly, a periodogram based on the residuals of a RV-only model as in Martioli et al. (2023) that includes a circular orbit at the period of the transiting planet shows two

peaks with false alarm probabilities (FAPs) below 0.1% at 5.46 and 60 days. We analyzed these two signals in detail in Sect. 3.3.

A GLS periodogram of all the available activity indicators is shown in Fig. A.3 of the Appendix. The GLS periodograms do not find any significant peaks (FAP below 10%) in any of the indicators. For $H\alpha$, a slight peak approaching FAP~10% can be found at 40 days, approximately double the inferred rotational period of TOI-2141 of $21 \pm 5$ d (Martioli et al. 2023) based on empirical activity-rotation relations (Mamajek & Hillenbrand 2008). This peak is also seen in SOPHIE's $H\alpha$ index (Martioli et al. 2023, see their Fig. C.6), but again below the significance threshold. Compatible with the reported rotational period of approximately 20 days, there is only a non-significant peak (FAP>10%) in the full width at half maximum (FWHM) of the CCF of the HARPS-N data at 17 d also seen in the bisector span of the SOPHIE CCFs (Martioli et al. 2023). The lack of long-term continuous photometric monitoring of the star prevents us from refining the stellar rotation period using an independent technique.

# 3. Analysis and results

## 3.1. Stellar characterization

The stellar spectroscopic parameters ($T_{eff}$, log g, microturbulence, [Fe/H]) were derived using the ARES+MOOG methodology (Sousa et al. 2021; Sousa 2014; Santos et al. 2013). The equivalent widths (EW) were automatically measured using the ARES code<fnref n="4" /> (Sousa et al. 2007, 2015). In this spectral analysis, we used a combined HARPS-N spectrum from individual exposures of TOI-2141, which was then used to measure the EWs for a list of lines presented in Sousa et al. (2008). We analyzed only the HARPS-N combined spectrum as it is of higher resolution (115 000 vs. 75 000) and S/N (630 vs. 507 on average) than the SOPHIE one. Besides, the SOPHIE spectrum was analyzed in detail in Martioli et al. (2023), for which we find below very consistent stellar parameters to ours. The best set of spectroscopic parameters was found by using a minimization process to find the ionization and excitation equilibrium. This process makes use of a grid of Kurucz model atmospheres (Kurucz 1993) and the latest version of the radiative transfer code MOOG (Sneden 1973). We also derived a more accurate trigonometric surface gravity using recent Gaia EDR3 data following the same procedure described in Sousa et al. (2021), which provided a consistent value compared to the spectroscopic surface gravity.

The stellar radius of TOI-2141 was computed using a MCMC modified infrared flux method (IRFM; Blackwell & Shallis 1977; Schanche et al. 2020). By constructing spectral energy distributions (SED) from two stellar atmospheric models catalogs (Kurucz 1993; Castelli & Kurucz 2003) that were constrained by

<footnote n="4">⁴ The last version, ARES v2, can be downloaded at https://github.com/sousasag/ARES</footnote>

<footnote n="3">³ https://www.ia2.inaf.it/</footnote>
our spectroscopically derived stellar parameters, we determined broadband synthetic photometry. Based on a comparison with observed fluxes in the following bandpasses: 2MASS J, H, and K, WISE W1 and W2, and Gaia G, G<sub>BP</sub>, and G<sub>RP</sub> (Skrutskie et al. 2006; Wright et al. 2010; Gaia Collaboration 2021), we computed the stellar bolometric flux that we converted into the effective temperature and angular diameter of TOI-2141. Lastly, we derived the stellar radius by combining the angular diameter with the offset-corrected Gaia parallax (Lindegren et al. 2021). To account for stellar atmosphere model uncertainties, we conducted a Bayesian model averaging of the stellar radius posteriors produced using each catalog. The weighted average value (0.950 ± 0.007 R<sub>☉</sub>) is reported in Table 2.

We used T<sub>eff</sub>, [Fe/H], and R<sub>⋆</sub> to derive the stellar mass M<sub>⋆</sub> and age t<sub>⋆</sub> accounting for two different stellar evolutionary models. In detail, we derived a first set of mass and age estimates thanks to the isochrone placement routine (Bonfanti et al. 2015, 2016) that interpolates the input parameters within pre-computed grids of PARSEC<sup>5</sup> v1.2S (Marigo et al. 2017) tracks and isochrones. Then, a second mass and age set was computed by the Code Liègeois d'Évolution Stellaire (CLES; Scuflaire et al. 2008) code that builds the best-fit evolutionary track based on the input parameters following a Levenberg-Marquardt optimisation scheme (Salmon et al. 2021). After carefully checking the consistency of the two respective pairs of outcomes via the χ<sup>2</sup>-based criterion outlined in Bonfanti et al. (2021a), we finally merged (i.e. we summed) the respective output distributions and we obtained M<sub>⋆</sub> = 0.896<sup>+0.059</sup><sub>-0.051</sub> M<sub>☉</sub> and t<sub>⋆</sub> = 9 ± 4 Gyr (see Bonfanti et al. 2021a, for further details).

Furthermore, we computed the right-handed, heliocentric Galactic space velocities (Johnson & Soderblom 1987), U, V, and W, using the Gaia EDR3 coordinates, proper motions, offset-corrected parallax (Lindegren et al. 2021), and RV (Gaia Collaboration 2021). These are reported in Table 2 and strongly indicate that the star is a member of the thin disk population.

## 3.2. CHEOPS analysis

We used the publicly available pycheops package (Maxted et al. 2022) to model systematic trends and fit transit light curves simultaneously. A 5-sigma clipping was first applied to the raw light curves to remove any outliers. Each visit was then detrended separately, where the set of detrending basis vectors used in respective detrending models were chosen based on the Bayes factor. A "glint" effect can be introduced into the light curve due to internal reflection of bright objects (e.g., moonlight) within 24° from the observed target (Maxted et al. 2022). We applied the glint function to remove glint effects using a spline fit with 8 to 11 knots. A Gaussian white noise component was also included in the model. Next, pycheops, implementing the LMFIT package (Newville et al. 2014), was run to perform a nonlinear least-squares minimisation and emcee (Foreman-Mackey et al. 2013) was used for the Bayesian analysis and posterior sampling. The final detrended light curves are shown in Figure A.1 and were used in subsequent joint fits.

## 3.3. Radial velocity analysis

We used juliet (Espinoza et al. 2019), a python library built on many publicly available tools for the modeling of transits (batman, Kreidberg 2015), RVs (radvel, Fulton et al. 2018),

<table><caption>**Table 2.** Stellar parameters of TOI-2141.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td colspan="3"><em>Name and identifiers</em></td></tr><tr><td>Name</td><td>BD+18 3330</td><td>BD</td></tr><tr><td>TIC</td><td>287256467</td><td>TIC</td></tr><tr><td>TOI</td><td>TOI-2141</td><td>Martioli et al. (2023)</td></tr><tr><td colspan="3"><em>Coordinates and spectral type</em></td></tr><tr><td>α</td><td>17:15:02.906</td><td>Gaia EDR3</td></tr><tr><td>δ</td><td>+18:20:26.76</td><td>Gaia EDR3</td></tr><tr><td>Epoch (ICRS)</td><td>J2000</td><td>Gaia EDR3</td></tr><tr><td colspan="3"><em>Magnitudes</em></td></tr><tr><td>V [mag]</td><td>9.46 ± 0.003</td><td>UCAC4</td></tr><tr><td>G [mag]</td><td>9.34408 ± 0.00021</td><td>Gaia EDR3</td></tr><tr><td>J [mag]</td><td>8.27 ± 0.02</td><td>2MASS</td></tr><tr><td>K<sub>s</sub> [mag]</td><td>7.87 ± 0.02</td><td>2MASS</td></tr><tr><td colspan="3"><em>Parallax and kinematics</em></td></tr><tr><td>π [mas]</td><td>12.957 ± 0.015</td><td>Gaia EDR3</td></tr><tr><td>d [pc]</td><td>77.7 ± 0.2</td><td>Gaia EDR3</td></tr><tr><td>μ<sub>α</sub> cos δ [mas yr<sup>-1</sup>]</td><td>52.229 ± 0.013</td><td>Gaia EDR3</td></tr><tr><td>μ<sub>δ</sub> [mas yr<sup>-1</sup>]</td><td>-98.286 ± 0.015</td><td>Gaia EDR3</td></tr><tr><td>U [km s<sup>-1</sup>]</td><td>15.56 ± 0.30</td><td>This work</td></tr><tr><td>V [km s<sup>-1</sup>]</td><td>-21.60 ± 0.25</td><td>This work</td></tr><tr><td>W [km s<sup>-1</sup>]</td><td>-36.53 ± 0.22</td><td>This work</td></tr><tr><td colspan="3"><em>Photospheric parameters</em></td></tr><tr><td>T<sub>eff</sub> [K]</td><td>5635 ± 61</td><td>This work</td></tr><tr><td>log(g [cm s<sup>-2</sup>])</td><td>4.40 ± 0.03</td><td>This work</td></tr><tr><td>[Fe/H]</td><td>-0.14 ± 0.04</td><td>This work</td></tr><tr><td>v<sub>micro</sub> [km s<sup>-1</sup>]</td><td>0.92 ± 0.01</td><td>This work</td></tr><tr><td colspan="3"><em>Physical parameters</em></td></tr><tr><td>M [M<sub>☉</sub>]</td><td>0.896<sup>+0.059</sup><sub>-0.051</sub></td><td>This work</td></tr><tr><td>R [R<sub>☉</sub>]</td><td>0.950 ± 0.007</td><td>This work</td></tr><tr><td>Age [Gyr]</td><td>9 ± 4</td><td>This work</td></tr></tbody></table>

**Notes.** BD: Schönfeld (1886); TIC: Stassun et al. (2018); Gaia EDR3: Gaia Collaboration (2021); UCAC4: Zacharias et al. (2013); 2MASS: Skrutskie et al. (2006).

and Gaussian Processes (george, Ambikasaran et al. 2015; celerite, Foreman-Mackey et al. 2017) to model the RV data first. We use the nested sampler dynasty (Speagle 2020) to explore the parameter space of a given prior volume and to compute efficiently the Bayesian model log evidence (ln Z). Thanks to this, we can compare models with different numbers of parameters accounting for the model complexity and the number of degrees of freedom with a sound statistical methodology. In our analysis, we considered that a model (M2) is greatly favored over another (M1) if Δln Z = ln Z<sub>M2</sub> - ln Z<sub>M1</sub> > 5 (Trotta 2008). If Δln Z ∈ [1, 2.5), we find that M2 is weakly favored over M1, so the simpler model with fewer degrees of freedom would be chosen. For intermediate cases, we consider that M2 is moderately favored over M1.

We carried out a Bayesian model comparison analysis to assess the statistical significance of the RV signals present in the SOPHIE and HARPS-N data (given in Fig. 2). We did not fit the photometry in this analysis. We followed an iterative approach where we add a Keplerian signal to the model at a time using wide, uninformative priors on all model parameters (see Hara & Ford 2023, for a review of statistical methods for analyzing RV datasets). Results are reported in Table 3.

First, we find that, as in Martioli et al. (2023), models with or without fitting for the eccentricity are statistically

<footnote n="5">5 PAOVA and TRIESTE Stellar Evolutionary Code: https://stev.oapd.inaf.it/cgi-bin/cmd</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** GLS periodograms of the SOPHIE and HARPS-N RV data. For each panel, the horizontal lines show the theoretical 10% (short-dashed line), 1% (long-dashed line), and 0.1% (dot-dashed line) false alarm probability levels. The vertical dashed lines mark the orbital frequencies of the transiting planet ($f_0 = 0.0547 \text{ d}^{-1}$) and of the signals detected in the RVs at 5.46 and 60 d. Panel a: Combined RV data corrected for an instrumental offset. Panel b: RV residuals following the subtraction of the signal of the transiting planet TOI-2141 b.</figcaption></figure>

<table><caption>**Table 3.** Model comparison of RV-only fits with juliet.</caption><thead><tr><th>Model</th><th>ln Z</th><th>K<sub>b</sub> (m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1p</td><td>-475.0</td><td>5.4 ± 0.5</td></tr><tr><td>1p\_ecc</td><td>-475.0</td><td>5.4 ± 0.5</td></tr><tr><td>1p + GP</td><td>-473.3</td><td>5.3 ± 0.4</td></tr><tr><td>2p (5-d)</td><td>-470.8</td><td>5.6 ± 0.4</td></tr><tr><td>2p (60-d)</td><td>-470.9</td><td>5.2 ± 0.4</td></tr><tr><td>2p + GP</td><td>-463.9</td><td>5.3 ± 0.3</td></tr><tr><td><b>3p</b></td><td><b>-459.3</b></td><td><b>5.3 ± 0.3</b></td></tr><tr><td>3p + GP</td><td>-465.5</td><td>5.4 ± 0.3</td></tr></tbody></table>

**Notes.** GP refers to a Gaussian process applied with an exponential sine-squared kernel from george. The model used for the final joint fit is marked in boldface.

indistinguishable ($\Delta \ln Z < 1$). Thus, we assumed circular orbits for all the signals in this analysis from here onward. Then, we find that adding a second signal to the model with an uninformative wide prior between 1 and 1000 days, results in a much better fit ($\Delta \ln Z \sim 4$). Different instances of that model result in a posterior distribution on the period centered at either 5.46 or 60.45 days, but there is no preference for including one over the other. However, the best fit to the data is with a model that includes all three signals ($\Delta \ln Z = 15.7$ over a 1p model, but >11 over 2p models).

This 3p model is preferred even in the presence of correlated noise, which we model using a Gaussian process (GP) exponential sine-squared kernel from george of the form

$$k_{i,j} = \sigma^2 \exp \left( -\alpha(t_i - t_j)^2 - \Gamma \sin^2 \left[ \frac{\pi|t_i - t_j|}{P_{\text{rot}}} \right] \right). \quad (1)$$

Testing with other GP kernels (exponential decay or quasi-periodic) give similar results ($\Delta \ln Z \sim 0.2$). We set log-uniform priors for $\sigma$ between 0.1 and 300 m s$^{-1}$, $\alpha$ between $10^{-8}$ and $10^{-3}$, $\Gamma$ between 0.01 and 10, and a uniform prior for $P_{\text{rot}}$ between 10 and 70 days. The choice of the prior on $P_{\text{rot}}$ is motivated to try to distinguish if the 60-day signal is of stellar origin, although the typical rotation period of Sun-like stars is between 10 and 30 days (e.g., McQuillan et al. 2014) and would be unexpected.

Furthermore, none of the activity indicators from the HARPS-N spectra show hints of stellar variability at that or any other period (Fig. A.3). The hyperparameters from these GP models show that $P_{\text{rot}} = 60.5 \pm 0.6 \text{ d}$ and $\alpha \sim 10^{-7}$, which indicates that the GP is modeling a signal close to be perfectly periodic with a very well-defined frequency unlikely to be of stellar origin (Stock et al. 2023).

We find that although a 2p model (with the transiting planet and the 5.46-day signal) and a GP explain the data well, assuming a circular orbit for the 60-day signal is not only a simpler model but statistically preferred ($\Delta \ln Z = 4.6$). A stacked Bayesian GLS periodogram analysis following Mortier & Collier Cameron (2017) reveals that all three significant RV signals (at 5.46, 18.26, and 60 days) are coherent and stable throughout the observing period, reinforcing the notion that they are of planetary nature (Fig. A.4). In addition, it is important to highlight that the derived semi-amplitude for the transiting planet ($K_b$) remains within $1\sigma$ for every model tested, ensuring a robust mass determination of the planet regardless of the chosen model.

To conclude, we find that the model that offers the highest statistical evidence is the three-planet model (assuming circular orbits), with no need to account for correlated noise. There is enough statistical evidence to support the claim that the 5.46- and 60-day signals as present (and significant) in the RV data and likely of planetary origin. This finding is in line with other works that follow-up transiting planets with RV instruments and find internal and/or external companions (e.g., Luque et al. 2019; Kemmer et al. 2020).

### 3.4. Joint fit

Using juliet, we jointly fit the TESS, CHEOPS, SOPHIE, and HARPS-N data. To speed up computational times, we used the TESS PDC photometry without additional detrending and the detrended CHEOPS data analyzed in Sect. 3.2. For the model, we assumed the configuration inferred from our Bayesian model comparison RV-only analysis presented in Sect. 3.3 and shown in Table 3. For the photometry part of our fit, we adopted a quadratic limb darkening law for TESS and CHEOPS data parameterized by the coefficients $q_1, q_2$ introduced by Kipping (2013) and fit them as free parameters with unconstrained priors. The transiting planet is modeled using the $r_1, r_2$ parameterization of the planet-to-star radius ratio $p = R_p/R_*$ and the impact parameter of the orbit $b$ introduced by Espinoza (2018). Additionally, rather than fitting the scaled planetary radius ($a/R_*$) for each planet, we use the stellar density ($\rho_*$) as a free parameter and set a normal prior based on the stellar properties described in Table 2. We added a jitter term $\sigma$ to the photometric uncertainties of TESS and CHEOPS data. For the radial velocity part of our fit, we assumed circular orbits for all planets. The priors on the period and phase of each planet were informed by our previous transit- and RV-only analyses, but wide enough to derive reliable uncertainties. We added an instrumental offset and a jitter term for each RV instrument and no GP. The model parameters, their priors, and posterior distributions are reported in Table 4. The best-fit models are shown in Figs. 1 and A.2 for the photometry and RV timeseries, and phase-folded to the periods of interest in Fig. 3.

## 4. Discussion

Table 5 shows the fundamental physical parameters derived for the three planets based on the model posteriors from Table 4
<table><caption>**Table 4.** Priors and posterior distributions (median and 68% credibility intervals) for each fit parameter of the final joint model obtained for the TOI-2141 system using juliet.</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Posterior</th><th>Units</th></tr></thead><tbody><tr><td>ρ<sub>★</sub></td><td>N(1473, 100)</td><td>1521<sup>+96</sup><sub>-95</sub></td><td>kg m<sup>-3</sup></td></tr><tr><td>P<sub>c</sub></td><td>U(4.8, 6.0)</td><td>5.4624<sup>+0.0026</sup><sub>-0.0027</sub></td><td>d</td></tr><tr><td>t<sub>0,c</sub></td><td>U(1976.0, 1982.0)</td><td>1979.19<sup>+0.58</sup><sub>-0.56</sub></td><td>BJD-2457000</td></tr><tr><td>K<sub>c</sub></td><td>U(0, 20)</td><td>2.40<sup>+0.33</sup><sub>-0.33</sub></td><td>m s<sup>-1</sup></td></tr><tr><td>P<sub>b</sub></td><td>N(18.26, 0.1)</td><td>18.261608<sup>+1.9e-05</sup><sub>-2e-05</sub></td><td>d</td></tr><tr><td>t<sub>0,b</sub></td><td>N(1992.5, 1.0)</td><td>1992.502<sup>+0.0014</sup><sub>-0.0013</sub></td><td>d</td></tr><tr><td>r<sub>1,b</sub></td><td>U(0, 1)</td><td>0.767<sup>+0.013</sup><sub>-0.014</sub></td><td>...</td></tr><tr><td>r<sub>2,b</sub></td><td>U(0, 1)</td><td>0.03037<sup>+0.00035</sup><sub>-0.00034</sub></td><td>...</td></tr><tr><td>K<sub>b</sub></td><td>U(0, 20)</td><td>5.25<sup>+0.35</sup><sub>-0.35</sub></td><td>m s<sup>-1</sup></td></tr><tr><td>P<sub>d</sub></td><td>U(50, 70)</td><td>60.45<sup>+0.24</sup><sub>-0.25</sub></td><td>d</td></tr><tr><td>t<sub>0,d</sub></td><td>U(1900.0, 1960.0)</td><td>1940.5<sup>+5.2</sup><sub>-5.2</sub></td><td>d</td></tr><tr><td>K<sub>d</sub></td><td>U(0, 20)</td><td>2.61<sup>+0.34</sup><sub>-0.34</sub></td><td>m s<sup>-1</sup></td></tr><tr><td>γ<sub>HARPS-N</sub></td><td>U(-19900.0, -19800.0)</td><td>-19829.97<sup>+0.26</sup><sub>-0.25</sub></td><td>m s<sup>-1</sup></td></tr><tr><td>σ<sub>HARPS-N</sub></td><td>J(0.1, 20)</td><td>1.70<sup>+0.23</sup><sub>-0.30</sub></td><td>m s<sup>-1</sup></td></tr><tr><td>γ<sub>SOPHIE</sub></td><td>U(-19900.0, -19800.0)</td><td>-19860.58<sup>+0.62</sup><sub>-0.60</sub></td><td>m s<sup>-1</sup></td></tr><tr><td>σ<sub>SOPHIE</sub></td><td>J(0.1, 20)</td><td>5.25<sup>+0.62</sup><sub>-0.60</sub></td><td>m s<sup>-1</sup></td></tr><tr><td>q<sub>1,TESS</sub></td><td>U(0, 1)</td><td>0.30<sup>+0.18</sup><sub>-0.13</sub></td><td>...</td></tr><tr><td>q<sub>2,TESS</sub></td><td>U(0, 1)</td><td>0.37<sup>+0.34</sup><sub>-0.25</sub></td><td>...</td></tr><tr><td>q<sub>1,CHEOPS</sub></td><td>U(0, 1)</td><td>0.32<sup>+0.08</sup><sub>-0.06</sub></td><td>...</td></tr><tr><td>q<sub>2,CHEOPS</sub></td><td>U(0, 1)</td><td>0.76<sup>+0.17</sup><sub>-0.24</sub></td><td>...</td></tr><tr><td>σ<sub>TESS</sub></td><td>J(0.1, 1000)</td><td>1.34<sup>+6.52</sup><sub>-1.11</sub></td><td>ppm</td></tr><tr><td>σ<sub>CHEOPS</sub></td><td>J(0.1, 1000)</td><td>230.85<sup>+6.65</sup><sub>-6.40</sub></td><td>ppm</td></tr></tbody></table>

**Notes.** The prior labels of *N*, *U*, and *J* represent normal, uniform, and Jeffrey's distributions, respectively.

<table><caption>**Table 5.** Derived planetary parameters obtained for the TOI 2141 system using the posterior values from the joint fit in Table 4 and stellar parameters from Table 2.</caption><thead><tr><th>Parameter</th><th>Planet c</th><th>Planet b</th><th>Planet d</th><th>Units</th></tr></thead><tbody><tr><td>p = R<sub>p</sub>/R<sub>★</sub></td><td>...</td><td>0.03037<sup>+0.00035</sup><sub>-0.00034</sub></td><td>...</td><td>...</td></tr><tr><td>b = (a<sub>p</sub>/R<sub>★</sub>) cos i<sub>p</sub></td><td>...</td><td>0.65<sup>+0.019</sup><sub>-0.02</sub></td><td>...</td><td>...</td></tr><tr><td>a<sub>p</sub>/R<sub>★</sub></td><td>13.39<sup>+0.27</sup><sub>-0.28</sub></td><td>29.93<sup>+0.62</sup><sub>-0.64</sub></td><td>66.47<sup>+1.37</sup><sub>-1.43</sub></td><td>...</td></tr><tr><td>i<sub>p</sub></td><td>...</td><td>88.755<sup>+0.063</sup><sub>-0.063</sub></td><td>...</td><td>deg</td></tr><tr><td>M<sub>p</sub></td><td>...</td><td>20.1<sup>+1.6</sup><sub>-1.5</sub></td><td>...</td><td>M<sub>⊕</sub></td></tr><tr><td>M<sub>p</sub> sin i</td><td>6.14<sup>+0.86</sup><sub>-0.86</sub></td><td>20.1<sup>+1.6</sup><sub>-1.5</sub></td><td>14.9<sup>+2.1</sup><sub>-2.1</sub></td><td>M<sub>⊕</sub></td></tr><tr><td>R<sub>p</sub></td><td>...</td><td>3.147<sup>+0.043</sup><sub>-0.042</sub></td><td>...</td><td>R<sub>⊕</sub></td></tr><tr><td>ρ<sub>p</sub></td><td>...</td><td>3.54<sup>+0.31</sup><sub>-0.3</sub></td><td>...</td><td>g cm<sup>-3</sup></td></tr><tr><td>g<sub>p</sub></td><td>...</td><td>19.9<sup>+1.7</sup><sub>-1.6</sub></td><td>...</td><td>m s<sup>-2</sup></td></tr><tr><td>a<sub>p</sub></td><td>0.0585<sup>+0.0011</sup><sub>-0.0012</sub></td><td>0.1322<sup>+0.003</sup><sub>-0.003</sub></td><td>0.2906<sup>+0.0055</sup><sub>-0.0057</sub></td><td>au</td></tr><tr><td>T<sub>eq, p</sub><sup>(a)</sup></td><td>1095.0<sup>+17.0</sup><sub>-17.0</sub></td><td>728.0<sup>+11.0</sup><sub>-11.0</sub></td><td>491.4<sup>+7.4</sup><sub>-7.2</sub></td><td>K</td></tr><tr><td>S</td><td>241.0<sup>+12.0</sup><sub>-11.0</sub></td><td>47.3<sup>+2.5</sup><sub>-2.4</sub></td><td>9.8<sup>+0.47</sup><sub>-0.43</sub></td><td>S<sub>⊕</sub></td></tr><tr><td>TSM</td><td>...</td><td>35.6<sup>+3.3</sup><sub>-2.9</sub></td><td>...</td><td>...</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>Assuming zero Bond albedo and perfect energy redistribution.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Top: TESS (left) and CHEOPS (right) photometry phase-folded to the 18.26-day period of the transiting planet TOI-2141 b. Bottom: RVs phase-folded to the period of all planets (TOI-2141 c, left; TOI-2141 b, center; TOI-2141 d, right). Residuals are shown in the bottom of each panel. In all panels, the black line is the best-fit model from the joint fit described in Sect. 3.4.</figcaption></figure>

and the stellar properties from Table 2. In summary, we refined the radius and mass determination of the known transiting planet TOI-2141 b (with 1.3% and 7.9% radius and mass uncertainties, respectively<fnref n="6" />) and reported the discovery of an inner and an outer non-transiting companion to the known planet.

Compared to Martioli et al. (2023), we find that TOI-2141 b is slightly larger (3.147 vs. 3.05 $R_⊕$) and less massive (20.1 vs. 24 $M_⊕$), making the planet puffier (3.5 vs. 4.6 g cm⁻³) although all differences are within the 1σ uncertainties reported in that work. Figure 4 shows the position of the planet in a mass-radius diagram with respect to other transiting planets orbiting Sun-like stars at relatively long-period orbits. We note that the precision on the bulk properties of TOI-2141 b is significantly higher than that of the comparison population and that the planet is still relatively dense compared to similar sub-Neptunes.

## 4.1. Interior structure and atmospheric evolution of TOI-2141 b

We used the publicly available tool ExoMDN (Baumeister & Tosi 2023), which uses machine learning techniques, particularly mixture density networks, to characterize planets quickly, to model the interior structure of TOI-2141 b in more detail. The

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Mass-radius diagram of transiting planets in the PlanetS catalog (Otegi et al. 2020; Parc et al. 2024) as of May 2025. The internal structure models are from Zeng et al. (2019). The location of TOI-2141 b is marked in red. For this planet, the marker size is larger than the uncertainty in the y-axis. Transiting planets with orbital periods longer than 15 days orbiting Sun-like stars (5500 K < T<sub>eff</sub> < 5900 K) are shown in blue. Compared to similar planets, the precision in the bulk properties of TOI-2141 b is remarkable.</figcaption></figure>

<footnote n="6">⁶ For comparison, the precision in radius and mass reported in Martioli et al. (2023) is of 7.5% and 16.7%, respectively.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** ExoMDN computed mass fractions of the core, mantle, water and atmosphere layers for the our derived planetary parameters in pink and the planetary parameters derived by Martioli et al. (2023) in purple.</figcaption></figure>

code is trained on ~5.6 million synthetic planets within mass ranges from 0.1–25 M⊕ and equilibrium temperatures from 100–1000 K. ExoMDN computes the interior structure of a planet using TATOOINE (Baumeister et al. 2020) and assumes that the planet consists of four distinct layers, which is a strong simplification for sub-Neptune planets (see e.g., Luo et al. 2024). The centre of the planet is found to be made up of a core consisting of solid iron. This core is surrounded by the second layer, a silicate mantle divided into an upper and lower mantle, assumed to have an Earth-like composition. On top of this is a water layer, which follows the equations of state from AQUA (Haldemann et al. 2020). Finally, the outer layer is an atmosphere assumed to be made up of hydrogen and helium with solar ratios. The atmosphere is treated as isothermal in the model at the equilibrium temperature of the planet.

We used the derived radius, mass, and equilibrium temperature from Table 5 as our input and use the default sample of 1000 uncertainty samples to predict the interior structure distribution for inputs within the error bars. We computed ten samples of each distribution to reach a total sample of 10 000. The resulting mass fractions are shown in pink in Fig. 5. ExoMDN predicts a core mass fraction of $0.103^{+0.144}_{-0.077}$ and a mantle mass fraction of $0.20^{+0.21}_{-0.14}$, while the water mass fraction lies at $0.67^{+0.14}_{-0.32}$ and the atmospheric mass fraction of $0.00039^{+0.01636}_{-0.00039}$. These errors may be underestimated given that the errors of an individual sample are from the search algorithm, not accounting for physical uncertainties.

We applied ExoMDN in the same way to model the interior structure with the planetary parameters derived by Martioli et al. (2023) as input. We obtained a core mass fraction of $0.15^{+0.22}_{-0.11}$, a mantle mass fraction of $0.26^{+0.30}_{-0.19}$, water mass fraction of $0.50^{+0.25}_{-0.29}$ and an atmospheric mass fraction of $0.00024^{+0.01743}_{-0.00024}$. These values are in agreement with the interior structure derived from the planetary parameters derived in this work. A comparison between the two is shown in Fig. 5, where the pink line represents the results from our work, while the purple line is the ExoMDN results of the planetary parameters from Martioli et al. (2023). This shows that our newly derived radius and mass allow us to better constrain the layers.

In addition, we modeled the atmospheric evolution of TOI-2141 b using the Planetary Atmosphere and Stellar RoTation RAte (PASTA; Bonfanti et al. 2021b) code, which works in an MCMC fashion within the MC3 Bayesian framework developed by Cubillos et al. (2017). The algorithm employs MESA evolutionary tracks (Choi et al. 2016), a gyrochronological relation in the form of a broken power-law (Bonfanti et al. 2021b),

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Posterior distributions (pdf) as derived from PASTA (blue histograms), with the violet area that highlights the 68.3% highest probability density interval. Left panel: rotation period pdf of TOI-2141 when the star was 150 Myr old ($P_{\text{rot,150}}$) compared with the $P_{\text{rot,150}}$ distribution of stars of comparable masses (black histogram) as extracted from the sample by Johnstone et al. (2015a). Right panel: initial atmospheric mass fraction ($f_{\text{atm}}^{\text{start}}$) pdf of TOI-2141 b compared with the model-inferred present-day atmospheric content (black line). The orange horizontal line marks the uniform prior that was imposed on $f_{\text{atm}}^{\text{start}}$.</figcaption></figure>

an empirical model to relate stellar rotation with the X-ray and extreme ultraviolet (XUV) flux (Bonfanti et al. 2021b, and references therein), a planetary structure model that links planetary observables to the atmospheric mass (Johnstone et al. 2015b), and a model for computing the atmospheric escape rate, which results from internal heating and XUV stellar irradiation (Kubyshkina et al. 2018a,b). PASTA neglects any migration that might have occurred after the proto-planetary disk dispersal and assumes that the atmospheric content of the planets (if any) is given by H-dominated atmospheres.

At each MCMC step, PASTA computes the evolutionary track of the exoplanet atmospheric mass fraction $f_{\text{atm}}(t)$ and the step is accepted if the present-day planet radius inferred from $f_{\text{atm}}(t_*)$ is consistent with its observational counterpart. As a result, we obtained the posterior distribution for the initial atmospheric mass fraction $f_{\text{atm}}^{\text{start}} = f_{\text{atm}}(t_{\text{disk}})$ of TOI-2141 b (blue histogram in the second panel of Fig. 6), where $t_{\text{disk}}$ is the time of dispersal of protoplanetary disk. According to our framework, TOI-2141 b has lost an atmospheric content roughly equal to 1% of its total mass during its lifetime. By modeling the stellar activity, PASTA also outputs the expected stellar rotation period at $t = 150$ Myr and compares it with the period distribution of coeval stars having masses comparable to TOI-2141 as extracted from Johnstone et al. (2015a) (thick blue and thin black histograms in the first
panel of Fig. 6, respectively). We used these results to infer that
the star was likely born as a slow rotator.

4.2. Considering whether the RV-only planets transit

This is a reasonable consideration given that the innermost
planet in the system, TOI-2141 c, has a transit probability of
7.5 ± 0.1%, using the values in Table 5 (and 100% if the orbits
would be perfectly coplanar). For comparison, planet b has a
3.3% transit probability and d just 1.5%. To investigate this ques-
tion, we carried out a joint fit as in Sect. 3.4, but assuming that
all three planets are transiting. Our priors were thus identical as
in Table 4, but adding the parameters r₁ and r₂ for both planets c
and d. The posteriors from this analysis are virtually identical to
those shown in Table 4 except for those parameters and are thus
not reported in this manuscript.

Figure A.5 shows the TESS photometry phase-folded to the periods of the c and d planets from that joint fit. There is no evidence that either of the planets were transiting during the time of TESS observations. The CHEOPS data, scheduled to coincide with transits of planet b, do not cover either of the two planets' predicted transit windows. For the innermost planet, this is surprising given the closeness to the star, but a model that assumes the planet to be transiting is significantly disfavored ($\Delta \ln Z = 351679.2 - 351693.6 = -14.4$). If planet c were to be transiting, our joint fit would constrain its size to be smaller than $0.34_{-0.23}^{+0.32} R_\oplus$, which combined with its inferred mass of ~6 $M_\oplus$ would imply a density of 850 g cm$^{-3}$, more than a hundred times denser than iron. For planet d, the inferred density would be even larger as the planet is more massive than c although its long period makes it more likely to be non-transiting from our line of sight. The minimum mutual inclination between planets b and c given these constrains $\Delta i_{\min} \approx \cos^{-1}\left(\frac{R_c}{a_b}\right) - \cos^{-1}\left(\frac{R_c}{a_c}\right) \ge 2.4$ deg.

4.3. System stability and dynamical constrains

The TOI-2141 system hosts three sub-Neptune-mass planets in relatively well-spaced orbits. Interestingly, the system seems to be close to an unusual 10/3 mean motion resonance chain, so the stability is not straightforward or warranted. To get a clear view on the stability of the system, we performed a dynamical analysis as in Correia et al. (2005, 2010). The system was integrated on a regular 2D mesh of initial conditions around the best fit (Table 4). Each initial condition was integrated for $10^4$ yr, using the symplectic integrator SABAC4 (Laskar & Robutel 2001), with a step size of $2.5 \times 10^{-4}$ yr and general relativity corrections. Then, we performed a frequency analysis (Laskar 1990, 1993) of the mean longitude of planet b over two consecutive time intervals of 5000 yr, and determined the main frequency, *n* and *n*', respectively. The stability is measured by $\Delta = |1 - n'/n|$, which estimates the chaotic diffusion of the orbits. The stability results are visually reported as color in Figs. 7 and 8. Orange and red represent chaotic, strongly unstable trajectories; yellow indicates the transition between stable and unstable regimes; green corresponds to moderately chaotic trajectories, but stable on Gyr timescales; cyan and blue give extremely stable quasi-periodic orbits.

In Fig. 7, we explore the stability of planet b by varying its
orbital period and eccentricity. We note that although the orbital
period of planet b is very well constrained, we opted to vary
this parameter, so that we could simultaneously view both 10/3
period ratios with planets c and d (marked with dashed lines). We

<figure><img src="image_8.png" /></figure>

**Fig. 7.** Stability analysis of the TOI-2141 planetary system in the plane (*P*<sub>b</sub>, *e*<sub>b</sub>), assuming coplanar orbits. For fixed initial conditions, the parameter space of the system is explored by varying the orbital period and the eccentricity of planet-b. The step size is 0.005 day in the orbital period and 0.0025 in the eccentricity. For each initial condition, the system is integrated over $10^4$ yr, and a stability indicator is calculated, which involves a frequency analysis of the mean longitude of planet b. The chaotic diffusion is measured by the variation in the frequency (see text). Red points correspond to highly unstable orbits, while blue points correspond to orbits that are likely to be stable on Gyr timescales. The vertical dashed lines correspond to both 10/3 period ratios with planets c and d. The black dot shows the best fit solution from Sect. 3.4.

observe that the best-fit solution (black dot) is close to, but still
outside both resonances, even if we include the observational
errors. Thus, the resonant motion can be discarded. Besides, the
best fit is embedded in a large stability region (blue); thus, we
can conclude that the orbits of all planets are trustworthy from a
stability point of view.

In Fig. 8, we explore the stability of the system by varying the eccentricities of planets b and d, which are the two most massive in the system. Although in the best-fit solution their eccentricities were fixed at zero, we observe that the system is resilient even for high values, up to about 0.3. We vary only $e_b$ and $e_d$, because planet c is the closest to the star, and so it is more likely that $e_c \sim 0$ owing to tidal damping.

Figure 8 shows the inclinations to the line-of-sight of planets c and d. We vary only $i_c$ and $i_d$, because planet b is transiting, and thus its inclination is well determined, $i_b = 88.76 \pm 0.06$ deg (Table 5). For guidance, we also show a dashed circle that corresponds to mutual inclinations with planet b of 60 deg. We observe that planet c can be stable for $i_c > 15$ deg, so its maximal mass is about 24 $M_\oplus$, while planet d can be stable for $i_d > 30$ deg, so its maximal mass is about 30 $M_\oplus$. This analysis shows all planets in the system, despite two being non-transiting, have masses below 30 $M_\oplus$ from stability arguments.

5. Conclusions

We refined the characterization of the TOI-2141 system using
new TESS and CHEOPS photometry, along with 61 HARPS-N
RVs. For the known transiting planet TOI-2141 b, we deter-
mined a radius of $3.15 \pm 0.04 R_\oplus$ and mass of $20.1 \pm 1.6 M_\oplus$,
corresponding to a bulk density of $3.54 \pm 0.31 \text{ g cm}^{-3}$. Interior
structure models suggest a substantial water-rich layer and a
thin hydrogen-dominated atmosphere, but stronger conclusions
cannot be made without future observations focusing on the
atmospheric composition. The high precision achieved in the
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** *Left:* stability analysis in the plane $(e_b, e_d)$, assuming coplanar orbits. For fixed initial conditions, the parameter space of the system is explored by varying the eccentricities of planets b and d with a step size of 0.0025. *Right:* stability analysis in the plane $(i_c, i_d)$, assuming $i_b = 88.76$ deg (Table 5) and the longitude of the node of all planets equal to zero. For fixed initial conditions, the parameter space of the system is explored by varying the inclinations of planets c and d with a step size of 1 deg. The dashed circle corresponds to mutual inclinations with planet b of 60 deg. The color codes are the same as in Fig. 7.</figcaption></figure>

planetary mass and radius of TOI-2141 b (7.9% and 1.3%, respectively) makes it a valuable benchmark for studies of sub-Neptune composition and evolution.

Bayesian model comparison analyses on the RV data strongly support a three-planet model, with additional Keplerian signals at 5.46 and 60.45 days. Despite the high transit probability of the inner companion, a joint analysis rules out transits of planets c and d with high confidence. This implies a minimum mutual inclination of ≥2.4° between planets b and c, suggesting a slightly misaligned but dynamically quiet system. This system adds to the growing body of evidence that many single-transiting planets are part of multi-planet systems with small, but significant inclinations.

## Data availability

The raw and detrended photometric time-series data are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/704/A174.

*Acknowledgements.* R.L. is funded by the European Union (ERC, THIRSTEE, 101164189). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. R.L. is supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51559.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. R.L. acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033. Data for this paper has been obtained under the International Time Programme of the CCI (International Scientific Committee of the Observatorios de Canarias of the IAC) with the Telescopio Nazionale Galileo operated on the island of La Palma by the Fundacion Galileo Galilei - INAF, Fundacion Canaria in the Observatorio Roque de los Muchachos. CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. The CHEOPS Consortium would like to gratefully acknowledge the support received by all the agencies, offices, universities, and industries involved. Their flexibility and willingness to explore

new approaches were essential to the success of this mission. CHEOPS data analysed in this article will be made available in the CHEOPS mission archive (https://cheops.unige.ch/archive browser/). K.W.F.L. was supported by Deutsche Forschungsgemeinschaft grants RA714/14-1 within the DFG Schwerpunkt SPP 1992, Exploring the Diversity of Extrasolar Planets. ABr was supported by the SNSA. ACMC acknowledges support from the FCT, Portugal, through the CFisUC projects UIDB/04564/2020 and UIDP/04564/2020, with DOI identifiers 10.54499/UIDB/04564/2020 and 10.54499/UIDP/04564/2020, respectively. A.C., A.D., B.E., K.G., and J.K. acknowledge their role as ESA-appointed CHEOPS Science Team Members. DG gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 "Gaseousor rocky? Unveiling the nature of small worlds". 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. S.G.S. acknowledge support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC). The Portuguese team thanks the Portuguese Space Agency for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number 4000142255. TWi acknowledges support from the UKSA and the University of Warwick. YAI acknowledges support from the Swiss National Science Foundation (SNSF) under grant 200020_192038. RAI, DBa, EPa, IRi, and EVi acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF "A way of making Europe" through projects PID2021-125627OB-C31, PID2021-125627OB-C32, PID2021-127289NB-I00, PID2023-150468NB-I00 and PID2023-149439NB-C41. SCCB acknowledges the support from Fundação para a Ciência e Tecnologia (FCT) in the form of work contract through the Scientific Employment Incentive program with reference 2023.06687.CEECIND and DOI 10.54499/2023.06687.CEECIND/CP2839/CT0002. LBo, VNa, IPa, GPi, RRa, and GSc acknowledge support from CHEOPS ASI-INAF agreement n. 2019-29-HH.0. CBr and ASI acknowledge support from the Swiss Space Office through the ESA PRODEX program. ACC acknowledges support from STFC consolidated grant number ST/V000861/1, and UKSA grant number ST/X002217/1. P.E.C. is funded by the Austrian Science Fund (FWF) Erwin Schroedinger Fellowship, program J4595-N. This project was supported by the CNES. A.De. This work was supported by FCT - Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 through the research grants UIDB/04434/2020, UIDP/04434/2020, 2022.06962.PTDC. O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. B.-O.D. acknowledges support from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00046. ADe, BEd, KGa, and JKo acknowledge their role as ESA-appointed CHEOPS Science Team Members. This project has received funding from the Swiss National Science Foundation for project 200021_200726. It has also been carried out within the
framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grant 51NF40_205606. The authors acknowledge the financial support of the SNSF. M.F. and C.M.P. gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19, 174/18). M.G. is an F.R.S.-FNRS Senior Research Associate. M.N.G. is the ESA CHEOPS Project Scientist and Mission Representative. BMM is the ESA CHEOPS Project Scientist. KGI was the ESA CHEOPS Project Scientist until the end of December 2022 and Mission Representative until the end of January 2023. All of them are/were responsible for the Guest Observers (GO) Programme. None of them relay/relayed proprietary information between the GO and Guaranteed Time Observation (GTO) Programmes, nor do/did they decide on the definition and target selection of the GTO Programme. CHe acknowledges financial support from the Österreichische Akademie 1158 der Wissenschaften and from the European Union H2020-MSCA-ITN-2019 1159 under Grant Agreement no. 860470 (CHAMELEON). Calculations were performed using supercomputer resources provided by the Vienna Scientific Cluster (VSC). This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. 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. AL and JK acknowledge support of the Swiss National Science Foundation under grant number TMSGI2_211697. ML acknowledges support of the Swiss National Science Foundation under grant number PCEFP2_194576. PM acknowledges support from STFC research grant number ST/R000638/1. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). NCSA acknowledges funding by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. GyMSz acknowledges the support of the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-125015, a PRODEX Experiment Agreement No. 4000137122, the Lendület LP2018-7/2021 grant of the Hungarian Academy of Science and the support of the city of Szombathely. V.V.G. is an F.R.S-FNRS Research Associate. J.V. acknowledges support from the Swiss National Science Foundation (SNSF) under grant PZ00P2_208945. N.A.W. acknowledges UKSA grant ST/R004838/1. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation.

## References

Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O'Neil, M. 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 252

Anglada-Escudé, G., & Butler, R. P. 2012, ApJS, 200, 15

Baumeister, P., & Tosi, N. 2023, A&A, 676, A106

Baumeister, P., Padovan, S., Tosi, N., et al. 2020, ApJ, 889, 42

Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron., 51, 109

Blackwell, D. E., & Shallis, M. J. 1977, MNRAS, 180, 177

Bonfanti, A., Ortolani, S., Piotto, G., & Nascimbeni, V. 2015, A&A, 575, A18

Bonfanti, A., Ortolani, S., & Nascimbeni, V. 2016, A&A, 585, A5

Bonfanti, A., Delrez, L., Hooton, M. J., et al. 2021a, A&A, 646, A157

Bonfanti, A., Fossati, L., Kubyshkina, D., & Cubillos, P. E. 2021b, A&A, 656, A157

Borsa, F., Scandariato, G., Rainer, M., et al. 2015, A&A, 578, A64

Bouchy, F., Díaz, R. F., Hébrard, G., et al. 2013, A&A, 549, A49

Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, in European Physical Journal Web of Conferences, 47, European Physical Journal Web of Conferences, 03005

Castelli, F., & Kurucz, R. L. 2003, in IAU Symposium, 210, Modelling of Stellar Atmospheres, eds. N. Piskunov, W. W. Weiss, & D. F. Gray, A20

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Correia, A. C. M., Udry, S., Mayor, M., et al. 2005, A&A, 440, 751

Correia, A. C. M., Couetdic, J., Laskar, J., et al. 2010, A&A, 511, A21

Cosentino, R., Lovis, C., Pepe, F., et al. 2012, SPIE Conf. Ser., 8446, 84461V

Cosentino, R., Lovis, C., Pepe, F., et al. 2014, SPIE Conf. Ser., 9147, 91478C

Cubillos, P., Harrington, J., Loredo, T. J., et al. 2017, AJ, 153, 3

Espinoza, N. 2018, RNAAS, 2, 209

Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, celerite: Scalable 1D Gaussian Processes in C++, Python, and Julia

Fortier, A., Simon, A. E., Broeg, C., et al. 2024, A&A, 687, A302

Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Haldemann, J., Alibert, Y., Mordasini, C., & Benz, W. 2020, A&A, 643, A105

Hara, N. C., & Ford, E. B. 2023, Annu. Rev. Statist. Appl., 10, 623

Hoyer, S., Guterman, P., Demangeon, O., et al. 2020, A&A, 635, A24

Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, Proc. SPIE, 9913, 99133E

Johnson, D. R. H., & Soderblom, D. R. 1987, AJ, 93, 864

Johnstone, C. P., Güdel, M., Brott, I., & Lüftinger, T. 2015a, A&A, 577, A28

Johnstone, C. P., Güdel, M., Stökl, A., et al. 2015b, ApJ, 815, L12

Kemmer, J., Stock, S., Kossakowski, D., et al. 2020, A&A, 642, A236

Kipping, D. M. 2013, MNRAS, 435, 2152

Kreidberg, L. 2015, PASP, 127, 1161

Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018a, ApJ, 866, L18

Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018b, A&A, 619, A151

Kurucz, R. L. 1993, SYNTHE spectrum synthesis programs and line data (Astrophysics Source Code Library)

Laskar, J. 1990, Icarus, 88, 266

Laskar, J. 1993, Physica D Nonlinear Phenomena, 67, 257

Laskar, J., & Robutel, P. 2001, Celest. Mech. Dyn. Astron., 80, 39

Lindgren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4

Luo, H., Dorn, C., & Deng, J. 2024, Nat. Astron., 8, 1399

Luque, R., Pallé, E., Kossakowski, D., et al. 2019, A&A, 628, A39

López, E. D., & Fortney, J. J. 2014, ApJ, 792, 1

Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264

Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77

Martioli, E., Hébrard, G., de Almeida, L., et al. 2023, A&A, 680, A84

Maxted, P. F. L., Ehrenreich, D., Wilson, T. G., et al. 2022, MNRAS, 514, 77

McQuillan, A., Mazeh, T., & Aigrain, S. 2014, ApJS, 211, 24

Mortier, A., & Collier Cameron, A. 2017, A&A, 601, A110

Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python

Nixon, M. C., & Madhusudhan, N. 2021, MNRAS, 505, 3414

Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43

Parc, L., Bouchy, F., Venturini, J., Dorn, C., & Helled, R. 2024, A&A, 688, A59

Perruchot, S., Kohler, D., Bouchy, F., et al. 2008, SPIE Conf. Ser., 7014, 70140J

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, Proc. SPIE, 9143, 914320

Rogers, L. A., & Seager, S. 2010, ApJ, 712, 974

Salmon, S. J. A. J., Van Grootel, V., Buldgen, G., Dupret, M. A., & Eggenberger, P. 2021, A&A, 646, A7

Santos, N. C., Sousa, S. G., Mortier, A., et al. 2013, A&A, 556, A150

Schanche, N., Hébrard, G., Collier Cameron, A., et al. 2020, MNRAS, 499, 428

Schönfeld, G. 1886, Eds Marcus and Weber's Verlag, 0

Scuflaire, R., Théado, S., Montalbán, J., et al. 2008, Ap&SS, 316, 83

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000

Sneden, C. A. 1973, PhD thesis, The University of Texas at Austin

Sousa, S. G. 2014, arXiv e-prints [arXiv:1407.5817]

Sousa, S. G., Santos, N. C., Israelian, G., Mayor, M., & Monteiro, M. J. P. F. G. 2007, A&A, 469, 783

Sousa, S. G., Santos, N. C., Mayor, M., et al. 2008, A&A, 487, 373

Sousa, S. G., Santos, N. C., Adibekyan, V., Delgado-Mena, E., & Israelian, G. 2015, A&A, 577, A67

Sousa, S. G., Adibekyan, V., Delgado-Mena, E., et al. 2021, A&A, 656, A53

Speagle, J. S. 2020, MNRAS, 493, 3132

Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018, AJ, 156, 102

Stock, S., Kemmer, J., Kossakowski, D., et al. 2023, A&A, 674, A108

Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985

Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100

Trotta, R. 2008, Contemp. Phys., 49, 71

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ, 145, 44

Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577

Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12

Zeng, L., Sasselov, D. D., & Jacobsen, S. B. 2016, ApJ, 819, 127

Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, PNAS, 116, 9723

<footnote n="1">¹ Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain</footnote>

<footnote n="2">² Institute of Space Research, German Aerospace Center (DLR), Rutherfordstrasse 2, 12489 Berlin, Germany</footnote>
³ Space Research Institute, Austrian Academy of Sciences, Schmiedl-Strasse 6, 8042 Graz, Austria

⁴ Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK

⁵ Department of Astronomy, Stockholm University, AlbaNova University Center, 10691 Stockholm, Sweden

⁶ Space Research and Planetary Sciences, Physics Institute, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland

⁷ Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland

⁸ Observatoire astronomique de l'Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland

⁹ CFISUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal

¹⁰ Dipartimento di Fisica, Università degli Studi di Torino, via Pietro Giuria 1, 10125 Torino, Italy

¹¹ ETH Zurich, Department of Physics, Wolfgang-Pauli-Strasse 2, 8093 Zurich, Switzerland

¹² Space sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

¹³ Instituto de Astrofísica e Ciencias do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal

¹⁴ Space Research and Planetary Sciences, Physics Institute, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland

¹⁵ Instituto de Astrofísica de Canarias, Vía Láctea s/n, 38200 La Laguna, Tenerife, Spain

¹⁶ Departamento de Astrofísica, Universidad de La Laguna, Astrofísico Francisco Sanchez s/n, 38206 La Laguna, Tenerife, Spain

¹⁷ European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

¹⁸ Admatis, 5. Kandó Kálmán Street, 3534 Miskolc, Hungary

¹⁹ Depto. de Astrofísica, Centro de Astrobiología (CSIC-INTA), ESAC campus, 28692 Villanueva de la Cañada (Madrid), Spain

²⁰ Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal

²¹ Max Planck Institute for Extraterrestrial Physics, Gießenbachstraße 1, 85748 Garching bei München

²² INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy

²³ Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK

²⁴ INAF, Osservatorio Astrofisico di Torino, Via Osservatorio, 20, 10025 Pino Torinese To, Italy

²⁵ Centre for Mathematical Sciences, Lund University, Box 118, 221 00 Lund, Sweden

²⁶ Aix Marseille Univ, CNRS, CNES, LAM, 38 rue Frédéric Joliot-Curie, 13388 Marseille, France

²⁷ ARTORG Center for Biomedical Engineering Research, University of Bern, Bern, Switzerland

²⁸ ELTE Gothard Astrophysical Observatory, 9700 Szombathely, Szent Imre h. u. 112, Hungary

²⁹ SRON Netherlands Institute for Space Research, Niels Bohrweg 4, 2333 CA Leiden, The Netherlands

³⁰ Centre Vie dans l'Univers, Faculté des sciences, Université de Genève, Quai Ernest-Ansermet 30, 1211 Genève 4, Switzerland

³¹ Leiden Observatory, University of Leiden, PO Box 9513, 2300 RA Leiden, The Netherlands

³² Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden

³³ National and Kapodistrian University of Athens, Department of Physics, University Campus, Zografos 157 84, Athens, Greece

³⁴ Astrobiology Research Unit, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

³⁵ Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17, 1180 Vienna, Austria

³⁶ Institute for Theoretical Physics and Computational Physics, Graz University of Technology, Petersgasse 16, 8010 Graz, Austria

³⁷ Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, 1121 Budapest, Konkoly Thege Miklós út 15-17, Hungary

³⁸ ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary

³⁹ IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ., 77 av. Denfert-Rochereau, 75014 Paris, France

⁴⁰ Institut d'astrophysique de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98bis blvd. Arago, 75014 Paris, France

⁴¹ Centro di Ateneo di Studi e Attività Spaziali G. Colombo – Università degli Studi di Padova, Via Venezia 15, 35131 Padova, Italy

⁴² Astrophysics Group, Lennard Jones Building, Keele University, Staffordshire ST5 5BG, UK

⁴³ European Space Agency, ESA – European Space Astronomy Centre, Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain

⁴⁴ INAF, Osservatorio Astrofisico di Catania, Via S. Sofia 78, 95123 Catania, Italy

⁴⁵ Weltraumforschung und Planetologie, Physikalisches Institut, Universität of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland

⁴⁶ Dipartimento di Fisica e Astronomia "Galileo Galilei", Università degli Studi di Padova, Vicolo dell'Osservatorio 3, 35122 Padova, Italy

⁴⁷ Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK

⁴⁸ German Aerospace Center (DLR), Markgrafenstrasse 37, 10117 Berlin, Germany

⁴⁹ Institut fuer Geologische Wissenschaften, Freie Universitaet Berlin, Malteserstrasse 74-100, 12249 Berlin, Germany

⁵⁰ Institut de Ciencies de l'Espai (ICE, CSIC), Campus UAB, Can Magrans s/n, 08193 Bellaterra, Spain

⁵¹ Institut d'Estudis Espacials de Catalunya (IEEC), 08860 Castelldefels (Barcelona), Spain

⁵² HUN-REN-ELTE Exoplanet Research Group, Szent Imre h. u. 112., Szombathely 9700, Hungary

⁵³ Leiden Observatory, University of Leiden, Einsteinweg 55, 2333 CA Leiden, The Netherlands

⁵⁴ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
**Appendix A: Additional figures**

<figure><img src="image_10.png" /><figcaption>**Fig. A.1.** CHEOPS photometry of TOI-2141 analysed using pycheops (see text). The transits of planet c were observed in five separate visits. Top panels show the normalised raw CHEOPS light curve (sky blue), the binned raw light curve (dark blue), and the best fitted GP+transit model (brown line). Bottom panel shows the detrended light curves (sky blue), the binned, detrended light curve and the best-fitted transit model (green).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. A.2.** RV time series from SOPHIE (pink) and HARPS-N (green). The error bars are broken into the internal measurement error and the added jitter in gray. The black line shows the best-fit model from Sect. 3.4.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

**Fig. A.3.** GLS periodograms of the spectral activity indicators from HARPS-N data. For each panel, the horizontal lines show the theoretical 10% (short-dashed line), 1% (long-dashed line), and 0.1% (dot-dashed line) false alarm probability levels. The vertical dashed lines mark the orbital frequencies of the transiting planet ($f_b$ = 0.0547 d⁻¹) and of the signals detected in the RVs at 5.46 and 60 d. **Panels a–c:** Cross-correlation full width half maximum (FWHM), bisector (BIS), and contrast (CTR) computed with YABI. **Panels d–h:** Chromatic index (CRX), differential line width (dLW), $H\alpha$, and Na doublet computed with *serval*. **Panels i–j:** S-index and $R'_{HK}$ index computed by YABI.
<figure><img src="image_13.png" /><figcaption>**Fig. A.4.** Stacked Bayesian GLS periodogram analysis of the HARPS-N dataset. **Top:** Periodogram of the original HARPS-N RV dataset (left) and the HARPS-N RV dataset after modeling of the transiting planet TOI-2141 b (right). The amount of observations is plotted against period, with the color scale indicating the logarithm of the power or probability of a coherent long-lived signal, where darker is more likely. **Bottom:** Power as a function of the number of observations for the three RV signals of interest in the system: the transiting planet at 18.26 d (top), the inner non-transiting planet at 5.46 d (middle), and the outer non-transiting planet at 60.45 d (bottom). The first panel is computed from the original RV dataset, the next two from the residuals of fitting the transiting planet. The monotonic increase of power as a function of number of observations is typical of bonafide planets and not of astrophysical or instrumental origin (Mortier & Collier Cameron 2017).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. A.5.** Phase-folded TESS photometry to the period of planet c (left) and d (right). White circles represent the binned photometry for readability. None of the planets appear to be transiting during the time of our observations.</figcaption></figure>