© The Authors 2026

1 Introduction

The first discovery of an exoplanet orbiting a main-sequence star (Mayor & Queloz 1995) was made with the radial velocity (RV) method anticipated by Struve (1952) (a review is available by Burt et al. 2025). Since then, more than 6000 planets have been discovered, as reported in the NASA Exoplanet Archive¹ (Christiansen et al. 2025). Even though detections via the transit method (Henry et al. 2000; Charbonneau et al. 2000) outnumbered those found via the RV method, the RV method remains the most fruitful for the mass determination of exoplanets.

A high-resolution spectrograph is necessary to achieve the required RV precision for detecting exoplanets. The High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003), installed in 2003 at 3.6 m telescope at the Observatory of La Silla, Chile, was the first instrument able to reach an RV precision below 1 m s⁻¹. HARPS-N (Cosentino et al. 2012, 2014), an improved version of HARPS, was installed at the Telescopio Nazionale Galileo at the Observatorio del Roque de Los Muchachos, Spain.

Long-term RV surveys enable the detection of planets on wide orbits. The ROcky Planets in Equatorial Stars (RoPES) program (Suárez Mascareño et al. 2018) is a long-term survey to investigate a sample of 17 G and K dwarf-type stars with HARPS and HARPS-N. The main aim of the project is to detect Earth-like and super-Earth-like planets in the habitable zone (HZ) of stars similar to the Sun, and to provide targets for future atmospheric characterization with facilities such as ArmazoNes high Dispersion Echelle Spectrograph (Marconi et al. 2022; Palle et al. 2025), Habitable World Observatory (Feinberg et al. 2026), and LIFE (Quanz et al. 2022). To achieve this goal, RoPES combines historical data and recent observations from the program itself.

For HD 161098, one of the stars in the program, we merged the HARPS and HARPS-N datasets with ESPRESSO observations from a recent observing campaign. ESPRESSO is the most precise spectrograph in the world, with a precision on RVs up to 10 cm s⁻¹ (Pepe et al. 2021). It is installed at the Very Large Telescope at the Paranal Observatory, Chile.

The paper is structured as follows. In Sect. 2, we discuss the observations we have used in our analysis. In Sect. 3, we describe the star HD 161098. In Sect. 4, we discuss the analysis on the system. In Sect. 5, we discuss the results of our analysis. In Sect. 6, we summarize our work.

2 Observations

2.1 HARPS and HARPS-N

HARPS and HARPS-N are fiber-fed, high-resolution spectrographs. They span wavelengths between 380 nm and 690 nm. The instruments have a resolving power of ~115 000. They are both pressure- and temperature-stabilized to minimize the drift in RV. Wavelength-calibrated, high-resolution spectra are provided by the data reduction software (DRS). For HARPS, we used version 3.3.6 of the DRS adapted from the ESPRESSO DRS², which corrects for various systematics, such as lamp aging and lamp changes, thus increasing the consistency of long-term surveys. We used version 3.0.1 of the DRS for HARPS-N. Both DRSs provide RVs calculated with the cross-correlation function (CCF) technique (Fellgett 1953). Other by-products of the CCF are the full width at half maximum (FWHM) and the bisector time span (BIS), among others.

HARPS underwent a fiber link update in June 2015 (Lo Curto et al. 2015), which introduced an RV offset. We considered observations collected before that date to be the H03 dataset, and observations collected after that date the H15. For each of these datasets, we fitted a zero-point and a jitter term. More details on the methods adopted in our analysis are reported in Appendix A. We refer to the HARPS-N dataset as HN.

We extracted HARPS and HARPS-N CCF RVs with YARARA (Cretignier et al. 2021, 2023). YARARA corrects for different systematics at the spectral level, such as cosmic rays, interference patterns, telluric lines, point spread function variability, ghosts, stitching of the detector, and thorium-argon lamp contamination. YARARA checks for anomalous CCF, anomalous S/N, or anomalous residuals in the spectra. Observations that do not surpass a certain threshold are not taken into account for weighting the recipes of the pipeline, and we discarded them from the analysis for this reason. We had a total of 339 HARPS and HARPS-N epochs, split between 128 H03 epochs, 101 H15 epochs, and 110 HN epochs. After the quality control of YARARA, we remained with 119 H03 epochs, 80 H15 epochs, and 102 HN epochs. We have a standard deviations of the different datasets of 2.17 m s⁻¹ for H03, 1.41 m s⁻¹ for H15, and 1.51 m s⁻¹ for HN. The mean error per dataset is 0.44 m s⁻¹ for H03, 0.51 m s⁻¹ for H15, and 0.33 m s⁻¹ for HN on nightly binned observations. For HO3, we collected single exposures of 900 s. We have a mean S/N on the single exposure of 176 for H03 at order 55. For H15, we collected three consecutive exposures of 300 s each, and we have a mean S/N of 93 at order 55 for a single exposure. For HN we collected three consecutive exposures of 300 s each with a mean S/N of 113 at order 55.

2.2 ESPRESSO

ESPRESSO (Pepe et al. 2021) is a high-resolution spectrograph. The instrument covers a wavelength range comprised between 378.2 nm and 788.7 nm. For our observations, we used the high-resolution 1-UT mode, which has a resolving power of ~140 000. ESPRESSO is contained in a vacuum vessel to prevent temperature and pressure shifts. We used DRS version 3.3.10 to extract CCF RVs and activity indicators. Hereafter, we refer to ESPRESSO observations as E19. We have 51 epochs of observations. The standard deviation of the E19 dataset is 1.28 m s⁻¹ for nightly binned observations. The mean error is 0.16 m s⁻¹ on nightly binned observations. We collected three consecutive exposures per night, with an exposure time of 300 s per exposure, with a mean S/N of 203 at order 145.

3 HD 161098: Stellar parameters

HD 161098 is a bright (V = 7.67 ± 0.01 mag; Høg et al. 2000) G8 V star, close to the Solar System (d = 29.75 ± 0.02 pc; Gaia Collaboration 2020). HD 161098 has an effective temperature of 5610 ± 50 K, with a mass of 0.837 ± 0.029 M_⊙, a radius of ${0.866}_{-0.020}^{+0.020}$ R_⊙, and a luminosity of 0.6769 ± 0.0015 L_⊙. In Table B.1, we present a summary of the main characteristics of the star.

We used the method described by Kopparapu et al. (2014) to calculate the boundaries of the HZ in the case of a 1 M_⊕ planet. We considered their recent Venus and early Mars regimes respectively as the inner and outer edges for the optimistic HZì. We considered the runaway and maximum greenhouse regimes as the conservative inner and outer edges of the HZ. We found the conservative HZ between 0.8155 ± 0.0025 au and ${1.4097}_{-0.0062}^{+0.0063}$ au, which corresponds to an orbital period between ${294.0}_{-5.1}^{+5.4}$ d and ${668}_{-12}^{+13}$ d for a circular orbit. We found the optimistic HZ comprised between 0.6216 ± 0.0014 au and 1.4808 ± 0.0066 au, which corresponds to an orbital period between ${195.6}_{-3.4}^{+3.5}$ d and ${719}_{-13}^{+14}$ d for a circular orbit.

4 Analysis

An analysis of stellar activity is presented in Sect. 4.1. We applied a multidimensional Gaussian process (GP) (Appendix B) analysis to simultaneously fit for activity indicators and RVs (Rajpaul et al. 2015; Barragán et al. 2023). Multidimensional GPs are less likely to overfit the time series, and especially preserve the long-period signals when compared to 1D GP (Nari et al. 2026). In Fig. 1 we show the time series used in this analysis and the corresponding GLS periodograms (Zechmeister & Kürster 2009).

4.1 Stellar activity

We used different activity proxies: the full width at half maximum of the CCF, directly provided by the DRS; the S index, related to the intensity of the chromospheric emission in the Ca II H & K lines (we used the method of Lovis et al. (2011) to calculate the S index); the Hα, related to the strength of the Hα emission line (we calculated the Hα index with the method described in Gomes da Silva et al. (2011)); the bisector time span, which measures the asymmetry in the profile of a CCF (Queloz et al. 2001); and the contrast (CON) of the CCF given by the relative ratio between the center of the CCF and the wings. We see in Fig. 1 the GLS periodogram of the different activity indicators points toward a magnetic cycle of period P ~ 3500 d. We modeled the magnetic cycle of the star with two sinusoids, at the period of the cycle and at its first harmonic. We tried different numbers of harmonics, and this solution was the favored one in terms of the Bayesian evidence (Appendix A). We used a uniform prior on the period of the cycle 𝒰(2000 d, 6000 d). We found a significant amplitude for the magnetic cycle in FWHM, S index, Hα, and CON. We determined slightly different magnetic cycle periods depending on the indicator: P_{S index} = ${3845}_{-78}^{+80}$ d, P_FWHM = 3590 ± 130 d, $P_{\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}}$ = ${3900}_{-97}^{+95}$ d, P_H–α = ${3810}_{-260}^{+280}$ d, and P_CON = 3740 ± 210 d. They are all compatible with each other within 1σ except for P_FWHM and P_{S index}, and P_FWHM and $P_{\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}}$, which are compatible with each other at 1.66σ and 1.89σ, respectively. We found an amplitude of the main component of the magnetic cycle defined at 14.6σ for S index, at 9.9σ for FWHM, at 11.3σ for $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$, at 5.4σ for Hα, and 8.6σ for CON.

We applied a 1D GP to model the rotation-induced effect of the star. We found P_rot = ${28.22}_{-0.34}^{+0.28}$ d in S index. This value is compatible with the values found in all the other activity indicators. We analyzed the activity indicators from the HARPS RVBank (Trifonov et al. 2020). Due to the limited number of nights available, it did not add relevant information to the analysis. We do not consider it for the analysis that follows.

Fig. 1 Time series of RV and activity indicators (left) and their respective GLS periodograms (right). The dashed red line is the 10% false alarm probability (FAP), the dashed-dotted green line is the 1% FAP, and the dotted blue line is the 0.1% FAP. We plot dashed lines at the periods of the candidate planets. Panel a: RV time series. Panel b: RV GLS periodogram. Panel c: FWHM time series. Panel d: FWHM GLS periodogram. Panel e: BIS time series. Panel f: BIS GLS periodogram. Panel g: S index time series. Panel h: S index GLS periodogram. Panel i: Hα time series. Panel j: Hα GLS periodogram. Panel k: CON time series. Panel 1: CON GLS periodogram. Panel m: $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ time series. Panel n: $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ GLS periodogram. The window function of the observations is shown in light gray. The vertical red lines in the periodograms indicate the period of the candidate planets we found in the analysis. The candidate at 72.5 d only appears in RV.

4.2 Planetary signal search

We searched for candidate planetary signals in our RV dataset by measuring a false inclusion probability (FIP) periodogram (Hara et al. 2022b). We fitted three sinusoids simultaneously, which shared a log-uniform prior between 2 d and 1000 d. We chose these priors to minimize the risk of fitting for harmonics of the magnetic cycle, which has a period of more than 3500 d. We modeled the stellar activity in a multidimensional GP framework (for details on the nested sampling setup, see Appendix A). We used FWHM and S index as ancillary activity indicators because they have the most significant determination of the cycle among independent activity indicators. We considered a magnetic cycle for each of the datasets. We shared the period and phase of the components of the cycle among the RVs and activity indicators. We show in Fig. 2 the results of the FIP analysis. We found two periods with FIP below 1%, one at 72.5 d and the other at 687.2 d. This is the threshold indicated in Hara et al. (2022b) as a detection criterion. We did not find other signals with FIP below 10%. The presence of only two signals of interest justifies the usage of a three-sinusoid model.

We also tried an iterative blind addition of planets to the model. First, we considered a one-sinusoidal model, along with the multidimensional GP framework with FWHM and S index. We used a log uniform prior on the period 𝒰(0.6, 6.9) corresponding to potential periods between ~2 d and ~1000 d. We found a sinusoidal signal of amplitude K = 0.68 ± 0.12 m s⁻¹ at a period of ${678}_{-13}^{+11}$ d. The one-sinusoidal model is favored in terms of Δ ln Z by + 7.0 over the flat model (Appendix A).

We added a second sinusoid using the same prior on the period of the one-sinusoidal model, 𝒰(0.6, 6.9), for both sinusoids. We found a Δ lnZ of +7.0 compared to the onesinusoidal model. We found a sinusoidal with amplitude K = ${0.63}_{-0.12}^{+0.11}$ m s⁻¹ and period P = ${676}_{-10}^{+11}$ d and a sinusoidal with amplitude K = 0.63 ± 0.10 m s⁻¹ and period P = ${72.580}_{-0.057}^{+0.058}$ d.

We performed a three-sinusoidal model to search for additional significant signals. The three-sinusoidal model is disfavored, compared to the two-sinusoidal model, by Δ lnZ = −0.2. We found the two-sinusoidal model to be the best model to explain the data. We will call the signal at ~72.5 d candidate HD 161098 b, and the signal at ~680 d candidate HD 161098 c. We performed different tests to verify the planetary nature of the candidates.

Fig. 2 FIP periodogram for HD 161098. The horizontal lines represent the different FIP level of 50%, 10%, 1%.

4.3 Stability of signals

We made a test to check the stability of the signals through the observing campaigns. We took inspiration from the apodized test (Hara et al. 2022a). We multiplied the amplitude of the sinusoid by a Gaussian term G: $$\mathrm{\Delta }RV=-K\cdot \sin (2\pi \cdot \frac{t-t_0}{P_{\mathrm{p}\mathrm{l}}})\cdot G(\mu ,\sigma )$$(1)

Here μ is the center of the Gaussian, while σ is the width of the Gaussian. The planetary signals are stable on a timescale larger than the time span of observations; it is not be possible to define the center of the Gaussian, and the width of the Gaussian should be longer than the baseline of the observations. We used normal priors, 𝒩(72.5d,0.3d) for candidate b and 𝒩(680d,10d) for candidate c.

We found for signal b μ_b = 2 455 900 ± 4100 BJD and lnσ_b = ${14.4}_{-4.0}^{+3.8}$. We found for signal c μ_c = $2454{500}_{-3100}^{+4400}$ BJD and lnσ_c = ${12.3}_{-3.8}^{+5.2}$. The center of the Gaussian is not well defined for both signals. The width of the Gaussian is above or comparable with the baseline of the observations within the error bars. The evidence of the model is worsened by ΔlnZ = −1.6 compared to a model with two sinusoids alone with the same priors on the periods. We show in Fig. 3 the product of the posterior distribution of the amplitudes of the signals multiplied by the Gaussian filters.

Candidate b is stable through the time series. We found an asymmetric shape of the signal for candidate c, with a preference for the epochs of the H03 dataset. We tested an informed model without apodization on the H15, HN, and E19 datasets without H03. We used as priors 𝒩(72.5 d, 5 d) for candidate b and 𝒩(680d,40d) for candidate c, to also test the period of the candidates. We found K_b = 0.69 ± 0.11 m s⁻¹ and P_b = 72.81 ± 0.086 d. We found K_c = 0.49 ± 0.13 m s⁻¹ and P_c = 692 ± 19 d. The amplitude of the two signals is within 1σ compatible with the amplitudes found in the full analysis. The period of candidate b is different at 3.1σ, while the period of candidate c is within 1σ, compatible with the period found in the full dataset. The time of inferior conjunction is BJD ${2460256}_{-85}^{+75}$ for candidate b and BJD 2460880.3 ± 2.6 for candidate c. Both are compatible within 1σ with the values found in the full analysis. The difference in the period of candidate b could be related to the different sampling of the signal, once we exclude the H03 dataset. We repeated the same test with a uniform prior on the period, 𝒰(50 d, 100 d) for candidate b and 𝒰(500 d, 1000 d) for candidate c. We found K_b = 0.67 ± 0.11 m s⁻¹ and P_b = ${72.81}_{-0.09}^{+0.013}$ d. We found K_c = ${0.53}_{-0.14}^{+0.16}$ m s⁻¹ and P_c = ${800}_{-115}^{+32}$ d. The absence of a large portion of the dataset and a baseline that ~3000 d shorter can explain the difficulties in detecting the long-period signal, but it could also be related to an artifact introduced by the H03 dataset. Only additional observations can solve the riddle.

Fig. 3 Apodized test for candidate planets of HD 161098.

4.4 Signals in activity indicators

We searched for the sinusoidal signals we found in RVs in the activity indicators. For simplicity, we used a multidimensional GP model with RVs and two activity indicators. Even so, we fitted the signal only in RV and one activity indicator at a time. We shared the period and phase of the sinusoidal signals. The priors on the periods were 𝒰(50 d, 100 d) for candidate HD 161098 b and 𝒰(500 d, 1000 d) for candidate HD 161098 c. We summarize the values obtained in Table 1. All the models with the candidate planets in the activity indicators are disfavored with respect to the model with candidates in RV only. None of the amplitudes of both signals in activity indicators are significant at 2σ.

We tried to use an independent phase for activity indicators. The results are in Table 1. We did not find signals in activity indicators at 3σ, and all the models worsen the lnZ with respect to the model with candidate planets in RVs only. Phase_b calculated in RVs differs from phase_b calculated in S index by 1.25σ. Phase_c calculated in RVs differs from phase_c calculated in S index by 1.4σ. Phase_b calculated in RVs differs from phase_c calculated in FWHM by 3.4σ. Phase_c calculated in RVs differs from phase_c calculated in FWHM by 1.3σ. Phase_b calculated in RVs differs from phase_c calculated in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ by 1.8σ. Phase_c calculated in RVs differs from phase_c calculated in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ by 1.3σ. Phase_c is poorly defined for both the signals. Phase_c RV = ${0.80}_{-0.32}^{+0.28}$ and Phase_c $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ = ${1.28}_{-0.36}^{+0.14}$.

We tried to search for the signal in the activity indicators alone without RVs. We made a multidimensional GP model with two activity indicators. We used uniform priors 𝒰(50 d, 100 d) and 𝒰(500 d, 1000 d) for signals b and c. We made three models for each activity indicator: only rotation and magnetic cycle, addition of signal b, and addition of signal c. We summarize the results in Table 2. We searched for the signals in the S index in a multidimensional GP analysis with S index and FWHM. In terms of evidence, the one-sinusoidal model is disfavored compared to the stellar activity-only model by ΔlnZ = −2.9. The two-sinusoidal model is disfavored compared to the stellar activity-only model by ΔlnZ = −0.7. The slight improvement of the two-sinusoidal model over the one-sinusoidal model is due to the 2σ detection of a signal at ${768}_{-16}^{+33}$ d.

For FWHM, we considered a multidimensional GP with FWHM and S index. We found that the one-sinusoidal model is disfavored by ΔlnZ = −2.4 and the two-sinusoidal model is disfavored by Δ lnZ = −1.6. For $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$, we used a multidimensional GP model with $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ and FWHM. We found that the one-sinusoidal model is disfavored by Δ lnZ = −2.6 and the two-sinusoidal model is disfavored by ΔlnZ = −3.0.

We plot in Figs. C.1 and C.2 the posterior distributions of the periods for signals b and c in activity indicators in the two-sinusoidal model, where we see some peaks in the posterior distribution of P_b close to the period of the signal found in RVs. We show in Table C.1 the significance of the different models we analyzed in activity indicators. Signal b is not well defined in period and is not significant in amplitude at 2σ. Signal c is defined in amplitude at 2σ in all the activity indicators we tested, but it points to a different period compared to the signal we found in RVs.

We do not conclude that the signals in RVs are related to activity. On the other hand, the presence of some hints of signals at the period of candidate b, and the presence of a signal probably related to a harmonic of the cycle at ~750 d, makes us cautious. We thus claim two signals we found as candidate planets. We will delegate future work to confirm or reject the planetary hypothesis.

4.5 Informed search for signals of interest

We made an informed search for the two candidates. We use the parameters found in the current analysis as the final results of our model in Sect. 5. We made a multidimensional GP analysis with FWHM and S index as ancillary indicators. We used uniform priors 𝒰(50 d, 100 d) and 𝒰(500 d, 1000 d) for candidates b and c. We ran a one-sinusoidal model with the short-period candidate. We found a Δ lnZ = +9.2 with respect to a flat model. We obtained K_b = 0.67 ± 0.10 m s⁻¹ and P = ${72.561}_{-0.059}^{+0.058}$ d. We ran a two-sinusoidal model. We found K_b = 0.63 ± 0.10 m s⁻¹ and P_b = ${72.578}_{-0.060}^{+0.059}$ d, K_c = 0.64 ± 0.12 m s⁻¹ and P_c = ${682.5}_{-9.9}^{+9.5}$ d. We have a time of inferior conjunction T_0b = ${2460877.5}_{-2.9}^{+2.8}$ BJD for the short-period signal and T_0c = $2460{236}_{-63}^{+57}$ BJD for the long-period signal. This model was favored compared to the one-sinusoidal model by Δ lnZ = + 8.6. We show in Fig. 4 the GLS periodogram of the residuals of the two-sinusoidal model in the case of informed search. We see only a nonsignificant peak at ~26 d with false alarm probability (FAP) below 10%. We show the best parameters for the candidates in Table 3. We show in Table C.2 the prior and posterior of the best model adopted.

We also tried a model with two Keplerians. We followed the parameterization of the eccentricity proposed by Anderson et al. (2011); Eastman et al. (2013), with a combination of e and ω: $\sqrt{e}\cdot \cos (\omega )$ and $\sqrt{e}\cdot \sin (\omega )$. We used a normal prior 𝒩(0, 0.3) for the parameters. We imposed the eccentricity to be less than 0.99. We found the two-Keplerian model to be disfavored by Δ lnZ = −1.4 We found an eccentricity on the short-period Keplerian of ${0.19}_{-0.13}^{+0.19}$ and an eccentricity for the long-period Keplerian of ${0.41}_{-0.21}^{+0.17}$. Both are not defined at 2σ.

Fig. 4 GLS periodogram of the RV residuals after subtracting the two-sinusoidal model in an informed search on the period.

5 Discussion

5.1 Planetary system

We found two signals of interest at periods of ~72.5d and ~680d. Due to some hints of the presence of the signals in activity indicators, and the doubt on the coherence of the long-period signal in all the datasets, we cannot claim the signals as confirmed planets. We claim the discovery of two candidate planets. Only additional observations will shed light on the nature of the two signals.

Candidate HD 161098 b has amplitude K_b = 0.63 ± 0.10 m s⁻¹ and period P_b = ${72.578}_{-0.060}^{+0.059}$ d. Such a planet would have a minimum mass M_b sini = 3.63 ± 0.59 M_⊕. The planet would orbit the star at a separation of 0.3209 ± 0.003 au. It would receive an insolation of ${6.48}_{-0.40}^{+0.42}$ S_⊕ and it would have an equilibrium temperature of 406.8 ± 6.5 K following Eq. (2) (Seager 2010): $$T_{\mathrm{e}\mathrm{q}}=T_{\star }\sqrt{\frac{R_{\star }}{2a}}(1-A{)}^{1/4}.$$(2)

We considered a Bond albedo of 0.3. Candidate HD 161098 c has an amplitude K_c = 0.64± 0.12 m s⁻¹ and an orbital period P_c = ${682.5}_{-9.9}^{+9.5}$ d. If it were a planet it would have a minimum mass of M_csini = ${7.8}_{-1.4}^{+1.5}$ M_⊕. The planet would orbit the star at a separation of 1.429 ± 0.021 au. It would receive an insolation of ${0.326}_{-0.021}^{+0.022}$ S_⊕ and it would have an equilibrium temperature of 192.7 ± 3.2 K. HD 161098 c would reside between the boundaries of the conservative and optimistic outer habitable zone. We show in Fig. 5 the position of the planets relative to the habitable zone. We measured v sin i = ${1.0}_{-0.4}^{+0.3}$ km s⁻¹. To measure vsini, we used the FWHM of the CCF in a way similar to that done in Rainer et al. (2023). We determined the inclination of the star with the method described by Masuda & Winn (2020). We found an angle I = 51 ± 20 deg between the line of sight and the rotation axis of the star. With coplanar planets we would have true masses M_b = ${4.7}_{-1.1}^{+2.4}$ M_⊕ and M_c = ${10.2}_{-2.5}^{+5.3}$ M_⊕. Candidate b would be a super-Earth, with different possible compositions. Candidate c would be a sub-Neptune-like planet. Fig. 6 shows the phase-folded plot of the two candidates. Fig. C.3 shows a zoomed-in image of the same plot. We show in Table C.3 the evolution of the root mean square (RMS) of the residuals of RVs after subtracting the different contributions: activity-related terms, HD 161098 b, and HD 161098 c. The RMS of the residuals always decreases for each dataset at every step. In Table 3 we show the parameters of the two candidate planets. In Table C.2, we show the priors used in our analysis with the results of the analysis. The candidate HD 161098 c would reside in the optimistic HZ of its parent star. Few discovered planets reside in the HZ of stars hotter than 4000 K. We show in Fig. C.4 the plot of all the planets with masses below 20 M_⊕ and insolation between 0.1 S_⊕ and 10 S_⊕ orbiting around stars hotter than 4000 K. We used this limit to circumscribe the research to G- and K-type stars.

Fig. 5 Positions of candidate planets of HD 161098 b and HD 161098 c relative to the HZ of the star.

Fig. 6 Phase-folded RV plot of candidate planets of HD 161098.

Fig. 7 Detection limits for HD 161098. Top panel: detection limits in amplitude. Bottom panel: detection limits in mass.

5.2 Detection limits

We followed Nari et al. (2026) to measure the detection limits. We took a grid of 1000 periods uniformly sampled in logarithmic space between ~1.5 d and ~5000 d. We fixed the period of the sinusoidal from the value taken from the grid, and we considered the 99th percentile of the posterior distribution of the amplitude as our detection limit. We show in Fig. 7 the results of our analysis. We found a mean detection limit in amplitude of ~38 cm s⁻¹. At a few periods, the 1st percentile of the posterior of the amplitude is above 10 cm s⁻¹. These outliers could be related to additional noise not fitted by our model, but could also hide the presence of low-amplitude signals of a different origin. In mass, we can detect planets with minimum mass inferior to 1 M_⊕ up to almost 10 d, while our sensitivity to planets with minimum mass below 10 M_⊕ spans until ~2000 d. We can detect planets with masses below 20 M_⊕ up to 5000 d orbital period.

5.3 Stellar activity

In the multidimensional GP analysis we found P_cycle = ${4090}_{-130}^{+140}$ d and P_rot = ${28.22}_{-0.35}^{+0.30}$ d. The magnetic cycle is longer than one-half the timespan of the observations. Future observations will refine the P_cycle. We found a timescale of evolution equal to ${107}_{-21}^{+26}$ d. This is ~3.8 times the rotation period of the star. This value is different from that of the Sun, and permits us to classify HD 161098 as a beater star, following the classification of Giles et al. (2017). A Sun-like star has a timescale of decay of the spots of approximately one rotation, while for a beater we have a lifetime of spots of a few rotations. A long decay time of the spots can suggest the presence of polar spots on the surface of the star (Giles et al. 2017).

We analyzed the dispersion in the residuals of the activity indicators and RVs at different phases of the cycle. First, we considered the dispersion between BJD 2 455 500 and BJD 2 456 500, at the maximum of the activity of the cycle, and then between BJD 2 458 000 and BJD 2 459 500, in a phase of a minimum of activity of the cycle. At the minimum of the cycle region, we have 41 epochs. At the maximum of the cycle, we have 59 epochs of observation. We derived these two regions from visual inspection of Fig. 1. We used the RMS of the residuals as our metric of reference. We report the results of the analysis in Table 4. In RVs we see a decrease in the dispersion of the residuals at the minimum of the cycle, while we see an increase in the dispersion for both FWHM and S index. We need to take into account multiple factors. The two regions do not have observations from overlapping instruments. The H15 dataset shows a larger dispersion in the S index compared to H03. Both the regions of minimum and maximum of the cycle show a standard deviation of the residuals comparable to or larger than the standard deviation of the full dataset, except for the S index at the maximum of the cycle. In the raw RVs, we see a lower dispersion at the minimum of the cycle, but this is not strongly confirmed in the activity indicators. We cannot conclude that observing a star at the minimum of its magnetic activity is related to a lower dispersion of the dataset. We cannot exclude this hypothesis, due to external factors such as the use of different instruments and the nonhomogeneity of the time series.

6 Conclusion

We made an in-depth analysis of the star HD 161098. We found evidence of two candidate planets. HD 161098 b has a period of ${72.578}_{-0.060}^{+0.059}$ d, an amplitude of 0.63 ± 0.10 m s⁻¹, a minimum mass of 3.63± 0.59 M_⊕; it orbits the star at a separation of 0.3207 ± 0.0037 au. HD 161098 c has a period of ${682.5}_{-9.9}^{+9.5}$ d, an amplitude of 0.64 ± 0.12 m s⁻¹, a minimum mass of ${7.8}_{-1.4}^{+1.5}$ M_⊕, and orbits the star at a separation of 1.429 ± 0.021 au. We found the inclination of the system to be I = 51± 20°. This implies, in the case that the two candidates are confirmed, that we are in the presence of a super-Earth in the inner orbit and a super-Earth or a mini-Neptune planet in the outer orbit. The outer candidate resides slightly outside the conservative HZ of the star, in the optimistic HZ.

HD 161098 b is stable through the observing seasons. For HD 161098 c, we see a preference for the H03 dataset. We tried to recover the signal without H03, and we could only find the signal at the correct period if we considered a normal prior on the period. The results of the apodized test do not lead us to claim the signal as a confirmed detection. We searched for the candidate signals in the activity indicators. For both HD 161098 b and HD 161098 c, we found the signal to be at ~2σ level in the activity indicators only when we used a different phase between the RV signal and the signal in the activity indicators. Due to this test, we point out the need for additional observations to confirm the planetary nature of the two signals.

We measured the detection limits of the dataset. We found a detection limit below 1 M_⊕ up to an orbital period of ~10 d. We can detect of planets with minimum mass below 20 M_⊕ up to 5000 d.

We checked the differences in the dispersion of the residuals of RVs and activity indicators between the maximum and the minimum of the magnetic cycle. We did not find significant evidence that the observations collected at the minimum of the cycle are better suited for the detection of exoplanets.

Data availability

A table for RVs and activity indicators is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/709/A123

Acknowledgements

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). JIGH, AKS, RR, CAP, NN, VMP, and ASM acknowledge financial support from the Spanish Ministry of Science and Innovation (MICINN) project PID2020-117493GB-I00. 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. This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extra-solar planets data visualization, exchange, and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch. XD acknowledges the support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement SCORE No. 851555) and from the Swiss National Science Foundation under the grant SPECTRE (No. 200021_215200). This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. This research has extensively used the SIMBAD database operated at CDS, Strasbourg, France, and NASA’s Astrophysics Data System. This research has made use of 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 Overleaf. Extensive use of numpy (van der Walt et al. 2011) and scipy (Virtanen et al. 2020). The main analysis was performed in Python 3 (Van Rossum & Drake 2009) running on a Ubuntu system (Sobell 2015). Programs IDs for the HARPS observations we used in the analysis are: 072.C-0488(E), 091.C-0936(A), 183.C-0972(A), 192.C-0852(A), Udry, 105.20PH.001, 109.2392.001, 112.25SF.001, 113.26U2.001, and 198.C-0836(A). Programs IDs for the HARPS-N observations we used in the analysis are: CAT14A_83, CAT15A_140, CAT16A_109, CAT17A_38, CAT17A_58, CAT18A_115, CAT19A_159, CAT20A_121, and ITP15_7. ESPRESSO observations we used in the analysis have program ID 115.28 BD.

References

Appendix A Methods

For the parameter estimation, we used the nested-sampling tool Dynesty (Speagle 2020). Dynesty gives an estimate of the natural logarithm of the evidence associated with a model, allowing an easy model comparison. We have considered as a criterion to accept a more complex model an improvement in lnZ of 5 or more. We use for the inference of parameters a number of livepoints equal to the maximum between 500 and 40 times the number of parameters of the model. For the FIP analysis of Sect. 4.2, we have used a number of live points equal to 5 times the difference between the maximum and minimum frequencies we consider for the model, divided by the frequency resolution measured as 2π/(max(RV)–min(RV)). In this way, we have 5 live points per element of frequency resolution. We have used it as a stopping criterion for convergence when the sampler remains with less than Δ lnZ = 0.01 to explore. To account for the different zero-points of different instruments, we always consider an offset term for each instrument in the analysis. We consider an offset specifically for the FWHM for the HN dataset after BJD 2459500. This is due to an inspection of the behavior of the FWHM of HN for multiple stars. A jitter term is added in quadrature to the nominal error of the different instruments for each time series. The jitter term takes into account all the sources of noise we are not modeling for, and the instrumental noise. To remove outliers from different datasets, we bin the observations nightly and we apply a cut on the dataset consisting of a 5σ clipping joined with the exclusion of measurements where the error is larger than three times the median error for each dataset. To implement Gaussian processes in our analysis we used S+LEAF (Delisle et al. 2022). S+LEAF allows only a certain kind of semi-separable matrix to be taken into account as a covariance matrix. In this way, the computational cost scales linearly with the dimension of the dataset, instead of with its cube, as it used to be in standard implementations. Before a fit for the jitter was available, exploratory GLS periodograms of our datasets were generated adding in quadrature a white noise term to the error on the measurements equal to the standard deviation of the dataset. For parameter determination we considered the median of the posterior as the final value of the parameter and the 16th and 84th percentile of the posterior as lower and upper error. We have used the package for the GLS periodogram described by Zechmeister & Kürster (2009). We split the FWHM HN dataset, applying a different zero-point for observations taken before 2459500 BJD and observations taken after this date. This is due to an offset present only in the FWHM of HARPS-N.

Appendix B Stellar activity

Extremely precise spectrographs at 1 m s⁻¹ precision made effects from the star become detectable in the RV time series. Stellar systems are nonequilibrium systems on a variety of timescales.

Stellar oscillations can originate short-term RV variations on the timescale of minutes with amplitudes up to tens of cm s⁻¹ (O’Toole et al. 2008). We can use tailored observing times to mitigate the stellar pulsation effect (Dumusque et al. 2011; Chaplin et al. 2019).

Granulation generates RV variations with a timescale of a few minutes to multiple days and amplitudes of a few m s⁻¹ (Dumusque et al. 2011; Mathur et al. 2011). Multiple observations of the same target through the same night can partially average out the effect of granulation (Dumusque et al. 2011). In large surveys, due to time constraints, is not always possible to implement this method, and the granulation effect remains a barrier toward the detection of tiny signals.

The Sun hosts an 11 yr magnetic cycle (Schwabe 1844). Lovis et al. (2011) shows this is common in other stars too. The amplitude of magnetic cycles signals in RV time series can reach tens of m s⁻¹ and mimic long-period planets. We modeled the cycle component in our time series with a sinusoidal and, eventually, additional harmonics.

The rotation of the star is an important source of stellar noise in RV time series. Inhomogeneities in stellar flux and suppression of convection associated with starspots create stellar-related RV variations (Saar & Donahue 1997). The amplitude of the effect of stellar rotation on RV can reach tens of m s⁻¹ (Saar & Donahue 1997). These variations can be quasi-periodic and mimic or hide planetary signals (Meunier & Lagrange 2021).

We use Gaussian processes (GP) to model the rotation signal in RVs (Rasmussen & Williams 2006; Aigrain & ForemanMackey 2023). Instead of parametrizing an analytical function to interpolate the signal of activity, GPs sample for the covariance between observations K = k(t_i,t_j,ϕ), where k(t_i, t_j, ϕ) is the covariance between observations made at time t_i and t_j, and ϕ are the hyper-parameters of the covariance function. The covariance function is the kernel of the GP.

A common covariance function used in the field of exoplanets to model rotation-related activity is the quasi-periodic kernel defined as: $$k(t,t^{\mathrm{\prime }})={\sigma }^2\exp (-\frac{{(t-t^{\mathrm{\prime }})}^2}{2{\lambda }^2})\exp (-\frac{{\sin }^2\left(\frac{\pi |t-t^{\mathrm{\prime }}|}{P}\right)}{2{\mathrm{\Gamma }}^2})$$(B.1)

Where σ is the amplitude of the kernel, λ is the timescale of evolution of the correlations, P is the rotation period, and Γ is the harmonic complexity. We used the MEP kernel, an approximation of the quasi-periodic kernel. The MEP kernel is representable as an S+LEAF matrix, so a matrix which is a sum of semi-separable and a LEAF matrix (Delisle et al. 2020). This decreases the computational cost from 𝒪(n³) to 𝒪(n), with n the number of measurements.

To mitigate the tendency of the GPs to overfit (Nari et al. 2026), we can train them together with activity indicators in a multidimensional GP fashion (Rajpaul et al. 2015; Barragán et al. 2023). This framework was first inspired by the FF’ method (Aigrain et al. 2012). The FF’ method models the RV variation with a linear combination of the squared flux obtained in simultaneous photometry and the flux multiplied by its time derivative: $$\mathrm{R}\mathrm{V}=V_{r}F(t)\dot{F}(t)+V_c{F}^2(t)$$(B.2)

The term with the derivative represents the flux modulation coming from the break of the flux symmetry and the term with F²(t) is linked to the suppression of convective blueshift in magnetized regions.

In the multidimensional GP framework, both RV and activity indicators are modeled with a linear combination of an underlying GP and its first derivative: $$\mathrm{R}\mathrm{V}=V_{c}G(t)+V_r\dot{G}(t),$$(B.3) $$A_1=A_{1c}G(t),$$(B.4) $$A_2=A_{2c}G(t)+A_{2r}\dot{G}(t).$$(B.5)

where V_c, V_r, A_1c, A_2c, and A_2r are the coefficients of RV and two hypothetical activity indicators.

Appendix C Figures and tables

Fig. C.1 Search for candidate planets signals in activity indicators, posterior distribution. The blue vertical lines represent the median value of the posterior distribution. Panel a) Posterior distribution of candidate signal b in FWHM. Panel b) Posterior distribution of the period of candidate b in FWHM. Panel c) Posterior distribution of the amplitude of candidate b in S index. Panel d) Posterior distribution of the period of the candidate b in S index. Panel e) Posterior distribution of the amplitude of candidate b in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$. Panel f) Posterior distribution of the period of candidate b in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$.

Fig. C.2 Search for candidate planet signals in activity indicators, posterior distribution. The blue vertical lines represent the median value of the posterior distribution. Figure a) Posterior distribution of candidate signal c in FWHM. Panel b) Posterior distribution of the period of candidate c in FWHM. Panel c) Posterior distribution of the amplitude of candidate c in S index. Panel d) Posterior distribution of the period of the candidate c in S index. Panel e) Posterior distribution of the amplitude of candidate c in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$. Panel f) Posterior distribution of the period of candidate c in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$.

Fig. C.3 Zoomed-in image of the phase-folded plot of candidate planets of HD 161098.

Fig. C.4 Comparison between the candidate system of HD 161098 and exoplanetary systems present in the literature for stars hotter than 4000K. We limit ourselves to the comparison with G and K-type stars because this is the scope of the RoPES survey. Only 3 exoplanets have been discovered orbiting the conservative HZ of their parent star (green region), and 13 planets have been discovered orbiting the optimistic HZ (yellow region). The parameters of the presented planets are from the NASA Exoplanet Archive (Christiansen et al. 2025).

All Tables

All Figures

Fig. 1 Time series of RV and activity indicators (left) and their respective GLS periodograms (right). The dashed red line is the 10% false alarm probability (FAP), the dashed-dotted green line is the 1% FAP, and the dotted blue line is the 0.1% FAP. We plot dashed lines at the periods of the candidate planets. Panel a: RV time series. Panel b: RV GLS periodogram. Panel c: FWHM time series. Panel d: FWHM GLS periodogram. Panel e: BIS time series. Panel f: BIS GLS periodogram. Panel g: S index time series. Panel h: S index GLS periodogram. Panel i: Hα time series. Panel j: Hα GLS periodogram. Panel k: CON time series. Panel 1: CON GLS periodogram. Panel m: $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ time series. Panel n: $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ GLS periodogram. The window function of the observations is shown in light gray. The vertical red lines in the periodograms indicate the period of the candidate planets we found in the analysis. The candidate at 72.5 d only appears in RV.

In the text

	Fig. 2 FIP periodogram for HD 161098. The horizontal lines represent the different FIP level of 50%, 10%, 1%.
In the text

	Fig. 3 Apodized test for candidate planets of HD 161098.
In the text

	Fig. 4 GLS periodogram of the RV residuals after subtracting the two-sinusoidal model in an informed search on the period.
In the text

	Fig. 5 Positions of candidate planets of HD 161098 b and HD 161098 c relative to the HZ of the star.
In the text

	Fig. 6 Phase-folded RV plot of candidate planets of HD 161098.
In the text

	Fig. 7 Detection limits for HD 161098. Top panel: detection limits in amplitude. Bottom panel: detection limits in mass.
In the text

Fig. C.1 Search for candidate planets signals in activity indicators, posterior distribution. The blue vertical lines represent the median value of the posterior distribution. Panel a) Posterior distribution of candidate signal b in FWHM. Panel b) Posterior distribution of the period of candidate b in FWHM. Panel c) Posterior distribution of the amplitude of candidate b in S index. Panel d) Posterior distribution of the period of the candidate b in S index. Panel e) Posterior distribution of the amplitude of candidate b in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$. Panel f) Posterior distribution of the period of candidate b in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$.

In the text

Fig. C.2 Search for candidate planet signals in activity indicators, posterior distribution. The blue vertical lines represent the median value of the posterior distribution. Figure a) Posterior distribution of candidate signal c in FWHM. Panel b) Posterior distribution of the period of candidate c in FWHM. Panel c) Posterior distribution of the amplitude of candidate c in S index. Panel d) Posterior distribution of the period of the candidate c in S index. Panel e) Posterior distribution of the amplitude of candidate c in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$. Panel f) Posterior distribution of the period of candidate c in $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$.

In the text

	Fig. C.3 Zoomed-in image of the phase-folded plot of candidate planets of HD 161098.
In the text

Fig. C.4 Comparison between the candidate system of HD 161098 and exoplanetary systems present in the literature for stars hotter than 4000K. We limit ourselves to the comparison with G and K-type stars because this is the scope of the RoPES survey. Only 3 exoplanets have been discovered orbiting the conservative HZ of their parent star (green region), and 13 planets have been discovered orbiting the optimistic HZ (yellow region). The parameters of the presented planets are from the NASA Exoplanet Archive (Christiansen et al. 2025).

In the text