© The Authors 2026

1 Introduction

M dwarf stars offer a favourable environment for exoplanet studies. Comprising approximately 70% of the stars in our galaxy (e.g. Reid & Gizis 1997; Bochanski et al. 2010) and representing the majority of nearby stellar systems (e.g. Henry et al. 2006), these low-mass low-luminosity stars amplify detectable signals through both radial velocity and transit methods. Their smaller radius makes transit detection more feasible, as planets occult a larger fraction of the stellar disk, while their lower mass enhances the gravitational reflex motion induced by orbiting planets. Additionally, M dwarfs have a higher occurrence rate of small planets (e.g. Bonfils et al. 2013; Sabotta et al. 2021), and this combined with the enhanced detectability of atmospheric features means they are ideal targets for studying their composition and evolutionary pathways (e.g. Benneke et al. 2019; Kreidberg et al. 2019).

Sub-Neptunes stand out as the most abundant exoplanets in our galaxy (Bonfils et al. 2013; Dressing & Charbonneau 2015; Fulton et al. 2017), yet their nature and formation mechanisms remain poorly understood. With radii between those of Earth and Neptune, they occupy a focus of parameter space where composition can vary considerably. This diversity is illustrated by the radius valley (Fulton et al. 2017; Cloutier & Menou 2020), a decrease in the number of observed planets between 1.5 and 2.0 R_⊕, suggesting distinct formation and evolution mechanisms that shape these planets’ structure and composition. Various scenarios have been proposed to explain this gap, including atmospheric loss through photoevaporation (Owen & Wu 2017), core-powered mass loss (Gupta & Schlichting 2020), intrinsically gas-poor formation (Lee & Chiang 2016; Venturini et al. 2020), or compositions dominated by volatiles (Zeng et al. 2019; Nixon & Madhusudhan 2021; Aguichine et al. 2021).

GJ 3090 b, first identified by Almenara et al. (2022) (hereafter A22) exemplifies the importance of the research on sub-Neptunes. With a reported radius of 2.171 ± 0.068 R_⊕ (Ahrer et al. 2025, hereafter A25), GJ 3090 b lies near the transition between super-Earths and sub-Neptunes, providing a unique opportunity to study atmospheric retention and interior composition (Rogers 2015). Sub-Neptunes in this regime exhibit a wide diversity in bulk densities, ranging from rocky cores with thin gaseous envelopes to volatile-rich water worlds (Zeng et al. 2019). Radial velocity (RV) measurements are crucial for constraining the masses of such planets, which in turn enable bulk density calculations and interior structure modelling (Borsato et al. 2019; Hoyer et al. 2021; Leleu et al. 2024). GJ 3090 b’s initial mass estimate of 3.34 ± 0.72 M_⊕ and equilibrium temperature of 693 ± 18 K, derived from High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003) measurements, placed it among the lower-density sub-Neptunes. However, precise mass constraints are critical for distinguishing between competing formation and evolutionary scenarios (Otegi et al. 2020).

Advances in instrumentation, including the new Near InfraRed Planet Searcher (NIRPS; Bouchy et al. 2025) at La Silla Observatory, have enhanced the precision of RV measurements in the near-infrared. NIRPS and HARPS have the capability to observe simultaneously, which can be useful for identifying stellar activity in RV data using its chromatic nature. In addition, new RV analysis techniques, including multidimensional Gaussian processes (GPs; Rajpaul et al. 2015), are now able to accurately separate planetary and stellar contributions in radial velocities, which enables us to precisely characterize known planets and find additional companions. As part of the NIRPS Guaranteed Time Observations (GTO) programme, GJ 3090 was monitored with the aim of improving the mass determination of its transiting planet GJ 3090 b. By combining these new NIRPS observations with an expanded HARPS dataset, which resulted in around 150 additional radial velocity measurements, we present a comprehensive reanalysis of the GJ 3090 system.

Section 2 details the observational datasets used in this work, including photometry, RVs, and differential temperatures (ΔT). In Section 3, we characterize the stellar properties of GJ 3090. Then we present the analysis of the photometry datasets in Section 4, from which we update the radius of GJ 3090 b and identify evidence for a possible long-term variability signal that may be associated with a magnetic cycle. Our RV modelling framework is introduced in Section 5, where we describe the treatment of instrumental offsets and jitter, the correction for this long-term signal, the modelling of Keplerian planetary signals, and the implementation of a multidimensional GP to mitigate stellar activity. The results are presented in Section 6, including updated mass and density constraints for GJ 3090 b and a search for additional planetary signatures in the radial velocities and photometry. Finally, we discuss the architecture of the GJ 3090 system in a broader exoplanetary context. We conclude with a summary of our findings in Section 7.

2 Observations

2.1 TESS photometry

GJ 3090 (TIC 262530407, TOI-177) was observed with TESS in four different sectors throughout its primary and extended missions: Sector 2 (2018-08-23 to 2018-09-20), Sector 3 (2018-09-20 to 2018-10-17), Sector 29 (2020-08-26 to 2020-09-21), and Sector 69 (2023-08-25 to 2023-09-20). For our study, we used both the Simple Aperture Photometry (SAP; Twicken et al. 2010; Morris et al. 2020) and the Presearch Data Conditioning SAP (PDCSAP; Stumpe et al. 2012, 2014; Smith et al. 2012) data products provided by the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016) at NASA Ames Research Center, available through the Mikulski Archive for Space Telescopes¹. On the one hand, as in A22, we chose the SAP photometry over the PDCSAP data to study the stellar rotation (Section 5), as the latter removes long-term trends and, in this case, was found to overcorrect this signal. For our stellar activity analysis, we removed in-transit photometric data points from the light curve to isolate the out-of-transit flux variations. This step was essential for using the light curve as a stellar activity indicator in our RV analysis. To reduce computation time while retaining adequate sampling to model the 18 d stellar rotation period found by A22, we binned the TESS light curves into 8-hour intervals.

Since TESS collected an additional sector (Sector 69) after the publication of A22, we updated the transit analysis accordingly, using the PDCSAP flux, which accounts for dilution caused by contaminating sources within the TESS aperture. This analysis is described in Sect. 2.1.

2.2 ASAS-SN photometry

To investigate the long-term variability of GJ 3090 and assess the possible presence of a magnetic cycle, we used archival photometric data from the All-Sky Automated Survey for Supernovae (ASAS-SN). ASAS-SN is a global network of 14 cm telescopes designed to monitor the entire visible sky on a nightly basis, originally developed to detect supernovae and other transient phenomena (Kochanek et al. 2017). The network consists of 24 telescopes distributed across six observatories worldwide, ensuring consistent coverage and minimizing weather-related data gaps. Sites are located at: Haleakalā Observatory, HI, USA; two at Cerro Tololo Observatory, Coquimbo, Chile; South African Astronomical Observatory, Sutherland, South Africa; McDonald Observatory, TX, USA; Ali Observatory, Tibet, China.

The ASAS-SN light curve was retrieved from the Sky Patrol V1 database², with a manual correction for the stellar proper motion applied by supplying the proper motion-corrected coordinates from SIMBAD in the query, and all measurements obtained when the Moon was within 90^° of the target were excluded. After a 3σ outlier removal, the dataset for GJ 3090 contains 2772 photometric measurements spanning May 2014 to September 2025, for a total time span of 4144 days (11.3 years). The data come from four different cameras labelled bf, bj, bF, and bn. The bf camera is distinct from the others in that it provides the only data obtained before 2018, the year it became out of service, and it uses a slightly redder effective bandpass (λ_eff=551 nm, V filter) compared to the other three cameras (λ_eff=480 nm, g filter). The sampling is irregular, with a mean time interval between measurements of about four days. In this work, the ASAS-SN data are used to identify the presence of a long-term signal on GJ 3090, which can motivate its modelling in other time series. This is done as an independent analysis described in Section 4.2.

2.3 NIRPS radial velocities

We acquired spectra of GJ 3090 on 98 nights with the NIRPS (Bouchy et al. 2017; Wildi et al. 2022), between 2022 Nov. 29 and 2025 Jan 10, as part of the NIRPS-GTO. NIRPS is an echelle spectrograph at the ESO 3.6−m telescope of the La Silla Observatory in Chile, with a wavelength range of 980−1800 nm. The instrument is equipped with a low-order adaptive optics system and two observing modes, High Accuracy (R ∼ 85000, 0.4^′′ fibre) and High Efficiency (R ∼ 70000, 0.9^′′ fibre) that can be utilized simultaneously with HARPS. All GJ 3090 observations used here are in the high-efficiency mode. Most data points were obtained by binning two 900s exposures done in the same night. We excluded from our analysis 14 night-binned spectra of shorter 300 s exposures acquired either during the NIRPS commissioning phase or for transit monitoring in another work package, as they exhibit higher dispersion and lower S/N than expected from photo-noise scaling. This conservative selection ensures a homogeneous dataset with consistent noise properties. The observations were reduced with the data reduction software APERO v0.7.292 (Cook et al. 2022), fully compatible with NIRPS, with RV extraction performed using the line-by-line method implemented in the open-source lbl package³, as described in Artigau et al. (2022). At the end, we are left with 84 data points and a median RV precision of 1.8 m s⁻¹ per point. We built a template spectrum of GJ 3090 using the telluric-corrected spectra from APERO in order to derive stellar parameters independently from the literature and to obtain the abundances of several elements with spectral features in the NIRPS spectral range. The top panel of Figure A.3 shows the last 59 RV measurements of NIRPS.

2.4 HARPS radial velocities

A22 obtained 55 radial velocities of GJ 3090 with HARPS in high-accuracy mode (HAM). We also obtained 26 new HAM measurements and 33 EGGS (high-efficiency mode) measurements simultaneous with NIRPS. All simultaneous observations with NIRPS were obtained with the same exposure times and binning procedure. The data were reduced with the HARPS DRS3.2.5 and radial velocities were extracted with version 0.65.003 of the lbl package. The template spectrum used to derive the radial velocities of the HAM mode was obtained by combining the archival and contemporaneous spectra in order to increase its S/N. HAM and EGGS were reduced with their own template because of the difference in their spectral resolution. We obtain 82 HAM points with a median RV precision of 1.9 m/s and 33 EGGS points with a median RV precision of 1.5 m/s for a total of 109 HARPS measurements. The top panel of Figure A.3 shows the last 45 RV measurements of HARPS.

2.5 HARPS (HAM) differential stellar temperatures

The ΔT time series was extracted from the HARPS (HAM) spectra following the method introduced in Artigau et al. (2024). This activity indicator measures the temperature difference between photospheric regions contributing to line asymmetries and the surrounding quiet stellar surface, providing a direct proxy for stellar heterogeneities such as spots and plages. Unlike traditional activity indicators, for example the H α index or bisector span, ΔT is grounded in physical temperature variations with a sub-Kelvin-precision. The ΔT dataset contains 78 measurements extracted from the HAM spectra with a median precision of 0.4 K (see A.3). Of all activity indicators available in lbl, it was the one with the strongest signal around the rotation period of GJ 3090 found in A22. We did not use the differential temperatures from either NIRPS or HARPS (EGGS), as both showed a weak correlation with the RVs, making the added complexity unjustified. In particular, including the NIRPS indicator would have required an additional GP dimension.

3 Spectroscopic stellar parameters

The effective temperature and chemical abundances of GJ 3090 are derived from the NIRPS template spectrum following the methodology of Jahandar et al. (2024, 2025). Using a grid of PHOENIX ACES stellar models (Husser et al. 2013) interpolated to log g=4.75 (for GJ 3090, log g=4.727 ± 0.029; A22) and convolved to the resolution of NIRPS, we fit a selection of strong spectral lines through χ² minimization, yielding T_eff=3709 ± 34 K. The abundances of specific chemical elements are then measured by fitting individual spectral lines for a fixed T_eff of 3700 K. We report the average abundances obtained from those fits in Table 1, including the overall metallicity of [M/H] = 0.03 ± 0.08, computed based on the average abundance of all the measured elements (Jahandar et al. 2024). We note the relative depletion of OH, which was observed in a previous sample of 31 M dwarfs observed in the near-infrared (Jahandar et al. 2025). It is still unknown whether this is due to an unusual depletion of oxygen in M dwarfs compared to FGK stars or a systematic underestimation of OH lines in the models. From the NIRPS spectrum, we measure the abundances of key refractory elements composing the bulk of planetary cores and mantles (Fe, Mg, and Si). The molar iron-to-magnesium ratio is Fe/Mg=0.69 ± 0.16, in agreement with the solar value (0.79 ± 0.10; Asplund et al. 2009).

Alternatively, we derived T_eff and metallicity ([Fe/H]) from the measurements of selected lines in the combined HARPS spectrum by using the machine learning tool ODUSSEAS⁴ (Antoniadis-Karnavas et al. 2020, 2024). This code applies a machine learning model trained with the same lines measured and calibrated in a reference sample of 47 M dwarfs observed with HARPS for which their [Fe/H] were obtained from photometric calibrations (Neves et al. 2012) and their T_eff from interferometric calibrations (Khata et al. 2021). With this method, we derived a T_eff=3691 ± 95 K and [Fe/H]=0.02 ± 0.11 dex, in excellent agreement with the values found with NIRPS spectrum.

The adopted T_eff and [Fe/H] are taken as the weighted average of the NIRPS and HARPS measurements, yielding T_eff = 3707 ± 32 K and [Fe/H] = 0.12 ± 0.06 dex. When compared to the measurements of A22, we find an agreement with the [Fe/H] value (−0.06 ± 0.12 dex) while the effective temperature has a 2σ discrepancy (3556 ± 70 K). Still, our measurements are both consistent with recent values obtained with ODUSSEAS using spectra from FEROS (T_eff=3737 ± 104 K and [Fe/H]= 0.04 ± 0.13 dex; Antoniadis-Karnavas et al. 2024).

For the stellar mass (M_*) and radius (R_*) of GJ 3090, we adopt the values of A22 which are computed from the empirical relations of Mann et al. (2015, 2019), calibrated on nearby M dwarfs with interferometric radii and dynamical masses, and based on the absolute K_s-band magnitude from 2MASS (Skrutskie et al. 2006).

We find that the effective temperature of GJ 3090 derived in this work (T_eff=3707 ± 32 K) is in better agreement with the empirical radius and mass (0.516 ± 0.016 R_* and 0.519 ± 0.013 M_*) than the temperature from A22 (T_eff=3556 ± 70 K). Using the Stefan-Boltzmann law with the luminosity from the SED fitting of A22 ($L_{*}=0.{0455}_{-0.0016}^{+0.0018}{L}_{\odot }$), the radius predicted from our T_eff differs from the Mann et al. (2015) empirical radius by only 0.006σ, whereas the radius predicted from the A22 temperature is off by 1.89σ. Converting these radii into masses using the radius-mass relations of Boyajian et al. (2012), our temperature yields a mass prediction that differs from the Mann et al. (2019) value by 1.59σ, compared to 2.76σ when using the A22 temperature.

4 Photometry analysis

4.1 Transit analysis with TESS

We performed the modelling of the TESS photometric data using the juliet framework (Espinoza et al. 2019), which combines several open-source tools for joint analysis of transits, radial velocities, and GPs. Transit modelling is carried out with batman (Kreidberg 2015). For Gaussian process regression, it employed either george (Ambikasaran et al. 2015a) or celerite (Foreman-Mackey et al. 2017a). juliet is using the nested sampling algorithm implemented in dynesty (Speagle 2020), it explores the parameter space by sampling from the priors rather than optimizing a likelihood function.

The transit model simultaneously fits for the stellar density ρ_*, planetary parameters, and photometric jitter. Instead of fitting directly for the planet-to-star radius ratio (p=R_p/R_*) and the orbital impact parameter (b = a/R_* cos i), we followed the parametrization proposed by Espinoza (2018), which involves fitting for the variables r₁ and r₂. This ensures a more complete and physically consistent exploration of the (p, b) parameter space. The explicit analytical expressions relating (r₁, r₂) to (p, b) are given in detail in Espinoza (2018) as well as in the documentation of the juliet package.

We used the power-2 limb-darkening law implemented within juliet by Parc et al. (2025), which has been shown to best reproduce intensity profiles of cool stars (Morello et al. 2017). The corresponding limb-darkening coefficients and their uncertainties for each photometric filter were computed using the LDCU tool⁵ (Deline et al. 2022). This tool, a refined version of the routine by Espinoza & Jordán (2015), calculates limb-darkening coefficients along with uncertainties by propagating the errors on stellar parameters through grids of stellar intensity profiles derived from ATLAS (Kurucz 1979) and PHOENIX (Husser et al. 2013) synthetic spectra. The resulting coefficients were used as Gaussian priors in our transit fit. Since the TESS PDCSAP flux is already corrected for contamination, we fixed the dilution factor to unity. To account for possible underestimation of the photometric uncertainties due to additional systematics, we added in quadrature a jitter term σ_i to the reported uncertainties. Additionally, to model the residual photometric variability, we included a GP component for each sector, using a MatÉrn 3/2 kernel. The GP is characterized by two hyperparameters: the amplitude (σ_GP) and the characteristic timescale (ρ_GP). We assumed a circular orbit so we fixed the eccentricity to zero and used normal priors around the ExoFOP values for the period and transit epoch. We used a normal prior for stellar density using the stellar mass and radii values from A22.

Fig. A.1 shows the four TESS sectors with the resulting model obtained with the juliet fit and Fig. 1 shows the phase-folded light curve of GJ 3090 b. The prior and posterior distributions are shown in Table A.2.

We derived a radius of 2.18 ± 0.06 R_⊕ for GJ 3090 b using the stellar radius from A22. This is in line with the previous values from A22 and A25. Moreover, we obtained an inclination of 86.9 ± 0.2^° and an impact parameter of 0.70 ± 0.05. The period and time of conjunction posteriors of the transit fit is used as priors for the RV fit.

4.2 Magnetic cycle with ASAS-SN

The long time span of the HARPS (HAM) dataset (∼ 6 years) allows us to scan for long-term signals indicative of long-period planets or magnetic cycles. The first hint of its presence is revealed by the ASAS-SN photometry dataset shown in Figure 2. Also, the second epoch of the archival HARPS data separated by 200 days shows a 6 m s⁻¹ offset relative to the first epoch which cannot be explained by a change in the instrument’s stability since no major modifications were made on the spectrograph between those epochs. Our multidimensional GP model fails to model this offset which then shows up in the residuals.

These hints encouraged us to look further into explanations for this long-term signal. To distinguish between the planetary and long-term stellar variability scenarios, we retrieved photometric data from the ASAS-SN survey, as photometry has been shown to exhibit signatures of magnetic activity cycles in stars (Suárez Mascareño et al. 2016). We perform Bayesian inference of a sinusoidal model (Eq. (4)) on the ASAS-SN dataset which is acceptable for long-term variability because few full cycles are observed. Each available camera from the ASAS-SN network was fitted with independent offset and jitter terms. The coloured points in Figure 2 show the magnitudes corrected using the fitted offset of each camera. Contrary to the best-fit model shown in the top panel, the periodogram displayed in the lower panel is computed from the light curve after subtracting only the median value of each camera, ensuring that the periodogram remains model-independent. We used the package UltraNest (Buchner 2021) with a reactive nested sampling algorithm, 440 live points, and a slice sampler performing 44 steps per iteration, distributed over 20 CPU cores. The best-fit curve is shown on top of the data in Figure 2. The posterior on the period is 3370 ± 120 days. There is another strong signal at 1300 d in the periodogram seen in Figure 2 which could also be a candidate for a long-term trend. In the case of the ASAS-SN fit, the 3370 d solution is overwhelmingly preferred by the FAP in the periodogram as well as by model comparison (Δ log z ≈ 90, obtained when comparing a model in which the uniform prior upper bound on the period is set to the RV timespan of ∼ 2000 d with another model in which the upper bound is set to the full ASAS-SN timespan of ∼ 4000 d), but we do not exclude the possibility that the lower 1300 d signal is present in the radial velocities. In order to account for this long-term signal, we decided to add a magnetic cycle model in our main multidimensional fit. We do not fit a magnetic cycle on the TESS photometry since the sectors are too sparse and discontinuous to detect signals on those timescales.

The magnetic cycle amplitudes of all modelled measurements converged to zero, except for the HARPS radial velocities, which yielded a best-fit amplitude of $K_H^{\mathrm{cyc}}=4.{14}_{-0.90}^{+1.07}\mathrm{m}/\mathrm{s}$. This result is not surprising since NIRPS does not yet have the time baseline required to be sensitive to such long-period signals. As for the ΔT indicator, it does not appear to be sufficiently sensitive or stable enough to capture long-term magnetic trends in our dataset.

As shown by the model comparison presented at the beginning of Section 6, a long-term sinusoidal component is statistically favoured. However, given the lack of a significant peak near ∼ 3370 d in the HARPS RV periodogram (see Figure A.2) and the absence of a corresponding signal in the ΔT indicator, we refrain from assigning a definitive physical interpretation. We therefore regard this variability as a candidate magnetic cycle or, more conservatively, an unclassified long-term stellar trend. It is nonetheless important to model it explicitly; otherwise, the GP would be forced to absorb the long-term variability, compromising the physical interpretability of its hyperparameters and biasing the stellar-activity component of the fit (Stock et al. 2023).

Fig. 1 Phase-folded TESS light curves of GJ 3090 b (grey points). The dark green circles are data binned to 10 min. The black lines represent the median model of each sector from the fit.

Fig. 2 Top panel: joint best-fit sinusoidal model to the ASAS-SN photometric data from all available cameras (bf, bj, bF, and bn). The solid lines represent the best-fit model. Bottom panel: Lomb-Scargle periodogram of the ASAS-SN photometry after camera offset corrections. The blue curve shows the power spectrum. The strongest peak corresponds to the best-fit magnetic cycle period of 3370 ± 120 days (red dashed line), while the orange dashed line indicates the stellar rotation period of 18.0 days and the green dashed line shows the sidereal lunar cycle that arises from the applied lunar correction. The grey dashed curve represents the window function.

5 Radial velocity modelling and inference

We denote all measurements by y_{s, q}, their associated uncertainties by Δ y_{s, q}, and their times of measurement by t_{s, q}, where the subscript s refers to the data source (either NIRPS, HARPS HAM, HARPS EGGS or TESS) and the subscript q refers to the type of quantity measured (either visible RVs, infrared RVs, ΔT or Flux).

The measurements are modelled as the sum of a systematic offset (γ), a stellar variability model, which include both rotational modulation (Rot) and a potential long−term magnetic cycle (Cyc), and a Keplerian model (Planet) which is only applied on radial velocities. The rotational modulation is modelled using a multidimensional GP on all physical quantities (q). The GP is trained simultaneously on the TESS photometry, the ΔT activity indicator obtained from the HARPS HAM spectra, and the radial velocities of NIRPS (infrared) and HARPS (visible) which sum up to a 4D GP. In addition to accounting for rotationally induced signals, we model a possible long−term cycle by including a single periodic component. To test whether this additional long−term component is statistically justified, we compared two models: one including all components except the long−term sinusoid (γ + Rot + Planet) and one including all components including the sinusoid (γ + Rot + Cyc + Planet). The model with the long-term component is favoured with Δ ln Z = 10.3. Equation (1) shows the model components for all physical quantities q and the next sections go through how each of the components are modelled: $$\begin{array}{rl}\mathrm{\Delta }{F}_{\mathrm{T}\mathrm{E}\mathrm{S}\mathrm{S}}& =\gamma +\text{\hspace{0.17em}Rot\hspace{0.17em}},\\ \mathrm{\Delta }{T}_{\mathrm{H}\mathrm{A}\mathrm{R}\mathrm{P}\mathrm{S}}& =\gamma +\text{\hspace{0.17em}Rot\hspace{0.17em}}+\text{\hspace{0.17em}Cyc\hspace{0.17em}},\\ \mathrm{\Delta }R{V}_{\mathrm{N}\mathrm{I}\mathrm{R}\mathrm{P}\mathrm{S}}& =\gamma +\text{\hspace{0.17em}Rot\hspace{0.17em}}+\text{\hspace{0.17em}Cyc\hspace{0.17em}}+\text{\hspace{0.17em}Planet\hspace{0.17em}},\\ \mathrm{\Delta }R{V}_{\mathrm{H}\mathrm{A}\mathrm{R}\mathrm{P}\mathrm{S}}& =\gamma +\text{\hspace{0.17em}Rot\hspace{0.17em}}+\text{\hspace{0.17em}Cyc\hspace{0.17em}}+\text{\hspace{0.17em}Planet.\hspace{0.17em}}\end{array}$$(1)

5.1 Source offsets and instrument jitters

Zero-point offsets (γ) and jitters (σ_jit) are included on all y_s,q and Δ y_s,q respectively: $$y_{s,q}\to y_{s,q}-{\gamma }_{s,q}$$(2) $$\Delta y_{s,q}\to \sqrt{\Delta y_{s,q}^2+{\sigma }_{\mathrm{jit};s,q}^2}\forall s$$(3)

Each data source (s) and observable (q) is assigned its own offset parameter and jitter term to align the datasets consistently and model excess noise or potential instrumental drifts not captured by the model, preventing underestimated uncertainties.

5.2 Magnetic cycle correction

The long-term magnetic cycle of period P_cyc is modelled by including a single periodic component, with each time series (q) having its own cycle epoch (τ^cyc) and amplitude (K^cyc) parameters to account for differences in how the cycle manifests across observables. $$y_{s,q}\to y_{s,q}-K_q^{cyc}sin\left(2\pi \frac{t_{s,q}-{\tau }_q^{cyc}}{P_{cyc}}\right)\forall q.$$(4)

This correction is crucial for disentangling planetary signals from long-term stellar activity trends. P_cyc is informed by the continuous ASAS-SN photometry data which shows a long-term signal.

5.3 Planetary signals

The radial velocity variations induced by planetary companions are modelled using the Keplerian equation, which describes the reflex motion of the host star under the gravitational influence of orbiting planets. This component of the model applies only to time series q belonging to the radial velocity subset, denoted R V={Δ RV_HARPS, Δ RV_NIRPS}. The full Keplerian model is given by: $$y_{s,q}\to y_{s,q}-K[cos(\omega +\nu \left(t\right))+ecos\left(\omega \right)]\forall q\in \mathbf{RV},$$(5) where K is the RV semi-amplitude, ω is the argument of periastron, e is the orbital eccentricity, and v(t) is the true anomaly at time t. This formulation captures the full Keplerian motion, allowing for eccentric orbits.

In cases where eccentricity is negligible or unconstrained by the data, a circular orbit approximation can be adopted. The model then simplifies to: $$y_{s,q}\to y_{s,q}+Ksin\left(2\pi \frac{t_{s,q}-t_0}{P}\right)\forall q\in \mathbf{RV},$$(6) where t₀ is the time of inferior conjunction and P is the orbital period. This approximation requires a change in sign relative to Eq. (5) (Stefanov et al. (2025b), hereafter S25). We note that the Keplerian model does not take planet-planet interactions into account.

5.4 Stellar-activity modelling

Gaussian processes are the basis of the stellar activity modelling in this work. GPs are non-parametric models that provide a flexible framework for capturing complex variability without assuming a specific functional form (Rasmussen & Williams 2005). Formally, a GP defines a distribution over functions such that any finite set of points drawn from it follows a joint multidimensional Gaussian distribution. GPs are particularly effective for modelling stellar variability in time-series data, as they do not require prior knowledge of the distribution or evolution of active regions on the stellar surface. However, their flexibility also poses a risk: without careful modelling, GPs can overfit the data and inadvertently absorb astrophysical signals such as planetary signatures (Damasso et al. 2019; Blunt et al. 2023). Stellar activity can generally be described as a quasi-periodic signal.

To model and correct the correlated noise introduced in the RV time series by stellar rotation, we use the multidimensional Gaussian process framework (Rajpaul et al. 2015). This state-of-the-art technique allows for the simultaneous fit of multiple time series using an underlying GP and its derivative following a scheme similar to the F F^′ formalism (Aigrain et al. 2012). In order to model this multidimensional GP, we use the S+LEAF package (Delisle et al. 2022, hereafter D22), which extends the semi-separable matrix formalism introduced by celerite (Foreman-Mackey et al. 2017b) to enable rapid evaluation of GP models, even for large datasets. This computational efficiency is crucial for handling extensive data while maintaining flexibility.

To maintain the linear scaling of the likelihood computation made possible by the semi-separable framework, S+LEAF provides a set of computationally efficient kernels that approximate the behaviour of more general GP kernels, such as the commonly used quasi-periodic kernel. These approximations preserve the key physical characteristics of stellar activity signals while allowing for fast evaluation. In this study, we tested three such kernels available in S+LEAF (D22): Double SHO, MatÉrn 3/2 exponential periodic (MEP), and exponential-sine periodic (ESP). These kernels differ in their mathematical complexity and approach to modelling quasi-periodic stellar activity. The Double SHO represents the simplest method, using two simple harmonic oscillators with basic trigonometric functions to capture the stellar rotation period and its first harmonic. The MEP kernel increases complexity by combining a MatÉrn-3/2 kernel with harmonic oscillators, employing exponential decay functions to better represent the quasi-periodic nature of stellar signals. The ESP kernel is the most sophisticated, utilizing modified Bessel functions to weight multiple harmonics, enabling it to model intricate periodic structures and complex quasi-periodic behaviour.

Equation (7) defines how the underlying GP is modelled on each physical quantity q: $$\begin{array}{rl}\mathrm{\Delta }R{V}_{\mathrm{H}\mathrm{A}\mathrm{R}\mathrm{P}\mathrm{S}}& ={\alpha }_{0}k+{\beta }_0{k}^{\mathrm{\prime }}\\ \mathrm{\Delta }R{V}_{\mathrm{N}\mathrm{I}\mathrm{R}\mathrm{P}\mathrm{S}}& ={\alpha }_{1}k+{\beta }_1{k}^{\mathrm{\prime }}\\ \mathrm{\Delta }{T}_{\mathrm{H}\mathrm{A}\mathrm{R}\mathrm{P}\mathrm{S}}& ={\alpha }_{2}k\\ \mathrm{\Delta }{F}_{\mathrm{T}\mathrm{E}\mathrm{S}\mathrm{S}}& ={\alpha }_{3}k\end{array}$$(7) where k is the underlying GP and k^′ is its derivative. The free parameters α and β scale the contribution of the GP and its derivative to each observable. The first order of the GP (k) represents activity-induced variations in photometry-like indicators, specifically the TESS light curve and the ΔT temperature indicator, while the stellar contribution to the radial velocities is modelled using both the first (k) and second order (k^′) of the same GP. The first-order term generally corresponds to brightness variations caused by the contrast between star spots and the stellar surface as well as suppression of convective blueshift in magnetized regions (spots), while the second-order term accounts for the Doppler shifts induced by these spots as they move across the stellar disk due to rotation (Rajpaul et al. 2015; Aigrain et al. 2012).

A key modelling choice in our analysis is to restrict the ΔT indicator to only the first-order GP component. This decision is motivated by both physical and empirical considerations. Physically, as in photometry, the temperature primarily traces the luminous contrast on the stellar surface rather than the Doppler shifts associated with moving regions, making it more appropriate to model ΔT with a first-order term alone. When allowing both the zeroth-order amplitude (α₂) and the derivative amplitude (β₂) to vary, the posteriors exhibited a symmetry between positive and negative values for these amplitudes, leading to a degeneracy in the solution without improving the fit quality. Restricting the model to only the zeroth-order amplitude removed this degeneracy and still provided an excellent fit to the data, further justifying its exclusion from the final model.

A key advantage of the multidimensional approach is that it mitigates the risk of misattributing planetary power to stellar activity compared to a serial GP approach, where a GP is first trained on activity indicators and then its posterior is used as a prior for the RV fit (Section 4.1.3; Tran et al. 2023). In the serial approach, the GP is optimized separately for each dataset, which can lead to an overestimation of the activity component in the RVs, as the GP has the flexibility to fit any residual variations, including planetary signals or magnetic cycles. This can result in the suppression or distortion of true planetary signals. Conversely, in the multidimensional approach, all datasets (RVs and activity indicators) are modelled jointly within the same covariance matrix, forcing the GP component to explain correlated variations consistently across all observables. This approach has been successfully used to confirm challenging Keplerian signals around highly active stars (Jones et al. 2022; Barragán et al. 2023; González Hernández et al. 2024).

5.5 Choice of priors

All priors can be found in Table A.3. For the planetary orbit parameters e and ω, we used uniform priors on $\sqrt{e}cos\omega$ and $\sqrt{e}sin\omega$ between −1 and 1, as this has been shown to better explore the parameter space (Lucy & Sweeney 1971). As the period of GJ 3090 b is more precisely constrained by the transit data than by the radial-velocity data, we chose to use the posteriors from the transit analysis (see Section 4.1) for the period and time of transit as priors in our analysis, since the additional TESS sector available since the system’s discovery has improved the constraints on these parameters. To search for additional planets, we initially ran a model with broad priors on the period (𝒰(3, 200)) to conduct a blind search. We then used the posterior from this fit to define a more constrained uniform prior for a second fit, allowing for improved posterior sampling and more robust Bayesian evidence measurements (Aitken & Akman 2013). We use a normal prior on P_cyc centered on the 3370 d periodicity identified by ASAS-SN (see Section 4.2) with an inflated uncertainty to include the 1300 d solution at 3σ. The rest of the parameters were given large uniform priors except for jitter parameters and GP or cycle amplitudes which were given log-uniform priors as prescribed in Stock et al. (2023).

5.6 Parameter inference

We optimized the model parameters using Bayesian inference via nested sampling (Skilling 2004, 2006) with the dynesty package (Speagle 2020). To ensure thorough exploration of the parameter space, we used a high number of live points, setting nlive = 200 × N_dim, where N_dim is the number of free parameters in the model. This high-density sampling improves the accuracy of posterior estimations, which are crucial for model comparison.

We employed the multi-ellipsoid bounding method (bound=‘multi’) to efficiently constrain the parameter space and reduce computational overhead. Sampling was performed using the random walk algorithm (sample=‘rwalk’), which is well-suited for handling correlated parameter distributions, with an increased number of MCMC walk steps (walks =10 × N_dim) to enhance convergence.

6 Results and discussion

6.1 Determining the number of planets in the GJ 3090 system

In addition to the known transiting planet GJ 3090 b, several significant signals appear in the periodogram of the combined RV dataset (top panel of Figure A.2). The most prominent peak occurs at 16 days, close to the estimated 18-day stellar rotation period (A22). Another noticeable signal appears near 13 days, which was first reported by A22 as a potential planetary candidate.

To assess the most plausible planetary configuration for the system, we conducted a rigorous Bayesian model comparison using the nested sampling log-evidence (ln Z). Figure 4 displays the Δ ln Z values for all tested combinations of stellar activity models (GP kernels) and planetary configurations similarly to S25. Each row corresponds to a planetary model with an increasing number of Keplerian signals, while each column represents a different GP kernel. We adopt the notation C for circular orbits (e=0) and K for Keplerian models with free eccentricity (e ≥ 0). The colour scale encodes the log-evidence relative to the best overall model, offering a joint view of stellar and planetary model complexity.

Our highest-evidence model (labelled b4) combines three circular Keplerian signals (C_b, C_c, and C_d) with periods of 2.85 days, 15.94 days, and 12.71 days respectively, plus a MEP kernel for stellar activity. The phase-folded radial velocity curves for all three planets are shown in Figure 3, illustrating their individual contributions. Figure A.3 displays the GP model prediction overlaid on both the RV datasets and activity indicators. The top-right panel shows the full model (Keplerian plus GP) on the most recent epoch of observations, demonstrating excellent visual agreement with the data in particular for NIRPS. The median residual values after subtracting the full model are 1.15 m s⁻¹ for NIRPS, 2.68 m s⁻¹ for HARPS (HAM), and 3.65 m s⁻¹ for HARPS (EGGS).

To quantify how the model evidence evolves with planetary multiplicity, we compute Δ ln Z in three different ways: let {ln Z_k}_i denote the set of log-evidence for planetary model i across all tested kernels k. First, we compute the average logevidence difference (see Eq. (8)) between two configurations i and j where N_k is the total number of tested GP kernels. This provides a global view of how favourable one planetary model is over another across all GP kernels. Second, we adopt a more conservative approach by computing the pessimistic (pes) log-evidence (see Eq. (9) which compares the worst-performing kernel of the better Keplerian model i with the best-performing kernel of the worse Keplerian model j. This illustrates the impact of the choice of kernel on planetary detection. Finally, we compute the best log-evidence difference (see Eq. (10)), which compares the log-evidence obtained with the overall best-performing kernel (MEP in our case) for both planetary architecture. This metric represents the most realistic assessment of the models’ relative performance, as it isolates the comparison to the kernel that best captures the stellar activity component. In practice, this best log-evidence is typically the one reported when confirming an exoplanet candidate, since once a kernel is shown to perform significantly better, the planetary configuration should be evaluated on that basis alone. $${(\Delta lnZ)}_{\mathrm{avg}}=\frac{1}{N_k}{\sum }_{k=1}^{N_k}(ln{Z}_{k,i})-\frac{1}{N_k}{\sum }_{k=1}^{N_k}(ln{Z}_{k,j}),$$(8) $${(\Delta lnZ)}_{\mathrm{pes}}={min}_k(ln{Z}_{k,i})-{max}_k(ln{Z}_{k,j}),$$(9) $${(\Delta lnZ)}_{\text{best}}=\underset{k}{best}(ln{Z}_{k,i})-\underset{k}{best}(ln{Z}_{k,j}).$$(10)

All Δ ln Z are then converted into equivalent sigma significances using Equations (8) and (9) of Benneke & Seager (2013).

The two-planet model, including GJ 3090 c at 16 days, is decisively favoured over the single-planet model, with σ_avg = 6.0, σ_pes=3.5, and σ_best=6.7. The second planet is detected with high statistical confidence (>5σ), justifying its classification as a confirmed planet. Its proximity to the stellar rotation period motivates a dedicated discussion in Section 6.4 to exclude stellar activity as its origin.

The three-planet model, which incorporates the 12.7-day signal as a candidate, receives modest additional support compared to the two-planet model with σ_avg=3.1 and σ_best=3.3 while the pessimistic metric prefers the two-planet model by σ_pes=6.0. If we add the 12.7-day planet as a second planet and compare it with the models containing only GJ 3090 b, we get σ_best=4.4. This warrants caution and indicates that additional data or improved stellar activity modelling may be necessary to establish the planetary nature of this signal.

We tested a four-planet model by adding a Keplerian signal with an open prior on the period (𝒰(3, 200)) even though there were no significant peaks in the periodogram of the residuals. We got no clear convergence to a planetary signal.

To quantify the improvement brought by NIRPS, we compare the semi-amplitude significance N_σ=K/σ_K of the planetary signals in the HARPS-only and combined NIRPS+HARPS fits (see Table A.1). For GJ 3090 b, N_σ increases from ∼ 5.5 (HARPS-only) to ∼ 9.8 in the combined fit, and for GJ 3090 c from ∼ 4.7 to ∼ 7.7, demonstrating a clear gain in precision. For the 12.7-day candidate, however, there is no gain: the HARPS-only and NIRPS-only semi-amplitudes differ significantly, and their combination does not improve the detection (N_σ ∼ 4.0).

Finally, as seen in Figure 4, eccentric models are consistently disfavoured by Bayesian evidence, suggesting that the eccentricities of the planets are likely low and not sufficiently constrained by the available data. As a result, we adopt model b4—which includes three circular Keplerian signals and employs the MEP kernel—as the baseline model for the remainder of the discussion.

Fig. 3 Phase-folded radial velocity measurements for the confirmed and candidate planets in the GJ 3090 system. The first panel shows the confirmed planet GJ 3090 b (2.85 days), the second panel presents the confirmed sub-Neptune GJ 3090 c (15.94 days), and the third panel displays the 12.7-day candidate planet. The data points from HARPS HAM, HARPS EGGS, and NIRPS are shown with their associated uncertainties, alongside the Keplerian model curves (black lines).

6.2 Explaining the absence of outer planet transits

Assuming coplanar orbits aligned with the inclination of GJ 3090 b(86.93 ± 0.17 deg) updated in this work, we computed the expected impact parameters of the outer planets using their semi-major axes and the stellar radius. The resulting values are b_c=2.22 ± 0.14 for GJ 3090 c and b_{cand 1}=1.91 ± 0.12 for the 12.7-day candidate. Since both values are significantly greater than 1, the planetary orbits do not intersect the stellar disk from our line of sight. This geometrical configuration naturally explains why only GJ 3090 b is observed in transit, while the outer planets remain undetected in photometry. The non-detection of GJ 3090 c and the 12.7-day candidate is shown in Figure A.6, which plots the phase-folded TESS light curves folded on the best-fit ephemeris of the two outer planets.

To further assess the detectability of potential transits, we used the spright package (Parviainen et al. 2023) to estimate the expected transit depths based on probabilistic mass–radius relations inferred from the distribution of known exoplanets. For GJ 3090 c, the predicted depth is ${2774}_{-804}^{+1135}\mathrm{ppm}$, well above the TESS detection threshold. For the 12.7-day candidate (GJ 3090 d), the estimated depth depends on its nature: if it is a super-Earth, the expected depth is ${679}_{-89}^{+97}\mathrm{ppm}$, whereas if it is a sub-Neptune, the depth increases to ${1171}_{-378}^{+459}\mathrm{ppm}$. Since its inferred mass is close to those of exoplanets lying in the radius valley, both compositional scenarios are plausible. However, in all cases, such depths would have been detectable in the available TESS photometry. This supports the interpretation that the absence of transit signals is due to geometric misalignment rather than insufficient signal strength.

Fig. 4 Heat map similar to that in S25 showing the comparison of Bayesian evidence (Δ ln Z) for different models. The x-axis displays different stellar-activity kernels in increasing order of complexity, while the y-axis shows different planetary configurations. Circular orbits (e= 0) are labelled ‘C’, while Keplerian orbits (e ≥ 0) are labelled ‘K’. The red rectangle highlights the best model based on the maximum Bayesian evidence.

6.3 New constraints on GJ 3090 b

Using the combined HARPS and NIRPS radial velocity datasets which include 152 new RV measurements since A22, we present updated constraints on the mass, radius and bulk density of GJ 3090 b. To get an absolute value for the mass, we use the inclination measured. Our analysis yields a mass of 4.52 ± 0.47 M_⊕, while A22 reported a value of 3.34 ± 0.72 M_⊕. The new mass is in good agreement with a reanalysis of the archival HARPS dataset made by Osborne et al. (2024), hereafter O24, that found $4.{59}_{-1.04}^{+1.13}{M}_{\oplus }$ by fitting a multidimensional GP model on RVs, FWHM (Full Width Half Maximum) and BIS (Bisector Inverse Slope). The mass obtained in this work is derived from the posterior of our highest-evidence model (b4), which includes three circular planets and uses the MEP kernel for activity modelling. The relative error on the mass of GJ 3090 b thus went from 22% (4.6σ, A 22), reevaluated at 24%(4.2σ, O 24) and finally to 10% (9.6σ, this work). The high precision on the mass of GJ 3090 b achieved in this work is primarily attributable to the GP modelling, the extended temporal baseline of the combined datasets, and the enhanced precision of the NIRPS measurements.

To ensure that the inferred mass of GJ 3090 b is robust against modelling choices, we compare its posterior distribution across five different model configurations. Three models vary the number of Keplerian signals (1,2, or 3 planets) while holding the kernel fixed (MEP), and two others vary the kernel type while keeping the three-planet configuration. As shown in Figure 5, the resulting mass posteriors are consistent well within 1σ, indicating that neither the inclusion of additional planets nor the choice of GP kernel significantly affects the mass constraint on GJ 3090 b.

A22 previously measured a radius of 2.13 ± 0.11 R_⊕ from 3 TESS sectors, which, combined with their mass estimate, implied a bulk density of $1.{89}_{-0.45}^{+0.52}\mathrm{g}{\mathrm{cm}}^{-3}$. More recently, a JWST transmission spectroscopy study of the planet (A25) refined the radius measurement to 2.171 ± 0.068 R_⊕ which is compatible with our re-analysis of the 4 available TESS sectors which yielded a radius of 2.18 ± 0.06 R_⊕. When combined with our updated mass, this yields a new bulk density of $2.{40}_{-0.30}^{+0.33}\mathrm{g}{\mathrm{cm}}^{-3}$.

This revised density places GJ 3090,b closer to the bulk population of known sub-Neptunes orbiting M dwarfs, which typically have mean densities of ∼ 2−3 g cm⁻³ (e.g. Parc et al. 2024). Furthermore, the updated mass and radius yield a lower transmission spectroscopy metric (TSM; Kempton et al. 2018). While A22 reported a TSM of ${221}_{-46}^{+66}$, we obtain a revised value of ${182}_{-23}^{+27}$ using the new parameters derived in this work. This decrease in TSM may partially account for the lack of detectable atmospheric features reported by A25 and Parker et al. (2025).

Fig. 5 Mass estimates of GJ 3090 b under different planetary configurations and different stellar-activity kernels. The models range from single-planet solutions (GJ 3090 b only) to multi-planet scenarios, including the 16-day confirmed planet (GJ 3090 c) and candidate planet at 13 days. The three-planet models are also shown with different kernels and a model with a non-zero eccentricity on GJ 3090 b. The adopted model is shown in black. The mass constraints from A22, O24, and an independent reanalysis of the 55 archival HARPS data points are also shown. A dashed horizontal line separates the models tested in this work from previously published mass estimates and a re-analysis of the HARPS archival data.

6.4 Assessing GJ 3090 c’s proximity to the stellar rotation period

The newly confirmed planet GJ 3090 c has an orbital period of 15.9407 ± 0.0058 days, which lies close to the star’s measured rotation period of 17.95 ± 0.15 days. Given that M dwarfs are known to exhibit differential rotation (Donati et al. 2008; Reinhold & Gizon 2015), we have to consider the possibility that the 16-day signal observed in the radial velocity (RV) data is an artefact of stellar activity rather than a planetary companion. To rigorously test this hypothesis, we conducted several analyses.

We begin by verifying the consistency of the 16-day signal across instruments. Unlike Keplerian signals, stellar activity-induced radial velocity (RV) variations are known to be chromatic (Reiners et al. 2010; Marchwinski et al. 2014; Trifonov et al. 2018), meaning their amplitude depends on the observed wavelength. Therefore, if a signal appears with consistent amplitude and phase in both the redder NIRPS and the bluer HARPS datasets, this supports a planetary interpretation.

To test this, we reran the RV analysis using only the NIRPS or only the HARPS data. In each case, we removed the other instrument’s RV time series from the model in Equation (1), while keeping all other components—including the photometry and activity indicators—unchanged. The same prior distributions were used as in the fit to the combined datasets.

As shown in Figure 6, both the HARPS-only and NIRPS-only fits yield highly consistent estimates for the orbital period, time of conjunction, and RV semi-amplitude of GJ 3090 c. To quantify this agreement, we compute the normalized difference between the two posteriors: $$\frac{|x_{\mathrm{N}}-x_{\mathrm{H}}|}{\sqrt{{\sigma }_{\mathrm{N}}^2+{\sigma }_{\mathrm{H}}^2}}.$$(11) Here x and σ represent the medians and 1σ uncertainties from the NIRPS (N) and HARPS (H) fits, respectively. This metric quantifies the level of agreement between two independent measurements in units of their combined uncertainty. A value close to zero indicates strong consistency, while values approaching or exceeding unity suggest significant tension. For the period, RV semi-amplitude, and time of conjunction of GJ 3090 c, we obtain normalized differences of 0.22,0.72, and 0.16, respectively. These values lie within the combined 1σ uncertainty, confirming the agreement between the two instruments. This indicates that the signal is achromatic and thus unlikely to originate from stellar activity. It provides strong support for the planetary nature of GJ 3090 c. Since there have been instances of achromatic activity signal in the literature (e.g. CortÉs-Zuleta et al. 2023; Larue et al. 2025), more tests are necessary.

A widely used diagnostic to assess the planetary nature of RV signals is to monitor the evolution of the false alarm probability (FAP) in a periodogram as additional observations are added chronologically (e.g. S25). This approach tests whether the signal persists consistently throughout the full time span of observations, rather than being concentrated in a specific epoch, an outcome that could be indicative of differential rotation.

Figure 7 shows the chronological evolution of the FAP for the 16-day signal in the RVs, compared to that in the ΔT and FWHM activity indicator, using all available HARPS (HAM and EGGS modes) and NIRPS data. As expected for a genuine Keplerian signal, the FAP associated with the RVs steadily decreases as more data are added, indicating a coherent and persistent signal across the full ∼ 6-year baseline. In contrast, the FAP of the same signal in the activity indicators remains negligible throughout the dataset, suggesting that the 16-day signal does not originate from photospheric activity. These results support the long-term dynamical stability of the GJ 3090 c signal.

To further assess the nature of the 16-day signal, we examine how its removal affects the correlation between HARPS (HAM) radial velocities and its differential temperature indicator ΔT. Because our RV model includes both first- and second-order activity terms, we analyse correlations with ΔT and its time derivative d ΔT/d t.

To compute the derivative of the differential temperatures, we employed a MEP kernel to model the ΔT time series from the HARPS (HAM) observations using the best-fit hyperparameters listed in Table A.3. We then calculated the time derivative of the GP prediction and interpolated these values onto the observation times while computing their uncertainties by sampling the posterior of the GP parameters and recomputing the predictions.

Before assessing these correlations, an important step is required. When plotting two out-of-phase periodic signals of the same period with different amplitudes on the y and x axis (in our case RVs and ΔT) of a graph, the emergent shape is an ellipse, similar to what was observed for the correlation between RV and FWHM for Proxima and GJ 526 (Suárez Mascareño et al. 2020; Stefanov et al. 2025a). To make the shapes linear, it is necessary to first subtract the best linear fit of each correlation (RVs vs ΔT and RVs vs d ΔT/d t) to the RVs of the opposing graph. This effectively isolates the first- and second-order contributions.

Subsequently, we repeated the linear fits after subtracting the best-fit Keplerian model for the 16-day signal from the RV data. If the signal is of stellar origin, the correlation between the RV residuals and the activity indicators (ΔT and d ΔT/d t) should weaken, as part of the activity signal which causes these correlations would have been removed from the RV data. Conversely, if the signal corresponds to a planetary origin, the correlations should strengthen, since we would have removed a signal to which the indicators are blind, isolating the remaining activity-driven variations.

To quantitatively assess the correlations, we used two metrics: the Pearson correlation coefficient (R) and the statistical significance of the linear slope, expressed as the number of standard deviations from zero slope (N_σ). The results are shown in Figure A.5. Notably, both correlations improve after removing the 16-day signal, with the second-order correlation (between RV and d ΔT/d t) showing a more significant improvement: Δ R_Pearson=0.10 and Δ N_σ=3.4. This increase in correlation strength indicates that the 16-day signal is unlikely to be caused by stellar activity and more likely originates from a planetary companion.

Fig. 6 Posterior distributions of the Keplerian parameters of GJ 3090 c derived from fits to the combined dataset (black), HARPS-only data (blue), and NIRPS-only data (red). The distributions show excellent agreement between HARPS and NIRPS, with tighter constraints in the combined dataset due to the increased number of RV points.

Fig. 7 Chronological evolution of the false alarm probability (FAP) for the 16-day signal. As new observations are added, the FAP in the RVs decreases steadily, consistent with a stable Keplerian signal. In contrast, the FAP in the ΔT and FWHM activity indicators remains high across the datasets. The orange and green horizontal lines represent the 3σ and 5σ detection significance respectively. This behaviour supports the long-term coherence of the GJ 3090 c signal and disfavours a stellar origin.

6.5 Comparison of GP frameworks

In this section, we aim to compare the performance of the serial, semi-shared, and multidimensional GP models and assess the sensitivity of each model to the 15.9-day planetary signal using an injection-recovery analysis. Specifically, which models are more prone to overfitting coherent signals, regardless of their origin or period, by absorbing them into the GP component and potentially misinterpreting them as stellar variability. To test this, we construct a synthetic dataset containing only stellar activity, simulated using a quasi-periodic kernel from the george package (Ambikasaran et al. 2015b), and add white noise matching the median noise level of the TESS, HARPS, and NIRPS datasets. The time sampling is identical to the real datasets, and we ensure that the root-mean-square (RMS) scatter matches that of the four observed time series.

We then inject circular Keplerian signals with varying periods and semi-amplitudes into this synthetic dataset. For each injection, we fit the data using our S+LEAF model and compute the Lomb-Scargle periodogram of the residuals. This process is repeated across a grid of periods from 1 to 1000 days and semi-amplitudes from 0.5 to 10 m s⁻¹. The resulting false alarm probabilities (FAPs) at each grid point are converted into detection significances using the standard inverse Gaussian cumulative distribution function (CDF), $$\sigma ={\Phi }^{-1}(1-\mathrm{FAP}),$$(12) where Φ⁻¹ is the quantile function (CDF⁻¹) of the standard normal distribution.

The resulting sensitivity map is shown in Figure 8, where the detection significance is encoded in the colour scale. Known planets from our dataset are overlaid with yellow markers, and vertical orange lines indicate the stellar rotation period (P_rot) and its first harmonic.

We also compare this multidimensional GP model to two other types of model which inform stellar activity from activity indicators in different ways: The serial model and the composite or semi-shared model. In the serial model, GP hyperparameters are inferred from activity indicators and their posteriors are used as priors to fit the RVs. In the composite or semi-shared model, each time serie has its own covariance matrix, but some GP hyperparameters are fitted simultaneously between all datasets. In our case the rotation period and the length scale of the GPs are shared.

As illustrated in Figure 8, the multidimensional GP model exhibits greater sensitivity than the serial GP across multiple regions of parameter space. It is the case near the stellar rotation period (∼ 18d), where activity-induced signals dominate. The multidimensional approach also shows a greater sensitivity at long periods where the detection region of the serial and composite GP models is restricted by overfitting. This means that the multidimensional model should be prioritized when looking for long period planets or magnetic cycles. In our case, the composite method seems to give better detection significance for small periods as well as near the first harmonic of the rotation period.

This injection test informs us that the 16-day signal corresponding to GJ 3090 c is recovered with a significance of 4.05σ using our multidimensional GP framework, whereas it is detected with lower significance with other methods: 1.10σ when using the composite GP approach and a null detection (0σ) for the serial GP. While the recovered significance is only ∼ 4σ, this detection map is based on a false alarm probability (FAP) statistic, which itself relies on the analytical approximation from Baluev (2008). As such, it should not be directly compared with the detection significance derived from Bayesian model comparison using the log-evidence difference (Δ ln Z). Some variation between the σ detection inferred from a FAP and that from Δ ln Z is therefore expected. In our study, this diagnostic primarily serves to illustrate the relative performance of the different GP approaches. Nevertheless, this exercise strengthens our conclusion that the multidimensional GP model is the most effective approach to prevent overfitting the planetary signal. Instead, it effectively distinguishes the planetary nature of GJ 3090 c from stellar activity. The lower sensitivity of the serial GP method, combined with the more limited dataset used in A22, may explain why this signal was not previously reported.

We note that the same analysis using the actual residuals from the data (after subtracting known planetary and magnetic cycle signals) has been tried, and the results were consistent with those from the synthetic dataset. We chose to show the results from the latter to avoid contamination from unmodelled signals in the real data (e.g. low-amplitude planets).

Fig. 8 Sensitivity map showing the detection significance (σ) of the multidimensional GP approach as a function of period and semi-amplitude K of injected planets. The colour map indicates the detection significance computed from the periodogram of the residuals. The 5σ detection contours are shown for the serial GP (dotted grey line), the composite or semi-shared GP (dashed grey line), and the multidimensional GP (solid grey line). Known planets are marked with yellow symbols. The vertical orange lines indicate the stellar rotation period and its first harmonic. The 16-day signal is recovered with a detection significance of 4.05σ using the multidimensional GP, 1.10σ using the composite GP and 0σ with the serial GP.

6.6 Stability and resonant analysis

We analyse the stability of the GJ 3090 system with the additions of both GJ 3090 c and the 12.7-day candidate and analyse the system’s resonant state using the planetary parameters derived in this work. We first perform an introductory simulation to determine whether or not this system is potentially stable. To assess the state of the system, we run a numerical simulation using Rebound and the default IAS15 integrator (Rein & Liu 2012). We initialize GJ 3090 and all three planets or planet candidates using the median values of planetary masses, periods, eccentricities, and epochs given in Table 2. As GJ 3090 c and the 12.7-day candidate are not transiting, we do not know their inclinations relative to GJ 3090 b; as such, we assume all three orbits are coplanar. We integrate the system for 1 Myr, saving planetary parameters for all three bodies every 1000 years.

We find that all three planets remain stable for the entirety of the 1 Myr integration, with periods for all three planets remaining nearly unchanged and eccentricities remaining very low (<0.02) for the whole integration. While further analysis of the stability of the GJ 3090 system is beyond the scope of this paper, it appears that the GJ 3090 system is not immediately unstable with the addition of this planet candidate.

In addition, we note that GJ 3090 c and the 12.7-day candidate have a period ratio very close to 5: 4, potentially suggestive of the planet being in a two-body mean motion resonance (MMR). As such, we attempt to identify if such a resonance exists between these two planets.

Two-body MMRs exist when a pair of planets have librating resonance angles φ₁, φ₂, given by $${\varphi }_1=p{\lambda }_i-q{\lambda }_j-(p-q){\omega }_i,$$(13) $${\varphi }_2=p{\lambda }_i-q{\lambda }_j-(p-q){\omega }_j,$$(14) where λ is the standard mean longitude of the planet and ω is the argument of periapsis of the planet, while the period ratio of the two planets is p/q for integers p, q (in our case, 5 and 4, respectively). If these two angles librate rather than circulate, the two planets are said to be in an MMR. Alternatively, one can use the composite two-body resonance angle φ₁₂, given by Laune et al. (2022) to be $$tan{\varphi }_{12}=\frac{f_1{e}_{1}sin{\varphi }_1+f_2{e}_{2}sin{\varphi }_2}{f_1{e}_{1}cos{\varphi }_1+f_2{e}_{2}cos{\varphi }_2},$$(15) where e is the eccentricity of a planet and f₁, f₂ are coefficients that depend on the values of p, q (see Table 1 of Deck et al. 2013). Planets can still be in resonance even if φ₁, φ₂ circulate, so long as the composite two-body resonance angle still librates.

We calculate all three of these angles every 1000 years across the entirety of the integration. We additionally increase the sampling rate drastically for the first 10 years of the integration, sampling 1000 times per year for only this period. The results are shown in Figures 9 and 10, respectively.

As can be seen in Figure 9, all three resonant angles appear to be librating around 180^°, 0^°, and 0^°, respectively, which would be indicative of GJ 3090 c and the 12.7-day candidate indeed being in a 5: 4 two-body MMR. However, as can be seen in Figure 10, a higher sampling frequency shows a very different result. The angles φ₁ and φ₂ appear to be circulating, just with a circulation frequency very near to 1 year. Meanwhile, φ₁₂ appears to be circulating much faster. While this integration shows that the system is not in resonance, this system appears to be very nearly resonant, with the slow circulation of φ₁ and φ₂ that takes many orbits to circulate potentially even indicating nodding in and out of resonance (Ketchum et al. 2013). As the resonant state of a system can depend on very slight changes in planet mass or period, however, more precise determination of planetary parameters is necessary to better determine the resonant state of the system.

Fig. 9 Plot of samples of φ1 (orang), φ2 (blue), and φ12 (green) between GJ 3090 c and the 12.7-day candidate. Samples are taken every 1000 years, across a 1 Myr integration of the GJ 3090 system. Horizontal dashed lines corresponding to 0° and 180° are shown as well. All three angles appear to librate, though this is merely an artefact of the sampling frequency.

Fig. 10 Plot of samples of φ1 (orange), φ2 (blue), and φ12 (green) between GJ 3090 c and the 12.7-day candidate. Samples are taken 1000 times a year, across the first 10 years of a 1 Myr integration of the GJ 3090 system. Horizontal dashed lines corresponding to 0° and 180° are shown as well. While Figure 9 appears to show all three angles librating, the higher sampling frequency shown here shows that all three angles circulate.

Fig. 11 Mass-radius diagram. GJ 3090 b (green) is shown alongside known M dwarf exoplanets with <30% uncertainty in both mass and radius (grey points). The red symbol indicates the mass–radius position of GJ 3090 b as previously measured by A22. Theoretical iso-composition lines include purely rocky compositions (black and red), H2O envelope (steam), planets with Earth-like rocky interior (dotted blue), and H/He envelope sub-Neptunes with Earth-like rocky interior (orange). These interiors are retrieved using the detailed modelling described in Sect. 6.8 and assume Te q=700 K with 1 Gyr system age.

6.7 Updated parameters of the GJ 3090 system

We update the stellar rotation period of the host star GJ 3090 to 17.95 ± 0.15 days, a small increase from the previous estimate of $17.{729}_{-0.036}^{+0.17}$ days reported by A22. This 1.3σ shift may be explained by the presence of the 16-day planetary signal associated with GJ 3090 c, which was most likely absorbed by the GP activity model used in the previous study. We recover both the 18-day rotational period and its first harmonic from the RV data after removing all planetary contributions, as shown in the periodogram in Figure A.2. The forest of peaks near the rotation period in the periodogram shows that the star most likely experienced differential rotation during the ∼ 6 yr timespan of our dataset.

We provide refined mass and radius measurements for GJ 3090 b of 4.52 ± 0.47 M_⊕ and 2.18 ± 0.06 R_⊕ respectively, and confirm the existence of a second planet, GJ 3090 c, orbiting the star with a period of 15.94 days and a mass of 10.0 ± 1.3 M_⊕. We also detect a 12.7-day planetary candidate originally reported by A22, now with an updated mass of 4.15 ± 1.04 M_⊕.

GJ 3090 b is placed in the mass–radius (M–R) diagram shown in Figure 11. This figure includes known exoplanets orbiting M dwarfs with better than 30% precision in both mass and radius (grey points). We overlay theoretical M-R relations from several interior models: from pure iron to pure silicate compositions (Zeng et al. 2019), water-rich envelopes (Aguichine et al. 2021), and H/He envelopes (Lopez & Fortney 2014). The revised mass of GJ 3090 b increases its density and brings it closer to the bulk population of sub-Neptunes around M dwarfs. Its position in the M-R diagram lies within a region of degeneracy between water-rich and H/He-dominated compositions, suggesting a range of plausible internal structures.

In terms of incident stellar flux, all three planets receive significantly more radiation than Earth, as summarized in Table 2. This places them well outside the habitable zone of the system, which is estimated to lie around orbital periods of ∼ 100 days. Finally, while we cannot completely rule out the existence of additional planetary companions in the system, no statistically significant peaks are observed in the periodograms of the RV residuals after modelling the three identified planets and stellar activity (Figure A.2).

6.8 Detailed interior modelling of GJ 3090 b

We investigate the composition of GJ 3090 b using detailed interior structure modelling coupled to the emcee affine-invariant ensemble sampler (Foreman-Mackey et al. 2013) to characterize the probability density of possible atmospheric mass fraction (AMF) and core mass fraction (CMF). The planetary interior is assumed to be composed of an Fe−Ni metallic core and a silicate mantle (SUPEREARTH Valencia et al. 2007), while the outermost layer consists of either a hydrogen-helium envelope or a water vapour atmosphere modelled using the CEPAM code (Guillot & Morel 1995), with equations of state from Saumon et al. (1995) for H/He and French et al. (2009) for H₂O. In all cases, we assume that the rocky interior contains no volatiles and follow the numerical set-up given in Plotnykov & Valencia (2020).

We consider two sub-Neptune composition scenarios: (1) the planet has H/He envelope (75% H₂ to .25% He) and (2) the planet has a pure H₂O envelope. For these scenarios, we impose stellar-informed priors on the rocky interior based on the host star’s refractory abundances taken from Table 1: $$\begin{array}{cc}\hfill \mathrm{Fe}/{\mathrm{Mg}}_{\text{planet}}& \sim \mathcal{N}(\mathrm{Fe}/{\mathrm{Mg}}_{\text{star}},{\sigma }_{\text{star}}^2)\hfill \\ \hfill \mathrm{Fe}/{\mathrm{Si}}_{\text{planet}}& \sim \mathcal{N}(\mathrm{Fe}/{\mathrm{Si}}_{\text{star}},{\sigma }_{\text{star}}^2)\hfill \end{array}$$ where all ratios are by weight. Additionally, the mantle mineralogy is allowed to vary in terms of the Bridgmanite to Wustite ratio (MgSiO₃ vs MgO, x_Wu). These assumptions effectively constrain the rocky core mass fraction (CMF $=\frac{\mathrm{rcmf}}{\mathrm{rcmf}+1}$, where rcmf is the core-to-mantle mass ratio, M_core/M_mantle) of the planet and mitigate the problem of compositional degeneracy.

Considering case (1), where GJ 3090 b has retained its primordial H/He envelope, we recover AMF=0.006 ± 0.002 and $\mathrm{CMF}=0.{35}_{-0.05}^{+0.04}$ when using stellar-informed priors (Appendix A.7). If we impose no prior on the rocky composition, the AMF increases slightly to $0.{007}_{-0.04}^{+0.06}$, while the CMF distribution becomes much broader, with CMF =0.5 ± 0.3 (also shown in Appendix A.7). For case (2), where the envelope is composed of pure water vapour, we find $\mathrm{WMF}=0.{43}_{-0.07}^{+0.08}$ and $\mathrm{CMF}=0.{35}_{-0.05}^{+0.04}$, with the posterior distributions shown in Appendix A.8.

6.9 Linking interior structure and atmospheric studies of GJ 3090b

The results of our interior structure modelling indicate that both the AMF and the WMF should be considered as upper limits. This is because degeneracies in the mass–radius relationship prevent a definitive separation between a purely H/He-dominated or water-dominated envelope, and the planet is more likely to possess a mixed composition with both components present. In particular, the constraint that AMF is below 1% is well under the predictions of atmosphere formation models, which typically expect initial H/He envelopes of at least a few percent by mass for sub-Neptunes at this size and insolation (Ginzburg et al. 2016; Venturini & Helled 2017; Rogers et al. 2023). Such a low H/He fraction would require either formation in a gaspoor environment (Lee & Connors 2021), rapid dissipation of the gas disk (Ercolano & Pascucci 2017), or significant loss through atmospheric escape or impacts (Inamdar & Schlichting 2016).

High-resolution CRIRES+ spectra of GJ 3090 b (Parker et al. 2025) revealed no signs of molecular species, consistent with either a highly metal-enriched atmosphere with mean molecular weight >7.1 g mol⁻¹ and H/He fraction <33%, or the presence of high-altitude aerosols at pressures below 3.3 × 10⁻⁵ bar that obscure molecular features. Similarly, JWST observations (A25) show highly muted spectral features and a lack of CH₄. A key result from the JWST study is the detection of metastable helium absorption at 10833Å at 5.5σ significance, representing the first such detection in a sub-Neptune with JWST. This demonstrates the presence of a H/He atmosphere, even if its overall mass fraction is small. The measured helium signal is an order of magnitude weaker than solar-metallicity models predict, and the most plausible explanation is that the atmosphere is highly enriched in heavy elements, which would reduce both the mass-loss rate and the amplitude of the helium feature. Based on hydrodynamic modelling, A25 estimate that photoevaporation would require ∼ 50 Gyr to strip the atmosphere entirely, implying that the current envelope is long-lived.

The combination of an H/He detection, the lack of molecular features in both CRIRES+ and JWST spectra, and the high inferred metallicity strongly supports the presence of a substantial water component mixed with H/He, making GJ 3090 b a compelling water-world candidate. This interpretation is reminiscent of the atmospheric properties inferred for K2-18 b (Benneke et al. 2019; Madhusudhan et al. 2020) and TOI-270 d (Holmberg & Madhusudhan 2024), both of which exhibit evidence for volatile-rich envelopes with high metallicities. In the case of GJ 3090 b, the atmospheric constraints thus provide a coherent picture with the interior structure modelling, pointing towards a volatile-rich planet where water likely constitutes a significant fraction of the envelope.

7 Summary and conclusions

We presented an updated characterization of the planetary system orbiting the nearby M2 dwarf GJ 3090 using radial velocity data from HARPS and NIRPS, combined with TESS photometry and the ΔT spectroscopic activity indicator. Using a multidimensional GP framework tailored to simultaneously model stellar activity across multiple datasets, we derived new mass constraints and provided a more robust interpretation of the system’s architecture.

Our analysis significantly improves the mass measurement of the transiting planet GJ 3090 b, now determined to be 4.52 ± 0.47 M_⊕. Combined with a radius of 2.18 ± 0.06 R_⊕ we derived from four TESS sectors, this yields a new density estimate of $2.{40}_{-0.30}^{+0.33}\mathrm{g}{\mathrm{cm}}^{-3}$, suggesting a volatile-rich composition. It places the planet in a highly degenerate section of the M-R diagram that can be reproduced by different ratios of H₂O and H/He, and brings it closer to the bulk distribution of sub-Neptunes around M dwarfs.

We confirm the existence of a second planet, GJ 3090 c, with a 15.94-day orbit and a minimum mass of 10.0 ± 1.3 M_⊕. Although the signal lies close to the stellar rotation period (P_rot ≈ 18 days), we provided multiple lines of evidence supporting its planetary nature. These include consistency between the HARPS and NIRPS RV signals, the inability of our GP activity model to absorb the 16-day signal in injection-recovery tests, its robust detection across model comparisons using different GP kernels, and the improved correlation between the HARPS RVs and the ΔT activity indicator once the 16-day planetary signal is subtracted.

We also report evidence for a third planet candidate at 12.7 days with a minimum mass of 4.15 ± 1.04 M_⊕, initially proposed by Almenara et al. (2022). While this signal does not reach the same level of statistical significance as GJ 3090 c, it remains a compelling target for continued monitoring as planets of this mass typically lie in the radius valley and could therefore provide insight into its origin.

Importantly, our study illustrates the advantages of multidimensional GP approaches in disentangling stellar activity from planetary signals, particularly for signals near the stellar rotation period. Compared to traditional sequential or semi-shared GP methods, our framework mitigates the risk of overfitting and increases sensitivity to both signals close to the rotation period and long-period signals.

The GJ 3090 system thus emerges as a compact multi-planetary system, possibly composed of three sub-Neptunes, with only the innermost planet observed in transit. Assuming coplanar orbits, we explain the absence of transits for the outer companions through geometric arguments based on their derived impact parameters. No additional significant signals are found in the residuals, but the system remains an interesting candidate for further RV monitoring.

Future observations, including further RV measurements and photometric campaigns, will help confirm the nature of the 12.7-day candidate and probe the presence of additional companions in the system. GJ 3090 b remains a particularly attractive target for atmospheric and interior studies, with its mass and radius now measured to 10% precision or lower.

Acknowledgements

PLam, CC, AL, LMo, RD, FBa, BB, LMa, RA, LB, ÉA, AB, AD-B, LD, NJC, OL, JS-A, PV, TV & JPW acknowledge the financial support of the FRQ-NT through the Centre de recherche en astrophysique du Québec as well as the support from the Trottier Family Foundation and the Trottier Institute for Research on Exoplanets. AL acknowledges support from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies under file #349961. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. AKS, JIGH, RR, ASM, FGT, NN, VMP & JLR acknowledge financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) projects PID2020-117493GB-I00 and PID2023-149982NB-I00. AKS acknowledges financial support from La Caixa Foundation (ID 100010434) under the grant LCF/BQ/DI23/11990071. LMo acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference number 589653]. RD, FBa, LMa, TA, ÉA, J-SM, MO, JS-A & PV acknowledges support from Canada Foundation for Innovation (CFI) program, the Université de Montréal and Université Laval, the Canada Economic Development (CED) program and the Ministere of Economy, Innovation and Energy (MEIE). SCB, ED-M, NCS, EC, ARCS & JGd acknowledge the support from FCT – Fundação para a Ciência e a Tecnologia through national funds by these grants: UIDB/04434/2020, UIDP/04434/2020. SCB acknowledges the support from Fundação para a Ciência e Tecnologia (FCT) in the form of a work contract through the Scientific Employment Incentive program with reference 2023.06687.CEECIND and DOI 10.54499/2023.06687.CEECIND/CP2839/CT0002. XB, XDe, ACar, TF & VY acknowledge funding from the French ANR under contract number ANR24CE493397 (ORVET), and the French National Research Agency in the framework of the Investissements d’Avenir program (ANR-15-IDEX-02), through the funding of the “Origin of Life” project of the Grenoble-Alpes University. The Board of Observational and Instrumental Astronomy (NAOS) at the Federal University of Rio Grande do Norte’s research activities are supported by continuous grants from the Brazilian funding agency CNPq. This study was partially funded by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brasil (CAPES) – Finance Code 001 and the CAPES-Print program. BLCM & AMM acknowledge CAPES postdoctoral fellowships. BLCM acknowledges CNPq research fellowships (Grant No. 305804/2022-7). NBC acknowledges support from an NSERC Discovery Grant, a Canada Research Chair, and an Arthur B. McDonald Fellowship, and thanks the Trottier Space Institute for its financial support and dynamic intellectual environment. DBF acknowledges financial support from the Brazilian agency CNPq-PQ (Grant No. 305566/2021-0). Continuous grants from the Brazilian agency CNPq support the STELLAR TEAM of the Federal University of Ceara’s research activities. JRM acknowledges CNPq research fellowships (Grant No. 308928/2019-9). ED-M further acknowledges the support from FCT through Stimulus FCT contract 2021.01294.CEECIND. ED-M acknowledges the support by the Ramón y Cajal contract RyC2022-035854-I funded by MICIU/AEI/10.13039/501100011033 and by ESF+. XDu 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). DE acknowledge support from the Swiss National Science Foundation for project 200021_200726. The authors acknowledge the financial support of the SNSF. ICL acknowledges CNPq research fellowships (Grant No. 313103/2022-4). CMo acknowledges the funding from the Swiss National Science Foundation under grant 200021_204847 “PlanetsInTime”. 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. GAW is supported by a Discovery Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada. 0 KAM acknowledges support from the Swiss National Science Foundation (SNSF) under the Postdoc Mobility grant P500PT_230225. RA acknowledges the Swiss National Science Foundation (SNSF) support under the Post-Doc Mobility grant P500PT_222212 and the support of the Institut Trottier de Recherche sur les Exoplanètes (IREx). We acknowledge funding from the European Research Council under the ERC Grant Agreement no. 337591-ExTrA. LB acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Spice Dune, grant agreement No. 947634). This material reflects only the authors’ views and the Commission is not liable for any use that may be made of the information contained therein. ARCS acknowledges the support from Fundaçao para a Ciência e a Tecnologia (FCT) through the fellowship 2021.07856.BD. LD acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) and from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies. FG acknowledges support from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies under file #350366. This work was supported by grants from eSSENCE (grant number eSSENCE@LU 9:3), the Swedish National Research Council (project number 2023 05307), The Crafoord foundation and the Royal Physiographic Society of Lund, through The Fund of the Walter Gyllenberg Foundation.” NN acknowledges financial support by Light Bridges S.L, Las Palmas de Gran Canaria. NN acknowledges funding from Light Bridges for the Doctoral Thesis “Habitable Earth-like planets with ESPRESSO and NIRPS”, in cooperation with the Instituto de Astrofísica de Canarias, and the use of Indefeasible Computer Rights (ICR) being commissioned at the ASTRO POC project in the Island of Tenerife, Canary Islands (Spain). The ICR-ASTRONOMY used for his research was provided by Light Bridges in cooperation with Hewlett Packard Enterprise (HPE). CPi acknowledges support from the NSERC Vanier scholarship, and the Trottier Family Foundation. CPi also acknowledges support from the E. Margaret Burbidge Prize Postdoctoral Fellowship from the Brinson Foundation. TV acknowledges support from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies under file #320056. We used the following Python packages for data analysis and visualization: numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), astropy (Astropy Collaboration 2022), juliet (Espinoza et al. 2019), spleaf (Delisle et al. 2022), and NIEVA (Stefanov et al., in prep.). This manuscript was written and compiled in OVERLEAF.

References

Appendix A Supplementary material

Fig. A.1 TESS PDCSAP flux light curves of the four different sectors with the best-fit juliet model for GJ 3090 b transits shown as a black line.

Fig. A.2 Periodograms of the RVs after having subtracted the best-fit offset of each instrument. The first panel shows the periodogram of the RVs of NIRPS and HARPS (HAM & EGGS) separately. The second panel shows the periodogram of the combined RVs, with significant peaks marked for known planets and signals. The subsequent panels show the periodograms after sequentially removing the signals of the magnetic cycle, GJ 3090c, GJ 3090b, the 12.7-day candidate, and finally the GP model. Analytical false alarm probability (FAP) levels of 10%, 1%, and 0.1% are indicated by dashed lines in red, green, and blue, respectively. The window function is shown as a grey dashed line.

Fig. A.3 Top panel: Subsets of the RV time series for HARPS and NIRPS. The left sub-panel shows the GP model fit to the RV data without the Keplerian signals, while the right sub-panel presents the combined Keplerian and GP model. Bottom panel: GP model applied to the activity indicators, specifically the ΔT and TESS data. The shaded regions represent the 1σ confidence intervals of the GP model predictions. We note that the top and bottom panels do not show the same epochs of observation. The top panel shows 135 days of data before the last spectrum and the bottom panel shows up to 195 days after the beginning of first TESS sector.

Fig. A.4 Lomb-Scargle periodogram of the HARPS (HAM) ΔT measurements. The solid blue line shows the power as a function of period, while the horizontal red, green, and blue dashed lines mark the 1%, 0.1%, and 0.01% analytical false alarm probability thresholds respectively. There are no significant peaks around the planets’ periods. There is a forest of statistically significant peaks between 13 and 20 days which shows a sign of potential differential rotation.

Fig. A.5 Correlations between HARPS (HAM) RVs and its ΔT indicator (left) and with the d ΔT/d t (right). The blue data points show the RVs that include the 16-day signal and the orange data points show the RVs after removing the best-fit Keplerian model of the 16-day signal. A linear function is fitted in each dataset. First-order correlation (top): R=−0.135 to R=−0.238 and Nσ=2.71 to Nσ=6.73. Second-order correlation (bottom): R=−0.386 to R=−0.479 and Nσ=17.2 to Nσ=20.6

Fig. A.6 Phase-folded light curves for the two outer planets. The left panel shows the light curve for GJ 3090 c, while the right panel corresponds to the 12.7-day candidate. The data points represent the TESS photometry corrected with the activity model. The red line is a binned version of the data (50 points per bin) so that we see the flux variations more easily. No transit is observed in either case.

Fig. A.7 Corner plot for the H2/He scenario of GJ 3090 b showing the posterior distributions for the planet mass (M), radius (R), atmospheric mass fraction (AMF), core mass fraction (CMF), and refractory element ratios (Fe/Mg and Fe/Si). The red contours correspond to the case with stellar-informed priors on the refractory abundances, while the black contours correspond to the case without such priors. The inclusion of stellar priors significantly reduces the spread in CMF and AMF values by constraining the rocky interior composition.

Fig. A.8 Corner plot for the pure H2O atmosphere scenario of GJ 3090 b showing the posterior distributions for the planet mass (M), radius (R), water mass fraction (WMF), core mass fraction (CMF), and refractory element ratios (Fe/Mg and Fe/Si).

All Tables

All Figures

	Fig. 1 Phase-folded TESS light curves of GJ 3090 b (grey points). The dark green circles are data binned to 10 min. The black lines represent the median model of each sector from the fit.
In the text

Fig. 2 Top panel: joint best-fit sinusoidal model to the ASAS-SN photometric data from all available cameras (bf, bj, bF, and bn). The solid lines represent the best-fit model. Bottom panel: Lomb-Scargle periodogram of the ASAS-SN photometry after camera offset corrections. The blue curve shows the power spectrum. The strongest peak corresponds to the best-fit magnetic cycle period of 3370 ± 120 days (red dashed line), while the orange dashed line indicates the stellar rotation period of 18.0 days and the green dashed line shows the sidereal lunar cycle that arises from the applied lunar correction. The grey dashed curve represents the window function.

In the text

Fig. 3 Phase-folded radial velocity measurements for the confirmed and candidate planets in the GJ 3090 system. The first panel shows the confirmed planet GJ 3090 b (2.85 days), the second panel presents the confirmed sub-Neptune GJ 3090 c (15.94 days), and the third panel displays the 12.7-day candidate planet. The data points from HARPS HAM, HARPS EGGS, and NIRPS are shown with their associated uncertainties, alongside the Keplerian model curves (black lines).

In the text

Fig. 4 Heat map similar to that in S25 showing the comparison of Bayesian evidence (Δ ln Z) for different models. The x-axis displays different stellar-activity kernels in increasing order of complexity, while the y-axis shows different planetary configurations. Circular orbits (e= 0) are labelled ‘C’, while Keplerian orbits (e ≥ 0) are labelled ‘K’. The red rectangle highlights the best model based on the maximum Bayesian evidence.

In the text

Fig. 5 Mass estimates of GJ 3090 b under different planetary configurations and different stellar-activity kernels. The models range from single-planet solutions (GJ 3090 b only) to multi-planet scenarios, including the 16-day confirmed planet (GJ 3090 c) and candidate planet at 13 days. The three-planet models are also shown with different kernels and a model with a non-zero eccentricity on GJ 3090 b. The adopted model is shown in black. The mass constraints from A22, O24, and an independent reanalysis of the 55 archival HARPS data points are also shown. A dashed horizontal line separates the models tested in this work from previously published mass estimates and a re-analysis of the HARPS archival data.

In the text

	Fig. 6 Posterior distributions of the Keplerian parameters of GJ 3090 c derived from fits to the combined dataset (black), HARPS-only data (blue), and NIRPS-only data (red). The distributions show excellent agreement between HARPS and NIRPS, with tighter constraints in the combined dataset due to the increased number of RV points.
In the text

Fig. 7 Chronological evolution of the false alarm probability (FAP) for the 16-day signal. As new observations are added, the FAP in the RVs decreases steadily, consistent with a stable Keplerian signal. In contrast, the FAP in the ΔT and FWHM activity indicators remains high across the datasets. The orange and green horizontal lines represent the 3σ and 5σ detection significance respectively. This behaviour supports the long-term coherence of the GJ 3090 c signal and disfavours a stellar origin.

In the text

Fig. 8 Sensitivity map showing the detection significance (σ) of the multidimensional GP approach as a function of period and semi-amplitude K of injected planets. The colour map indicates the detection significance computed from the periodogram of the residuals. The 5σ detection contours are shown for the serial GP (dotted grey line), the composite or semi-shared GP (dashed grey line), and the multidimensional GP (solid grey line). Known planets are marked with yellow symbols. The vertical orange lines indicate the stellar rotation period and its first harmonic. The 16-day signal is recovered with a detection significance of 4.05σ using the multidimensional GP, 1.10σ using the composite GP and 0σ with the serial GP.

In the text

	Fig. 9 Plot of samples of φ1 (orang), φ2 (blue), and φ12 (green) between GJ 3090 c and the 12.7-day candidate. Samples are taken every 1000 years, across a 1 Myr integration of the GJ 3090 system. Horizontal dashed lines corresponding to 0° and 180° are shown as well. All three angles appear to librate, though this is merely an artefact of the sampling frequency.
In the text

Fig. 10 Plot of samples of φ1 (orange), φ2 (blue), and φ12 (green) between GJ 3090 c and the 12.7-day candidate. Samples are taken 1000 times a year, across the first 10 years of a 1 Myr integration of the GJ 3090 system. Horizontal dashed lines corresponding to 0° and 180° are shown as well. While Figure 9 appears to show all three angles librating, the higher sampling frequency shown here shows that all three angles circulate.

In the text

Fig. 11 Mass-radius diagram. GJ 3090 b (green) is shown alongside known M dwarf exoplanets with <30% uncertainty in both mass and radius (grey points). The red symbol indicates the mass–radius position of GJ 3090 b as previously measured by A22. Theoretical iso-composition lines include purely rocky compositions (black and red), H2O envelope (steam), planets with Earth-like rocky interior (dotted blue), and H/He envelope sub-Neptunes with Earth-like rocky interior (orange). These interiors are retrieved using the detailed modelling described in Sect. 6.8 and assume Te q=700 K with 1 Gyr system age.

In the text

	Fig. A.1 TESS PDCSAP flux light curves of the four different sectors with the best-fit juliet model for GJ 3090 b transits shown as a black line.
In the text

Fig. A.2 Periodograms of the RVs after having subtracted the best-fit offset of each instrument. The first panel shows the periodogram of the RVs of NIRPS and HARPS (HAM & EGGS) separately. The second panel shows the periodogram of the combined RVs, with significant peaks marked for known planets and signals. The subsequent panels show the periodograms after sequentially removing the signals of the magnetic cycle, GJ 3090c, GJ 3090b, the 12.7-day candidate, and finally the GP model. Analytical false alarm probability (FAP) levels of 10%, 1%, and 0.1% are indicated by dashed lines in red, green, and blue, respectively. The window function is shown as a grey dashed line.

In the text

Fig. A.3 Top panel: Subsets of the RV time series for HARPS and NIRPS. The left sub-panel shows the GP model fit to the RV data without the Keplerian signals, while the right sub-panel presents the combined Keplerian and GP model. Bottom panel: GP model applied to the activity indicators, specifically the ΔT and TESS data. The shaded regions represent the 1σ confidence intervals of the GP model predictions. We note that the top and bottom panels do not show the same epochs of observation. The top panel shows 135 days of data before the last spectrum and the bottom panel shows up to 195 days after the beginning of first TESS sector.

In the text

Fig. A.4 Lomb-Scargle periodogram of the HARPS (HAM) ΔT measurements. The solid blue line shows the power as a function of period, while the horizontal red, green, and blue dashed lines mark the 1%, 0.1%, and 0.01% analytical false alarm probability thresholds respectively. There are no significant peaks around the planets’ periods. There is a forest of statistically significant peaks between 13 and 20 days which shows a sign of potential differential rotation.

In the text

Fig. A.5 Correlations between HARPS (HAM) RVs and its ΔT indicator (left) and with the d ΔT/d t (right). The blue data points show the RVs that include the 16-day signal and the orange data points show the RVs after removing the best-fit Keplerian model of the 16-day signal. A linear function is fitted in each dataset. First-order correlation (top): R=−0.135 to R=−0.238 and Nσ=2.71 to Nσ=6.73. Second-order correlation (bottom): R=−0.386 to R=−0.479 and Nσ=17.2 to Nσ=20.6

In the text

Fig. A.6 Phase-folded light curves for the two outer planets. The left panel shows the light curve for GJ 3090 c, while the right panel corresponds to the 12.7-day candidate. The data points represent the TESS photometry corrected with the activity model. The red line is a binned version of the data (50 points per bin) so that we see the flux variations more easily. No transit is observed in either case.

In the text

Fig. A.7 Corner plot for the H2/He scenario of GJ 3090 b showing the posterior distributions for the planet mass (M), radius (R), atmospheric mass fraction (AMF), core mass fraction (CMF), and refractory element ratios (Fe/Mg and Fe/Si). The red contours correspond to the case with stellar-informed priors on the refractory abundances, while the black contours correspond to the case without such priors. The inclusion of stellar priors significantly reduces the spread in CMF and AMF values by constraining the rocky interior composition.

In the text

	Fig. A.8 Corner plot for the pure H2O atmosphere scenario of GJ 3090 b showing the posterior distributions for the planet mass (M), radius (R), water mass fraction (WMF), core mass fraction (CMF), and refractory element ratios (Fe/Mg and Fe/Si).
In the text