© The Authors 2025

1 Introduction

The detection of low-mass exoplanets has been fueled by technological and methodological advancements. The arrival of high-precision, ultra-stabilised spectrographs, opened the possibility of detecting low-mass exoplanets (Santos et al. 2004; Mayor et al. 2009). Later, improvements in RV determination yielded better RV precision than what was originally expected to be achievable (Anglada-Escudé et al. 2016). Stellar activity is nowadays the main limiting factor when attempting to detect planets that induce m s⁻¹ RV signals (Robertson et al. 2014; Suárez Mascareño et al. 2017b). In recent years, the techniques used to correct stellar activity have progressed significantly, enabling the detection of planets inducing signals smaller than those of stellar activity (Haywood et al. 2014; Faria et al. 2022). Simultaneously, new statistical methods have been developed to optimally incorporate the stellar variability model in the computation of the statistical significance of candidate planets (Hara et al. 2022b). With all these tools at hand, it is now easier to identify the presence of very low-mass planets, and asses their planetary nature. A consequence of these improvements in radial velocity extraction and modelling tools is that it is now possible to find planets hidden in archival data.

GJ 536 is a nearby M1V dwarf that hosts a close-in low-mass planet (Suárez Mascareño et al. 2017a). The planet was identified using a traditional periodogram analysis on the HARPS and HARPS-N RV data, derived from cross-correlation functions (CCFs). The rotation period of the star was measured as 43 days, and modelled using a double sinusoidal model. At the time, no other signals were identified in the data.

In this paper, we report the detection of a second low-mass planet orbiting inner to the edge of the optimistic habitable zone of the star. We perform a re-extraction of the RVs using a Line-By-Line algorithm, we complete the RV data using the publicly available CARMENES and HIRES data, perform a global model using multidimensional Gaussian processes (GP) regression, and asses the detections using state-of-the-art statistical models.

2 Observations and data

This work is based on archival observations from the High Accuracy Radial velocity Planet Searcher (HARPS, Mayor et al. 2003), HARPS-N Cosentino et al. (2012), the Calar Alto high-Resolution search for M-dwarfs with Exoearths with Near-infrared and optical Échelle Spectrographs (CARMENES, Quirrenbach et al. 2014), and The High Resolution Echelle Spectrometer (HIRES, Vogt et al. 1994).

HARPS and HARPS-N are fibre-fed high-resolution echelle spectrographs with similar resolving power and spectral range, R ~ 115000 between ~380 and ~690 nm and have been designed to attain very high long-term radial-velocity precision. HARPS is installed at the 3.6-m ESO telescope in La Silla Observatory (Atacama, Chile), while HARPS-N is installed at the 3.6-m Telescopio Nazionale Galileo of the Roque de los Muchachos Observatory (La Palma, Spain). Both are contained in temperature- and pressure-controlled vacuum vessels to avoid spectral drifts due to temperature and air pressure variations, thus ensuring its stability. They are equipped with their own pipeline providing extracted and wavelength-calibrated spectra, as well as RV measurements and other data products such as cross-correlation functions and their bisector profiles.

CARMENES consists of visual (VIS) and near-infrared (NIR) vacuum-stabilised spectrographs covering 520–960 nm and 960–1710 nm with a spectral resolution of 94 600 and 80 400, respectively. It is located at the 3.5-m Zeiss telescope at the Centro Astronómico Hispano Alemán (Almería, Spain). The spectra are extracted with the CARACAL pipeline, based on flat-relative optimal extraction (Zechmeister et al. 2014). The wavelength calibration was performed by combining hollow cathode (U-Ar, U-Ne, and Th-Ne) and Fabry-Pérot etalons exposures. The instrument drift during the nights is tracked with the Fabry-Pérot in the simultaneous calibration fibre.

HIRES is a slit-fed, grating cross-dispersed echelle spectrograph with a spectral resolution of 55–80 000 depending on the input slit decker (Vogt et al. 1994). HIRES is installed at the Keck I telescope on Mauna Kea (Hawaii, USA) and uses an iodine cell placed in the converging beam of the telescope as a wavelength calibrator (Marcy & Butler 1992). This cell imprints incoming stellar light with a dense forest of I2 lines from 500–620 nm, providing a precise wavelength solution in this region. The I2 lines also act as a proxy for the point spread function (PSF) of the instrument, and so changes to the instrument profile caused by temperature or pressure variations will be reflected in the observed absorption line profiles (Butler et al. 1996). While HIRES has a spectral range of 370–800 nm, only this 120 nm iodine region is used when measuring precise radial velocities.

In addition to the spectroscopic data, we used the photometric time series obtained by the All Sky Automated Survey (ASAS; Pojmanski 1997). The ASAS data is composed of V-band observations taken by the ASAS-3 survey, obtained from the ASAS website¹. The light curves include the photometric measurements obtained using four different apertures. We averaged over the four apertures to obtain the final photometric series, reducing the short-term scatter of the data.

Last, GJ 536 has been observed by the All Sky Automated Survey for SuperNovae (ASAS-SN; Kochanek et al. 2017) and the Transiting Exoplanet Survey Satellite (TESS) (Ricker et al. 2015).The default aperture photometry data of ASAS-SN were strongly contaminated by the moon. We attempted to clean them by using the saturated stars photometry method (Winecki & Kochanek 2024) and removing epochs with lunar separations of less than 90 degrees. The resulting time series showed no evidence of astrophysical variations. TESS observed GJ 536 in sector 50. We found no evidence of transits in the TESS data (see Appendix F).

2.1 Radial velocities

We extracted the HARPS and HARPS-N RVs using the Line-By-Line (LBL²) code developed by Artigau et al. (2022), and based on Dumusque (2018). The LBL algorithm performs an outlier-resistant template matching to each individual line in the spectra. For non-telluric corrected spectra, it produces its own telluric correction. In addition to the velocity, it derives other quantities such as a differential line-width (dLW, similar to Zechmeister et al. 2018). The dLW is obtained from the second derivative of the template and can be understood as a change in the line full width at half maximum (FWHM), assuming a Gaussian profile. Using the dLW variations, and an estimation of the average FWHM of the lines, the LBL algorithm estimates a variation of the FWHM of the lines. We used these measured changes in the FWHM as our main activity indicator. In addition, the LBL algorithm derives a chromatic RV slope over the full spectral range (chromatic index, CRX, similar to Zechmeister et al. 2018). In addition, we used the Barycentric Earth Radial Velocity (BERV) to track potential leftover effects from the telluric correction.

The time series data of CARMENES were taken from the CARMENES DR1 (Ribas et al. 2023). These RVs were obtained using SERVAL. This software builds a high signal-to-noise template by co-adding all the existing observations, and then performs a maximum likelihood fit of each observed spectrum against the template, yielding a measure of the Doppler shifts and their uncertainty. The time series provided includes additional indicators, such as the FWHM of a CCF built using the RACOON pipeline (Lafarga et al. 2020).

We used the HIRES RV time series provided by Rosenthal et al. (2021), which, in addition, included time series of Mount Wilson S -index. The HIRES RVs were obtained by measuring the Doppler shift of starlight relative to a reference spectrum of molecular iodine, which is at rest in the observatory frame (Butler et al. 1996).

2.2 Activity proxies

The presence of active regions on the stellar surface affects the flux emitted by the star and its velocity field, distorting the shape of the lines and the point spread function measured by the instrument (Dravins 1982). To measure these variations, we use the full width at half maximum (FWHM), provided by the LBL code and provided by the CARMENES DR1.

The intensity of the emission of the cores of the Ca II H&K lines is linked to the strength of the magnetic field of the star, which in turn is well correlated with the stellar rotation period for FGKM stars (Noyes et al. 1984). The measured emission intensity of the line cores also changes when active regions move across the stellar disc, helping us trace the rotation of the star. We calculate the Mount Wilson S index (S _MW) for the HARPS and HARPS-N data following Vaughan et al. (1978). We define two triangular-shaped passbands with a FWHM of 1.09 Å centred at 3968.470 Å and at 3933.664 Å for the Ca II H&K line cores. For the continuum, we use two bands 20 Å in width centred at 3901.070 Å (V) and 4001.070 Å (R). Then the Snindex is defined as: $$S=\alpha \frac{{\tilde{N}}_H+{\tilde{N}}_K}{{\tilde{N}}_R+{\tilde{N}}_V}+\beta ,$$(1)

where ${\tilde{N}}_H$, ${\tilde{N}}_K$, ${\tilde{N}}_R$, and ${\tilde{N}}_V$ are the mean fluxes per wavelength unit in each passband, while α and β are calibration constants fixed at α = 1.111 and β = 0.0153 (Lovis et al. 2011). The S index serves as a measurement of the Ca II H&K core flux normalised to the neighbour continuum.

The HIRES S_MW measurements from Rosenthal et al. (2021) are provided without estimates for their uncertainties. We downloaded the original spectra and computed the photon noise-uncertainties by propagating the photon noise in each band, as with the HARPS and HARPS-N data (e.g. Laliotis et al. 2023).

We complement the data with the photometric time series coming from the All Sky Automated Survey (ASAS; Pojmanski 1997). We used the observations obtained by the ASAS-3. It consists of V-band data taken over a long baseline. We obtained the ASAS-3 data from the ASAS website³. The light curves include the photometric measurements obtained using four different apertures. We averaged over the four apertures to obtain the final photometric series, which reduced the short-term scatter of the data.

2.3 Summary of time series

Overall, we analysed 235 nightly binned RV measurements, distributed as 147 from HARPS, 5 from HARPS-N, 25 from CARMENES, and 58 from HIRES. The combined RV time series shows an RMS of 3.55 m s⁻¹, and a median uncertainty of 81 cm s⁻¹. We analysed 177 FWHM measurements from HARPS, HARPS-N, and CARMENES. The combined FWHM time series shows an RMS of 6.68 m s⁻¹ with a median uncertainty of 81 cm s⁻¹. We analysed 210 S_MW measurements from HARPS, HARPS-N, and HIRES. The data shows an RMS of 0.983 with a median uncertainty of 0.038. Last, we analysed 325 photometric epochs, with an RMS of 13.1 parts-per-thousand (ppt) and a median uncertainty of 8.3 ppt. Figure 1 shows the time series of RV, S_MW index, FWHM, and V-fluxes. Figure 1 shows the time series of RV, S_MW index, FWHM, and V-fluxes. Table 1 shows a summary of the spectroscopic data considered for the analysis, with details for every individual instrument. We found the RV data to show weak levels of core-lationship with the S_MW and FWHM and data. More details can be found in Appendix A.

3 Fundamental stellar parameters of GJ 536

We derived the fundamental stellar parameters of GJ 536 using the Bayesian inference code of del Burgo & Allende Prieto (2016, 2018) on the PARSEC v1.2S library of stellar evolution models (Bressan et al. 2012; Chen et al. 2014; Tang et al. 2014; Chen et al. 2015), which shows good statistical agreement with the dynamical masses of main-sequence detached eclipsing binaries (within 4%, del Burgo & Allende Prieto 2018).

We arranged a grid of stellar evolution models with ages ranging from 2 to 13 800 million years and steps of 5%, and [M/H] from −2.18 to 0.51 with steps of 0.02 dex, adopting the photometric passband calibration of Riello et al. (2021), with the zero points of the VEGAMAG system. We then applied the Bayesian method described in del Burgo & Allende Prieto (2016), taking as imputs the absoltute G magnitude and the colour G_BP – G_RP, the absolute G magnitude M_G from Gaia DR3 (Gaia Collaboration 2021), assuming a null extinction, and the abundance ratio of iron to hydrogen [Fe/H]=−0.08 ± 0.09 (Maldonado et al. 2020). With this method, we inferred the stellar mass, radius, luminosity, effective temperature, and surface gravity. It is important to note that the uncertainties of derived stellar masses and radii are typically smaller than the dispersion between different models. As these parameters can affect the computation of planetary parameters, we enlarged the error bars of the stellar parameters in all our calculations following Tayar et al. (2022). The final uncertainty in the stellar mass of 4% is consistent with the measured statistical accuracy of the code. Table 2 shows the inferred parameters of GJ 536.

We measured a rotational period of 43.63 ± 0.85 d, which can be used for dating GJ 536 from the gyrochronology-based relation of Mamajek & Hillenbrand (2008), taking into account the colour B-V= 1.47 of this red dwarf (Koen et al. 2010). The resulting gyro-kinematic age is 4.2 ± 0.7 Ga, which is lower but nearly within 1-σ when compared to the isochronal age, which matches that of Maldonado et al. (2020) (8±4 Ga).

GJ 536 has a magnetic cycle originally reported by Suárez Mascareño et al. (2017a), with a period of ~825 days, and recently updated by Ibañez Bustos et al. (2025) to a period of ~1600 days, with a second component at longer period (~4000 days).

Using the reported stellar parameters, we estimate the limits of the habitable zone (HZ) of the star to be 0.176–0.372 au, following Kopparapu et al. (2014, 2017), for the runaway green-house to early-Mars limits. These limits correspond to orbital periods of 37.1 ± 1.1 days to 114.2 ± 3.4 days.

Fig. 1 Spectroscopic and photometric data. Time series of RV, SMW, FWHM, and FluxV used in this study.

4 Analysis

We performed a global analysis of the data, always including the full set of RVs, S_MW, FWHM, and photometric data. Every time series includes a zero point per instrument, a model for the cycle, and a model for the stellar rotation. The cycle model is defined as a set of sinusoidal signals, with a common period and phase for all time series, and independent amplitudes. The stellar rotation is modelled as a multi-dimensional Gaussian Process (Rajpaul et al. 2015), using the S+LEAF code (Delisle et al. 2022) ⁴. The RV data of HARPS and HARPS-N includes polynomial terms against BERV. In addition, the RV data includes a circular/Keplerian model for each planet in the model. A more detailed description of the model can be found in Appendix B.

We optimised the parameters of the models using Bayesian inference through the nested sampling (Skilling 2004, 2006) code Dynesty (Speagle 2020; Koposov et al. 2023). We sampled the parameter space using random slice sampling, which is well suited for the high-dimensional spaces (Handley et al. 2015a,b) resulting from modelling several time series at once. We used a number of live points of 20 × N_free in models with narrow priors in period, and 100 × N_free in models with wide priors in frequency, to ensure the discoverability of the narrow frequency posteriors.

We assessed the significance of the detection of the signals using the False inclusion probability (FIP) framework (Hara et al. 2022b). This method uses the posterior distribution of the nested sampling run, and computes the probability of having a planet within a certain orbital frequency interval based on the relative weight of the samples within each frequency interval. The FIP framework has been demonstrated to minimise both false detections and missed detections, compared to other detection criteria such as the false alarm probability, or the Bayesian evidence. Based on an extensive set of simulations, Hara et al. (2022b, 2024) suggested that a FIP threshold of 50% maximises the number of detections while thresholds between 10% and 1%, minimise the number of mistakes (sum of false positives and false negatives). We adopt a threshold of 1% to consider the significant detection of a signal. This method has already been shown to be effective in the detection and confirmation of exo-planets from radial velocity time series (Suárez Mascareño et al. 2023).

5 Results

We performed a blind search for planetary signals under the FIP formalism. We run a model with up to three planetary signals with wide priors and assessed the significance of the detection based on the posterior distribution of the frequency parameters. Fig. 2 shows the FIP periodogram of the data. The periodogram shows two signals with a FIP <0.1%. The first signal, with a period of 8.70877 ± 0.00059 days and an RV semi-amplitude of 3.03 ± 0.15 m s⁻¹, corresponds to the planet GJ 536 b. The second, with a period of 32.761 ± 0.015 days and an RV semi-amplitude of 1.79 + 0.22 m s⁻¹, is a new signal not previously known. There are no hints of additional signals. We cross-checked the results measuring the difference in Bayesian evidence (Ln Z) between the models with no planets, 1 planet, and 2 planets. We measure a Δ Ln Z of 53.9 favouring the 1-planet model against the no-planet model, and Δ Ln Z of 30 favouring the 2-planet model over the 1-planet model. These results strongly favour the 2-planet model over the 1- or 0-planet models.

As the FIP analysis already marginalises the probability to have a planet over the activity model, it is assumed that the detected signals are likely not of stellar origin. At the very least, they are signals not accounted for the activity model, or any of the remaining components defined in the model. We measure the stellar rotation period to be 43.63 ± 0.90 days, with a timescale of coherence of the signal of ${78}_{-24}^{+38}$ days. We measure the stellar cycle period to be ${3390}_{-61}^{+120}$ days.

Fig. 2 FIP periodogram of the full dataset. The highlighted peaks corresponds to the period of 8.7d (GJ 536 b), and 32.8 d.

5.1 Presence of the signals in activity indicators

Despite having implemented an activity model in the planet search, it is possible that some activity signals might remain in the residuals of the activity model, which could lead to them being misidentified as planets. In this case, we would expect these signals to also show up in the activity indicators data. While some of their properties (amplitude, phase) would not need to match those of the RV signals, we would expect the signals to show the same periodicity. To test whether some activity signals remain at the periods of the candidate RV signals, we run a model that included two sinusoidal functions in each activity indicator with priors in common with the RV signals. To specifically test the regions around the detected signals, we run a guided search with normal (𝒩) priors centred at the determined periods, with a width of 10% of the period. These values correspond to periods 𝒩[8.7,0.9] days and 𝒩[32.8, 3.3] days. The priors for the phase and amplitude remain the same.

Figure 3 shows the posterior distributions obtained for the RV, S_MW, FWHM, and Flux data. We only obtain well-defined posteriors for the RV data. In the case of the three activity proxies, we obtain amplitudes that are 1–1.5σ consistent with zero, and the posteriors of the period and phase are very similar to the priors. No signal appears to remain in the activity indicators with similar properties as those found in the RV data.

5.2 Aliasing of activity signals

The gaps in observability of the target stars, along with the difficulties of maintaining continuous surveys over years, typically lead to poorly sampled time series prone to aliasing problems. Poorly sampled periodic signals will not only appear at their natural frequency, but at alias frequencies due to the combination of the sampling frequency and the signal frequency. To double-check the extent in which the stellar rotation signal could cause signals that might be mistaken for any of the other two RV signals detected, we computed the window function of the data. We then produced the periodograms of the rotation signal in RV, and evaluated its main periodogram peaks and their aliases.

The window function of our dataset has three main features: a peak at 1 day, another at 361 days, and one at 2791 days. The periodogram of the RV signal has its main activity periodicity at 43.7 days, another at 45.9 days, and then two harmonics of the 43.7d signal (21.85d and 14.55d). All those structures are presented as a forest of peaks, with the 361d aliases explaining the majority of those peaks. Figure C.1 shows the window function, the periodogram of the rotation-induced signals, and the position of their aliases. The periods of the detected RV signals do not coincide with any of the rotation-induced peaks, their alises, or the structures in the window function.

5.3 Stability of the detected signals

We test the stability of the signals over time using apodised signals (Gregory 2016; Hara et al. 2022a). In this test, the suspected planetary signals are multiplied by Gaussian functions with potential values of their centre (µ) over the full baseline of observations, and a large range of potential widths (σ). As planetary signals are stable over time, we would expect to retrieve an undefined µ and a large σ. Retrieving very defined parameters for the Gaussian model would imply the signal is concentrated in a specific part of the baseline. We ran this test using the same narrow priors for the periods described before (𝒩[P, 10% P])

Figure 4, top panels, shows the measured evolution of the amplitude, using the posterior of the apodised semi-amplitude at each point in time. The distribution of amplitudes of both signals is consistent with the values obtained in the static model. The range is wider at the beginning, coinciding with the data with the lowest observing cadence. In both cases we obtain that the centre is consistent with any point of the baseline, and the σ is > 120 000 days.

We confirm this stability by measuring the parameters of both signals independently on the first and second half of the dataset. We split the data in two halves with equal amount of RV measurements, and ran a guided search with the same priors for the periods. Figure 4, bottom panels, shows the RV semi-amplitude versus the time of inferior conjunction of the two signals for the first and second half of the dataset, using the 95% confidence region of the parameters. We obtain consistent parameters (within 1 σ) for both datasets. The parameters recovered in the first half of the data are less precise, which is consistent with the behaviour seen in the apodised test, due to the early data (HIRES) being less precise and having a lower cadence of observations.

5.4 Adopted model

Following the results of Sections 5.1, 5.2 and 5.3, we conclude that the signal with period 32.8 days is of planetary origin (GJ 536 c). We adopt a model that includes two sinusoidal signals to account for planets GJ 536 b and c. To establish the final parameters of these signals, we ran a model using 𝒰 priors for the periods, centred at the known period and with a range of ± 40% of the period. The rest of the priors for all parameters remain the same. We compared circular and Keplerian models and the results were consistent with each other. We measure upper limits to the eccentricities of 0.1 for GJ 536 b and 0.35 for GJ 536 c (95% confidence). The rest of the parameters were consistent within 1σ. The circular model was statistically favoured, with a Δ In Z = + 3.3 in favour of the simpler model. We adopt the circular model. The full set of priors and parameters of the adopted model is available in Appendix E.

GJ 536 b has an orbital period of 8.70874 ± 0.00056 days and induce an RV semi-amplitude of 3.03 ± 0.15 m s⁻¹. We measure a phase of 0.208 ± 0.014, corresponding to T_c = 10553.10 ± 0.12 (BJD − 2450 000 days). GJ 536 c has an orbital period of 32.761 ± 0.015 days, and induces an RV semi-amplitude of 1.80 ± 0.20 m s⁻¹. We measure a phase of 0.614 ± 0.029, which corresponds to T_c = 10547.35 ± 0.96 (BJD − 2450 000 days).

We measure the stellar cycle to have a period of ${3387}_{-62}^{+110}$ days, with most of its RV-induced signal manifesting at one quarter the cycle period. It induces an RMS of 0.9 m s⁻¹ in RV. We measured the stellar rotation period to be 43.63 ± 0.85 days, with a timescale of evolution of ${80}_{-25}^{+39}$ days. We measure RV jitters of 0.69 ± 0.20 m s⁻¹ for the HARPS-03 data, 1.42 ± 0.40 m s⁻¹ for the HARPS-15 data, 1.57 ± 0.95 m s⁻¹ for the HARPS-N data, 0.75 ± 0.44 m s⁻¹ for the CARMENES data, and 2.00 ± 0.43 for the HIRES data.

Figure 5 shows the full RV dataset, along with the best-fit model, the residuals, and their respective periodograms. The model appropriately describes the data, leaving behind a residual RMS of 1.30 m s⁻¹. The periodogram of the raw dataset has a single dominant peak, at the period of GJ 536 b. No peaks at the periods of the stellar rotation, or GJ 536 c, appear. Both signals having been significantly detected in RV in our global model. This highlights the limitations of the use of periodograms to either identify signals in the data, or asses their significance. Methods that explore the full likelihood space, accounting for all potential planetary signals simultaneously, along with the variations due to stellar activity, are more efficient and robust at identifying low-amplitude planetary signals (Dumusque et al. 2017; Hara et al. 2022b; Hara & Ford 2023). Figure 6 shows zoomed-in views of specific observing campaigns and Figure 7 shows the phase-folded plots of the planetary-induced RV signals. Table E.1 shows the full set of priors, and the measured parameters of the adopted model.

Figure D.1 shows the data, periodograms, and the best-fit model, of the activity proxies used in the analyses. The cycle variations are easy to visually identify in the S_MW and FWHM time series. The periodograms both show peaks at different components of the cycle model; S_MW at P_cyc and P_cyc/2, and FWHM at P_cyc/2. The periodograms of both feature significant peaks at the rotation period. The periodogram of the photometry features a non-significant peak at the rotation period. Figure D.2 shows zoomed-in views of specific observing campaigns.

Fig. 3 Comparison of the posterior of RV and activity indicators. Posterior distributions of the signal parameters of the RV signals compared with the result of performing the same model in the activity indicators.

Fig. 4 Stability of the signals. The top panels show the evolution of the amplitude of the apodised signals. The shaded red area shows the confidence intervals of the amplitudes of the apodised signals over time. The horizontal black lines show the median value and 1-σ range of the static model. The bottom panels show the RV semi-amplitude vs time of inferior conjunction for the solutions obtained with the first half (H1) and second half (H2) of the dataset. The lines encapsulate the 95% confidence interval.

Fig. 5 Adopted RV model. The top panels show the RV data (detrended from BERV), with the best model fit, along with the periodogram of the data. The bottom panels shows the residuals after the fit, along with their periodogram. Figure D.1 shows the model of the activity proxies.

5.5 Estimating the inclination of the system

In addition to tradditional data driven GP kernels, the S+LEAF package includes a physics-based GP model to represent stellar activity (FENRIR, Hara & Delisle 2025), which is able to perform statistical Doppler imaging. As opposed to traditional Doppler imaging, which can retrieve a specific intensity map of the stellar surface at a given time, statistical Doppler imaging consists in retrieving the average properties of the star, in particular the projected obliquity of the stellar rotation axis and the average latitude of the active regions. The models take into account the photometric RV effect, as well as the inhibition of convective blueshift. Hara & Delisle (2025) shows that the FENRIR model is capable of estimating the inclination of the rotation axis, and the latitudes of active regions, of the Sun, using photometry and RV data.

We ran a model that included the two planets and a FEN-RIR model for photometry (ASAS) and radial velocity, with the same priors for the planetary signals in RV as in the adopted model. We run spot-only, faculae-only, and combined models. In all cases, we used the SpotsOppLatMag version of the model, which assumes that spots appear at the two same latitutdes, symmetric with respect to the equator of the star. We also assumed that spots appear preferentially at opposite latitudes. As we do not find a significant contribution of the magnetic cycle in either photometry or RV, we did not include a magnetic cycle model.

The spot-only model was the only one that produced consistent run-to-run results. Both faculae-only and combined models resulted in large degeneracies between the parameters. This is consistent with the expectations that activity variations in M-dwarfs is mostly spot-dominated. Using this formalism, we derive a rotation period (43.66 ± 0.69 d) and timescale of evolution (${66}_{-26}^{+42}$ d) consistent with the adopted model. The resulted ativity-induced RV signal is modelled as a combination of the photometric effect (5.4 ± 1.2 m s⁻¹) and the inhibition of convection (1.75 ± 0.40 m s⁻¹). The inclination of the stellar rotation axis is estimated to be ${58}_{-19}^{+16}$ degrees with respect to our line of sight, with active regions consistent with being located at latitudes of 36 ± 18 degrees. This is consistent with the spot latitude being roughly aligned with our line of sight. These results must be taken with caution, as statistical Doppler imaging is a new technique which needs further validation.

We analysed the TESS light curve of GJ 536 looking for signs of transits, finding no such features. The lack of transits is not surprising, if we assume that the measured inclination is not too far from the real inclination of the rotaton axis. The TESS light curve, and methods used, can be found in Appendix F.

6 Discussion

We obtained significant detections of two planetary signals in RV. These signals correspond to the already known planet GJ 536 b (Suárez Mascareño et al. 2017a), and a new confirmed planet GJ 536 c. GJ 536 b is a planet with a minimum mass of 6.37 ± 0.38 M_⊕, an orbital period of 8.70874 ± 0.00056 days, and an eccentricity lower than 0.1 (95% confidence). It orbits at a distance of 0.0668 ± 0.0012 au and receives around 10 times the Earth’s insolation. GJ 536 c has a minimum mass of 5.89 ± 0.70 M_⊕, an orbital period of 32.761 ± 0.015 days, and an eccentricity lower than 0.35 (95% confidence). It orbits at a distance of 0.1617 ± 0.0028 au, and receives an insolation slightly lower than that received by Venus, making it most likely too hot to be habitable. Table 3 shows the parameters of the planets of the system of GJ 536, using the adopted (circular) model. Figure 8 shows the position of the planets in the context of known low-mass exoplanets, as a function of orbital period and stellar insolation. GJ 536 c is among the group of lowest-mass exoplanets known at moderate orbital periods, and one with comparatively low insolation.

Planets GJ 536 b and c have very similar minimum masses (6.37 ± 0.38 M_⊕ and 5.89 ± 0.70 M_⊕, respectively). Their nature would be extremely dependent of the inclination of their orbits with respect to the plane of the sky. For large inclinations (i > 45°), they will most likely be massive super-Earths, water-worlds, or mini-Neptunes. With inclinations <45°, they will most likely be mini-Neptunes or even Neptune-like planets. For them to be gas giants, it would be necessary that the inclination of the orbital plane is <5°. If we assume the inclination of the rotation axis derived using the FENRIR model, and the orbital plane to be coplanar with the rotation axis, their masses would be ${7.6}_{-1.0}^{+2.3}{\text{M}}_{\oplus }$ and ${7.1}_{-1.2}^{+2.2}{\text{M}}_{\oplus }$, respectively, making them most likely super-Earths or water-worlds. This estimation, however, needs to be taken with caution, given the number of assumptions required for the computation.

Our current measure of the mass of GJ 536 b (6.37 ± 0.38 M_⊕) is slightly larger than the measurement of Suárez Mascareño et al. (2017a) (5.36 ± 0.69 M_⊕). The difference comes from a larger measured RV amplitude of 3.03 ± 0.15 m s⁻¹, versus the original 2.60 ± 0.33 m s⁻¹. For the current analysis, we used a larger dataset, an RV extraction method for the HARPS data more suited for M-dwarf spectra, and a more sophisticated model for stellar activity. The limitations of the original study likely lead to a slight underestimation of the RV amplitude.

Fig. 6 Zoom to selected observing campaigns. RV data of with the best model fit.

Fig. 7 Phase-folded plots of the planetary-induced RV signals. RV variations induced by GJ 536 b and c with the best model fit (top panels), and the residuals after the fit (bottom panels).

Fig. 8 GJ 536 b and c in context. The upper panel shows the masses of known planets, for masses below 100 M⊕, and periods shorter than 1000 days, with the planets of the system of GJ 536 highlighted, and the solar system planets as reference. The middle panel shows the same, but as a function of insolation. The lower panel shows the contrast in reflected light of low-mass planets orbiting bright stars (my < 12), compared to their angular separation to their parent star. The dotted vertical lines indicate the limits in which ANDES will be able to observe (Palle et al. 2025). The red symbols for the planets of GJ 536 show the potential range of contrasts. The green squares show the position of the ANDES golden sample.

6.1 Direct imaging characterisation

The angular separation of GJ 536 b and c is of 6.41 ± 0.11 mas, and 15.51 ± 0.26 mas, respectively. Assuming their radii are between 1.5 and 2.5 R_⊕ (as most planets of similar masses), and albedos of 0.3, we estimate a F_p/F_S contrast in reflected light between 8.7 ⋅ 10⁻⁸ and 2.4 ⋅ 10⁻⁷ for planet b, and between 1.8 ⋅ 10⁻⁸ and 1.5 ⋅ 10⁻⁸ for planet c. While GJ 536 b is too close to its star, GJ 536 c is a good candidate for direct imaging with future giant telescopes, such as the ELT with ANDES Marconi et al. (2022), and missions such as LIFE Quanz et al. (2022). Figure 8, lower panel, shows their position in comparison with other low-mass planets (m_p < 100 M_⊕) orbiting bright stars (V < 12). GJ 536 c is among the 20–30 most favourable targets for atmosphere characterization via direct imaging, depending on the specific radius assumption and magnitude cut.

6.2 Limits on additional companions

Using the results from the adopted model, we measured the compatibility limits for an additional planetary signal at a wide range of orbital periods. We froze most of the parameters of the model (trends, cycle, GP, planets) and left free only the white noise and zero point RV components. We included a third sinusoid in the model. We built a grid of 1000 bins over periods 1-1000 days with equal width in log₁₀ space. We ran this model with a narrow prior for the period within each bin. From the posterior distribution of each try, we computed the 1% and 99% limits in RV amplitude. Figure 9 shows the result of this exercise. The shaded area encapsulates this range. Any potential signal left in the data should show as a deviation from zero in RV amplitude (or mass). We found that, for periods shorter than 10 days, we can exclude the presence of any additional signal with an amplitude larger than 50 cm s⁻¹ (minimum masses 0.6–1.5 M_⊕). Within the habitable zone, we can exclude, for the most part, the presence of signals with amplitudes between 0.5 and 1.2 m s⁻¹ (m sin i 2–4 M_⊕). At periods between the edge of the habitable zone (114 d) and 1000 days, we can exclude the presence of additional signals beyond 0.5–1.2 m s⁻¹ (m sin i of 2–10 M_⊕).

There is one period bin in which the measured RV amplitude is not consistent with zero. At 44 days we measure an RV semi-amplitude of 76 ± 22 cm s⁻¹. This signal did not appear in our global model optimisation and the period is consistent with the measured stellar rotation. We believe it to be a leftover of stellar rotation due to having fixed the activity parameters rather than optimising them.

Fig. 9 Compatibility limits. The upper panel shows the RV amplitude limit (1–99% area) as a function of orbital period. The lower panels show the same, but for planetary masses. The blue rectangles show the 3σ amplitude/mass range of the detected planets.

6.3 Stellar activity of GJ 536

Along with the planetary system, we characterised the activity variations of GJ 536. We analysed the time series of the V-mag flux, Ca II H&K lines, the RV, and the FWHM of the stellar lines simultaneously with a model that combined a cycle component and Gaussian processes.

We measured a cycle period of ${3387}_{-62}^{+110}$ days (${9.27}_{-0.17}^{+0.30}$ years). The period is in the ballpark of the longer period reported by Ibañez Bustos et al. (2025). The shape of the magnetic cycle is significantly different from a sinusoidal. Instead, it is better described by a combination of several, at the natural period and its harmonics. This shape of the cycle explains the period measurements at ~1700 days and ~850 days estimated in previous works (Suárez Mascareño et al. 2017a; Ibañez Bustos et al. 2025). We found that the main variability of the cycle did not manifest at the same period for all time series. The S_MW data showed variations at the main period, half, and one quarter. The FWHM data only at half the period, and the RV data only at one quarter. Figure 10 shows the comparison between the cycle models in the three time series. The cycle creates a variation in SMW of 0.05 RMS (total RMS 0.1), of 3.2 m s⁻¹ RMS in FWHM (total RMS 6.9 m s⁻¹), and of 1.0 m s⁻¹ RMS in RV (total RMS 3.5 m s⁻¹).

Fully disentangling whether we see a complex cycle, which different proxies track in different ways, or the real cycle is shorter (~1700 days) and the S_MW index data is affected by some additional source of noise is difficult with the data at hand. We compared models with longer and complex cycles, against shorter and simpler cycles, and tried to evaluate the best using the difference in Bayesian evidence (∆ lnZ). The model featured in the analysis (long period with power at P_cyc/2 and P_cyc/4) was favoured if we used all time series simultaneously (RV, S_MW, FWHM, and Flux_V; ∆ lnZ > 3.5 over the second-best model). Removing the RVs from the equation, the favoured period remains consistent, although simpler shapes were favoured (P_cyc and P_cyc/2 only).

Using GP regression, we measured the stellar rotation period to be 43.63 ± 0.85 days, with a timescale of evolution of ${80}_{-26}^{+39}$ days. The measured period is consistent with previous measurements (Suárez Mascareño et al. 2015, 2017a), with a timescale of coherence of about two times that period. The rotation signal creates a variation in photometry of 3.27 ppt RMS, of 0.05 RMS in S_MW, of 4.2 m s⁻¹ RMS in FWHM, and of 1.8 m s⁻¹ RMS in RV.

Analysing the posterior distribution of the GP parameters, we find that the rotation variations in Ca II H&K and photometry are anti-correlated. The variations in Ca II H&K and FWHM are positively correlated. The variations in RV are positively correlated with the variations in Ca II H&K and with their gradient, and negatively with those of photometry and their gradient. The RV component linked to the gradient of the flux dominates over that linear with the flux. This behaviour is consistent with a spot-dominated stellar surface in which plages surround the spots. Lower fluxes (more spots) come together with higher Ca II H&K fluxes. The stellar RV variations are consistent with being produced by both inhibition of convection and the photometric effect. This is reinforced by the results coming from the FENRIR model.

Fig. 10 Comparison of the cycle variations. Cycle model in SMW, FWHM, and RV, showing the relationship between the components. The models have been scaled as if all original time series had the same RMS.

7 Conclusions

We revisited the system GJ 536 using all available spectroscopic data from HARPS, HARPS-N, CARMENES, and HIRES, and updated RV-extraction, and analysis methods. We performed a joint model that combined RV with information of different activity proxies, into the multi-dimensional Gaussian process framework.

Our analysis provided updated parameters of the known planet GJ 536 b, and the discovery of a new planet (GJ 536 c). GJ 536 b has an orbital period of 8.70874 ± 0.00056 days, a minimum mass of 6.37 ± 0.38 M_⊕, and receives ~10 times the Earth’s insolation. GJ 536 c has an orbital period of 32.761 ± 0.015 days, a minimum mass of 5.89 ± 0.70 M_⊕, and receives an insolation slightly lower than that received by Venus. Baed on statistical Doppler imaging, we estimated an inclination of the rotation axis of ${58}_{-19}^{+16}$ degrees. If we assume the orbital plane to be coplanar with the rotation axis, their masses would be ${7.6}_{-1.0}^{+2.3}$ M_⊕ and ${7.1}_{-1.2}^{+2.2}$ M_⊕, respectively. Given its projected angular separation of 15.5 mas, expected planet to star contrast of 1.8 ⋅ 10⁻⁸ to 3.0 ⋅ 10⁻⁸, and the brightness of the star, GJ 536 c is among the select group of low-mass exoplanets amenable for atmospheric characterization using its reflected light.

In addition to the analysis of the planetary system, we analysed the activity of the star. GJ 536 has a ${9.27}_{-0.17}^{+0.30}$ years activity cycle, detectable in the variations of the Ca II H&K fluxes, the FWHM of the stellar lines (at half the period), and the RV (at one quarter the period). The star rotates every 43.63 ± 0.81 days, and the induced signals are consistent with a spot-dominated surface in which plages and spots coexist.

The analysis of GJ 536 highlights the advances in data extraction and activity mitigation that have been made over the past years. These advances open the door to detecting planets that may have previously escaped detection, revealing new worlds hidden within archival data.

Data availability

The time series data used in this paper are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/702/A121.

Acknowledgements

A.S.M., J.I.G.H., and R.R. acknowledge financial support from the Spanish Ministry of Science and Innovation (MICINN) project PID2020-117493GB-I00 and from the Government of the Canary Islands project ProID2020010129. CdB acknowledges support from the Agencia Estatal de Investigación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) under grant WEAVE: EXPLORING THE COSMIC ORIGINAL SYMPHONY, FROM STARS TO GALAXY CLUSTERS and the European Regional Development Fund (ERDF) with reference PID2023-153342NB-I00/10.13039/501100011033, as well as from a Beatriz Galindo Senior Fellowship (BG22/00166) from the MICIU. The University of La Laguna (ULL) and the Department of Economy, Knowledge, and Employment of the Government of the Canary Islands are also gratefully acknowledged for the support provided to CdB (2024/347). The project that gave rise to these results received the support of a fellowship from the “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/DI23/11990071. 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). Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (NASA) This work is based on data obtained via the HARPS public database at the European Southern Observatory (ESO). We are grateful to all the observers of the following ESO projects, whose data we are using: 072.C-0488, 085.C-0019, 183.C-0972, and 191.C-087. We are grateful to the crews at the ESO observatories of Paranal and La Silla. This research has made extensive use of the SIMBAD database, operated at CDS, Strasbourg, France, and NASA’s Astrophysics Data System. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. The manuscript was written using VS Code. Main analysis performed in Python3 (Van Rossum & Drake 2009) running on Ubuntu (Sobell 2015) systems and MS. Windows running the Windows subsystem for Linux (WLS). Extensive usage of Numpy (van der Walt et al. 2011). Extensive usage of Scipy (Virtanen et al. 2020). All figures built with Matplotlib (Hunter 2007). The bulk of the analysis was performed on desktop PC with an AMD Ryzen^TM 9 9950X (16 cores, 32 threads, 3.5–4.7 GHz), provided by ASM, and a server hosting 2x AMD Epyc^TM 7663 (56 cores, 112 threads, per cpu), provided by the SUBSTELLAR ERC AdG.

Appendix A Correlations between RV and activity indicators

We studied the level of correlation between RV and the used spectroscopic activity indicators. Fig. A.1 shows th RV measurements against S index and FWHM. We computed the level of correlation using Spearman’s rank correlation coefficient (Spearman 1904). As the correlation coefficient does not take uncertainties into account, we performed 10 000 independent measures by taking random values from the distribution given by the error bars of the data, and computed the median result and its standard deviation. We measured a correlation coefficient of 0.216 ± 0.032 for the RV vs S index and 0.147 ± 0.042 for the RV vs FWHM. These values imply a weak correlation between the RV and the activity indicators. In addition, we modelled the relationship between these variables as a first-order polynomial. We measured a RV vs S index slope of 0.57 ± 0.25 m s⁻¹ and a RV vs FWHM slope of 0.098 ± 0.043. Both slopes are 3σ consistent with zero, supporting the weak level of correlation.

Fig. A.1 RV against S index ×10 (left panel) and FWHM (right panel) for GJ 536. The green line, and shaded areas, show the best fit to the data and the confidence intervals of the models.

Appendix B Model description

We performed global analyses of the data that included photometry, FWHM, S_MW index, and RV. Every time series includes a zero point per instrument, and models for the stellar cycle and rotation. The stellar rotation is modelled as a Gaussian Process. The HARPS FWHM data included a linear term to account for the focus drift of the instrument. The HARPS RV data, in addition, includes a polynomial term against the BERV. On top of them, the RV data includes a circular/Keplerian model for each planet in the model. To avoid numerical issues, we scaled up the S_MW index as S_MW index ×10.

The full model is defined as: $$\begin{array}{c}\Delta Flus=V0+Cycle+Rot,\\ \Delta S_{MW}\times 10=V0+Cycle+Rot,\\ \Delta FWHM=V0+f_{Focus}+Cycle+Rot,\\ \Delta RV=V0+f_{BERV}+Cycle+Rot+Planets,\end{array}$$(B.1)

where V0 represents the zero-point of each time series, with priors 𝒩[0,10] ppt for the photometric data, 𝒩[0,10] for the S_MW × 10, 𝒩[0,10] m s⁻¹ for the FWHM, and 𝒩[0,10] m s⁻¹ for the RVs. The rest of the components are described in the following subsections.

B.1 Cycle model

Based on the results of a S_MW-only model, the cycle component is defined as: $$\begin{array}{c}\Delta y=-A_1\cdot sin\left(2\pi \left(t-t_1\right)/P_{\text{cyc}}\right)\\ -A_2\cdot sin\left(4\pi \left(t-t_2\right)/P_{\text{cyc}}\right)\\ -A_3\cdot sin\left(8\pi \left(t-t_3\right)/P_{\text{cyc}}\right)\end{array}$$(B.2)

with: $$\begin{array}{l}{t}_1=t_0+P_{\text{cyc}}\cdot \left({\varphi }_1-1\right)\hfill \\ t_2=t_0+P_{\text{cyc}}\cdot \left({\varphi }_2-1\right)/2\hfill \\ t_3=t_0+P_{\text{cyc}}\cdot \left({\varphi }_3-1\right)/4\hfill \end{array}$$(B.3)

with t₀ = 7880 (BJD − 2450 000) being the integer date of the last RV observation. The cycle has independent amplitudes for photometry, S_MW × 10, FWHM, and RV. The period of the cycle has a prior Ln P 𝒰[7,8.7] (1100–6000 days).The phases have priors 𝒰[−0.5, 0.5]. The amplitudes in S_MW × 10 are restricted to be positive (𝒰[0, 2]), while the amplitudes in Photometry, FWHM and RV can be either positive or negative, to account for opposition of phase (𝒩[0, 10] in their respective units for all of them).

B.2 Stellar rotation model

To model the rotation we opted to work within the multidimensional Gaussian Processes (GP) framework (Rajpaul et al. 2015), which is based on the assumption that there exists an underlying function governing the behaviour of the stellar activity in all time series, which we denote G(t). G(t) manifests in each time series (Photometry, RV, etc.) as a linear combination of itself and its gradient, G′(t), with a set of amplitudes for each time series, following the idea of the FF′ formalism (Aigrain et al. 2012). We used the S+LEAF code (Delisle et al. 2022)⁵, which extends the formalism of semi-separable matrices introduced with celerite (Foreman-Mackey et al. 2017) to allow for fast evaluation of GP models even in the case of large datasets. The S+LEAF code supports a wide variety of GP kernels with fairly different properties. We opted for a combination of two stochastic harmonic oscillators (SHO), one with a period equal to the rotation period, and the second at half the rotation period. The kernel is defined as: $$\begin{array}{l}k\left(\tau ,P_{\text{rot}},L\right)=k_{SHO1}\left(\tau ,{\alpha }_1,P_1,Q_1\right)+k_{SHO2}\left(\tau ,{\alpha }_2,P_2,Q_2\right)\\ +\left({\sigma }^2\left(t\right)+{\sigma }_j^2\right)\cdot {\delta }_{\tau }\end{array}$$(B.4)

where k denotes the GP term, α_i is the standard deviation of the process, P_i is the period of the oscillator, related to the rotation period, and Q_i is the quality factor of the oscillator, related to the timescale of evolution L. These parameters are defined as: $$\begin{array}{l}{\alpha }_1=1;P_1=P_{\text{rot}};Q_1=\frac{\pi L}{P_{\text{rot}}}\hfill \\ {\alpha }_2=\beta ;P_2=0.5\cdot P_{\text{rot}};Q_2=\frac{2\pi L}{P_{\text{rot}}}\hfill \end{array}$$(B.5)

Equation (B.4) also includes a term of uncorrelated noise (σ), independent for every instrument, added quadratically to the diagonal of the covariance matrix to account for all unmodelled sources of variation, such as uncorrected activity, instrumental instabilities, or additional planets. δ_τ is the Kronecker delta function, and τ represents a time interval between two measurements, t − t′. The white noise model uses log-normal priors, centred around the dispersion of the data and with a sigma similar to the dispersion of the data. This parametrization makes it difficult for the model to converge to very low white noise values, while not forbidding them, which helps prevent overfitting.

The amplitudes α_i are related with the amplitude of the underlying function, not to any of the specific time series. These amplitudes are degenerate with the amplitudes of the component at each time series. We fixed α₁ to 1, and let α₂ free to account for scaling of the second oscillator, with a prior β = 𝒰[0, 1].

We used a prior 𝒰[35, 50] days for the rotation period, and a prior ln L 𝒩[4.5, 0.8] days for the timescale of evolution (${90}_{-50}^{+100}$ days), following Giles et al. (2017).

We used the S_MW × 10 data, and photometry, as main dataset guiding the GP, which meant not including a gradient amplitude for them. The GP model, over all time series, is described as: $$\begin{array}{c}\Delta S{MW}_{}\times 10=A_1{k}_1+\beta \cdot A_1{k}_2,\\ \Delta Flux=A_2{k}_1+\beta \cdot A_2{k}_2,\\ \Delta FWHM=A_3{k}_1+B_3{k^{\prime }}_1+\beta \cdot A_3{k}_2+\beta \cdot B_3{k^{\prime }}_2,\\ \Delta RV=A_4{k}_1+B_4{k^{\prime }}_1+\beta \cdot A_4{k}_2+\beta \cdot B_4{k^{\prime }}_2,\end{array}$$(B.6)

The scaling factor of the S_MW × 10 (A₁) has a prior 𝒰[0, 10], the scalling factor of the photometryc data (A₂) has a prior 𝒰[−20, 20] ppt. The scaling factors of the FWHM and RV linear components (A₃ and A₄) have priors 𝒰[−50, 50], while the scaling factors of the gradient components (B₃ and B₄) have priors 𝒰[−100, 100].

B.3 Correlation with the BERV

There are several effects potentially affecting the spectra, such as telluric contamination, ghosts, or CCD stitching, anchored at the detector reference frame (Dumusque et al. 2015; Coffinet et al. 2019; Cretignier et al. 2021). Combined with the movement of the Earth around the Sun, this results in periodic shifts with 1-year periodicity. To account for potential contamination, we included a polynomial term against the BERV (see Fig. 3 of Dumusque et al. 2015). This term is defined as shown in eq. B.7, with different parameters for each instrument. We only found the BERV term to be significant for the HARPS data.

The correlation takes the form of: $$\Delta RV=\alpha \cdot BERV+b\cdot BER{V}^2$$(B.7)

with a and b being the parameters of the polynomial, with priors 𝒩[0, 0.1].

B.4 HARPS focus drift

HARPS had a known focus drift that manifested as a linear trend in FWHM data. The focus drift was corrected in 2015, when the fibers were upgraded. To account for this, we include a polynomial against time in the model of the FWHM of the HARPS-03. $$\Delta FWHM=\alpha \cdot \left(t-t_0\right)$$(B.8)

with a being the slope, with prior 𝒩 [0, 0.1].

B.5 Planetary model

All previously stated components of the model define our null-hypothesis model, i.e. a model aimed to account for stellar activity and instrumental effect. On top of them, we include a planetary model where the number of planets will vary. The planetary signals are defined as circular for the detection tests, and later as full Keplerians to characterise their eccentricities. This choice avoids the potential pitfall of having a poorly sampled eccentric signal mimicking two circular signals and significantly cuts down on computing time.

RV variations due to circular planetary orbits are defined as: $$y\left(t\right)=-K\cdot \sin \left(2\pi \cdot \left(t-t_0\right)/P_{pl}\right)$$(B.9)

where t₀ = 7880 + P_pl ⋅ (ϕ_pl − 1). This parametrization ensures that our t₀ coincides to the inferior conjunction of the planets, and is within the baseline of observations.

When conducting a blind search for planets, the orbital period is parametrised as angular frequency ω = 2π/P_pl. We use a prior 𝒰[2π/1000, 2π/2]. When performing a guided search using the published solutions, we directly sample the periods and use 𝒩 priors centred around the published solution. The phases ϕ_pl are parametrised as 𝒰[−0.25, 0.75].

We parameterise the planetary amplitude K as ln K with a prior 𝒰[−5,2] m s⁻¹ (which keeps 50% of the prior space below 22 cm s⁻¹). This parametrization expands the parameter space around amplitudes consistent with zero, reducing potential biases towards large posteriors in noise dominated data, as described in Rajpaul et al. (2024).

RV variations due to planetary elliptical orbits are defined as: $$y\left(t\right)=-K\left(\cos \left(\eta +\omega \right)+e\cos \omega \right)$$(B.10)

where the true anomaly η is related to the solution of the Kepler equation, which depends on the orbital period of the planet P_orb and the orbital phase ϕ. This phase corresponds to the periastron time, which depends on the mid-point transit time T₀, the argument of periastron ω, and the eccentricity of the orbit e.

We parametrise $e={\left(\sqrt{e}cos\left(\omega \right)\right)}^2+{\left(\sqrt{e}sin\left(\omega \right)\right)}^2$ and $\omega ={\arctan }^2\left(\sqrt{e}sin\left(\omega \right),\sqrt{e}cos\left(\omega \right)\right)$. We then sample $\sqrt{e}cos\left(\omega \right)$ and $\sqrt{e}sin\left(\omega \right)$ with priors 𝒩[0, 0.3]. This parametrization favours low eccentricities, which are expected for short-period signals.

Appendix C Aliasing

Fig. C.1 Aliasing of activity signals. The top panel shows the window function of the RV data, with the main peakis suspected of creating aliases hightlighted. The rest of the panels show the periodogram of the activity-only RVs, with the suspected peaks created by activity highlighted, along with their main aliases.

Appendix D Adopted model – additional figures

Fig. D.1 Adopted activity model. Activity proxies’ data, with the best model fit, along with the periodogram of the data.

Fig. D.2 Zoom to selected observing campaigns. Activity proxies’ data of with the best model fit.

Appendix E Parameters of the model

Appendix F TESS photometry

GJ 536 (TIC 119147875) was observed by the Transiting Exoplanet Survey Satellite (TESS) (Ricker et al. 2015) with the 2 min cadence in sector 50. We analysed the pre-search data conditioned simple aperture photometry (PDCSAP) fluxes (Smith et al. 2012; Stumpe et al. 2012). The PDCSAP light curve shows an RMS of 400 ppm (155 ppm in 30-minute bins). Figure F.1 shows the TESS PDCSAP light curve. No aparent features consistent with transits appear in the data. We applied the box least squares (BLS) periodogram (Kovács et al. 2002; Hartman & Bakos 2016) to the TESS time series data to search for transit features. No transits have been found in the BLS periodogram.

Fig. F.1 TESS light curve. Blue symbols show the 2-minute cadence data. The red symbols represent 30-minute bins.

References

All Tables

All Figures

	Fig. 1 Spectroscopic and photometric data. Time series of RV, SMW, FWHM, and FluxV used in this study.
In the text

	Fig. 2 FIP periodogram of the full dataset. The highlighted peaks corresponds to the period of 8.7d (GJ 536 b), and 32.8 d.
In the text

	Fig. 3 Comparison of the posterior of RV and activity indicators. Posterior distributions of the signal parameters of the RV signals compared with the result of performing the same model in the activity indicators.
In the text

Fig. 4 Stability of the signals. The top panels show the evolution of the amplitude of the apodised signals. The shaded red area shows the confidence intervals of the amplitudes of the apodised signals over time. The horizontal black lines show the median value and 1-σ range of the static model. The bottom panels show the RV semi-amplitude vs time of inferior conjunction for the solutions obtained with the first half (H1) and second half (H2) of the dataset. The lines encapsulate the 95% confidence interval.

In the text

	Fig. 5 Adopted RV model. The top panels show the RV data (detrended from BERV), with the best model fit, along with the periodogram of the data. The bottom panels shows the residuals after the fit, along with their periodogram. Figure D.1 shows the model of the activity proxies.
In the text

	Fig. 6 Zoom to selected observing campaigns. RV data of with the best model fit.
In the text

	Fig. 7 Phase-folded plots of the planetary-induced RV signals. RV variations induced by GJ 536 b and c with the best model fit (top panels), and the residuals after the fit (bottom panels).
In the text

Fig. 8 GJ 536 b and c in context. The upper panel shows the masses of known planets, for masses below 100 M⊕, and periods shorter than 1000 days, with the planets of the system of GJ 536 highlighted, and the solar system planets as reference. The middle panel shows the same, but as a function of insolation. The lower panel shows the contrast in reflected light of low-mass planets orbiting bright stars (my < 12), compared to their angular separation to their parent star. The dotted vertical lines indicate the limits in which ANDES will be able to observe (Palle et al. 2025). The red symbols for the planets of GJ 536 show the potential range of contrasts. The green squares show the position of the ANDES golden sample.

In the text

	Fig. 9 Compatibility limits. The upper panel shows the RV amplitude limit (1–99% area) as a function of orbital period. The lower panels show the same, but for planetary masses. The blue rectangles show the 3σ amplitude/mass range of the detected planets.
In the text

	Fig. 10 Comparison of the cycle variations. Cycle model in SMW, FWHM, and RV, showing the relationship between the components. The models have been scaled as if all original time series had the same RMS.
In the text

	Fig. A.1 RV against S index ×10 (left panel) and FWHM (right panel) for GJ 536. The green line, and shaded areas, show the best fit to the data and the confidence intervals of the models.
In the text

	Fig. C.1 Aliasing of activity signals. The top panel shows the window function of the RV data, with the main peakis suspected of creating aliases hightlighted. The rest of the panels show the periodogram of the activity-only RVs, with the suspected peaks created by activity highlighted, along with their main aliases.
In the text

	Fig. D.1 Adopted activity model. Activity proxies’ data, with the best model fit, along with the periodogram of the data.
In the text

	Fig. D.2 Zoom to selected observing campaigns. Activity proxies’ data of with the best model fit.
In the text

	Fig. F.1 TESS light curve. Blue symbols show the 2-minute cadence data. The red symbols represent 30-minute bins.
In the text