© The Authors 2026

1 Introduction

Multi-planetary systems are of high interest in planetology. Beyond the question of whether the Solar System is unique or not, multiple extra-solar systems exhibit a large diversity (e.g. Lissauer et al. 2014; Muresan et al. 2024), providing us with new insights into the processes that govern planetary formation and evolution (e.g. Ford 2014; Mishra et al. 2023). Moreover, they enable direct comparative studies both within a single system (e.g. Ciardi et al. 2013; Millholland et al. 2017; Weiss et al. 2018; Otegi et al. 2022) and between systems (e.g. Alibert 2019; Gilbert & Fabrycky 2020; Mishra et al. 2023). Within a given system, the planets formed contemporaneously from a common proto-planetary disc, although at different locations. Their present-day properties may have preserved imprints of their common dynamical history and formation pathways, albeit with a different evolution.

As of February 2026, more than 1026 planetary systems with at least two confirmed planets have been detected¹. These systems provide valuable constraints on a wide range of planetary parameters such as masses, radii, orbital spacing, eccentricities (e.g. Van Eylen & Albrecht 2015), and dynamical stability (e.g. Laskar & Petit 2017; Stalport et al. 2022). However, the current observational picture of transiting planetary system architectures remains strongly affected by observational biases. Transit surveys such as Kepler (Borucki et al. 2010) are more sensitive to compact, short-period systems, while planets on wider orbits remain poorly explored due to their low geometric probability. Combining observational results with population synthesis models has led to important theoretical predictions that now require confirmation by upcoming observations (e.g. Mordasini et al. 2009; Emsenhuber et al. 2021). The Transiting Exoplanet Survey Satellite (TESS) survey expands on the Kepler sample, but still primarily detects short-period planets (P ≤ 50 days). As of today, no theoretical model is able to accurately reproduce the observed population of multi-planetary systems (e.g. Turtelboom et al. 2025).

One intriguing pattern identified in many Kepler systems is the so-called peas-in-a-pod configuration, in which planets are found to be similar in size and evenly spaced in terms of orbital distance (Weiss et al. 2018). This architecture has been interpreted as the outcome of smooth and synchronous formation, followed by convergent migration and resonance trapping within the protoplanetary disc (Izidoro et al. 2017; Goldberg & Batygin 2022). Systems such as TRAPPIST-1 (Gillon et al. 2017) and Kepler-223 (Mills et al. 2016) are typical of this architecture. However, recent studies have questioned the universality of this trend (e.g. Murchikova & Tremaine 2020; Zhu 2020), suggesting that additional observations are required to determine how common it truly is. Notably, Millholland et al. (2022) reported a possible breakdown of the pattern for planets with orbital periods of between 100 and 300 days. Expanding the sample of systems with longer period exoplanets is therefore crucial to test theories of formation and evolution.

In this context, HIP 41378 is a unique planetary system; this bright star hosts five transiting planets with long orbital periods, as revealed by the K2 mission (Vanderburg et al. 2016). The system is composed of three inner planets named b, c, and g (P_b = 15 days, P_c = 31 days, P_g ≈ 64 days) (Leonardi et al. 2025); and three outer planets named d, e, and f, with orbital periods of up to 1.5 years. While planets b and c are mini-Neptunes, planets d, e, and f are larger, reaching up to the size of Saturn. Planet g is not transiting. In this system, HIP 41378 f is of particular interest: with a period of 542 days (T_eq ≈ 294 K) and a radius of 9.2 ± 0.1 R_⊕, it stands out as an extremely low-density ‘super-puff’ planet with a 95% credible upper limit on its density of 0.13 g cm⁻³ (Santerne et al. 2019). The orbital periods and masses of HIP 41378 d and e remain poorly constrained (Grouffal et al. 2022; Sulis et al. 2024). Determining the architecture and compositions of these three outer planets is essential for a complete understanding of the system.

In this paper, we present a global characterisation of the HIP 41378 planetary system based on radial-velocity (RV) observations obtained with the HARPS, HARPS-N, HIRES, and ESPRESSO spectrographs between 2016 and 2025. The goals of this work are to measure the masses of all planets, refine the orbital periods of the long-period planets, and place the system into context with other multi-planetary architectures. HIP 41378 thus serves as an ideal testbed for future investigations of temperate, multi-planet systems. An overview of the system and a summary of previous results are presented in Section 2. Then, the observations are described in Section 3. In Section 4, we analyse the RVs and describe our investigation into the planet’s orbital periods and stellar activity. In Section 5, we present a joint analysis of the photometry and RVs. Finally, we discuss the possible orbital periods of the outer planets, the nature of HIP 41378 f, and the dynamics of the system in Section 6. The conclusions are given in Section 7.

2 System overview

The system around HIP 41378 (V-mag = 8.93) has been the subject of numerous studies, which have progressively refined our understanding of its architecture. Five transiting planets were discovered in the 70-day campaign C5 of the K2 space mission (Vanderburg et al. 2016), and four of them were re-observed in transit during the C18 campaign in 2018 (Berardo et al. 2019; Becker et al. 2019). The two inner sub-Neptunes (planets b and c, with orbital periods of 15.6 and 31.7 days) have been observed in transits during the two K2 campaigns, as well as in eight sectors of TESS and subsequent observations with CHEOPS (Leonardi et al. 2025). The outer part of the system hosts three transiting planets, named HIP 41378 d, e, and f. Planets d and f were observed once in each K2 campaign, leading to a range of possible orbital periods. Since planet e was observed only in the first campaign, its orbital period is unknown. The characterisation of the system improved through the asteroseismology study of the star. Lund et al. (2019) used the K2 C18 short cadence data for the asterosismic study of the star and derived precise stellar parameters: a stellar radius of 1.300 ± 0.009 R_⊙, a stellar mass of ${1.22}_{-0.02}^{+0.03}$ M_⊙, and an age of ${2.07}_{-0.27}^{+0.36}$ Gyr (see Table 1). These precise stellar parameters, combined with transit durations of planets d, e, and f, allowed for improved estimates of their orbital periods, all exceeding 200 days. The period of planet f was further confirmed at 542 days by the first RV analysis Santerne et al. (2019).

The planet HIP 41378 f has drawn considerable attention due to its extremely low density, which remains unexplained by current theories (e.g. Alam et al. 2022; Belkovski et al. 2022). Six transits of this planet have been observed to date: two with K2, one from the ground with NGTS (Bryant et al. 2021), one in transmission spectroscopy with the Hubble Space Telescope (HST) (Alam et al. 2022), one via the Rossiter–McLaughlin effect (Grouffal et al. 2025), and an additional event from the ground with Tierras, TRAPPIST-North, and LCOGT (García-Mejía et al. 2026). Its orbital period has been measured at 542.07 days, and its radius is R_f = 9.2 ± 0.1 R_⊕. RV monitoring by Santerne et al. (2019) constrained its mass to M_p = 12 ± 3 M_⊕, corresponding to an extraordinarily low bulk density of ρ_f = 0.09 ± 0.02 g cm⁻³ for such a temperate, long-period, and mature planet. Several hypotheses have been proposed to explain this super-puff (Lee & Chiang 2016) nature, including the presence of circumplanetary rings (Akinsanmi et al. 2020) or high-altitude photochemical hazes (Wang & Dai 2019; Gao & Zhang 2020). Belkovski et al. (2022) attempted to model HIP 41378 f’s structure with H/He envelope models (Lopez et al. 2012) and hydrostatic equilibrium calculations (Howe et al. 2014; Howe & Burrows 2015), but found that no ring-free scenario was compatible with a stable planetary configuration. Transmission spectroscopy with the HST could not distinguish between the scenarios (Alam et al. 2022). Further atmospheric characterisation with, for instance, the James Webb Space Telescope (JWST) would provide insights into this unusual temperate exoplanet, which differs markedly from the highly irradiated hot Jupiters. Confirming its mass and comparing it with the other planets in the system are crucial steps towards understanding its true nature. The two other long-period planets, d and e, remain incompletely characterised. Only upper limits on the masses of planets d and e were derived in Santerne et al. (2019), and their orbital periods are still uncertain. With only a single observed transit, HIP 41378 e currently has the least information available. HIP 41378 d was tentatively detected via the Rossiter–McLaughlin effect by Grouffal et al. (2022), consistent with a period of 278 days. However, a non-detection with CHEOPS reported by Sulis et al. (2024) casts doubt on this estimate, suggesting either that large transit timing variations (TTVs) led to a missed transit or a different orbital period. Based on TESS observations, which never detected a transit of planet d, Sulis et al. (2024) proposed three possible orbital periods: 278, 371, or 1113 days. Determining the true orbital periods and masses of planets d and e is essential for a complete understanding of the HIP 41378 system. An additional non-transiting planet called HIP 41378 g, with a period of nearly 60 days, was proposed as a candidate by Santerne et al. (2019) based on RV analysis. Recently, thanks to TESS and CHEOPS data, Leonardi et al. (2025) measured TTVs longer than three hours for planet c and confirmed the existence of HIP 41378 g.

Figure 1 places HIP 41378 in the broader context of confirmed multi-planetary systems with at least five known planets. Only a small fraction of multi-planetary system have precise mass and radius measurements, and very few include long-period giants (P > 100 days). HIP 41378 provides a unique opportunity to probe planetary architectures beyond the compact, short-period systems that dominate current populations.

Fig. 1 Confirmed planetary system with more than five planets detected and at least two with determination of both mass and radius. Colours represent planetary masses and circle sizes represent the radius. Planets represented as a cross have no radius measurement, and unfilled circles mark cases where no mass measurement was available. The planetary parameters are from Nasa Exoplanet Archive.

3 Observations

3.1 Photometry

3.1.1 K2

The target HIP 41378 was first observed in photometry during the K2 extended mission (Howell et al. 2014) of the Kepler space telescope (Borucki et al. 2010). It was first monitored continuously over ~75 days from 27 April 2015 to 10 July 2015 as part of the C5 campaign, in long-cadence mode (30 minutes); it was published in Vanderburg et al. (2016) and re-observed over ~52 days from 12 May 2018 to 2 July 2018 as part of the C18 campaign, in short cadence (1 minute) (Berardo et al. 2019; Becker et al. 2019).

3.1.2 TESS

The Transiting Exoplanet Survey Satellite (TESS) observed HIP 41378 during eight sectors: sector 7 (from January 7 to February 1 2020); sector 34 (from January 14 to February 8 2021); sectors 44, 45, and 46 (from October 12 to December 30 2021); sector 61 (from January 18 to February 12 2023); sector 72 (from November 11 to December 7 2023); and sector 88 (from January 14 to February 11 2025). The data from each sector were processed by the TESS Science Processing Operation Center Pipeline (SPOC, Jenkins et al. 2016).

3.1.3 Other photometric observations

Transits of HIP 41378’s planets were also observed by the Spitzer space telescope (Berardo et al. 2019), revealing large TTVs for planet c. These TTVs were further observed by the CHEOPS space telescope and analysed in Leonardi et al. (2025). The transit of planet f was also observed from the ground by NGTS in 2019 (Bryant et al. 2021); by HST in 2021 (Alam et al. 2022); and by the Tierras, LCOGT, and TRAPPIST-N observatories in 2024 (García-Mejía et al. 2026), also revealing large TTVs for this planet.

3.2 Radial velocities

3.2.1 HARPS

We observed HIP 41378 with the HARPS spectrograph mounted on the ESO 3.6 m telescope at La Silla Observatory in Chile (Pepe et al. 2002). Observations were carried out between 2017-01-27 and 2019-04-18 under ESO programme IDs 198.C-0169 and 0102.C-0171. In total, 370 spectra were obtained over 218 nights, spanning a temporal baseline of 811 days. The median photon-noise uncertainty of the nightly binned RVs is 1.88 m s⁻¹, with a mean of 1.95 m s⁻¹, and the RV root-mean-square (rms) amounts to 3.26 m s⁻¹. These were reduced using the latest data reduction software (v3.2.5), and RVs were binned nightly. A calibration noise of 0.4 m s⁻¹ was added in quadrature to the nightly binned uncertainties.

3.2.2 HARPS-N

The target HIP 41378 was also observed with the HARPS-N spectrograph (Cosentino et al. 2012), mounted on the Telescopio Nazionale Galileo at the Roque de Los Muchachos observatory, in La Palma, Spain. A total of 168 spectra from 2016-02-13 to 2019-05-05 across 161 nights, as part of the GTO programme. The nightly binned RVs have a median photon-noise uncertainty of 1.89 m s⁻¹, a mean of 2.01 m s⁻¹, and an RV rms of 3.43 m s⁻¹. The data were reduced with the latest version of the data reduction software (v3.0.1), and RVs were binned nightly. Similarly to HARPS, a calibration noise of 0.4 m s⁻¹ was added in quadrature to the nightly binned uncertainties.

3.2.3 HIRES

The HIRES spectrograph collected a total of 75 nightly binned (218 RV measurements in total) epochs between 2016-10-16 and 2019-05-07, spanning a temporal baseline of 942 days. The nightly binned RVs have a median photon-noise uncertainty of 1.19 m s⁻¹, a mean of 1.36 m s⁻¹, and an RV rms of 4.81 m s⁻¹. These observations were conducted using the C2 decker with a typical signal to noise ratio (S/N) of 200 per pixel (approximately 250 k on the exposure meter, 5-minute exposures). Wavelength calibration was performed with an iodine cell (Butler et al. 1996). To average over short-timescale stellar oscillations, most HIRES observations were obtained as three consecutive exposures. A high-resolution template spectrum was acquired using the B3 decker on 2016 October 10 under 1.1j seeing. This template consisted of a triple exposure with a combined S/N of 340 per pixel (each exposure reaching 250k on the exposure meter) and was taken without the iodine cell. The HIRES data collection, reduction, and analysis were done according to the California Planet Search methodology (Howard et al. 2010).

3.2.4 ESPRESSO

We monitored HIP 41378 with the ESPRESSO spectrograph (Pepe et al. 2021) mounted on the 8.2-m ESO-VLT telescope at the Cerro Paranal observatory, in Chile. Observations were carried out between 2020-12-15 and 2024-02-17 as part of programme ID 105.20 G 0. A total of 272 nightly binned spectra were obtained over 272 nights, spanning a temporal baseline of 1159 days. The nightly binned RVs have a median photon-noise uncertainty of 0.63 m s⁻¹, a mean of 0.68 m s⁻¹, and an RV rms of 2.78 m s⁻¹. The spectra were also reduced with the most up-to-date data-reduction software v3.3.1.

3.2.5 Other spectroscopic observations

The host star was also observed by the PFS and the SOPHIE spectrographs as reported in Santerne et al. (2019). However, only 20 and 11 (respectively) epochs were secured with these instruments over a relatively short time span. Given that HIP 41378 hosts five transiting exoplanets at long orbital periods, as well as extra non-transiting planets and a low-amplitude stellar variability (see Sections 4 and 5), these time series are not long and dense enough to significantly constrain the system parameters. However, they add complexity in the modelling to account for their potential systematics. As a consequence, we did not include them in the subsequent analysis, but we reported them in this manuscript for completeness and transparency.

The star was also observed in spectroscopy during the transit of planet d in 2019 and 2022 (Grouffal et al. 2022) and planet f in 2022 (Grouffal et al. 2025). These in-transit data were not used in the analysis below.

Fig. 2 GLS periodograms of HARPS, HARPS-N, HIRES, and ESPRESSO data. Planet periods b, c, and f are highlighted in pink, as well as the candidate period of planet d (278). Subsequent panels show periodograms after subtracting b, c, and f, and then including h and g. False-alarm probability levels at 1% and 10% are indicated as horizontal dotted black lines.

4 Radial-velocity analysis

We first analysed the RV data from HARPS, HARPS-N, HIRES, and ESPRESSO independently of the photometry, while adopting constraints from the transiting planets with well-determined orbital periods, namely planets b, c, and f. The RV modelling was performed using the recently developed jaxoplanet package (Hattori et al. 2024), which extends the functionality of exoplanet (Foreman-Mackey et al. 2021) and starry (Luger et al. 2019) within the JAX framework. JAX enables high-performance numerical computation with automatic differentiation and just-in-time compilation, significantly improving the efficiency of Markov chain Monte Carlo (MCMC) inference.

4.1 Planet detection

Preliminary exploration of the RV time series using the generalised Lomb–Scargle (GLS) periodograms (Zechmeister & Kürster 2009) revealed significant signals corresponding to planets b (15.5 days), c (31.7 days), and f (542 days), consistent with their transit ephemerides (upper panel of Fig. 2). We modelled these three planets using the full RV dataset, imposing Gaussian priors on their orbital periods and transit times centred on the observed values. RV semi-amplitudes were assigned uniform priors between 0 and 10 m s⁻¹. Eccentricities followed a truncated normal distribution with σ = 0.083 (Van Eylen et al. 2019), and the argument of periastron was uniform over [−π, π], with the parametrisation in e cos ω and e sin ω for numerical stability. Instrumental offsets were given broad normal priors centred on the median RV of each dataset, and jitter terms followed half-normal distributions. Posterior sampling was performed with a Hamiltonian Monte Carlo method using the No-U-Turn Sampler (NUTS; Duane et al. 1987; Hoffman et al. 2014) as implemented in NumPyro (Phan et al. 2019; Bingham et al. 2019). Four chains were run with 2000 warm-up steps and 3000 draws each (target acceptance 0.99; maximum tree depth 12).

After subtracting planets b, c, and f, the GLS periodograms show additional signals at 63.3 days and near 278 days, and a long-period variation beyond 1000 days (second panel of Fig. 2). The latter was modelled with a Keplerian with a uniform prior between 800 and 4000 days, converging to a candidate solution at 2828 ± 287 days (planet candidate h). Including planets b, c, f, and h, the 63-day signal remains highly significant and matches the candidate planet g reported by Santerne et al. (2019) (P = 62.06 ± 0.32 days).

With planets b, c, f, g, and h included, no significant peaks remain in the periodogram at the expected periods of planet d (bottom panel of Fig. 2), and no additional signals are detected. However, the existence of planets d and e is proven by their transits, implying RV semi-amplitudes below the current detection threshold and suggesting that both planets are low mass.

4.2 Model comparison for planets d and e

The orbital periods of HIP 41378 d and e cannot be constrained from photometry alone, and no significant signals are detected in the periodograms that would allow an unambiguous identification in the RV data. Planet d has three candidate periods: 278, 371, and 1113 days (Sulis et al. 2024), the longest of which requires a high eccentricity to reproduce the observed transit. To investigate which of these configurations is most consistent with the RV data, we performed a series of fits and model-comparison tests, progressively increasing the model complexity.

4.2.1 Markov chain Monte Carlo analysis

We first modelled the full RV dataset, including all confirmed planets except planet e, to isolate the contribution of planet d. A broad uniform prior was adopted for the period of planet d between 250 and 400 days, encompassing the two most plausible solutions (278 and 371 days). The resulting posterior distribution (left panel of Fig. 3) is strongly concentrated near P_d ≃ 278 days, while secondary modes at longer periods (e.g. ~385 days) have significantly lower probability densities. This preliminary analysis disfavours the 371-day solution.

However, the 278-day period is inconsistent with the non-detection of a transit by CHEOPS in 2023 (Sulis et al. 2024). We therefore explored three discrete orbital periods for planet d (278, 371, and 1113 days), again excluding planet e. The right panel of Fig. 3 shows the median and 99.5% highest density interval (HDI) for the RV semi-amplitude of planet d under these assumptions. Only the 278-day model yields a significant (> 3σ) non-zero semi-amplitude.

To assess whether the preference for the 278-day signal could instead be attributed to planet e, we extended the analysis by explicitly including planet e and performing model comparison using leave-one-out cross-validation (LOO CV). This approach provides a robust estimate of out-of-sample predictive performance (Vehtari et al. 2015, 2024) and is well suited to high-dimensional RV models. The LOO estimates were computed from the posterior samples obtained with NumPyro – without rerunning the MCMC for each data subset – using Paretosmoothed importance sampling (Vehtari et al. 2015) provided in the Python package ArviZ.

For each assumed value of P_d, we explored different uniform period ranges for P_e, reflecting the possible relative orbital ordering of the two planets: (1) if P_d = 278 days, we tested three possible ranges for P_e: first, assuming planet e lies interior to planet d (100–278 days), second, assuming planet e lies between planets d and f (300–540 days), and third, assuming planet e lies outside planet f (600- 1000 days); (2) if P_d = 371 days, we tested two ranges: 100–370 and 600–1000 days; (3) if P_d = 1113 days, we adopted a broad range: 100–540 days. Only models with well-converged posteriors (i.e. $\hat{R}\approx 1$ and minimal divergences) were included. Here, $\hat{R}$ denotes the Gelman–Rubin statistic, which measures MCMC convergence by comparing within and between-chain variances. The results, ranked by expected log point-wise predictive density (elpd LOO), are shown in Fig. 4.

The best-performing model includes all seven planets with P_d = 278 days and P_e between 300 and 540 days. Models without planet e or with P_d = 371 days remain statistically competitive, but they generally favour a signal near P_e ~ 278 days, which strongly suggests the presence of a planet at this period.

Fig. 3 Tests for the orbital period of HIP 41378 d. Left panel: posterior density of the period of planet d, considering a uniform prior between 250 and 400 days. The vertical dashed lines indicate the possible periods from transit observations. Right panel: resulting RV semi-amplitude and 99% HDI for the three periods tested: 278, 371, and 1113 days.

4.2.2 Nested sampling

Given the low RV semi-amplitude and the limited transit constraints for planets d and e, we performed an independent nested sampling analysis using RV-only data first, with HARPS, HARPS-N, and ESPRESSO data. Orbital periods for planets b, c, and f were constrained with Gaussian priors, while broad uniform priors were adopted for planets g (P_g ∈ [61, 65] days) and the long-period candidate h (P_h ∈ [800, 4000] days), allowing us to focus on the poorly constrained periods of planets d and e.

We used jaxns (Albert 2020) with 900 live points. For planet d, we tested the two possible periods of 278 and 371 days using a uniform prior between 200 and 400. The nested sampling analysis strongly favours the 278-day solution.

For planet e, we adopted a broad uniform prior (P_e ∈ [150, 450] days). The posterior distribution is dominated by the ≈390 day solution, with less significant secondary peaks at shorter periods (see Fig. C.1).

To better constrain the orbital period of planet e, we performed three tests combining photometric and RV data, considering configurations in which planet e is located inside planet d, between planets d and f, or outside of planet f. The configuration in which planet e lies between planets d and f is strongly favoured, with a difference in Bayesian evidence of Δ log Z ≃ 15 compared to the scenario in which planet e is exterior to planet f, and Δ log Z ≃ 6.5 compared to the configuration in which planet e is interior to planet d. These results are consistent with the MCMC analysis and model comparison with LOO-CV. The complete nested sampling analysis can be found in Appendix C.

Fig. 4 Comparison of different RV models, including subsets of the HIP 41378 planets, ranked by their expected log point-wise predictive density (elpd LOO). Each model label indicates the included planets and, for planet d, the tested orbital period in days. Higher elpd LOO values correspond to better out-of-sample predictive performance. The horizontal line segments represent the standard error of the elpd. The lighter line segment with triangles indicates the elpd difference of one model relative to the best model.

4.3 RV analysis and instrument-specific datasets

4.3.1 Global fitting

For the global data fitting, we adopted Gaussian priors on the orbital periods of planets b, c, and f, centred on the transit-derived values with narrow standard deviations. Uniform priors were used for the remaining planets: 60–65 days for g, 250–400 days for d, 360–450 days for e, and 700–4000 days for h. These priors are based on two assumptions based on knowledge from photometry and model comparisons from Section 4.2: the period of planet d is either 278 or 371 days, and we consider that planet e is before planet f in the system, considering its shorter transit duration (≈13 hours) and assuming low eccentricities in the multi-planetary system. The analysis was carried out sequentially, first using HARPS, HARPS-N, and HIRES data as in Santerne et al. (2019), then ESPRESSO data alone, and finally the full dataset combining all four instruments. Again, we used NumPyro to run four MCMC chains in parallel, each with 2000 warm-up steps and 3000 draws. The posteriors for periods and RV semi-amplitudes are shown in Fig. H.1. The RV signal of planet e becomes clearly detectable only when ESPRESSO data are included, reflecting its low semi-amplitude and incomplete phase coverage in the other datasets. The period of planet d consistently converges to 278 days across all instruments. The candidate planet h shows the largest variation between instruments, being most visible in ESPRESSO data.

The preferred model, with all planets included, still shows some correlated residuals, with an rms of 2.60 m s⁻¹. The HIRES data (with 75 measurements over a total of 713) show the largest errors and the largest residual rms, with a weighted rms of 4.22 m s⁻¹. Excluding the HIRES data improves the overall fit, reducing the weighted rms of the residuals to 2.38 m s⁻¹. Nevertheless, some correlations persist and are likely related to stellar activity. We therefore investigate the impact of stellar activity in Section 4.4.

Fig. 5 Upper panel: ℓ1 periodogram of the RV data, including offsets in the base model. Lower panel: ℓ1 periodogram of the RV data, including offsets, transiting planets, and order 11 polynomial in the base model. Planets b, c, g, and f are highlighted by their letters.

4.3.2 Blind search for RV signals

We performed a fully blind search for periodic signals using the ℓ₁ periodogram (Hara et al. 2017), which reconstructs the RV time series as a sparse sum of sinusoids on a predefined frequency grid. The method takes the RV data as input, a frequency grid, a noise covariance model, and a linear base model, and returns a periodogram-like output with reduced aliasing and associated FAP estimates.

We analysed the combined HIRES, HARPS, HARPS-N, and ESPRESSO data over a frequency range of 0–0.95 cycles/day with a step of 1/(10T_obs). The noise was assumed white, with instrument-specific jitter terms fixed to the fitted values (Table I.1). As a linear base model, we first only used the offsets of each instrument. The results are shown in Fig. 5; the known planets are clearly visible, as is a long-period signal at 3250 d. The log₁₀ FAPs of the 15, 31, 63, 562, and 3250 d signal are −15.5, −14.8, −5.2, −22.1, and −15.9. We note that the only other significant peak is at 835 d, with a log₁₀ FAP of −3.6; however, as discussed below, the interpretation of this signal is difficult.

We next refined the analysis in two ways. First, we note that the ℓ₁ periodogram seeks to find a Fourier spectrum that reproduces the data but has as many vanishing values as possible. This was done by minimising the sum of the absolute values of the Fourier amplitudes, that is, the ℓ₁ norm of the Fourier spectrum. As a consequence, amplitudes were underestimated, which might have impacted our ability to detect the smallest signals. To circumvent this issue, following (Hara et al. 2020), we included sinusoids at the periods of the confirmed transiting planets (15.57209, 31.70591, 63.19, 542.07975 d) in the base model. The asymmetric 3250 d peak suggests that a purely periodic model is inadequate. We therefore added a low-amplitude polynomial trend. Testing polynomial degrees from 5 to 15 and selecting via cross-validation (Hara et al. 2020), we adopted degree 11 (lower panel of Fig. 5). In this configuration, no additional signals reach formal significance (FAP > 5%). However, there are signals at 6.5 and 8 days (1.13 d is an alias of 8.48 d), at the expected stellar rotation period (see Section 4.4). In Appendix E, we show that this signal is not stable over time, which supports the hypothesis that it is of stellar origin.

Additionally, the ℓ₁ periodogram exhibits a peak at 279 d, which is fully compatible with the 278 d period of planet d. We finally note a 146 d signal. A planet at this period would continue the resonance chain, and could be due to a low-mass, non-transiting planet (see Section 6.1.3). Although the 279 signal is not statistically significant, finding it at this particular period, without any prior on the position of putative signals, is a strong hint of its planetary origin.

4.4 Stellar activity analysis

4.4.1 Stellar rotation period

The HIP 41378 object is a relatively quiet F-type star. An analysis of the K2 short-cadence photometry by Grouffal et al. (2025) yielded a stellar rotation period of P_rot = 6.8 ± 0.4 days, which was derived using a quasi-periodic Gaussian process (GP) model of the light curve.

The average Ca II chromospheric activity index $\log \left(R_{HK}^{\mathrm{\prime }}\right)$ is measured as −4.88 ± 0.03 for ESPRESSO, −4.94 ± 0.02 for HARPS, and −4.95 ± 0.02 for HARPS-N, indicating a low and stable level of stellar activity over time. Using the empirical relations from Noyes et al. (1984) and Mamajek & Hillenbrand (2008), these values translate into an estimated stellar rotation period of 7.75 ± 1.42 days according to ESPRESSO and 8.29 ± 1.52 days for HARPS and HARPS-N. These estimates are consistent with the photometric rotation period within uncertainties.

4.4.2 Frequency analysis of the activity indicators

We investigated several stellar activity indicators using HARPS, HARPS-N, and ESPRESSO data, including the full width at half maximum (FWHM) of the cross-correlation function (CCF), the bisector span (BIS; Queloz et al. 2001), the CCF contrast, and absorption-line diagnostics such as the Ca II H&K S-index and the Hα line. A detailed description of this analysis is provided in Appendix A. Of particular interest is the long-period RV signal attributed to candidate planet h. If this signal were present in the activity indicators, it could point to a stellar origin, such as a magnetic cycle. All indicators except the BIS exhibit significant long-period power (periods ≳ 900 days), suggesting the presence of long-term stellar activity. However, the orbital period of the candidate planet is longer than the signal observed in the activity indicators (~1000–1700 days).

After subtracting the long-period component from the FWHM, contrast, and S-index using sinusoidal fits, no significant power remains at shorter periods associated with stellar rotation. The BIS only shows a very low-amplitude signal at ~8 days, but a quasi-periodic GP fit does not converge. Overall, no robust rotational signal is detected in any activity indicator.

In summary, HIP 41378 exhibits evidence of long-period stellar activity variability, likely related to a magnetic cycle, which may partially contribute to the RV signal. In contrast, the absence of a clear rotational signature in the activity indicators confirms the star’s low activity level. The star’s inclination of i_* = ${42}_{-7}^{+6}$ deg measured through the Rossiter–McLaughlin effect and considering the rotation period (Grouffal et al. 2025) may suppress the observed rotational modulation, consistent with the low amplitude in both photometry and activity indicators.

4.5 1D GP analysis

After modelling the six confirmed planets and the long-period signal of the candidate planet h, correlated residuals remain in the RVs. The periodogram of the residuals (Fig. B.1) reveals excess power during short periods of ~8 days, similar to the stellar rotation period, suggesting a stellar activity origin. To account for this correlated noise, we added the fit of a GP to the fit of the seven planets in the RV data of the four instruments. GPs are a popular mathematical tool for modelling stellar activity-induced signals that describe stochastic variations (see Aigrain & Foreman-Mackey 2023 for a complete review).

We adopted a quasi-periodic (QP) kernel of the form $$k(t,t^{\mathrm{\prime }})=A^2\exp [-\frac{{(t-t^{\mathrm{\prime }})}^2}{2{\ell }^2}]\exp [-\mathrm{\Gamma }{\sin }^2\left(\frac{\pi (t-t^{\mathrm{\prime }})}{P_{\mathrm{r}\mathrm{o}\mathrm{t}}}\right)],$$(1)

where A is the GP amplitude, ℓ is the exponential decay timescale, P_rot is the QP timescale associated with stellar rotation, and Γ controls the coherence of the periodic component. The GP kernel was implemented with tinyGP (Foreman-Mackey et al. 2024). We adopted uniform priors on all the GP hyperparameters: P_rot ∈ [5, 10] days, A ∈ [0.01, 5.0] m s⁻¹, ℓ ∈ [1, 500] days, and Γ ∈ [0.1, 10]. Posterior sampling was performed using four parallel chains with 2000 warm-up steps and 3000 draws each. The planetary system parameters were given the same priors as in the RV-only analysis. The fit converged to a stellar rotation period of 8.12 ± 0.09 days, consistent with the $\log \left(R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}\right)$ estimate and the signal observed in the residuals’ periodograms. The characteristic evolution timescale of active regions converged towards l = ${100}_{-70}^{+60}$, indicating long-lived surface structures persisting over ≈ 10–15 stellar rotations and compatible with a star with low activity. The Γ parameter converged towards Γ = ${0.94}_{-0.83}^{+0.52}$, which could indicate the signature of a small number of dominant active regions.

We note a modest discrepancy between the photometric rotation period (~6–7 days) and the RV-derived period (~8 days). This difference may be explained by a combination of stellar inclination, differential rotation (Grouffal et al. 2025), and the latitude distribution of active regions. Photometric variability is dominated by equatorial spots, which rotate faster under differential rotation, whereas RVs can be more sensitive to higher latitude structures that may remain visible throughout the stellar rotation.

Including the GP in addition to the planetary system parameters significantly improves the RV fit. The rms of the residuals decreases from 2.38 to 2.16 m s⁻¹, and the GP model is strongly favoured by LOO CV, with Δelpd_loo ≈ +24 relative to the no-GP model. Figure 6 compares the resulting posterior distributions of the RV semi-amplitudes.

Including a GP increases uncertainties on the semi-amplitudes of the long-period planets, as expected. The GP absorbs part of the long-timescale stellar variability, preventing the orbital model from artificially fitting activity-induced trends. This leads to more conservative and realistic estimates of the RV semi-amplitudes: K. For these reasons, we adopted the GP model as the final RV model for subsequent analyses combined with photometric data.

We also tested a multi-dimensional GP framework (Rajpaul et al. 2015) incorporating activity indicators and their derivatives. However, given the absence of a clear rotational signal in the activity indicators, this approach did not yield a significant improvement over the 1D GP model.

Fig. 6 Comparison of the semi-amplitude of the seven planets in the HIP 41378 system for a model with GPs (in blue) and without GPs (in pink). The 68% HDI is represented as a box and the 99% HDI as a horizontal bar.

5 Joint analysis of photometry and RV

To jointly constrain the orbital and physical parameters of the HIP 41378 planetary system, we performed a combined analysis of K2 photometry and RV data from HARPS, HARPS-N, HIRES, and ESPRESSO using jaxoplanet (Hattori et al. 2024). Only K2 light curves were included, as no transits of planets d, e, or f have been detected with TESS or CHEOPS. Independent fits of TESS and CHEOPS data for planets b and c, including transit-timing variations of planet c, are presented in Leonardi et al. (2025). Including these datasets did not improve constraints on the outer planets.

The priors are listed in Table I.1. They assume the best model found with model comparison from the RV-only analysis.

For transiting planets, we imposed informative priors on the orbital inclination and sampled uniformly in cos i. Planetary masses were derived from the inclination-dependent RV semi-amplitudes. Transit depths were parametrised in log-space with Gaussian priors centred on the values from Vanderburg et al. (2016), and planetary radii were computed from these depths. Bulk densities were derived from the sampled mass and radius. Quadratic limb-darkening coefficients were sampled using the parametrisation from Kipping (2013).

For non-transiting planets, the time of inferior conjunction, T₀, was not directly observed. We therefore parametrised the orbital phase using a dimensionless parameter, ϕ, which represents the phase offset relative to the mean observation time. The reference epoch was then defined as t₀ = t_ref + ϕ × P, where P is the orbital period and t_ref is the mean of the observation times. This parametrisation centres the epoch near the dataset, reducing parameter degeneracies and numerical correlations during sampling.

The sampling was performed using NUTS (Duane et al. 1987; Hoffman et al. 2014) as implemented in NumPyro (Phan et al. 2019; Bingham et al. 2019). Three chains were run in parallel, each with 2000 warm-up steps and 3000 draws, with a target acceptance probability of 0.99. Using NUTS in JAX allows for fast gradient-based sampling, even for high-dimensional problems; the full RV analysis with four instruments and 43 free parameters completes in minutes. Achieving stable convergence required non-centred re-parametrisation of latent variables. In hierarchical models, rather than sampling directly from θ ~ 𝒩(μ, σ), we sampled z ~ 𝒩(0,1) and set θ = μ + σ · z. This avoided funnel-shaped posterior geometries and improved sampler efficiency. This approach was applied to both Gaussian and uniform priors. A quasi-periodic GP was included in the RV model to account for stellar activity using the same priors as the 1D-GP analysis.

Figure 7 shows the RV time series with the median posterior model, the GP component, and residuals. Phase-folded RVs for each planet are shown in Figure 9. The main results of the joint fit are summarised in Table 2, with full posteriors and priors in Table I.1. Figure J.1 shows the two K2 light curves with the median of the best-fit model.

All planetary masses are below 10 M_⊕, except for planets f and h. We derived m_b = 6.9 ± 0.5 M_⊕, m_c = ${5.97}_{-0.86}^{+0.56}$ M_⊕, m_g = ${5.75}_{-0.83}^{+0.76}$ M_⊕, m_d = ${6.53}_{-3.15}^{+2.65}$ M_⊕, m_e = ${7.62}_{-4.63}^{+3.20}$ M_⊕, m_f = 25 ± 5 M_⊕, and m_h = ${43}_{-13}^{+16}$ M_⊕. For planets g and h, the reported values are minimum masses due to the unknown inclination. The mass of planet f is significantly higher than the 12 ± 3 M_⊕ previously reported by Santerne et al. (2019), likely due to the inclusion of ESPRESSO data and improved modelling of long-period signals. These mass estimates lead to the following bulk densities for the five transiting planets: ρ_b = 2.07 ± 0.16 g cm⁻³, ρ_c = 1.60 ± 0.21 g cm⁻³, ρ_d = 0.74 ± 0.33 g cm⁻³, ρ_e = ${0.30}_{-0.18}^{+0.13}$ g cm⁻³, and ρ_f = ${0.16}_{-0.04}^{+0.03}$ g cm⁻³. Densities decrease with increasing orbital distance. For reference, Jupiter and Saturn have densities of 1.33 and 0.68 g cm⁻³, respectively. Planets e and f are notably less dense than Saturn. All planets exhibit low eccentricities, except candidate planet h. The uncertainty on the orbital period of the candidate planet h remains large, despite its high RV amplitude; further observations are needed to confirm its planetary nature.

We assessed the completeness of our RV analysis by computing the detection limit on the residuals after subtracting the median best fit of planets b, c, g, d, e, f, and the candidate h, as well as the quasi-periodic GP. We used the diffusive nested sampler kima (Faria et al. 2018) following the approach of Standing et al. (2022), Standing et al. (2023), John et al. (2023), and Balsalobre-Ruza et al. (2025). This method maps all the Keplerian signals compatible with the residuals and defines a detection threshold, which is shown in Figure 8. Signals lying above this threshold would have been detected in our analysis. We find that no additional planets with orbital periods shorter than ~7 years and RV semi-amplitudes larger than 1 m s⁻¹ can be present in the system. Earth-mass planets remain below the sensitivity threshold of our dataset, while long-period companions such as Jupiter- or Saturn-mass planets at periods ≳ 30 years cannot be ruled out. Complementary constraints are provided by astrometric acceleration measurements. HIP 41378 appears in the catalogues of Kervella et al. (2022) and Brandt (2021). No statistically significant proper-motion anomaly is reported. In the separation range of 3–10 au, where the technique reaches its maximum sensitivity, companions more massive than ~2 Jupiter masses can be excluded. These astrometric limits therefore reinforce our RV detection limits at intermediate separations. The details of the method are provided in Appendix D.

Fig. 7 RV timeseries of HIP 41378 with HARPS (pink), HARPS-N (green), HIRES (purple) and ESPRESSO (blue), along with best fit model (median of the posterior) in pink and the GP component in black. The upper panel is a zoom on ESPRESSO data. The lower panel represents the residuals of the fit (data – model).

Fig. 8 Greyscale hex-bin plot representing the posterior samples obtained from kima on the full RV dataset residuals with the number of planets fixed to one. The blue line represents the 99th percentile detection limit. The RV semi-amplitudes of Jupiter, Saturn, Neptune, and Earth are represented by dashed pink lines for comparison. The top axis shows the corresponding angular separation.

6 Discussion

6.1 Architecture of the HIP 41378’s system

The star HIP 41378 hosts a total of five transiting planets, a confirmed non-transiting one, and a candidate planet at very long orbital period. The system appears to be composed of two dynamically distinct sub-systems: an inner system with planets b, c, and g close to a 1:2:4 MMR and an outer system with planets d, e, and f. This separation is also reflected in the planetary inclinations. Planet b has an inclination of i_b = 88.847 ± 0.048 degrees, and planet c has an inclination of i_c = 88.475 ± 0.013 degrees. Assuming planet g is coplanar to planets b and c with an inclination close to 88.5°, we can compute its expected impact parameter knowing its semi-major axis of 0.3307 ± 0.0023. Using b = a cos i/R_*, we find b ≈ 1.05 > 1, implying that planet g does not transit the stellar disc and is not compatible with any transits observed for this planet. The outer system, in contrast, appears to be tilted by roughly one degree relative to the inner system, with inclinations for planets d, e, and f between 89.8 and 89.9 degrees. We find no significant RV signal between the planet g (P_g=63 d) and planet d (P_d=278 d), resulting in an upper-mass limit of about 4–5M_⊕(see Fig. 8). The TTVs of planets b and c were fully explained by the presence of planet g (Leonardi et al. 2025). Hence, the existence of two slightly non-coplanar sub-systems could account for the observed gap in planetary periods between 63 and 278 days.

Fig. 9 HARPS-N, HARPS, HIRES, and ESPRESSO RV phase-folded on the results of the joint RV+transit fit for the seven planets in the system. The pink line represents the median of the best-fit model RV curve with its 68% HDI in the shaded area.

6.1.1 The orbital period of HIP 41378 d

Our RV analysis favours an orbital period of P_d = 278.3618 ± 0.0003 days for HIP 41378 d, consistent with one of the three transit-based solutions proposed to date (Sulis et al. 2024). The period is recovered independently in all RV data subsets and is statistically preferred by cross-validation model comparison. However, the non-detection of a transit by CHEOPS in 2023 (Sulis et al. 2024) remains puzzling. A possible explanation is the presence of large TTVs, as already observed for HIP 41378 f (Bryant et al. 2021) and suggested for HIP 41378 d (Sulis et al. 2024). Dynamical interactions with planets e and f could shift transit times by several hours and up to several days, which is enough to move the expected transit outside the CHEOPS observation window. The TTV analysis remains limited by the small number of observed transits: only seven transits of planet f were detected (García-Mejía et al. 2026, Leonardi et al., in prep.). Only two photometric epochs have been secured for planet d, as well as a third tentative detection in spectroscopy (Grouffal et al. 2022). Without a new transit, P_d cannot be definitively confirmed, although periods near 371 days are strongly disfavoured in the RV analysis. Future transit monitoring with CHEOPS could detect additional transits – enabling a full TTV analysis and a definitive orbital solution allowing for follow-up observation – including transmission spectroscopy to compare with planet f (Alam et al. 2022).

Fig. 10 Possible orbital periods for HIP 41378 e. The posterior distribution of the orbital period of planet e, obtained from the RV-only analysis with a broad uniform prior, is shown as a blue histogram. Period ranges excluded by K2 (pink), TESS (purple), and CHEOPS (orange) photometric observations are indicated as pink shaded areas.

6.1.2 The orbital period of HIP 41378 e

Combining the RV data with K2 photometry yields a most probable orbital period of P_e = ${393}_{-3}^{+4}$ days. As shown in Fig. 10, most periods between 370 and 410 days are excluded by K2, TESS, and CHEOPS, leaving a narrow viable range at around 386–390 days, implying a missed transit shortly after K2 Campaign 18.

Alternative aliases appear in the nested sampling analysis (Appendix C), notably near ~190 and ~760 days, representing half and twice the favoured period, respectively. These likely arise from harmonics or sampling effects of the ~390-day signal. Periods of nearly 200 days were additionally ruled out by TESS.

Another signal at ~146 days was detected in the ℓ₁ periodogram analysis (Section 4.3.2) and has been proposed to complete a near-three-body resonance configuration for all triplets of planets (Section 6.1.3). While the RV phase coverage for a 146-day period is complete, attempting to place planet e at this period using photometric timing did not yield a converging solution. Using a uniform prior of between 100 and 250 days for planet e in the nested sampling analysis, we did not recover a signal near 146 days. The origin of this 146-day signal thus remains unclear and may correspond to another non-transiting planet. Model comparison continues to favour P_e = ${393}_{-3}^{+4}$ days, although incomplete phase coverage and possible yearly systematics warrant further RV monitoring.

This estimate neglects potential TTVs. Given the large TTVs observed for planets d and f, similar variations for e could shift transit times by several hours, reopening some photometrically excluded windows. A definitive confirmation of the orbital period will therefore require either a dedicated transit search or a future joint RV–photometry analysis that includes TTVs, once additional transits of planets d and f become available.

The favoured period lies close to a 5:7 commensurability, with both d and f and exact resonance implying P_e ~ 388.45 days. Within the 386d–390 day interval, the next two opportunities to observe a transit of planet e should occur between BJD=2 461 002 and BJD=2 461 042 (or near BJD=2461414.9 assuming exact resonance), and between BJD=2461388 and BJD=2 461 432 (near BJD=2461803.4 in the resonance case).

Given that HIP 41378 e has a transit depth of 1.2 mmag and duration of 13.2 hours, any ground-based observatory that can secure nightly observations with sub-millimag precision, such as NGTS (Bryant et al. 2021), has a >50% probability of detecting one in-transit observation. Such a detection would significantly refine the ephemeris and enable a targeted space-based follow-up with, for example, CHEOPS. This would constitute only the second observed transit of planet e since the event detected by K 2 in 2015. Only a future transit will unambiguously confirm the precise orbital period of HIP 41378 e.

6.1.3 Analysis of the resonances and possible additional planets

Using the RV posteriors, we examined the proximity of the system to two- and three-body MMR (Fig. 11). The inner c–g pair is consistent with a 2:1 commensurability, while b–c lies near to but not inside resonance (Leonardi et al. 2025).

The compact d–e–f triplet is compatible with several threebody resonances, notably a chain of second-order 7:5 MMRs. In contrast, the gap in period ratios between the pairs g–d and f–h results in three sub-systems (c–b–g, d–e–f, and h) that do not have strong (i.e. resonant) gravitational interactions. We note, however, that the current posterior of h makes a 3:1 MMR with planet f possible.

Given the relatively loose constraints on h, we now focus on the rest of the system and see how missed planets could dynamically reconnect the inner and outer triplets of known planets considering that the period ratios between g and d is ≈4.4. Following Leleu et al. (2021) and Luque et al. (2023), we minimised the distance to two and three-body MMRs when adding planets. We discuss three cases, which are listed below: (i) one additional planet, minimising the distance to two-body MMRs; (ii) one additional planet, minimising the distance to three-body MMRs; and (iii) two additional planets, minimising distance to two-body MMRs.

(i) The gap between g and d could be filled by an additional planet, resulting in a 3:2 MMR followed by a 3:1 MMR, with a planet x at ~95 days. Equivalently, 3:1 followed by 3:2 works, with a planet x at ~188 days; this is the case shown in Figure F.1.

(ii) A slightly more dynamically compact solution can be found by putting a planet x at ~148 days, resulting in a near-three-body resonance for all triplets involving said planet x; however, this results in pairs g–x and x–d at period ratios of 2.34 and 1.88, respectively, which are hence far from significant two-body MMRs. This is the case shown in Figure F.2.

(iii) Adding two planets allowed us to fully reconnect the system through a chain of two-body MMRs (5:3, 4:3, 2:1), leading to six equivalent solutions (Table 3); one is shown in Fig. 11.

The near-resonant nature of b–c–g and d–e–f motivated this search. Consecutive second-order 7:5 MMRs, as suggested for d–e–f, are rare (only TOI-1136 shows a similar case; Dai et al. 2023). While cases (i) and (ii) do not strictly lie inside two-body resonances, this is not necessarily problematic, as b–c itself is slightly offset, and HIP 41378 may represent a partially disrupted chain. Case (iii) produces compact resonant configurations comparable to other known chains, consistent with the system’s puffy and nearly coplanar architecture (Leleu et al. 2024).

Fig. 11 Subsequent triplet of HIP 41378 system with respect to MMRs. The solid purple lines show first-order two-planet MMRs, the dashed purple lines show second-order two-body MMRs, the solid red line shows the 2:1 MMR between planet n and n + 2, and the grey lines show zeroth-order three-body MMRs of the form k/Pn − (k + q)/Pn + 1 + q/Pn+2 = 0 (e.g. Cerioni et al. 2022). The top panel shows the RV posterior of the known planets. The bottom panel shows the example of case (iii) – solution 1, for which we added two planets, x and y, in order to reconnect planets c–f in a resonant chain.

6.1.4 A signal at very long orbital period

We detected a low-amplitude, long-period signal (P_h = ${2602}_{-433}^{+468}$) in the RV data, which is also present, albeit less clearly and at different orbital periods, in some stellar activity indicators. This signal is primarily driven by the ESPRESSO dataset (see Figure H.1), particularly the third observing campaign (November 2022–May 2023). The phase of this signal does not match any known transiting planet in the system, and it rules out planets d and e.

We considered the possibility of an instrumental origin. However, no such offset is seen in ESPRESSO constant stars (Figueira et al. 2025), nor in the activity indicators at the relevant level. This suggests a physical origin for the signal.

An alternative explanation is that this signal may be related to a stellar magnetic cycle. Some activity indicators (e.g. FWHM, CCF contrast, S-index) show low-amplitude variations on timescales of 900–1700 days, which is roughly compatible with the period of the RV signal. However, modelling the RVs jointly with these indicators using a squared-exponential GP did not converge, and the signal in the RVs remains significant, independent of activity trends.

Its origin therefore remains uncertain: it may correspond to a long-period planet candidate, a stellar magnetic cycle, or, most plausibly, a combination of both. We tentatively modelled it as a Keplerian signal but stress that its planetary nature is not established. Continued long-term RV monitoring and activity characterisation are required. Future Gaia DR4/DR5 time-series astrometry may detect an acceleration or orbital motion and, combined with the RV dataset, help constrain the mass and nature of this companion.

6.1.5 Total mass of the planetary system

The HIP 41378 system hosts six planets up to ~1.5 AU. Given the observed resonances in the inner and outer sub-systems (see Sections 6.1 and 6.1.3), the architecture is likely complete up to the orbit of planet f. The most likely missed planet would be between the planets d and g, but would be less massive than a few Earth masses (see Fig. 8). The total planetary mass of the system up to 1.5 AU is of ~58 M_⊕. This value is substantially greater than that of the Solar System (~2 M_⊕), which might be explained by a more massive disc of a more massive host star. However, if we integrate the total planetary mass of the Solar System, up to the orbit of Neptune (~450 M_⊕), HIP 41378 is missing several hundreds of M_⊕, even including planet h. We thus speculate that HIP 41378 hosts some massive planets beyond 1.5 AU. Given the detection limits (see Fig. 8), Jupiter-mass planets could be compatible with our decade-long RV time series if their orbital periods exceed ~20 000 days (55 years). Detecting these potential outermost and massive planets would require substantially longer RV time baselines. These companions might also be probed by next-generation high-contrast imaging instruments such as the Habitable Worlds Observatory (HWO, National Academies of Sciences, Engineering, and Medicine 2021) and the Planetary Camera and Spectrograph (PCS, Kasper et al. 2021) on the Extremely Large Telescope (ELT), although likely at the limit of detectability. A Jupiter-like planet around HIP 41378 with a 55-year orbital period would have an angular separation of about 0.13″ and would be easily resolved. The expected reflected-light contrast, of the order of a few ×10⁻¹⁰, could be accessible with HWO considering the predicted performances (Feinberg et al. 2026). Alternatively, HIP 41378 and the Solar System might have formation mechanisms of different efficiency to form massive planets, especially as the former is a slightly metal-poor star.

Fig. 12 Mass–radius diagram for HIP 41378 transiting planets. Confirmed planets from the NASA Exoplanet Archive with 20% precision on mass and radius are represented as grey dots. Planets with densities under 0.2 g/cm−3 are highlighted. The colours show the importance of mass loss through αevap. Theoretical composition models for rocky planets and water worlds from Zeng et al. (2019) are represented as dashed lines, and rocky planets with H/He atmospheres from Lopez & Fortney (2014) are represented as dotted lines. The isodensity curves are show straight lines. Parts of this plot were prepared using mr-plotter (Castro-González et al. 2023).

6.2 The nature of HIP 41378 f

The five transiting planets of HIP 41378 span a wide range of orbital distances, yet they show a trend of decreasing density with radius and with distance from the star. The mass-radius diagram is represented in Figure 12. Planets b and c are consistent with sub-Neptunes with moderate H/He envelopes, while planets d and e are less dense.

For HIP 41378 f, our revised mass of 25 ± 5M_⊕ is twice as large as the 12 ± 3M_⊕ reported in Santerne et al. (2019). Its measured density is now estimated as ρ_f = ${0.16}_{-0.04}^{+0.03}$ g cm⁻³. The planet still belongs to the class of super-puff planets.

6.2.1 The planet’s low density

The low density of HIP 41378 f (ρ_f = ${0.16}_{-0.04}^{+0.03}$ g cm⁻³) is among the lowest measured for exoplanets and may at first appear difficult to reconcile with standard planet formation models. Indeed, conventional theories indicate that dense Earth-like cores in protoplanetary discs start accreting efficiently when the core mass is in excess of about 10 M_⊕ (Mizuno 1980; Pollack et al. 1996), leading to planetary densities much in excess of the measured value. However, when both accretion rates and opacities are reduced compared to traditional values, much smaller cores, even down to about ~1–5 M_⊕, can become supercritical and accrete a significant hydrogen–helium envelope within the lifetime of the discs (Ikoma et al. 2000).

The first example of planets with extremely low densities, also known as super-puffs, are the planets in the Kepler-51 system that all range between 0.03 and 0.05 g/cm³ (Masuda 2014; Masuda et al. 2024). As shown by Lee & Chiang (2016) these planets may form at large orbital distances in quiescent discs and migrate, potentially retaining their large radius and small densities. Their low density implies that the most stringent limitation for their survival may be resisting to the XUV stellar flux that they receive and the subsequent mass loss (Lopez & Fortney 2014).

Figure 12 shows a mass-radius diagram focusing on low-density planets with precise masses and radii. The colours correspond to the values of a mass-loss parameter α_evap defined as the ratio of the time integral of the UV flux received by the planet to its binding energy: $${\alpha }_{\text{evap\hspace{0.17em}}}\equiv E_{\mathrm{U}\mathrm{V}}/E_{\text{bind\hspace{0.17em}}},$$(2)

where $$E_{\mathrm{U}\mathrm{V}}=\frac{\pi R_{\mathrm{p}}^2}{a^2}{\int }_0^t{F}_{\mathrm{U}\mathrm{V}}dt,$$(3)

R_p being the planetary radius, a its orbital semi-major axis, and F_UV(t) the UV flux of the star, which was taken as equal to 29.7 (age/Gyr)^−1.23 erg/s/cm² following Ribas et al. (2005). We also accounted for a saturation of that relation for ages of < 100 Myr. The relation is formally valid only for solar-type stars, and we neglected the mass dependence (stars highlighted in Fig. 12 range from 0.69 to 1.25 M_⊙). The binding energy is approximated as $E_{\text{bind}}=\frac{G{M}_{\mathrm{p}}^2}{R_{\mathrm{p}}}$, where G is the gravitational constant and M_p the planetary mass.

As discussed first in a slightly different form by Lecavelier Des Etangs (2007) and by Lopez & Fortney (2014), α_evap indicates the magnitude of mass loss due to XUV radiation. We expect α_evap < 1 and this is mostly what we find for the planets shown in Fig. 12. For the Kepler-51 planets, we obtain α_evap values of 1.5, 0.4, and 0.08 for planets b, c, and d. Planet b appears most vulnerable to atmospheric escape, although uncertainties in its mass may bring α_evap below unity. Given the young age of the system (0.7 ± 0.5 Gyr; Masuda et al. 2024), continued mass loss may still alter their bulk properties.

In contrast, the planets orbiting HIP 41378 exhibit significantly smaller values of α_evap. The largest value in the system is 0.04 (planet b), while HIP 41378 f has α_evap = 0.001. The large orbital separation of planet f and the system’s age of ~2 Gyr imply that atmospheric escape has likely played only a minor role in shaping its present-day structure.

Among the low-density planets occupying a similar region of the mass–radius diagram as HIP 41378 f are planets b and e around the young star K2-309 A (V1298 Tau), a 20–30 Myr, 1.1 M_⊙ star. They are characterised by α_evap = 0.02 and 0.013, implying that they should suffer more mass loss over the lifetime of the star. Karalis et al. (2025) showed that the planets’ densities may be accounted for with a gas-to-core ratio that approaches unity. They probably migrated from beyond 1 au before reaching their present-day orbits at 0.17 and 0.27 au, respectively.

Overall, while the low density of HIP 41378 f is unusual, our analysis indicates that its envelope is energetically stable against XUV-driven escape over the system’s lifetime. Its current structure, therefore, does not require extreme ongoing mass loss to be explained.

As we detail in Section 6.2.2, the presence of planetary rings remains a possible interpretation of the large measured radius of HIP 41378 f. However, alternative formation pathways must be studied in detail. The possibility that incoming planetesimals and pebbles may be filtered by other planetary companions at larger orbital distances (e.g. Guillot et al. 2014) may help us to understand their formation. Detailed studies of these kinds of planets and a characterisation of their atmospheres would allow a much better understanding of the growth of planetary cores and the importance of neighbouring planetary companions during the process.

6.2.2 The ring hypothesis

A full understanding of the formation pathway of HIP 41378 f is essential to explain its low density. We therefore revisited the alternative ring scenario proposed by Akinsanmi et al. (2020) and assessed its compatibility with the updated system parameters. This hypothesis remains viable, although it requires a somewhat higher intrinsic density than initially assumed. Adopting R_true = 3.9 R_⊕ as the true planet’s radius yields a true density for planet f of 2.3 g cm⁻³ and an observed radius of R_f ≈ 2.4 R_true, which is close to the Roche limit, as expected for a Saturn-like tidal ring. Using Eq. (1) of Akinsanmi et al. (2020), we find the required ring particle density is 2.47 g cm⁻³, which is consistent with a rocky composition and long-term survival at the planet’s equilibrium temperature (294 K).

As shown by Akinsanmi et al. (2020), if a ring is present, it must lie close to the sky plane, implying an obliquity of about 90°, which is similar to that of Uranus. Such a tilt requires a dynamical mechanism, as proposed by Saillenfest et al. (2023) and Lu et al. (2025), but it remains to be verified that these mechanisms are still compatible with the updated system parameters in Table 2. These mechanisms rely on the existence of resonances between the nodal precession modes of the planet’s orbital plane (produced by the planet-planet gravitational interactions) and its spin-axis precession (produced by the stellar torque applied on the planet’s equatorial bulge). The precession modes of planet f’s orbital plane depend on the masses and semi-major axes of all planets in the system. A good estimate for the dominant modes is given by the Lagrange–Laplace theory. Figure 13 shows the probability density of these modes considering the parameters from Table 2 and their uncertainties. By comparing Fig. 13 with Fig. 3 of Saillenfest et al. (2023), we note that the updated parameters decrease the uncertainties of all frequencies, especially frequency s₄, without producing any major rearrangement in the system.

The inclusion of the distant candidate planet h only adds a very low-frequency term (s₇), corresponding to a slow precession of the inner system around the total angular momentum.

In the scenario of Saillenfest et al. (2023), a migrating moon around planet f gradually changes its spin-axis precession frequency until it matches one of the modes (s) shown in Fig. 13. At that point, a resonance is triggered, and the planet’s spin axis is tilted on a gigayear timescale up to an obliquity of about 90°, and it eventually leads to tidal disruption of the moon and ring formation. This mechanism has been shown to be generic (see e.g. Saillenfest et al. 2020, 2021, 2022). The top horizontal axis in Fig. 13 shows the minimum mass, m_min, that the moon must have had in order to tilt planet f through a resonance with a frequency of s along the bottom horizontal axis. The values of m_min were computed using Eq. (5) of Saillenfest et al. (2023), the assumed radius R_true = 3.9 R_⊕, and the spin rate and oblateness of Uranus. If R_true is taken as R_f instead, then the values obtained for m_min are multiplied by roughly a factor of three. Figure 13 shows that in order to be tilted through a resonance with mode s₆, planet f must have had a moon with mass ratio of m/M_f equal to a few 10⁻⁴. This is consistent with what is expected for regular moons around giant planets (see Canup & Ward 2006; Saillenfest et al. 2023). For a 25 M_⊕ planet, this would correspond to a moon of ~0.005–0.01 M_⊕, which is below the maximum satellite mass of ~0.15 M_⊕ found to be stable over long timescales for HIP 41378 by Harada et al. (2023). Although such a moon would be undetectable, it would be sufficient to drive the mechanism.

The mechanism of Lu et al. (2025) is similar, except that it relies on the early disc-driven planetary migration in order to trigger the resonance, and it does not explain the existence of the ring. If we consider that planet f did not have a moon during the planetary migration phase, its spin-axis precession frequency would have been very slow (see Eqs. (6) and (7) of Lu et al. 2025); yet, Fig. 13 shows that frequencies s₃ and s₇ extend to very low absolute values, still allowing for resonances with planet f’s slow spin-axis precession. If instead planet f had a moon during planetary migration, Fig. 13 shows that a resonance with frequency s₆ would also be reachable. Therefore, the mechanism proposed by Lu et al. (2025) also remains a viable mechanism when considering the updated parameters for the HIP 41378 system given in Table 2.

Fig. 13 Probability density of the inclination proper modes of the HIP 41378 system with and without the candidate planet h. Histograms are built from 106 realisations of the Lagrange-Laplace system with the mass and semi-major axis uncertainties listed in Table 2. For planets g and h, the minimum mass is taken as the true mass. The histograms for frequencies s3 and s7 peak above the top border of the graph. Frequency s5 is identically equal to zero from the conservation of angular momentum. The upper axis shows the minimum moon mass needed for planet f to be fully tilted through a resonance with a given frequency value s.

7 Conclusions

Following the first set of transits detected by K2 in 2015, we monitored the HIP 41378 system from 2016 to 2025 in RV. Over 1000 spectra were obtained in one decade over 621 individual nights, mainly with the HARPS, HARPS-N, HIRES, and ESPRESSO spectrographs. The analysis of these data, combined with the K2 photometry and constraints from TESS and CHEOPS, allowed us to determine the orbital periods and masses of six planets in the system (ordered by increasing distance from the host star): b, c, g, d, e, f. The RV amplitudes range from 0.5 m s⁻¹ to 1.8 m s⁻¹ with typical precision levels between 6% and 50%.

Compared to previous studies on this system, the present work provides several improvements. These are listed below:

• We find strong evidence that planets d and e have orbital periods of P_d ~ 278 days and P_e ~ 386–390 days;

• We refined the masses for all planets, including planets d and e, and doubled the mass estimate for the super-puff planet f;

• We detected a long-period planet candidate, HIP 41378 h, of m sin i ~ 43 M_⊕ and a period of six to eight years;

• We observe a clear trend of decreasing density with semi-major axis for the transiting planets of the system;

• We investigated the system architecture, which appears to consist of two dynamically distinct sub-systems. Additional undetected planets could exist and potentially connect the two sub-systems in a resonant chain.

However, there are still open questions concerning this system that are not fully addressed here and require further observations:

• A third (second) photometric transit detection of planet d (e, respectively) is needed to fully secure their orbital period and derive an accurate ephemeris. The CHEOPS space telescope is perfectly suited for that, although the observation time span should last long enough to account for the large TTVs (Sulis et al. 2024);

• The nature of the low-density planet f is still uncertain. Its properties are consistent with either an inflated temperate planet or the presence of rings. Observing the transit in the near-infrared with, for example, the JWST or ARIEL, should allow us to characterise the planet’s atmosphere and radius more accurately, thereby testing the ‘inflating planet’ scenario and confirm or rule out the ring hypothesis (Alam et al. 2022);

• The candidate planet h needs to be secured via long-term spectroscopic observations or a detailed characterisation of the stellar variability across one decade;

• The architecture of the system appears almost complete up to 1.5 AU and down to a few Earth masses, although one or two planets might still be undetected between planets g and d. Further observations could deeply explore this rich system. Comparison with the Solar System indicates that more massive planets might still be discovered beyond ~50 AU, which would require next-generation high-contrast imaging instruments to be detected.

The diversity of planets’ properties in this bright system, from warm sub-Neptunes to temperate inflated planets, makes HIP 41378 a unique laboratory for comparative planetology within the system and with the Solar System. However, this requires maintaining a precise ephemeris in the future to schedule follow-up observations. More transit observations of the outermost planets (Leonardi et al., in prep.) are thus needed in this regard.

Data availability

Radial velocity data are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/710/A275

Acknowledgements

We acknowledge the contribution from the PFS team with their early observations of the system. We made use of the SIMBAD database, CDS, Strasbourg Astronomical Observatory, France and of the VizieR catalogue access tool, CDS, Strasbourg Astronomical Observatory, France (DOI: 10.26093/cds/vizier). The project leading to this publication has received funding from the Excellence Initiative of Aix-Marseille University – A*Midex, a French “Investissements d’Avenir programme” AMX-19-IET-013. This work was supported by the “Programme National de Planétologie” (PNP) of CNRS/INSU co-funded by CNES. This work was supported by the ‘Action Thématique Physique Stellaire” of CNRS/INSU PN ASTRO co-funded by CEA and CNES. This work is based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. 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 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. The work of HPO 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. NCS is co-funded 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. This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds by these grants: UIDB/04434/2020 DOI: 10.54499/UIDB/04434/2020, UIDP/04434/2020 DOI: 10.54499/UIDP/04434/2020. We acknowledge financial contribution from the INAF Large Grant 2023 “EXODEMO”. This research was funded in part by the UKRI (Grants ST/X001121/1, EP/X027562/1). Some of the data presented herein were obtained at Keck Observatory, which is a private 501 (c) 3 non-profit organisation operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. We recognise and acknowledge the cultural role and reverence that the summit of Maunakea has within the indigenous Hawaiian community. We are deeply grateful to have the opportunity to conduct observations from this mountain. 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. S.C.C.B.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. P.F. acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033. P.F. is also funded by the European Union (ERC, THIRSTEE, 101164189). BA acknowledges the support of the Swiss National Science Foundation under grant number PCEFP2 194576. C.A.W. would like to acknowledge support from the UK Science and Technology Facilities Council (STFC, grant number ST/X00094X/1). 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). AL acknowledges support of the Swiss National Science Foundation under grant number TMSGI2_211697. ACC acknowledges support from STFC consolidated grant number ST/V000861/1, UKRI/ERC Synergy Grant EP/Z000181/1 (REVEAL), and UKSA grant number ST/X002217/1. O.D.S.D. acknowledges support from e-CHEOPS (PEA No. 4000142255.

References

Appendix A Stellar activity

Figure G.1 presents the time series and GLS periodograms of several stellar activity indicators: the CCF FWHM, BIS, CCF contrast, and the Ca II H&K S-index. These diagnostics are used to search for signatures of stellar rotation and long-term magnetic variability, and to assess their potential impact on the RV measurements.

Prior to computing the periodograms, we removed linear trends and inter-instrument offsets from each indicator. In the periodograms, we indicate the orbital periods of the transiting planets, as well as the stellar rotation periods inferred from $\log \left(R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}\right)$ and photometry. The highest peaks and their one-year aliases are also marked. We additionally computed Spearman rank correlation coefficients between each activity indicator and the RV time series, both globally and per instrument, and found little to no correlation in all cases.

Below, we analyse each indicator in detail.

Appendix A.1 FWHM

The FWHM periodogram displays a dominant peak at ~1735 days, with a secondary peak at ~463 days corresponding to a one-year alias. The long-period signal can be well modelled by a sinusoid with a period of 1714 days. After subtracting this signal, no significant power remains at shorter periods that could be associated with stellar rotation.

This long-term variability is primarily driven by a quadratic trend visible in the ESPRESSO data. To illustrate this, we computed GLS periodograms separately for ESPRESSO and for HARPS/HARPS-N, as shown in Fig. A.1. Once the quadratic trend is removed from the ESPRESSO FWHM time series, no significant periodic signal is detected.

Therefore, the FWHM does not provide clear evidence for rotationally modulated stellar activity in this system. Moreover, the long-period signal detected in the FWHM (1714 days) is significantly different from the long-term signal identified in the RV time series, suggesting that it is unlikely to be responsible for the RV variability, or at least not in its entirety.

Appendix A.2 BIS

After removing a linear trend, the BIS periodogram does not exhibit any statistically significant periodic signal. It is the only activity indicator that does not display long-period signals. Low-amplitude peaks are present near 8 days, close to the expected stellar rotation period, and at one-day aliases, but their power remains well below significance thresholds.

To further investigate the possible presence of rotational modulation, we modelled the BIS time series using a GP with a quasi-periodic kernel (Eq. (1)), implemented with tinyGP. We adopted uniform priors on the rotation period (2–20 days), amplitude (0–5 m s⁻¹), exponential decay timescale (1–500 days), and Γ (0.1–10 days), and included instrument-dependent offsets and jitter terms.

Despite exploring several prior configurations on the rotation period, the GP fit failed to converge to a stable solution. This lack of convergence is consistent with the absence of a coherent and sufficiently strong quasi-periodic signal in the BIS time series, and suggests that rotationally induced activity is not significantly detected in this indicator.

Fig. A.1 Frequency analysis of the FWHM. Upper panel: Periodograms of the FWHM timeserie with only HARPS and HARPS-N data in pink, only ESPRESSO data in black and all instruments together in grey. The two rotation periods derived from $\log \left(R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}\right)$ and photometry are represented as, respectively, solid and dotted green lines. Vertical dotted black lines represent the period of the known planets in the system. Lower panel: Periodogram of ESPRESSO data after removing a quadratic-trend. The horizontal dashed lines in the two panels represents the FAP levels at 0.1, 1 and 10%.

Appendix A.3 Contrast

The CCF contrast shows the strongest correlation with the RVs among the considered indicators, particularly in the ESPRESSO data. Its GLS periodogram exhibits a dominant peak at ~1013 days, which may be indicative of a low-amplitude stellar magnetic cycle. This long-period signal is well reproduced by a sinusoidal model at the corresponding period.

After subtracting this long-term variability, several low-power peaks remain at shorter periods, but none reach statistical significance or coincide robustly with the expected rotation period. We attempted to model the contrast time series with a quasi-periodic GP, both before and after removal of the long-period sinusoidal component, using the same priors as for the BIS.

In all cases, the GP modelling suffered from poor convergence, preventing a reliable determination of the stellar rotation period. This again suggests that rotational modulation in the contrast is weak or incoherent over the observational baseline.

Appendix A.4 S-index

The Ca II H&K S-index periodogram shows a significant peak near 948 days, close to the long-period signal detected in the CCF contrast, and likely associated with the same magnetic activity cycle. A sinusoidal fit yields a period of 969 ± 27 days for this long-term variability.

Once this signal is removed, no significant power remains at shorter periods that could be attributed to stellar rotation. As for the other indicators, quasi-periodic GP modelling did not converge, further indicating the absence of a detectable rotationally modulated signal in the S-index.

Overall, none of the activity indicators show a robust periodic signal at the expected stellar rotation period, nor do they exhibit strong correlations with the RV time series. The long-term variability detected in the FWHM, contrast, and S-index may be attributed to a low-amplitude magnetic cycle, however, its characteristic period does not fully coincide with the long-period signal observed in the RVs, which is significantly longer.

Appendix B Residuals analysis

Fig. B.1 Periodogram of the residuals of the RV analysis without GP (in black) and with GP (in red). The horizontal dashed lines represents the FAP levels at 0.1, 1 and 10%. The shaded area correspond to the stellar rotation period.

Appendix C Nested sampling analysis of planets d and e

We performed a nested sampling analysis on the HARPS, HARPS-N and ESPRESSO RV data to evaluate the orbital periods of planets d and e and to compare different models in a Bayesian framework. Nested sampling is a Monte Carlo algorithm that simultaneously explores the posterior distribution and computes the Bayesian evidence, 𝒵, for model comparison. We used the jaxns library (Albert 2020), which implements a fully JAX-based nested sampling framework optimised for high-dimensional and correlated parameter spaces.

We considered RV-only models, including the known planets b, c, g and f and the candidate planet h. We first included only planet d in order to test its orbital period. A uniform prior between 200 and 400 days was adopted to encompass the two possible periods (278 and 371 days). All other priors were identical to those used in the MCMC analysis presented in Section 4.2.1. We ran jaxns with 900 live points and resampled the weighted posterior to obtain 5000 effective samples. Figure C.1 shows the resulting posterior distribution, which exhibits a clear preference for an orbital period of 278 days for HIP 41378 d.

We then included planet e, adopting a uniform prior between 100 and 500 days, since longer periods are unlikely given the observed transit duration. An orbital period of P_e ≈ 393 days is preferred. An alias of this period is also present in the posterior distribution, near 190 days.

In order to constrain the orbital period of planet e, which is the most challenging to determine, we performed a nested sampling analysis combining transit data (the two K2 light curves) and RV data from the four spectrographs. We included the fit of all the other planets of the system with gaussian prior on their orbital periods. We tested three scenarios: (1) planet e located between planets d and f (uniform prior between 300 and 500 days), (2) planet e located exterior to planet f (uniform prior between 550 and 900 days), and (3) planet e located interior to planet d (uniform prior between 100 and 200 days). For each scenario, we ran jaxns with 900 live points.

Fig. C.1 Posterior distributions of the orbital periods of planets d (left) and e (right) obtained from nested sampling with uniform priors of the ESPRESSO data. The two main possible orbital periods for planet d are highlighted in pink (278 days) and purple (371 days). For planet e, the signal at ≈ 393 days and its alias at ≈ 200 days are highlighted in blue dashed lines.

The Bayesian evidence allows us to quantify the relative support for models with different assumptions on the orbital period of planet e. Table C.1 summarises the results. The quantity Δ log Z is computed relative to the model with the strongest Bayesian evidence; differences larger than 5 are generally interpreted as strong evidence against the disfavored model. The configuration in which planet e lies between planets d and f is strongly favoured, with a difference in Bayesian evidence of Δ log Z ≈ 15 compared to the alternative configuration in which planet e is located exterior to planet f.

Appendix D Detection limit

We computed the RV detection limit using the method introduced in Standing et al. (2022) and implemented with the diffusive nested sampler kima. The analysis was performed on the RV residuals obtained after subtracting the median best-fit signals of the seven known planets and the quasi-periodic GP.

We modelled the residuals assuming a single additional Keplerian signal using the known object mode of kima, which samples the posterior distribution of all signals compatible with the data but not yet detected. The orbital period was assigned a log-uniform prior between 3 and 10⁵ days, while the RV semi-amplitude was drawn from a uniform prior between 0.07 and 100 m s⁻¹. The eccentricity followed a truncated normal distribution centred on zero with a dispersion of 0.083, and uniform priors were adopted for the argument of periastron and orbital phase.

To ensure a well-sampled posterior, we performed three independent runs with different random seeds, each producing more than 20 000 effective posterior samples. The resulting (P_p, K_p) posterior densities were combined and displayed as a greyscale hexbin map.

The detection limit was computed by dividing the period range into logarithmically spaced bins and, within each bin, evaluating the maximum 99th percentile of the K_p distribution. To mitigate small-number statistics, we required at least two posterior samples to lie above the chosen percentile in each bin. The detection limit was not computed beyond the first period bin, where the posterior became significantly affected by the upper boundary of the K_p prior.

Appendix E Tests of periodicity

Fig. E.1 Apodized sine periodogram of the RV data, including offsets, transiting planets, and order 11 polynomial in the base model. Left: apodized periodogram for different apodization timescales, that is Eq. (E.1) for τ = 10Tobs (blue), Tobs/3 (orange), Tobs/9 (green) and Tobs/27 (red), denoting by Tobs the observation timespan of the data. Horizontal dashed lines represent the value of the ASP at the period corresponding to the global maximum (here attained for τ = Tobs/3 Middle: zoom around the maximum peak. The right panel represents a statistical test: assuming the true value of the timescale given in the x axis, the circles represent the expected value of the peaks corresponding to other timescales, the error bar represent the standard deviation. For instance, assuming τ = 10Tobs is the true timescale, we would expect an apodized sine with τ = Tobs/3 to have a χ2 difference of ≈ 12 ±5, which is 5.4 σ away from the measured value (39, orange dashe line). The error bars are horizontally shifted to improve the readability of the figure, but they correspond to the same assumed timescale.

To assess the origin of the signals, we apply an apodized sine periodogram (ASP), which tests whether the signals have a consistent amplitude, phase, and frequency over time. It consists in computing the difference between the χ² of a linear base model H, and the χ² of a model K(ω, t₀, τ) defined as the linear model of H plus an apodized sinusoid $e^{-\frac{{(t-t_0)}^2}{2{\tau }^2}}(A\cos \omega t+B\sin \omega t)$ (Hara et al. 2022). The ASP is defined as the χ² difference between the two models, maximized over the center of the time window. $$z(\omega ,t_0,\tau )=\underset{t_0}{max}({\chi }_H^2-{\chi }_{K(\omega ,t_0,\tau )}^2).$$(E.1)

We use the same base model as in Section 4.3.2, and use four values of the signal time scale, τ = 10T_obs, T_obs/3, T_obs/9, and T_obs/27. The ASP is shown in Fig. E.1 (left), and exhibits a strong 8.2 signal with a time scale of τ = T_obs/3, the centre panel is simply a zoom in.

We want to test whether the signal is statistically compatible with a constant one (i.e. with τ = 10T_obs). Denoting by (t_{τ, ω}) the value of t₀ maximising the value of the periodogram (E.1) for a given ω and τ, we compute the distribution of: $$D_z=z(\omega ,t_{(\tau ,\omega )},\tau )-z(\omega ,t_{({\tau }^{\mathrm{\prime }},\omega )},{\tau }^{\mathrm{\prime }})$$(E.2)

with the hypothesis that the model K(ω,t_(τ,ω), τ, A^⋆, B^⋆) is correct, where the fitted cosine and sine amplitudes A^⋆, B^⋆ are obtained by fitting model K to the data. D_z can easily be expressed as a generalised χ² distribution, with mean and variance is given by an analytical expression (Hara et al. 2022). We find that a stable, periodic signal is excluded (Fig. E.1, right panel).

As a final remark, we note that the ASP is dominated by short-term signals, possibly due to systematics. Nonetheless, the two peaks at 278 and 146 days exhibited in Section 4.3.2 also appear as stable signals in Fig. E.1, left panel.

Appendix F Additional resonant maps

Fig. F.1 Position of the posterior of HIP 41378 with respect to MMR, with the addition of a planet x between g and d while minimising distance to 2-body MMRs (case i). See Figure 11 for more details.

Fig. F.2 Position of the posterior of HIP 41378 with respect to mean motion resonance, with the addition of a planet x between g and d while minimising distance to 3-body MMRs (case ii). See Figure 11 for more details.

Appendix G Activity Indicators periodograms

Fig. G.1 Time series and GLS periodograms of RVs and stellar activity indicators for HIP 41378 from HARPS, HARPS-N, and ESPRESSO. From top to bottom: RVs, FWHM, BIS, CCF contrast, and S-index. Periods of planets b, c, g, d, e and f are shown as grey dashed lines. Stellar rotation periods from photometry and $\log \left(R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}\right)$ are indicated as green dotted and solid lines, respectively. Maximum power periods are shown as dark pink lines, with aliases in light pink. False-alarm probabilities at 1% and 10% are indicated with horizontal lines.

Appendix H RV analysis priors and posteriors distributions

Fig. H.1 Prior and posterior distributions of the orbital periods (left panels) and RV semi-amplitudes (right panels) for the seven planets. The priors are shown as grey distributions. The dashed orange lines indicate the posterior distributions obtained from the RV analysis using the HARPS, HARPS-N, and HIRES instruments. The dashed blue lines correspond to the posterior distributions derived from the ESPRESSO data alone. The pink distributions represent the posterior results from the combined analysis of all four instruments.

Appendix I Global analysis results

Appendix J K2 transit light curves

Fig. J.1 Light curves of HIP 41378. Upper panel: K2 C5 campaign (long cadence). Bottom panel: K2 C18 campaign (short cadence). The pink curve represents the median of the best-fit model from photometry and RV analysis. The names of the planets are highlighted with their letters.

All Tables

All Figures

	Fig. 1 Confirmed planetary system with more than five planets detected and at least two with determination of both mass and radius. Colours represent planetary masses and circle sizes represent the radius. Planets represented as a cross have no radius measurement, and unfilled circles mark cases where no mass measurement was available. The planetary parameters are from Nasa Exoplanet Archive.
In the text

	Fig. 2 GLS periodograms of HARPS, HARPS-N, HIRES, and ESPRESSO data. Planet periods b, c, and f are highlighted in pink, as well as the candidate period of planet d (278). Subsequent panels show periodograms after subtracting b, c, and f, and then including h and g. False-alarm probability levels at 1% and 10% are indicated as horizontal dotted black lines.
In the text

	Fig. 3 Tests for the orbital period of HIP 41378 d. Left panel: posterior density of the period of planet d, considering a uniform prior between 250 and 400 days. The vertical dashed lines indicate the possible periods from transit observations. Right panel: resulting RV semi-amplitude and 99% HDI for the three periods tested: 278, 371, and 1113 days.
In the text

Fig. 4 Comparison of different RV models, including subsets of the HIP 41378 planets, ranked by their expected log point-wise predictive density (elpd LOO). Each model label indicates the included planets and, for planet d, the tested orbital period in days. Higher elpd LOO values correspond to better out-of-sample predictive performance. The horizontal line segments represent the standard error of the elpd. The lighter line segment with triangles indicates the elpd difference of one model relative to the best model.

In the text

	Fig. 5 Upper panel: ℓ1 periodogram of the RV data, including offsets in the base model. Lower panel: ℓ1 periodogram of the RV data, including offsets, transiting planets, and order 11 polynomial in the base model. Planets b, c, g, and f are highlighted by their letters.
In the text

	Fig. 6 Comparison of the semi-amplitude of the seven planets in the HIP 41378 system for a model with GPs (in blue) and without GPs (in pink). The 68% HDI is represented as a box and the 99% HDI as a horizontal bar.
In the text

	Fig. 7 RV timeseries of HIP 41378 with HARPS (pink), HARPS-N (green), HIRES (purple) and ESPRESSO (blue), along with best fit model (median of the posterior) in pink and the GP component in black. The upper panel is a zoom on ESPRESSO data. The lower panel represents the residuals of the fit (data – model).
In the text

	Fig. 8 Greyscale hex-bin plot representing the posterior samples obtained from kima on the full RV dataset residuals with the number of planets fixed to one. The blue line represents the 99th percentile detection limit. The RV semi-amplitudes of Jupiter, Saturn, Neptune, and Earth are represented by dashed pink lines for comparison. The top axis shows the corresponding angular separation.
In the text

	Fig. 9 HARPS-N, HARPS, HIRES, and ESPRESSO RV phase-folded on the results of the joint RV+transit fit for the seven planets in the system. The pink line represents the median of the best-fit model RV curve with its 68% HDI in the shaded area.
In the text

	Fig. 10 Possible orbital periods for HIP 41378 e. The posterior distribution of the orbital period of planet e, obtained from the RV-only analysis with a broad uniform prior, is shown as a blue histogram. Period ranges excluded by K2 (pink), TESS (purple), and CHEOPS (orange) photometric observations are indicated as pink shaded areas.
In the text

Fig. 11 Subsequent triplet of HIP 41378 system with respect to MMRs. The solid purple lines show first-order two-planet MMRs, the dashed purple lines show second-order two-body MMRs, the solid red line shows the 2:1 MMR between planet n and n + 2, and the grey lines show zeroth-order three-body MMRs of the form k/Pn − (k + q)/Pn + 1 + q/Pn+2 = 0 (e.g. Cerioni et al. 2022). The top panel shows the RV posterior of the known planets. The bottom panel shows the example of case (iii) – solution 1, for which we added two planets, x and y, in order to reconnect planets c–f in a resonant chain.

In the text

Fig. 12 Mass–radius diagram for HIP 41378 transiting planets. Confirmed planets from the NASA Exoplanet Archive with 20% precision on mass and radius are represented as grey dots. Planets with densities under 0.2 g/cm−3 are highlighted. The colours show the importance of mass loss through αevap. Theoretical composition models for rocky planets and water worlds from Zeng et al. (2019) are represented as dashed lines, and rocky planets with H/He atmospheres from Lopez & Fortney (2014) are represented as dotted lines. The isodensity curves are show straight lines. Parts of this plot were prepared using mr-plotter (Castro-González et al. 2023).

In the text

Fig. 13 Probability density of the inclination proper modes of the HIP 41378 system with and without the candidate planet h. Histograms are built from 106 realisations of the Lagrange-Laplace system with the mass and semi-major axis uncertainties listed in Table 2. For planets g and h, the minimum mass is taken as the true mass. The histograms for frequencies s3 and s7 peak above the top border of the graph. Frequency s5 is identically equal to zero from the conservation of angular momentum. The upper axis shows the minimum moon mass needed for planet f to be fully tilted through a resonance with a given frequency value s.

In the text

Fig. A.1 Frequency analysis of the FWHM. Upper panel: Periodograms of the FWHM timeserie with only HARPS and HARPS-N data in pink, only ESPRESSO data in black and all instruments together in grey. The two rotation periods derived from $\log \left(R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}\right)$ and photometry are represented as, respectively, solid and dotted green lines. Vertical dotted black lines represent the period of the known planets in the system. Lower panel: Periodogram of ESPRESSO data after removing a quadratic-trend. The horizontal dashed lines in the two panels represents the FAP levels at 0.1, 1 and 10%.

In the text

	Fig. B.1 Periodogram of the residuals of the RV analysis without GP (in black) and with GP (in red). The horizontal dashed lines represents the FAP levels at 0.1, 1 and 10%. The shaded area correspond to the stellar rotation period.
In the text

	Fig. C.1 Posterior distributions of the orbital periods of planets d (left) and e (right) obtained from nested sampling with uniform priors of the ESPRESSO data. The two main possible orbital periods for planet d are highlighted in pink (278 days) and purple (371 days). For planet e, the signal at ≈ 393 days and its alias at ≈ 200 days are highlighted in blue dashed lines.
In the text

Fig. E.1 Apodized sine periodogram of the RV data, including offsets, transiting planets, and order 11 polynomial in the base model. Left: apodized periodogram for different apodization timescales, that is Eq. (E.1) for τ = 10Tobs (blue), Tobs/3 (orange), Tobs/9 (green) and Tobs/27 (red), denoting by Tobs the observation timespan of the data. Horizontal dashed lines represent the value of the ASP at the period corresponding to the global maximum (here attained for τ = Tobs/3 Middle: zoom around the maximum peak. The right panel represents a statistical test: assuming the true value of the timescale given in the x axis, the circles represent the expected value of the peaks corresponding to other timescales, the error bar represent the standard deviation. For instance, assuming τ = 10Tobs is the true timescale, we would expect an apodized sine with τ = Tobs/3 to have a χ2 difference of ≈ 12 ±5, which is 5.4 σ away from the measured value (39, orange dashe line). The error bars are horizontally shifted to improve the readability of the figure, but they correspond to the same assumed timescale.

In the text

	Fig. F.1 Position of the posterior of HIP 41378 with respect to MMR, with the addition of a planet x between g and d while minimising distance to 2-body MMRs (case i). See Figure 11 for more details.
In the text

	Fig. F.2 Position of the posterior of HIP 41378 with respect to mean motion resonance, with the addition of a planet x between g and d while minimising distance to 3-body MMRs (case ii). See Figure 11 for more details.
In the text

Fig. G.1 Time series and GLS periodograms of RVs and stellar activity indicators for HIP 41378 from HARPS, HARPS-N, and ESPRESSO. From top to bottom: RVs, FWHM, BIS, CCF contrast, and S-index. Periods of planets b, c, g, d, e and f are shown as grey dashed lines. Stellar rotation periods from photometry and $\log \left(R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}\right)$ are indicated as green dotted and solid lines, respectively. Maximum power periods are shown as dark pink lines, with aliases in light pink. False-alarm probabilities at 1% and 10% are indicated with horizontal lines.

In the text

Fig. H.1 Prior and posterior distributions of the orbital periods (left panels) and RV semi-amplitudes (right panels) for the seven planets. The priors are shown as grey distributions. The dashed orange lines indicate the posterior distributions obtained from the RV analysis using the HARPS, HARPS-N, and HIRES instruments. The dashed blue lines correspond to the posterior distributions derived from the ESPRESSO data alone. The pink distributions represent the posterior results from the combined analysis of all four instruments.

In the text

	Fig. J.1 Light curves of HIP 41378. Upper panel: K2 C5 campaign (long cadence). Bottom panel: K2 C18 campaign (short cadence). The pink curve represents the median of the best-fit model from photometry and RV analysis. The names of the planets are highlighted with their letters.
In the text