# The SOPHIE search for northern extrasolar planets

N. Heidari, I. Boisse, N C Hara, T G Wilson, F. Kiefer, G. Hébrard, F.

Philipot, S. Hoyer, K G Stassun, G W Henry, et al.

► To cite this version:

N. Heidari, I. Boisse, N C Hara, T G Wilson, F. Kiefer, et al.. The SOPHIE search for northern extrasolar planets. Astronomy and Astrophysics - A&A, 2024, 681, pp.A55. 10.1051/0004-6361/202347897. insu-04433639

HAL Id: insu-04433639

https://insu.hal.science/insu-04433639

Submitted on 2 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The SOPHIE search for northern extrasolar planets

N. Heidari, I. Boisse, N C Hara, T G Wilson, F. Kiefer, G. Hébrard, F.

Philipot, S. Hoyer, K G Stassun, G W Henry, et al.

► To cite this version:

N. Heidari, I. Boisse, N C Hara, T G Wilson, F. Kiefer, et al.. The SOPHIE search for northern extra-solar planets. Astronomy & Astrophysics - A&A, 2024, 681, pp.A55. 10.1051/0004-6361/202347897. insu-04433639

HAL Id: insu-04433639

https://insu.hal.science/insu-04433639v1

Submitted on 2 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The SOPHIE search for northern extrasolar planets

## XIX. A system including a cold sub-Neptune potentially transiting a $V = 6.5$ star HD 88986*

N. Heidari<sup>1,2,3,4,**</sup> ◦, I. Boisse<sup>1</sup>, N. C. Hara<sup>5</sup>, T. G. Wilson<sup>6</sup>, F. Kiefer<sup>7</sup>, G. Hébrard<sup>2,8</sup>, F. Philipot<sup>7</sup>, S. Hoyer<sup>1</sup>, K. G. Stassun<sup>9</sup> ◦, G. W. Henry<sup>10</sup>, N. C. Santos<sup>11,12</sup> ◦, L. Acuña<sup>13</sup> ◦, D. Almasian<sup>4</sup>, L. Arnold<sup>8,14</sup>, N. Astudillo-Defru<sup>15</sup>, O. Attia<sup>5</sup> ◦, X. Bonfils<sup>16</sup>, F. Bouchy<sup>5</sup>, V. Bourrier<sup>5</sup>, B. Collet<sup>1</sup> ◦, P. Cortés-Zuleta<sup>1</sup>, A. Carmona<sup>16</sup> ◦, X. Delfosse<sup>16</sup>, S. Dalal<sup>17</sup> ◦, M. Deleuil<sup>1</sup> ◦, O. D. S. Demangeon<sup>11,12</sup> ◦, R. F. Díaz<sup>20</sup>, X. Dumusque<sup>5</sup>, D. Ehrenreich<sup>5</sup>, T. Forveille<sup>16</sup> ◦, M. J. Hobson<sup>13,19</sup>, J. S. Jenkins<sup>24,31</sup> ◦, J. M. Jenkins<sup>25</sup> ◦, A. M. Lagrange<sup>7</sup>, D. W. Latham<sup>21</sup> ◦, P. Larue<sup>16</sup> ◦, J. Liu<sup>1</sup>, C. Moutou<sup>18</sup>, L. Mignon<sup>5</sup>, H. P. Osborn<sup>22,23</sup>, F. Pepe<sup>5</sup>, D. Rapetti<sup>25,26</sup>, J. Rodrigues<sup>11,12,27</sup>, A. Santerne<sup>1</sup>, D. Segransan<sup>5</sup>, A. Shporer<sup>28</sup>, S. Sulis<sup>1</sup> ◦, G. Torres<sup>29</sup>, S. Udry<sup>5</sup>, F. Vakili<sup>3</sup>, A. Vanderburg<sup>23</sup> ◦, O. Venot<sup>30</sup>, H. G. Vivien<sup>1</sup>, and J. I. Vines<sup>32</sup>

*(Affiliations can be found after the references)*

Received 7 September 2023 / Accepted 14 November 2023

## ABSTRACT

Transiting planets with orbital periods longer than 40 d are extremely rare among the 5000+ planets discovered so far. The lack of discoveries of this population poses a challenge to research into planetary demographics, formation, and evolution. Here, we present the detection and characterization of HD 88986 b, a potentially transiting sub-Neptune, possessing the longest orbital period among known transiting small planets (<4 $R_⊕$) with a precise mass measurement ($\sigma_M/M > 25\%$). Additionally, we identified the presence of a massive companion in a wider orbit around HD 88986. To validate this discovery, we used a combination of more than 25 yr of extensive radial velocity (RV) measurements (441 SOPHIE data points, 31 ELODIE data points, and 34 HIRES data points), Gaia DR3 data, 21 yr of photometric observations with the automatic photoelectric telescope (APT), two sectors of TESS data, and a 7-day observation of CHEOPS. Our analysis reveals that HD 88986 b, based on two potential single transits on sector 21 and sector 48 which are both consistent with the predicted transit time from the RV model, is potentially transiting. The joint analysis of RV and photometric data show that HD 88986 b has a radius of $2.49 \pm 0.18 R_⊕$, a mass of $17.2_{-3.8}^{+4.0} M_⊕$, and it orbits every $146.05_{-0.40}^{+0.43}$ d around a subgiant HD 88986 which is one of the closest and brightest exoplanet host stars (G2V type, $R = 1.543 \pm 0.065 R_⊕$, $V = 6.47 \pm 0.01$ mag, distance = 33.37 $\pm$ 0.04 pc). The nature of the outer, massive companion is still to be confirmed; a joint analysis of RVs, HIPPARCOS, and Gaia astrometric data shows that with a $3\sigma$ confidence interval, its semi-major axis is between 16.7 and 38.8 au and its mass is between 68 and 284 $M_{Jup}$. HD 88986 b's wide orbit suggests the planet did not undergo significant mass loss due to extreme-ultraviolet radiation from its host star. Therefore, it probably maintained its primordial composition, allowing us to probe its formation scenario. Furthermore, the cold nature of HD 88986 b ($460 \pm 8$ K), thanks to its long orbital period, will open up exciting opportunities for future studies of cold atmosphere composition characterization. Moreover, the existence of a massive companion alongside HD 88986 b makes this system an interesting case study for understanding planetary formation and evolution.

**Key words.** planets and satellites: detection – techniques: photometric – techniques: radial velocities

## 1. Introduction

Among the 5000+$^1$ planets discovered so far (and many more candidates), transiting planets (~4000) have a considerable impact on our understanding of the formation and evolution of planetary systems. Such planets, when orbiting a bright host star that allows radial velocity (RV) follow-up, can be accurately characterized in terms of fundamental parameters such as mass and density, allowing for their internal structure to be modeled (e.g., Heidari et al. 2022; Delrez et al. 2021). Moreover, these objects provide us with a great opportunity to gather

information about the composition and temperature of their atmospheres by transmission and/or emission spectroscopy (e.g., Tabernero et al. 2020).

Among all the known transiting exoplanets with precise mass and radius measurements (Otegi et al. 2020), those with orbital periods exceeding 40 d are extremely rare, representing only ~1% of the total population (as of June 7, 2023). This scarcity poses a significant challenge to our understanding of planet demographics, formation, evolution, and the potential for habitability. The limited number of such long-period exoplanets has compelled many studies on exoplanet occurrence rates to focus primarily on planets with relatively short periods (e.g., Silburt et al. 2015; Petigura et al. 2013; Fulton et al. 2017). Moreover, Kopparapu et al. (2013) showed that the inner boundary of the habitable zone, encircling main-sequence stars with spectral types earlier than approximately M4 ($T_{eff} > 2800$ K), is

<footnote>* Appendix D is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/681/A55</footnote>

<footnote>** CNES postdoctoral fellowship.</footnote>

<footnote n="1">¹ https://exoplanetarchive.ipac.caltech.edu</footnote>
# The SOPHIE search for northern extrasolar planets

## XIX. A system including a cold sub-Neptune potentially transiting a $V = 6.5$ star HD 88986 ⋆

N. Heidari<sup>1,2,3,4,**</sup> ◦, I. Boisse<sup>1</sup>, N. C. Hara<sup>5</sup>, T. G. Wilson<sup>6</sup>, F. Kiefer<sup>7</sup>, G. Hébrard<sup>2,8</sup>, F. Philipot<sup>7</sup>, S. Hoyer<sup>1</sup>, K. G. Stassun<sup>9</sup> ◦, G. W. Henry<sup>10</sup>, N. C. Santos<sup>11,12</sup> ◦, L. Acuña<sup>13</sup> ◦, D. Almasian<sup>4</sup>, L. Arnold<sup>8,14</sup>, N. Astudillo-Defru<sup>15</sup>, O. Attia<sup>5</sup> ◦, X. Bonfils<sup>16</sup>, F. Bouchy<sup>5</sup>, V. Bourrier<sup>5</sup>, B. Collet<sup>1</sup> ◦, P. Cortés-Zuleta<sup>1</sup>, A. Carmona<sup>16</sup> ◦, X. Delfosse<sup>16</sup>, S. Dalal<sup>17</sup> ◦, M. Deleuil<sup>1</sup> ◦, O. D. S. Demangeon<sup>11,12</sup> ◦, R. F. Díaz<sup>20</sup>, X. Dumusque<sup>5</sup>, D. Ehrenreich<sup>5</sup>, T. Forveille<sup>16</sup> ◦, M. J. Hobson<sup>13,19</sup>, J. S. Jenkins<sup>24,31</sup> ◦, J. M. Jenkins<sup>25</sup> ◦, A. M. Lagrange<sup>7</sup>, D. W. Latham<sup>21</sup> ◦, P. Larue<sup>16</sup> ◦, J. Liu<sup>1</sup>, C. Moutou<sup>18</sup>, L. Mignon<sup>5</sup>, H. P. Osborn<sup>22,23</sup>, F. Pepe<sup>5</sup>, D. Rapetti<sup>25,26</sup>, J. Rodrigues<sup>11,12,27</sup>, A. Santerne<sup>1</sup>, D. Segransan<sup>5</sup>, A. Shporer<sup>28</sup>, S. Sulis<sup>1</sup> ◦, G. Torres<sup>29</sup>, S. Udry<sup>5</sup>, F. Vakili<sup>3</sup>, A. Vanderburg<sup>23</sup> ◦, O. Venot<sup>30</sup>, H. G. Vivien<sup>1</sup>, and J. I. Vines<sup>32</sup>

*(Affiliations can be found after the references)*

Received 7 September 2023 / Accepted 14 November 2023

## ABSTRACT

Transiting planets with orbital periods longer than 40 d are extremely rare among the 5000+ planets discovered so far. The lack of discoveries of this population poses a challenge to research into planetary demographics, formation, and evolution. Here, we present the detection and characterization of HD 88986 b, a potentially transiting sub-Neptune, possessing the longest orbital period among known transiting small planets (<4 $R_⊕$) with a precise mass measurement ($\sigma_M/M > 25\%$). Additionally, we identified the presence of a massive companion in a wider orbit around HD 88986. To validate this discovery, we used a combination of more than 25 yr of extensive radial velocity (RV) measurements (441 SOPHIE data points, 31 ELODIE data points, and 34 HIRES data points), Gaia DR3 data, 21 yr of photometric observations with the automatic photoelectric telescope (APT), two sectors of TESS data, and a 7-day observation of CHEOPS. Our analysis reveals that HD 88986 b, based on two potential single transits on sector 21 and sector 48 which are both consistent with the predicted transit time from the RV model, is potentially transiting. The joint analysis of RV and photometric data show that HD 88986 b has a radius of $2.49 \pm 0.18 R_⊕$, a mass of $17.2_{-3.8}^{+4.0} M_⊕$, and it orbits every $146.05_{-0.40}^{+0.43}$ d around a subgiant HD 88986 which is one of the closest and brightest exoplanet host stars (G2V type, $R = 1.543 \pm 0.065 R_⊕$, $V = 6.47 \pm 0.01$ mag, distance = 33.37 $\pm$ 0.04 pc). The nature of the outer, massive companion is still to be confirmed; a joint analysis of RVs, HIPPARCOS, and Gaia astrometric data shows that with a $3\sigma$ confidence interval, its semi-major axis is between 16.7 and 38.8 au and its mass is between 68 and 284 $M_{Jup}$. HD 88986 b's wide orbit suggests the planet did not undergo significant mass loss due to extreme-ultraviolet radiation from its host star. Therefore, it probably maintained its primordial composition, allowing us to probe its formation scenario. Furthermore, the cold nature of HD 88986 b ($460 \pm 8$ K), thanks to its long orbital period, will open up exciting opportunities for future studies of cold atmosphere composition characterization. Moreover, the existence of a massive companion alongside HD 88986 b makes this system an interesting case study for understanding planetary formation and evolution.

**Key words.** planets and satellites: detection – techniques: photometric – techniques: radial velocities

## 1. Introduction

Among the 5000+<sup>1</sup> planets discovered so far (and many more candidates), transiting planets (~4000) have a considerable impact on our understanding of the formation and evolution of planetary systems. Such planets, when orbiting a bright host star that allows radial velocity (RV) follow-up, can be accurately characterized in terms of fundamental parameters such as mass and density, allowing for their internal structure to be modeled (e.g., Heidari et al. 2022; Delrez et al. 2021). Moreover, these objects provide us with a great opportunity to gather

information about the composition and temperature of their atmospheres by transmission and/or emission spectroscopy (e.g., Tabernero et al. 2020).

Among all the known transiting exoplanets with precise mass and radius measurements (Otegi et al. 2020), those with orbital periods exceeding 40 d are extremely rare, representing only ~1% of the total population (as of June 7, 2023). This scarcity poses a significant challenge to our understanding of planet demographics, formation, evolution, and the potential for habitability. The limited number of such long-period exoplanets has compelled many studies on exoplanet occurrence rates to focus primarily on planets with relatively short periods (e.g., Silburt et al. 2015; Petigura et al. 2013; Fulton et al. 2017). Moreover, Kopparapu et al. (2013) showed that the inner boundary of the habitable zone, encircling main-sequence stars with spectral types earlier than approximately M4 ($T_{eff} > 2800$ K), is

<footnote>* Appendix D is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/681/A55</footnote>

<footnote>**CNES postdoctoral fellowship.</footnote>

<footnote n="1">¹ https://exoplanetarchive.ipac.caltech.edu</footnote>
longer than ~11 d. This emphasizes the importance of exploring planetary systems with longer orbital periods in our pursuit of habitable planets. Additionally, new scientific exploration such as detecting exomoons has yet to be achieved. The importance of exomoon discovery is highlighted by our Moon's influence on Earth's spin dynamics (Li & Batygin 2014) as well as the prospective habitability of icy moons inside our own Solar System (Reynolds et al. 1983). The lack of successful exomoon discoveries could be linked to the dearth of long-period planets, as planets with longer orbital periods are more likely to harbor moons (Dobos et al. 2021). These examples highlight the importance of detecting and studying this “missing population”, thereby advancing our comprehension of exoplanet populations.

However, long-period transiting planets are particularly difficult to detect. The two primary methods used for detecting exoplanets, RV and transit photometry, each have their own challenges when it comes to detecting these elusive long-period planets. The RV method requires high-precision measurements over extended time spans to capture a full orbital period, while the transit method faces challenges due to both a lower transit probability and the limited baseline observations in most photometric surveys. These challenges are further compounded when dealing with smaller planets with a shallow transit depth and small RV semi-amplitude.

Here, we present the detection and characterization of a planetary system orbiting a $V = 6.5$ mag G2-type star: HD 88986 b, the longest-period transiting sub-Neptune among the known small planets ($<4 R_⊕$), accompanied by an outer massive companion. To accomplish this discovery, we employed a variety of observations, including photometric, spectroscopic, and astrometric techniques. The structure of this paper is organized as follows: Sect. 2 provides a description of the spectroscopic observations utilized in our analysis. In Sect. 3, the stellar properties of the host star are discussed. To identify the stellar rotational period, in Sect. 4 we present an analysis of various stellar activity indicators. In Sect. 5, we analyze RV data that led to the discovery of the sub-Neptune planet and a long-term curvature. In Sects. 6 and 7, we present the detection of a single transit event in TESS data sector 21 and perform a joint analysis, respectively. In Sect. 8, the result of our investigations into detecting HD 88986 b's second transit event in additional photometric data is presented. Finally, in Sects. 9 and 10, we discuss the origin of the long-term trend and conclude on this system, respectively. The paper also includes appendices presenting updated and new procedures for spectroscopic extractions on SOPHIE data.

## 2. Spectroscopic observations

We have intensive spectroscopic observations of HD 88986 spanning a remarkable duration of 25 yr. These observations comprise a total of 506 data points obtained through the utilization of three high-resolution spectrographs: SOPHIE (Bouchy et al. 2013), HIRES (Vogt & Penrod 1988), and ELODIE (Baranne et al. 1996). The whole dataset is shown in Fig. 1. It displays a clear long-term curvature on a time scale of at least 25 yr, as well as other variations on shorter time scales.

### 2.1. High-resolution spectroscopy with SOPHIE

HD 88986 has been monitored by the SOPHIE high-precision spectrograph mounted on the 1.93-m telescope at the Haute-Provence Observatory (OHP, France). The observations were carried out as part of Recherche de Planètes Extrasolaires (RPE)

<figure><img src="image_2.png" /><figcaption>Fig. 1. Radial velocity measurements of HD 88986 from ELODIE, HIRES, HIRES+, SOPHIE and SOPHIE+.</figcaption></figure>

program 1, also known as SP1, which is a high-precision program to search for Neptunes and Super-Earths orbiting bright stars in the solar neighborhood (Courcol et al. 2015; Hara et al. 2020; Heidari et al. 2022). Over a period of more than 15 yr, spanning from 2007 December 7th to 2023 March 12th, 441 high-resolution spectra were collected for this star (see Fig. 1). The observations were conducted in SOPHIE high resolution (HR) mode (resolving power of $\lambda/\Delta\lambda \approx 75\,000$ at 550 nm), with simultaneous thorium-argon (Th-Ar) or Fabry-Perot (FP) calibration light measurements. The latter allows us to track instrumental drift, ensuring precise and accurate RV measurements.

In June 2011, hexagonal fibers were installed in the SOPHIE spectrograph. This led to achieving an RV precision of $1-2\ m\ s^{-1}$ but also about $50\ m\ s^{-1}$ jump in the measured RVs of the standard stars (Bouchy et al. 2013). Hence, we separated the data before June 2011 (12 data points with the name SOPHIE) from the data after (429 data points with the name SOPHIE+). The exposure time was set for both data sets at 600–900 s to average the stellar oscillations, achieving a median signal-to-noise ratio (S/N) of 158 per pixel at 550 nm.

The RV was derived using the SOPHIE data reduction system (DRS, Bouchy et al. 2009). The DRS encompasses several crucial steps, including spectrum extraction, removal of telluric lines, correction for CCD charge transfer inefficiency (CTI), computation of the cross-correlation function (CCF) between the spectra and a binary mask, barycentric RV correction, and ultimately fitting Gaussian profiles to the CCFs to extract the RVs (Baranne et al. 1996; Pepe et al. 2002). We note that to extract HD 88986 RVs, we used a G2 mask. Additionally, prior to RV extraction, we corrected the spectra for the atmospheric dispersion effect (see Appendix A).

Once the RVs were extracted, we performed a correction for the nightly instrumental drift, which was measured through simultaneous calibration light observations. For this purpose, SOPHIE fiber A was dedicated to star observations, while fiber B was used to monitor a calibration lamp. In this configuration, sky background observations were not possible. Thus, it is crucial to identify and flag spectra contaminated by moonlight to ensure the accuracy of the RV analysis. By considering the phase and position of the Moon at the time of observation (see Appendix B for more details), we identified 31 spectra contaminated by the Moon. Consequently, we conservatively excluded these spectra from our analysis. We note that the inclusion or exclusion of these data points has no significant effect on our final results. Furthermore, we removed 1 data point that was identified as a $3\sigma$ RV outlier, along with an additional 17 measurements due to low signal-to-noise ratio ($S/N < 50$) and invalid calibration lamp flux. Consequently, a total of 50 data points were excluded from the SOPHIE+ measurements, and our final analysis incorporated
longer than ~11 d. This emphasizes the importance of exploring planetary systems with longer orbital periods in our pursuit of habitable planets. Additionally, new scientific exploration such as detecting exomoons has yet to be achieved. The importance of exomoon discovery is highlighted by our Moon's influence on Earth's spin dynamics (Li & Batygin 2014) as well as the prospective habitability of icy moons inside our own Solar System (Reynolds et al. 1983). The lack of successful exomoon discoveries could be linked to the dearth of long-period planets, as planets with longer orbital periods are more likely to harbor moons (Dobos et al. 2021). These examples highlight the importance of detecting and studying this “missing population”, thereby advancing our comprehension of exoplanet populations.

However, long-period transiting planets are particularly difficult to detect. The two primary methods used for detecting exoplanets, RV and transit photometry, each have their own challenges when it comes to detecting these elusive long-period planets. The RV method requires high-precision measurements over extended time spans to capture a full orbital period, while the transit method faces challenges due to both a lower transit probability and the limited baseline observations in most photometric surveys. These challenges are further compounded when dealing with smaller planets with a shallow transit depth and small RV semi-amplitude.

Here, we present the detection and characterization of a planetary system orbiting a $V = 6.5$ mag G2-type star: HD 88986 b, the longest-period transiting sub-Neptune among the known small planets ($<4 R_⊕$), accompanied by an outer massive companion. To accomplish this discovery, we employed a variety of observations, including photometric, spectroscopic, and astrometric techniques. The structure of this paper is organized as follows: Sect. 2 provides a description of the spectroscopic observations utilized in our analysis. In Sect. 3, the stellar properties of the host star are discussed. To identify the stellar rotational period, in Sect. 4 we present an analysis of various stellar activity indicators. In Sect. 5, we analyze RV data that led to the discovery of the sub-Neptune planet and a long-term curvature. In Sects. 6 and 7, we present the detection of a single transit event in TESS data sector 21 and perform a joint analysis, respectively. In Sect. 8, the result of our investigations into detecting HD 88986 b's second transit event in additional photometric data is presented. Finally, in Sects. 9 and 10, we discuss the origin of the long-term trend and conclude on this system, respectively. The paper also includes appendices presenting updated and new procedures for spectroscopic extractions on SOPHIE data.

## 2. Spectroscopic observations

We have intensive spectroscopic observations of HD 88986 spanning a remarkable duration of 25 yr. These observations comprise a total of 506 data points obtained through the utilization of three high-resolution spectrographs: SOPHIE (Bouchy et al. 2013), HIRES (Vogt & Penrod 1988), and ELODIE (Baranne et al. 1996). The whole dataset is shown in Fig. 1. It displays a clear long-term curvature on a time scale of at least 25 yr, as well as other variations on shorter time scales.

### 2.1. High-resolution spectroscopy with SOPHIE

HD 88986 has been monitored by the SOPHIE high-precision spectrograph mounted on the 1.93-m telescope at the Haute-Provence Observatory (OHP, France). The observations were carried out as part of Recherche de Planètes Extrasolaires (RPE)

<figure><img src="image_2.png" /><figcaption>Fig. 1. Radial velocity measurements of HD 88986 from ELODIE, HIRES, HIRES+, SOPHIE and SOPHIE+.</figcaption></figure>

program 1, also known as SP1, which is a high-precision program to search for Neptunes and Super-Earths orbiting bright stars in the solar neighborhood (Courcol et al. 2015; Hara et al. 2020; Heidari et al. 2022). Over a period of more than 15 yr, spanning from 2007 December 7th to 2023 March 12th, 441 high-resolution spectra were collected for this star (see Fig. 1). The observations were conducted in SOPHIE high resolution (HR) mode (resolving power of $\lambda/\Delta\lambda \approx 75\,000$ at 550 nm), with simultaneous thorium-argon (Th-Ar) or Fabry-Perot (FP) calibration light measurements. The latter allows us to track instrumental drift, ensuring precise and accurate RV measurements.

In June 2011, hexagonal fibers were installed in the SOPHIE spectrograph. This led to achieving an RV precision of $1-2\,m\,s^{-1}$ but also about $50\,m\,s^{-1}$ jump in the measured RVs of the standard stars (Bouchy et al. 2013). Hence, we separated the data before June 2011 (12 data points with the name SOPHIE) from the data after (429 data points with the name SOPHIE+). The exposure time was set for both data sets at 600–900 s to average the stellar oscillations, achieving a median signal-to-noise ratio (S/N) of 158 per pixel at 550 nm.

The RV was derived using the SOPHIE data reduction system (DRS, Bouchy et al. 2009). The DRS encompasses several crucial steps, including spectrum extraction, removal of telluric lines, correction for CCD charge transfer inefficiency (CTI), computation of the cross-correlation function (CCF) between the spectra and a binary mask, barycentric RV correction, and ultimately fitting Gaussian profiles to the CCFs to extract the RVs (Baranne et al. 1996; Pepe et al. 2002). We note that to extract HD 88986 RVs, we used a G2 mask. Additionally, prior to RV extraction, we corrected the spectra for the atmospheric dispersion effect (see Appendix A).

Once the RVs were extracted, we performed a correction for the nightly instrumental drift, which was measured through simultaneous calibration light observations. For this purpose, SOPHIE fiber A was dedicated to star observations, while fiber B was used to monitor a calibration lamp. In this configuration, sky background observations were not possible. Thus, it is crucial to identify and flag spectra contaminated by moonlight to ensure the accuracy of the RV analysis. By considering the phase and position of the Moon at the time of observation (see Appendix B for more details), we identified 31 spectra contaminated by the Moon. Consequently, we conservatively excluded these spectra from our analysis. We note that the inclusion or exclusion of these data points has no significant effect on our final results. Furthermore, we removed 1 data point that was identified as a $3\sigma$ RV outlier, along with an additional 17 measurements due to low signal-to-noise ratio (S/N < 50) and invalid calibration lamp flux. Consequently, a total of 50 data points were excluded from the SOPHIE+ measurements, and our final analysis incorporated
379 SOPHIE+ spectra, ensuring the reliability of our dataset for subsequent analysis.

SOPHIE experiences a long-term instrumental variation (Courcol et al. 2015), which is tracked by observing so-called “constant stars” every night. To account for this variation, following Courcol et al. (2015), we constructed a master time series using the RVs of these constant stars, which we then subtracted from the HD 88986 data. A detailed description of our update on the construction of the RV master constant time series can be found in Appendix C. The mean uncertainty of our final RV SOPHIE+ data is 1.2 m s⁻¹, with a root mean square (RMS) of 15.30 m s⁻¹. The final SOPHIE dataset is provided in Table D.1.

## 2.2. High-resolution spectroscopy with HIRES

The star was observed using the HIRES spectrograph from 1996 December 2nd to 2014 January 19th, spanning 17 yr during which 51 high-resolution spectra were obtained (see Fig. 1). Detailed information regarding data reduction and observation can be found in Butler et al. (2017). HIRES data experience a small jump of 1.5 ± 0.1 m s⁻¹ resulting from a CCD change in August 2004, as well as a long-term drift (≤1 m s⁻¹) and a small intra-night drift, as identified by Tal-Or et al. (2019). To account for these effects, we utilized the systematically corrected HIRES data obtained from the Vizier catalog access tool<fnref n="2" /> following the methodology outlined in Tal-Or et al. (2019). The mean uncertainty of the HIRES data is 1.2 m s⁻¹, with a RMS of 11.1 m s⁻¹. To address any residual offsets from the CCD change in the HIRES data, we fit an offset between the data obtained before (HIRES) and after (HIRES+) the CCD change. The HIRES RVs are presented in Table D.2. We note that our moon contamination criteria (see Appendix B) were not applied to HIRES data since these data are only utilized to investigate the origin of the long-term curvature and are not intended for high-precision detection of the sub-Neptune.

## 2.3. High-resolution spectroscopy with ELODIE

ELODIE was a high-resolution spectrograph, mounted on the 1.93-m telescope at OHP, which was in particular used to discover the first exoplanet in 1995 (Mayor & Queloz 1995). HD 88986 was observed by ELODIE from 1997 February 28 to 2004 January 29, gathering 31 high-resolution spectra (see Fig. 1). The K0 numerical mask is used to extract the RVs (Baranne et al. 1996). The exposure time varied from 600 to 900 s, resulting in a mean uncertainty of 9.0 m s⁻¹ and RMS of 13.0 m s⁻¹, which is close to the intrinsic stability of ELODIE. We note that 3 data points were removed due to their low S/N (<50). They are listed in Table D.3. Similar to the HIRES data, the moon contamination criteria were not applied, as this data was not utilized for the detection of the sub-Neptune but to investigate the origin of the long-term curvature.

# 3. Stellar parameters

HD 88986 is a G2V star with a G magnitude of 6.3. To obtain the stellar atmospheric parameters, we co-added all SOPHIE+ spectra (428) after correcting for the RV variation of the star, barycentric Earth radial velocity, and background pollution due to the calibration lamp (Heidari 2022; Hobson 2019). This results in a high S/N per pixel spectrum of 3174 at 550 nm. Then we calculated the effective temperature ($T_{eff}$), metallicity ([Fe/H]), and

<table><caption>**Table 1.** Stellar properties of HD 88986.</caption><thead><tr><th>Identifiers</th><th></th><th></th></tr></thead><tbody><tr><td colspan="3">TIC 1042868<br>HD 88986<br>HIP 50316<br>Gaia DR3 741184091114529792<br>2MASS J10162809+2840571</td></tr><tr><th>Parameter</th><th>HD 88986</th><th>References</th></tr><tr><td colspan="3">Astrometric properties</td></tr><tr><td>Parallax (mas)</td><td>29.9864 ± 0.0205</td><td>Gaia DR3</td></tr><tr><td>AEN ε (mas)</td><td>0.135</td><td>Gaia EDR3</td></tr><tr><td>Significance of ε</td><td>20.6</td><td>Gaia EDR3</td></tr><tr><td>Distance (pc)</td><td>33.37 ± 0.04</td><td>Gaia DR3</td></tr><tr><td>α (h m s)</td><td>10:16:28.08</td><td>Gaia DR3</td></tr><tr><td>δ (d m s)</td><td>28:40:56.94</td><td>Gaia DR3</td></tr><tr><td colspan="3">Photometric properties</td></tr><tr><td>B - V</td><td>0.635 ± 0.006</td><td>HIP</td></tr><tr><td>V(mag)</td><td>6.47 ± 0.01</td><td>HIP</td></tr><tr><td>Gaia (mag)</td><td>6.315 ± 0.003</td><td>Gaia DR3</td></tr><tr><td>Gaia<sub>BP</sub> (mag)</td><td>6.628 ± 0.003</td><td>Gaia DR3</td></tr><tr><td>Gaia<sub>RP</sub> (mag)</td><td>5.822 ± 0.004</td><td>Gaia DR3</td></tr><tr><td>TESS (mag)</td><td>5.8706 ± 0.0061</td><td>TESS</td></tr><tr><td>J (mag)</td><td>5.247 ± 0.024</td><td>2MASS</td></tr><tr><td>H (mag)</td><td>4.946 ± 0.023</td><td>2MASS</td></tr><tr><td>K<sub>s</sub> (mag)</td><td>4.884 ± 0.020</td><td>2MASS</td></tr><tr><td>W<sub>1</sub> (mag)</td><td>4.895 ± 0.239</td><td>WISE</td></tr><tr><td>W<sub>2</sub> (mag)</td><td>4.762 ± 0.085</td><td>WISE</td></tr><tr><td>W<sub>3</sub> (mag)</td><td>4.933 ± 0.014</td><td>WISE</td></tr><tr><td>W<sub>4</sub> (mag)</td><td>4.873 ± 0.029</td><td>WISE</td></tr><tr><td colspan="3">Spectroscopic properties</td></tr><tr><td>Spectral type</td><td>G2V</td><td>MK classification</td></tr><tr><td>log g (cm s<sup>-2</sup>)</td><td>4.16 ± 0.03</td><td>Sect. 3</td></tr><tr><td>ξ<sub>f</sub> (km s<sup>-1</sup>)</td><td>1.11 ± 0.02</td><td>Sect. 3</td></tr><tr><td>log(R'<sub>HK</sub>)</td><td>-5.04 ± 0.10</td><td>Sect. 4</td></tr><tr><td>ν sin i (km s<sup>-1</sup>)</td><td>3.3 ± 1.0</td><td>SOPHIE DRS</td></tr><tr><td>[Fe/H]</td><td>0.06 ± 0.02</td><td>Sect. 3</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>5861 ± 17</td><td>Sect. 3</td></tr><tr><td colspan="3">Bulk properties</td></tr><tr><td>Mass (M<sub>☉</sub>)</td><td>1.25 ± 0.05[0.06]<sup>(*)</sup></td><td>SED Sect. 3</td></tr><tr><td>Radius (R<sub>☉</sub>)</td><td>1.543 ± 0.010[0.065]<sup>(*)</sup></td><td>SED Sect. 3</td></tr><tr><td>P<sub>rot</sub> (d)</td><td>25<sup>+8</sup><sub>-6</sub></td><td>Sect. 4</td></tr><tr><td>Age (Gyr)</td><td>7.9 ± 1.3[1.58]<sup>(*)</sup></td><td>Sect. 3</td></tr></tbody></table>

**Notes.** (*) Adopted the systematic uncertainty floor as suggested by Tayar et al. (2022) throughout this study.

surface gravity (log g), using the procedure described in Santos et al. (2013) and Sousa et al. (2018). The resulting $T_{eff}$, [Fe/H], and log g together with other stellar parameters of HD 88986, are presented in Table 1.

We performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia EDR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius, following the procedures described in

<footnote n="2">² https://vizier.u-strasbg.fr/viz-bin/VizieR</footnote>
379 SOPHIE+ spectra, ensuring the reliability of our dataset for subsequent analysis.

SOPHIE experiences a long-term instrumental variation (Courcol et al. 2015), which is tracked by observing so-called “constant stars” every night. To account for this variation, following Courcol et al. (2015), we constructed a master time series using the RVs of these constant stars, which we then subtracted from the HD 88986 data. A detailed description of our update on the construction of the RV master constant time series can be found in Appendix C. The mean uncertainty of our final RV SOPHIE+ data is 1.2 m s⁻¹, with a root mean square (RMS) of 15.30 m s⁻¹. The final SOPHIE dataset is provided in Table D.1.

## 2.2. High-resolution spectroscopy with HIRES

The star was observed using the HIRES spectrograph from 1996 December 2nd to 2014 January 19th, spanning 17 yr during which 51 high-resolution spectra were obtained (see Fig. 1). Detailed information regarding data reduction and observation can be found in Butler et al. (2017). HIRES data experience a small jump of 1.5 ± 0.1 m s⁻¹ resulting from a CCD change in August 2004, as well as a long-term drift (≤1 m s⁻¹) and a small intra-night drift, as identified by Tal-Or et al. (2019). To account for these effects, we utilized the systematically corrected HIRES data obtained from the Vizier catalog access tool<fnref n="2" /> following the methodology outlined in Tal-Or et al. (2019). The mean uncertainty of the HIRES data is 1.2 m s⁻¹, with a RMS of 11.1 m s⁻¹. To address any residual offsets from the CCD change in the HIRES data, we fit an offset between the data obtained before (HIRES) and after (HIRES+) the CCD change. The HIRES RVs are presented in Table D.2. We note that our moon contamination criteria (see Appendix B) were not applied to HIRES data since these data are only utilized to investigate the origin of the long-term curvature and are not intended for high-precision detection of the sub-Neptune.

## 2.3. High-resolution spectroscopy with ELODIE

ELODIE was a high-resolution spectrograph, mounted on the 1.93-m telescope at OHP, which was in particular used to discover the first exoplanet in 1995 (Mayor & Queloz 1995). HD 88986 was observed by ELODIE from 1997 February 28 to 2004 January 29, gathering 31 high-resolution spectra (see Fig. 1). The K0 numerical mask is used to extract the RVs (Baranne et al. 1996). The exposure time varied from 600 to 900 s, resulting in a mean uncertainty of 9.0 m s⁻¹ and RMS of 13.0 m s⁻¹, which is close to the intrinsic stability of ELODIE. We note that 3 data points were removed due to their low S/N (<50). They are listed in Table D.3. Similar to the HIRES data, the moon contamination criteria were not applied, as this data was not utilized for the detection of the sub-Neptune but to investigate the origin of the long-term curvature.

## 3. Stellar parameters

HD 88986 is a G2V star with a G magnitude of 6.3. To obtain the stellar atmospheric parameters, we co-added all SOPHIE+ spectra (428) after correcting for the RV variation of the star, barycentric Earth radial velocity, and background pollution due to the calibration lamp (Heidari 2022; Hobson 2019). This results in a high S/N per pixel spectrum of 3174 at 550 nm. Then we calculated the effective temperature ($T_{eff}$), metallicity ([Fe/H]), and

<table><caption>**Table 1.** Stellar properties of HD 88986.</caption><thead><tr><th>Identifiers</th><th></th><th></th></tr></thead><tbody><tr><td colspan="3">TIC 1042868<br>HD 88986<br>HIP 50316<br>Gaia DR3 741184091114529792<br>2MASS J10162809+2840571</td></tr><tr><th>Parameter</th><th>HD 88986</th><th>References</th></tr><tr><td colspan="3">Astrometric properties</td></tr><tr><td>Parallax (mas)</td><td>29.9864 ± 0.0205</td><td>Gaia DR3</td></tr><tr><td>AEN ε (mas)</td><td>0.135</td><td>Gaia EDR3</td></tr><tr><td>Significance of ε</td><td>20.6</td><td>Gaia EDR3</td></tr><tr><td>Distance (pc)</td><td>33.37 ± 0.04</td><td>Gaia DR3</td></tr><tr><td>α (h m s)</td><td>10:16:28.08</td><td>Gaia DR3</td></tr><tr><td>δ (d m s)</td><td>28:40:56.94</td><td>Gaia DR3</td></tr><tr><td colspan="3">Photometric properties</td></tr><tr><td>B - V</td><td>0.635 ± 0.006</td><td>HIP</td></tr><tr><td>V(mag)</td><td>6.47 ± 0.01</td><td>HIP</td></tr><tr><td>Gaia (mag)</td><td>6.315 ± 0.003</td><td>Gaia DR3</td></tr><tr><td>Gaia<sub>BP</sub> (mag)</td><td>6.628 ± 0.003</td><td>Gaia DR3</td></tr><tr><td>Gaia<sub>RP</sub> (mag)</td><td>5.822 ± 0.004</td><td>Gaia DR3</td></tr><tr><td>TESS (mag)</td><td>5.8706 ± 0.0061</td><td>TESS</td></tr><tr><td>J (mag)</td><td>5.247 ± 0.024</td><td>2MASS</td></tr><tr><td>H (mag)</td><td>4.946 ± 0.023</td><td>2MASS</td></tr><tr><td>K<sub>s</sub> (mag)</td><td>4.884 ± 0.020</td><td>2MASS</td></tr><tr><td>W<sub>1</sub> (mag)</td><td>4.895 ± 0.239</td><td>WISE</td></tr><tr><td>W<sub>2</sub> (mag)</td><td>4.762 ± 0.085</td><td>WISE</td></tr><tr><td>W<sub>3</sub> (mag)</td><td>4.933 ± 0.014</td><td>WISE</td></tr><tr><td>W<sub>4</sub> (mag)</td><td>4.873 ± 0.029</td><td>WISE</td></tr><tr><td colspan="3">Spectroscopic properties</td></tr><tr><td>Spectral type</td><td>G2V</td><td>MK classification</td></tr><tr><td>log g (cm s<sup>-2</sup>)</td><td>4.16 ± 0.03</td><td>Sect. 3</td></tr><tr><td>ξ<sub>f</sub> (km s<sup>-1</sup>)</td><td>1.11 ± 0.02</td><td>Sect. 3</td></tr><tr><td>log(R'<sub>HK</sub>)</td><td>-5.04 ± 0.10</td><td>Sect. 4</td></tr><tr><td>ν sin i (km s<sup>-1</sup>)</td><td>3.3 ± 1.0</td><td>SOPHIE DRS</td></tr><tr><td>[Fe/H]</td><td>0.06 ± 0.02</td><td>Sect. 3</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>5861 ± 17</td><td>Sect. 3</td></tr><tr><td colspan="3">Bulk properties</td></tr><tr><td>Mass (M<sub>☉</sub>)</td><td>1.25 ± 0.05[0.06]<sup>(*)</sup></td><td>SED Sect. 3</td></tr><tr><td>Radius (R<sub>☉</sub>)</td><td>1.543 ± 0.010[0.065]<sup>(*)</sup></td><td>SED Sect. 3</td></tr><tr><td>P<sub>rot</sub> (d)</td><td>25<sup>+8</sup><sub>-6</sub></td><td>Sect. 4</td></tr><tr><td>Age (Gyr)</td><td>7.9 ± 1.3[1.58]<sup>(*)</sup></td><td>Sect. 3</td></tr></tbody></table>

**Notes.** (*) Adopted the systematic uncertainty floor as suggested by Tayar et al. (2022) throughout this study.

surface gravity (log g), using the procedure described in Santos et al. (2013) and Sousa et al. (2018). The resulting $T_{eff}$, [Fe/H], and log g together with other stellar parameters of HD 88986, are presented in Table 1.

We performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia EDR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius, following the procedures described in

<footnote n="2">² https://vizier.u-strasbg.fr/viz-bin/VizieR</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Spectral energy distribution of HD 88986. Red symbols represent the observed photometric measurements, whereas the horizontal bars represent the effective width of the passband. Blue symbols are the model fluxes from the best-fit Kurucz atmosphere model (black).</figcaption></figure>

Stassun & Torres (2016); Stassun et al. (2017); Stassun & Torres (2018). We obtained the $B_T V_T$ magnitudes from Tycho-2 (Høg 2001), the JHK$_S$ magnitudes from 2MASS (Cutri et al. 2003), the W1–W4 magnitudes from WISE (Wright et al. 2010), the uvby Strömgren magnitudes from Paunzen (2015), and the $G G_{BP} G_{RP}$ magnitudes from Gaia (Gaia Collaboration 2016, 2021). We also used the UV measurement at 274 nm from the TDI UV satellite (Thompson et al. 1978). Together, the available photometry spans the full stellar SED over the wavelength range 0.2–22 µm (see Fig. 2). Then, we performed a fit using Kurucz’s stellar atmosphere models with the $T_{\text{eff}}$, log g, and [Fe/H] adopted from the spectroscopic analysis. The remaining free parameter is the extinction $A_V$, which we fixed at zero due to the proximity of the system to Earth.

The resulting fit (Fig. 2) has a reduced $\chi^2$ of 1.4. Integrating the model SED gives the bolometric flux at Earth, $F_{\text{bol}} = 7.28 \pm 0.17 \times 10^{-8} \text{ erg s}^{-1} \text{ cm}^{-2}$. Taking the $F_{\text{bol}}$ and $T_{\text{eff}}$ together with the Gaia parallax, gives the stellar radius, $R_\star = 1.543 \pm 0.010 R_\odot$, placing the star within the subgiant range (1.5–3 $R_\odot$; Huber et al. 2017; Berger et al. 2018). In addition, we can estimate the stellar mass from the empirical relations of Torres et al. (2010), giving $M_\star = 1.19 \pm 0.07 M_\odot$, which is consistent with the value of $1.25 \pm 0.05 M_\odot$ determined empirically via $R_\star$ and log g. We acknowledge the possibility that our formal error budget for radius, mass, and stellar age could be underestimated, as suggested by Tayar et al. (2022). For stars such as HD 88986, the systematic uncertainty floor could rather be up to ≈4.2% for radius, ≈5% for mass, and ≈20% for age. Hence, we conservatively adopt these relative uncertainties, reported in brackets in Table 1, to obtain more realistic stellar parameter errors throughout this study.

Finally, we can use the observed log $R'_{\text{HK}}$ (see Sect. 4) to estimate the stellar age via empirical activity-age relations. We obtain an age of $\tau_\star = 7.9 \pm 1.3$ Gyr via the empirical relations of Mamajek & Hillenbrand (2008).

## 4. Stellar rotation and activity

To study the host star’s activity and rotational period, we used log $R'_{\text{HK}}$, S-index, Na index, and CCF bisector measurements from the SOPHIE+ spectrograph, along with the S-index values from the HIRES and HIRES+ spectrographs. The bisector spans

are derived from the SOPHIE DRS using the method described by Queloz et al. (2001). We extracted the Na index as introduced in Da Silva et al. (2011). The HIRES and HIRES+ S-index values were acquired from Butler et al. (2017). For the extraction of log $R'_{\text{HK}}$ and S-index from the SOPHIE+ spectra, we followed the procedure outlined by Noyes et al. (1984) and Boisse et al. (2010), respectively. The key step before deriving the log $R'_{\text{HK}}$ and S-index is subtracting the background contamination light from the Th-Ar or FP calibration lamp from the stellar spectra. To correct this, we used the direct measurement method (see Appendix E).

To estimate the rotational period of the star, we summed 176 HD 88986 SOPHIE+ spectra which fulfilled two criteria. First, the spectra with $S/N > 50$ in the first (bluest, $\lambda \sim 3955 \AA$) order of spectra where CaII H&K lines are located. Second, the spectra with minimal contamination due to the background light. This led to the value of log $R'_{\text{HK}} = -5.04 \pm 0.10$. Furthermore, we investigated potential changes in magnetic activity over an approximately 11-yr SOPHIE+ observation span. For this purpose, the dataset was divided into three distinct observational seasons spanning from 2012 to 2023: 2012–2016, 2016–2019, and 2019–2023. We observed a gradual increase in the log $R'_{\text{HK}}$ parameter, with values transitioning from $-5.12 \pm 0.10$ during the initial subset to $-5.05 \pm 0.10$ in the second subset, and ultimately stabilizing at $-5.00 \pm 0.10$ in the final subset. All these values are in good agreement with the value of log $R'_{\text{HK}} = -5.22$ and log $R'_{\text{HK}} = -5.07$ reported by Radick et al. (2018) and Hall et al. (2007), respectively. Finally, we estimated a rotational period of $25^{+8}_{-6}$ d following Noyes et al. (1984), a value consistent with the Mamajek & Hillenbrand (2008) recipe, which yields $26.3 \pm 3.1$ d.

We searched for rotational modulation in the Simple Aperтуre Photometric (SAP; see Twicken et al. 2010; Morris et al. 2020) TESS data (see Sects. 6 and 8.2 for details of the observations) using the Lomb-Scargle periodogram (Lomb 1976; Scargle 1982; VanderPlas 2018). No convincing signal was found. This was expected, given the star’s quiet nature and also the limited ~27 d observation window of TESS, which made the clear visibility of a 25 d signal difficult. We also note that additional photometric data from the T8 automatic photoelectric telescope (APT) did not show any photometric variability related to the stellar rotational period (see Sect. 9.1).

To constrain the stellar rotation, we conducted an analysis of the activity indicator periodogram using the Data and Analysis Center for Exoplanets (DACE; Delisle et al. 2016)$^3$. For this analysis, we utilized the SOPHIE+ S-index to ensure comparability with the S-index measurements obtained from HIRES and HIRES+ instruments. We excluded a total of 68 S-index data points due to their dependency on $S/N$ (<30 in SOPHIE+ order 1, $\lambda \sim 3955 \AA$) and significant contamination caused by the calibration lamp. Additionally, one data point was excluded due to its identification as a $5\sigma$ outlier. In a similar vein, 102 data points from the Na index were omitted due to their reliance on $S/N$ (<70 in SOPHIE+ order 30, $\lambda \sim 5931 \AA$) and their susceptibility to contamination by the telluric lines. Figure 3 displays the periodogram of activity indicators. The HIRES and HIRES+ S-index measurements, along with the SOPHIE+ S-index, reveal periodic signals at 29.6 d (false alarm probability (FAP) <1%, Baluev 2008) and 31.5 d (FAP <10%), respectively. These results are consistent with the estimated rotational period of the star at $25^{+8}_{-6}$ d. Furthermore, there is an activity signal at 141.2 d in the

<footnote>$^3$ Available at https://dace.unige.ch</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Spectral energy distribution of HD 88986. Red symbols represent the observed photometric measurements, whereas the horizontal bars represent the effective width of the passband. Blue symbols are the model fluxes from the best-fit Kurucz atmosphere model (black).</figcaption></figure>

Stassun & Torres (2016); Stassun et al. (2017); Stassun & Torres (2018). We obtained the $B_T V_T$ magnitudes from Tycho-2 (Høg 2001), the JHK$_S$ magnitudes from 2MASS (Cutri et al. 2003), the W1–W4 magnitudes from WISE (Wright et al. 2010), the uvby Strömgren magnitudes from Paunzen (2015), and the $G G_{BP} G_{RP}$ magnitudes from Gaia (Gaia Collaboration 2016, 2021). We also used the UV measurement at 274 nm from the TDI UV satellite (Thompson et al. 1978). Together, the available photometry spans the full stellar SED over the wavelength range 0.2–22 µm (see Fig. 2). Then, we performed a fit using Kurucz’s stellar atmosphere models with the $T_{\text{eff}}$, log g, and [Fe/H] adopted from the spectroscopic analysis. The remaining free parameter is the extinction $A_V$, which we fixed at zero due to the proximity of the system to Earth.

The resulting fit (Fig. 2) has a reduced $\chi^2$ of 1.4. Integrating the model SED gives the bolometric flux at Earth, $F_{\text{bol}} = 7.28 \pm 0.17 \times 10^{-8} \text{ erg s}^{-1} \text{ cm}^{-2}$. Taking the $F_{\text{bol}}$ and $T_{\text{eff}}$ together with the Gaia parallax, gives the stellar radius, $R_\star = 1.543 \pm 0.010 R_\odot$, placing the star within the subgiant range (1.5–3 $R_\odot$; Huber et al. 2017; Berger et al. 2018). In addition, we can estimate the stellar mass from the empirical relations of Torres et al. (2010), giving $M_\star = 1.19 \pm 0.07 M_\odot$, which is consistent with the value of $1.25 \pm 0.05 M_\odot$ determined empirically via $R_\star$ and log g. We acknowledge the possibility that our formal error budget for radius, mass, and stellar age could be underestimated, as suggested by Tayar et al. (2022). For stars such as HD 88986, the systematic uncertainty floor could rather be up to ≈4.2% for radius, ≈5% for mass, and ≈20% for age. Hence, we conservatively adopt these relative uncertainties, reported in brackets in Table 1, to obtain more realistic stellar parameter errors throughout this study.

Finally, we can use the observed log $R'_{\text{HK}}$ (see Sect. 4) to estimate the stellar age via empirical activity-age relations. We obtain an age of $\tau_\star = 7.9 \pm 1.3$ Gyr via the empirical relations of Mamajek & Hillenbrand (2008).

## 4. Stellar rotation and activity

To study the host star’s activity and rotational period, we used log $R'_{\text{HK}}$, S-index, Na index, and CCF bisector measurements from the SOPHIE+ spectrograph, along with the S-index values from the HIRES and HIRES+ spectrographs. The bisector spans

are derived from the SOPHIE DRS using the method described by Queloz et al. (2001). We extracted the Na index as introduced in Da Silva et al. (2011). The HIRES and HIRES+ S-index values were acquired from Butler et al. (2017). For the extraction of log $R'_{\text{HK}}$ and S-index from the SOPHIE+ spectra, we followed the procedure outlined by Noyes et al. (1984) and Boisse et al. (2010), respectively. The key step before deriving the log $R'_{\text{HK}}$ and S-index is subtracting the background contamination light from the Th-Ar or FP calibration lamp from the stellar spectra. To correct this, we used the direct measurement method (see Appendix E).

To estimate the rotational period of the star, we summed 176 HD 88986 SOPHIE+ spectra which fulfilled two criteria. First, the spectra with $S/N > 50$ in the first (bluest, $\lambda \sim 3955$ Å) order of spectra where CaII H&K lines are located. Second, the spectra with minimal contamination due to the background light. This led to the value of log $R'_{\text{HK}} = -5.04 \pm 0.10$. Furthermore, we investigated potential changes in magnetic activity over an approximately 11-yr SOPHIE+ observation span. For this purpose, the dataset was divided into three distinct observational seasons spanning from 2012 to 2023: 2012–2016, 2016–2019, and 2019–2023. We observed a gradual increase in the log $R'_{\text{HK}}$ parameter, with values transitioning from $-5.12 \pm 0.10$ during the initial subset to $-5.05 \pm 0.10$ in the second subset, and ultimately stabilizing at $-5.00 \pm 0.10$ in the final subset. All these values are in good agreement with the value of log $R'_{\text{HK}} = -5.22$ and log $R'_{\text{HK}} = -5.07$ reported by Radick et al. (2018) and Hall et al. (2007), respectively. Finally, we estimated a rotational period of $25^{+8}_{-6}$ d following Noyes et al. (1984), a value consistent with the Mamajek & Hillenbrand (2008) recipe, which yields $26.3 \pm 3.1$ d.

We searched for rotational modulation in the Simple Aperтуre Photometric (SAP; see Twicken et al. 2010; Morris et al. 2020) TESS data (see Sects. 6 and 8.2 for details of the observations) using the Lomb-Scargle periodogram (Lomb 1976; Scargle 1982; VanderPlas 2018). No convincing signal was found. This was expected, given the star’s quiet nature and also the limited ~27 d observation window of TESS, which made the clear visibility of a 25 d signal difficult. We also note that additional photometric data from the T8 automatic photoelectric telescope (APT) did not show any photometric variability related to the stellar rotational period (see Sect. 9.1).

To constrain the stellar rotation, we conducted an analysis of the activity indicator periodogram using the Data and Analysis Center for Exoplanets (DACE; Delisle et al. 2016)$^3$. For this analysis, we utilized the SOPHIE+ S-index to ensure comparability with the S-index measurements obtained from HIRES and HIRES+ instruments. We excluded a total of 68 S-index data points due to their dependency on $S/N$ (<30 in SOPHIE+ order 1, $\lambda \sim 3955$ Å) and significant contamination caused by the calibration lamp. Additionally, one data point was excluded due to its identification as a $5\sigma$ outlier. In a similar vein, 102 data points from the Na index were omitted due to their reliance on $S/N$ (<70 in SOPHIE+ order 30, $\lambda \sim 5931$ Å) and their susceptibility to contamination by the telluric lines. Figure 3 displays the periodogram of activity indicators. The HIRES and HIRES+ S-index measurements, along with the SOPHIE+ S-index, reveal periodic signals at 29.6 d (false alarm probability (FAP) <1%, Baluev 2008) and 31.5 d (FAP <10%), respectively. These results are consistent with the estimated rotational period of the star at $25^{+8}_{-6}$ d. Furthermore, there is an activity signal at 141.2 d in the

<footnote>$^3$ Available at https://dace.unige.ch</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Periodogram of RVs and activity indicators of HD 88986. From top to bottom: HIRES and HIRES+ S-index, SOPHIE+ S-index, bisector, RVs, and residuals of RVs after Keplerian fit on the 146.1 d. The vertical red line illustrates the planet candidates on 146.1 d, which have no corresponding peak in activity indicators. The vertical gray strip marks the estimated rotational period of the star. Also, the horizontal lines show the FAP level of 10, 1, and 0.1%, respectively (Baluev 2008).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** $ℓ_1$ periodogram of SOPHIE+ data. The identified periods are shown in red.</figcaption></figure>

periodogram of the HIRES and HIRES+ S-index, which will be further discussed in the following section.

# 5. Detection of the sub-Neptune HD 88986 b in the SOPHIE+ RVs

## 5.1. RV periodogram

To perform our RV periodogram analysis, we only used SOPHIE+ RV data because it contains an extensive number of data points (378), a long baseline, a higher RV accuracy, and superior sampling compared to other instruments. The analysis employed the default noise model in DACE, assuming an additional Gaussian white noise with nominal error bars of 1.5 m s⁻¹ on the SOPHIE+ RVs. After removing the long-term curvature using a second-order polynomial model, the periodogram of SOPHIE+ RVs showed significant peaks at 104.8 and 146.1 d, falling below the analytical FAP (Baluev 2008) of 0.1%

(see Fig. 3, fourth panel). These two signals are yearly aliases of one another.

To determine the favored alias between two signals, we compute the $ℓ_1$ periodogram. This one takes in a frequency grid and an assumed covariance matrix of the noise as input. It aims to find a representation of the RV time series as a sum of a small number of sinusoids whose frequencies are in the input grid. It outputs a figure that has a similar aspect as a regular periodogram but with fewer peaks due to aliasing. The peaks can be assigned a FAP, whose interpretation is close to the FAP of a regular periodogram peak. The signals found to be statistically significant might vary from one noise model to another. To explore this aspect, as in Hara et al. (2020), we considered several candidate noise models based on the periodicities found in the ancillary indicators. We tried 1200 noise models, all Gaussian, such that the covariance is the sum of a white noise term of amplitude $σ_W$, a red noise term with Gaussian decay of amplitude $σ_R$ and timescale $τ_R$, and a quasi-periodic component (Haywood et al. 2014) with amplitude $σ_{QP}$, timescale $τ_{QP}$ and period $P^*$ and harmonic complexity equal to 1. We tried all combinations with $σ_W, σ_R, σ_{QP} = 0, 0.3, 0.6, 0.9, 1.2, 1.5$ m s⁻¹, $τ_R = 0, 3, 6$ day, $P^* = 29$ d and $τ_{QP} = 20, 40, 60$ d or $P^* = 40$ d and $τ_{QP} = 20, 50, 80$ d. We ranked the models with cross-validation as well as Laplace approximation of the Bayesian evidence. We find that for the 20% highest ranked models, a peak with a period between 145 and 149 d consistently appears. The model with the highest Laplace approximation of the evidence is obtained with $σ_W = σ_{QP} = 1.2$ m s⁻¹, $σ_R = 1.5$ m s⁻¹, $τ_R = 0$ d, $τ_{QP} = 20$ d. The corresponding $ℓ_1$ periodogram is shown in Fig. 4. The highest peak appears at 146.5 d, which is compatible with the 146.1 d signal given the frequency resolution set by the timespan of observations. This signal presents a FAP of $2.45 × 10^{-6}$, which is clearly statistically significant. We, therefore, conclude that the true signal is at 146 d, while the signal at 104.8 d represents its yearly aliases. Other signals appearing on the $ℓ_1$ periodogram are not statistically significant.

Subsequently, we performed a circular Keplerian fit (see Sects. 5.2 and 5.3) to remove the signal at 146.1 d and investigated the resulting RV residuals. The periodogram displayed two signals at 29.2 d and 40.0 d, with FAP values below 10% (see Fig. 3, bottom panel).

To conclude our analysis, we checked that the 146 d signal has a constant phase and amplitude following the methodology of Hara et al. (2022a). It simply consists of computing *g* the phase and amplitude of a signal at a given period with a moving time window. To perform this analysis, it is crucial to have realistic error bars. As visible in Fig. 3, the SOPHIE+ S-index exhibits low-frequency variations that are most likely due to a magnetic cycle. We expect a higher dispersion of RVs when S-index values are high (Borgniet et al. 2015; Meunier 2021;
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Periodogram of RVs and activity indicators of HD 88986. From top to bottom: HIRES and HIRES+ S-index, SOPHIE+ S-index, bisector, RVs, and residuals of RVs after Keplerian fit on the 146.1 d. The vertical red line illustrates the planet candidates on 146.1 d, which have no corresponding peak in activity indicators. The vertical gray strip marks the estimated rotational period of the star. Also, the horizontal lines show the FAP level of 10, 1, and 0.1%, respectively (Baluev 2008).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** $\ell_1$ periodogram of SOPHIE+ data. The identified periods are shown in red.</figcaption></figure>

periodogram of the HIRES and HIRES+ S-index, which will be further discussed in the following section.

# 5. Detection of the sub-Neptune HD 88986 b in the SOPHIE+ RVs

## 5.1. RV periodogram

To perform our RV periodogram analysis, we only used SOPHIE+ RV data because it contains an extensive number of data points (378), a long baseline, a higher RV accuracy, and superior sampling compared to other instruments. The analysis employed the default noise model in DACE, assuming an additional Gaussian white noise with nominal error bars of 1.5 m s⁻¹ on the SOPHIE+ RVs. After removing the long-term curvature using a second-order polynomial model, the periodogram of SOPHIE+ RVs showed significant peaks at 104.8 and 146.1 d, falling below the analytical FAP (Baluev 2008) of 0.1%

(see Fig. 3, fourth panel). These two signals are yearly aliases of one another.

To determine the favored alias between two signals, we compute the $\ell_1$ periodogram. This one takes in a frequency grid and an assumed covariance matrix of the noise as input. It aims to find a representation of the RV time series as a sum of a small number of sinusoids whose frequencies are in the input grid. It outputs a figure that has a similar aspect as a regular periodogram but with fewer peaks due to aliasing. The peaks can be assigned a FAP, whose interpretation is close to the FAP of a regular periodogram peak. The signals found to be statistically significant might vary from one noise model to another. To explore this aspect, as in Hara et al. (2020), we considered several candidate noise models based on the periodicities found in the ancillary indicators. We tried 1200 noise models, all Gaussian, such that the covariance is the sum of a white noise term of amplitude $\sigma_W$, a red noise term with Gaussian decay of amplitude $\sigma_R$ and timescale $\tau_R$, and a quasi-periodic component (Haywood et al. 2014) with amplitude $\sigma_{QP}$, timescale $\tau_{QP}$ and period $P^*$ and harmonic complexity equal to 1. We tried all combinations with $\sigma_W, \sigma_R, \sigma_{QP} = 0, 0.3, 0.6, 0.9, 1.2, 1.5$ m s⁻¹, $\tau_R = 0, 3, 6$ day, $P^* = 29$ d and $\tau_{QP} = 20, 40, 60$ d or $P^* = 40$ d and $\tau_{QP} = 20, 50, 80$ d. We ranked the models with cross-validation as well as Laplace approximation of the Bayesian evidence. We find that for the 20% highest ranked models, a peak with a period between 145 and 149 d consistently appears. The model with the highest Laplace approximation of the evidence is obtained with $\sigma_W = \sigma_{QP} = 1.2$ m s⁻¹, $\sigma_R = 1.5$ m s⁻¹, $\tau_R = 0$ d, $\tau_{QP} = 20$ d. The corresponding $\ell_1$ periodogram is shown in Fig. 4. The highest peak appears at 146.5 d, which is compatible with the 146.1 d signal given the frequency resolution set by the timespan of observations. This signal presents a FAP of $2.45 \times 10^{-6}$, which is clearly statistically significant. We, therefore, conclude that the true signal is at 146 d, while the signal at 104.8 d represents its yearly aliases. Other signals appearing on the $\ell_1$ periodogram are not statistically significant.

Subsequently, we performed a circular Keplerian fit (see Sects. 5.2 and 5.3) to remove the signal at 146.1 d and investigated the resulting RV residuals. The periodogram displayed two signals at 29.2 d and 40.0 d, with FAP values below 10% (see Fig. 3, bottom panel).

To conclude our analysis, we checked that the 146 d signal has a constant phase and amplitude following the methodology of Hara et al. (2022a). It simply consists of computing $g$ the phase and amplitude of a signal at a given period with a moving time window. To perform this analysis, it is crucial to have realistic error bars. As visible in Fig. 3, the SOPHIE+ S-index exhibits low-frequency variations that are most likely due to a magnetic cycle. We expect a higher dispersion of RVs when S-index values are high (Borgniet et al. 2015; Meunier 2021;
<figure><img src="image_6.png" /><figcaption>Fig. 5. Amplitude (red) and phase (green) of a 146.1 d signal as a function of time for different sizes of time windows. Solid lines correspond to estimate and shaded areas to ± 1 σ uncertainties. Denoting by $T_{\text{obs}}$ the total time span of observations, a) and b) are obtained with windows of size $T_{\text{obs}}/3 = 1234.2$ d and $T_{\text{obs}}/9 = 411.42$ d, respectively.</figcaption></figure>

Hara & Delisle 2023). Following Díaz et al. (2016), we added a white noise jitter term, as well as a jitter scaled with the value of the log $R'_{HK}$, and fit those along with a polynomial line and the 146 d signal. We used those values to compute the amplitude and phase as a function of time, shown in Fig. 5. We here consider windows of size $T_{\text{obs}}/3 = 1211$ days and $T_{\text{obs}}/9 = 411$ days where $T_{\text{obs}}$ is the total time span of observations. The most notable feature is the drop in amplitude in Fig. 5b at BJD 2 459 000–2 460 000. We performed the quantitative analysis presented in Hara et al. (2022a) and determined that the phase and amplitudes are consistent with being constant. The drop in amplitude is likely due to the fact that the epoch BJD 2 459 000–2 460 000 corresponds to an activity maximum, and it is possible that the activity pattern changes and masks the planetary signal in the corresponding observational seasons.

To summarize, both periodograms exhibit a strong periodicity of about 146 d in RVs, which falls within a different period than the estimated star’s rotation period ($25^{+8}_{-6}$ d; see Sect. 4). Notably, the periodogram of the HIRES and HIRES+ S-index reveals an activity signal at 141.2 d, which differs from the RV signal by five days and has a relatively low strength (FAP ~ 10%). None of the SOPHIE+ activity indicators showed a periodicity at 146.1 d. Moreover, no correlation was found between RV residuals after removing the trend and both S-index (Pearson’s coefficient $R = 0.03$) and bisector (Pearson’s coefficient $R = 0.13$). Furthermore, our analysis above demonstrates consistent phase and amplitude of the 146.1 d signal. Consequently, it is unlikely that RV periodicity at 146.1 d has an activity origin. Additionally, it is noteworthy that a survey encompassing all SOPHIE constant stars, as well as other stars observed by SOPHIE, did not reveal the same RV periodicity, confirming that the observed periodicity is not due to instrumental artifacts. Therefore, the RV signal at 146.1 d is likely to have a planetary origin. Throughout the rest of the paper, we attribute this periodicity to the planet HD 88986 b.

Furthermore, given the two periodic signals of 31.5 d in the SOPHIE+ S-index and 29.6 d in the HIRES and HIRES+ S-index, along with the estimated star rotation period ($25^{+8}_{-6}$ d), the RV residuals signal at 29.2 d is likely the result of stellar rotational modulation. We take into account this signal in our analysis for the rest of the paper, testing three different

methods to model it (see below in Sect. 5.2). Finally, since the RV residuals signal at 40.0 d is statistically insignificant, further observations are required to determine whether this signal has an astrophysical origin or is simply noise.

We note that among the RVs presented in this paper, only those obtained with SOPHIE+ allow the low-amplitude signal of planet HD 88986 b to be detected. This is due both to the large number of measurements and to their high accuracy. ELODIE, SOPHIE, HIRES, and HIRES+ RVs do not assist here in the discovery of that planet. Therefore, to avoid any offsets or potential systematics among instruments, we only use the SOPHIE+ data in the fits of HD 88986 b presented below in the continuation of this section and Sect. 7. The other RVs datasets are only included in Sect. 9.2 for constraining the outer companion.

## 5.2. RV models

To adequately describe the data, we took into account several SOPHIE+ RV-only models and performed model comparisons. The RV analysis was carried out by juliet (Espinoza et al. 2019), which employs radvel (Fulton et al. 2018) to model RVs and george (Ambikasaran et al. 2015) and celerite (Foreman-Mackey et al. 2017) to model potential activity effects on the data through Gaussian process regression (GPs). For each tested model, juliet computes the Bayesian log evidence (ln Z). If the Bayesian log-evidence difference ($\Delta \ln Z$) between a model and another exceeds two, it is moderately favored over the latter; a difference greater than five indicates strong favorability (Trotta 2008). The models are indistinguishable when the difference in Bayesian log evidence is $\Delta \ln Z \le 2$. In this case, the model with the fewest free parameters would be chosen.

We defined our RV model as follows:

$$M(t) = K(t) + \epsilon(t) + \mu + Q(t^2) + A(t), \quad (1)$$

where $K(t)$ is the Keplerian model and the $\epsilon(t) \sim N(0, \sigma(t)^2 + \sigma_w^2)$, is a normal distribution ($N$) of white-Gaussian noise where $\sigma(t)$ is the uncertainty of each RV point at time $t$, and $\sigma_w$ is a jitter term. Additionally, $\mu$ is a systematic RV offset of the instrument, and $Q$ and $A$ are defined as quadratic and linear terms, respectively, to model the long-term curvature. The model
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Amplitude (red) and phase (green) of a 146.1 d signal as a function of time for different sizes of time windows. Solid lines correspond to estimate and shaded areas to ± 1 σ uncertainties. Denoting by $T_{\text{obs}}$ the total time span of observations, a) and b) are obtained with windows of size $T_{\text{obs}}/3 = 1234.2$ d and $T_{\text{obs}}/9 = 411.42$ d, respectively.</figcaption></figure>

Hara & Delisle 2023). Following Díaz et al. (2016), we added a white noise jitter term, as well as a jitter scaled with the value of the log $R'_{HK}$, and fit those along with a polynomial line and the 146 d signal. We used those values to compute the amplitude and phase as a function of time, shown in Fig. 5. We here consider windows of size $T_{\text{obs}}/3 = 1211$ days and $T_{\text{obs}}/9 = 411$ days where $T_{\text{obs}}$ is the total time span of observations. The most notable feature is the drop in amplitude in Fig. 5b at BJD 2 459 000–2 460 000. We performed the quantitative analysis presented in Hara et al. (2022a) and determined that the phase and amplitudes are consistent with being constant. The drop in amplitude is likely due to the fact that the epoch BJD 2 459 000–2 460 000 corresponds to an activity maximum, and it is possible that the activity pattern changes and masks the planetary signal in the corresponding observational seasons.

To summarize, both periodograms exhibit a strong periodicity of about 146 d in RVs, which falls within a different period than the estimated star’s rotation period ($25^{+8}_{-6}$ d; see Sect. 4). Notably, the periodogram of the HIRES and HIRES+ S-index reveals an activity signal at 141.2 d, which differs from the RV signal by five days and has a relatively low strength (FAP ~ 10%). None of the SOPHIE+ activity indicators showed a periodicity at 146.1 d. Moreover, no correlation was found between RV residuals after removing the trend and both S-index (Pearson’s coefficient $R = 0.03$) and bisector (Pearson’s coefficient $R = 0.13$). Furthermore, our analysis above demonstrates consistent phase and amplitude of the 146.1 d signal. Consequently, it is unlikely that RV periodicity at 146.1 d has an activity origin. Additionally, it is noteworthy that a survey encompassing all SOPHIE constant stars, as well as other stars observed by SOPHIE, did not reveal the same RV periodicity, confirming that the observed periodicity is not due to instrumental artifacts. Therefore, the RV signal at 146.1 d is likely to have a planetary origin. Throughout the rest of the paper, we attribute this periodicity to the planet HD 88986 b.

Furthermore, given the two periodic signals of 31.5 d in the SOPHIE+ S-index and 29.6 d in the HIRES and HIRES+ S-index, along with the estimated star rotation period ($25^{+8}_{-6}$ d), the RV residuals signal at 29.2 d is likely the result of stellar rotational modulation. We take into account this signal in our analysis for the rest of the paper, testing three different

methods to model it (see below in Sect. 5.2). Finally, since the RV residuals signal at 40.0 d is statistically insignificant, further observations are required to determine whether this signal has an astrophysical origin or is simply noise.

We note that among the RVs presented in this paper, only those obtained with SOPHIE+ allow the low-amplitude signal of planet HD 88986 b to be detected. This is due both to the large number of measurements and to their high accuracy. ELODIE, SOPHIE, HIRES, and HIRES+ RVs do not assist here in the discovery of that planet. Therefore, to avoid any offsets or potential systematics among instruments, we only use the SOPHIE+ data in the fits of HD 88986 b presented below in the continuation of this section and Sect. 7. The other RVs datasets are only included in Sect. 9.2 for constraining the outer companion.

## 5.2. RV models

To adequately describe the data, we took into account several SOPHIE+ RV-only models and performed model comparisons. The RV analysis was carried out by juliet (Espinoza et al. 2019), which employs radvel (Fulton et al. 2018) to model RVs and george (Ambikasaran et al. 2015) and celerite (Foreman-Mackey et al. 2017) to model potential activity effects on the data through Gaussian process regression (GPs). For each tested model, juliet computes the Bayesian log evidence (ln Z). If the Bayesian log-evidence difference ($\Delta \ln Z$) between a model and another exceeds two, it is moderately favored over the latter; a difference greater than five indicates strong favorability (Trotta 2008). The models are indistinguishable when the difference in Bayesian log evidence is $\Delta \ln Z \le 2$. In this case, the model with the fewest free parameters would be chosen.

We defined our RV model as follows:

$$M(t) = K(t) + \epsilon(t) + \mu + Q(t^2) + A(t), \quad (1)$$

where $K(t)$ is the Keplerian model and the $\epsilon(t) \sim N(0, \sigma(t)^2 + \sigma_w^2)$, is a normal distribution ($N$) of white-Gaussian noise where $\sigma(t)$ is the uncertainty of each RV point at time $t$, and $\sigma_w$ is a jitter term. Additionally, $\mu$ is a systematic RV offset of the instrument, and $Q$ and $A$ are defined as quadratic and linear terms, respectively, to model the long-term curvature. The model
<table><caption>**Table 2.** Different tested models on the SOPHIE+ RV-only data along with model comparisons with juliet.</caption><thead><tr><th>Model</th><th>Δ ln Z</th><th>period (d)</th><th>K (m s<sup>-1</sup>)</th><th>T<sub>c</sub> (BJD-2400000 d)</th></tr></thead><tbody><tr><td>No planet</td><td>-97.0</td><td>-</td><td>-</td><td>-</td></tr><tr><td>No planet+EXP-GP</td><td>10.2</td><td>-</td><td>-</td><td>-</td></tr><tr><td>No planet+QP-GP</td><td>17.5</td><td>-</td><td>-</td><td>-</td></tr><tr><td>1Co</td><td>0</td><td>146.0 ± 0.5</td><td>1.7 ± 0.2</td><td>58 897 ± 3</td></tr><tr><td>1eccentricity-free</td><td>1.6</td><td>145.8<sup>+0.4</sup><sub>-0.3</sub></td><td>1.9 ± 0.3</td><td>58 884<sup>+8</sup><sub>-7</sub></td></tr><tr><td>1Co+sinusoidal</td><td>5.3</td><td>146.2 ± 0.4</td><td>1.6 ± 0.2</td><td>58 899 ± 3</td></tr><tr><td>1Co+EXP-GP</td><td>20.0</td><td>146.3 ± 0.6</td><td>1.6 ± 0.4</td><td>58 897<sup>+5</sup><sub>-6</sub></td></tr><tr><td>1Co+QP-GP</td><td>22.3</td><td>146.3 ± 0.6</td><td>1.7 ± 0.4</td><td>58 897 ± 6</td></tr></tbody></table>

**Notes.** K refers to the RV semi-amplitude and T_c is the time of mid-transit from RV fit. Additionally, in the model names, “Co” represents a circular orbit. The final choice model is indicated in bold.

is tested for both scenarios of the eccentricity-free and circular orbit (Co). To explore the possible effect of stellar activity on the planet's parameters, we model the stellar activity in three different ways:

- a sinusoidal orbit;

- an exponential GP kernel (EXP-GP) with the form of $k_{i,j} = \sigma_{GP}^2 \exp(-|t_i - t_j|/T_{GP})$, where $\sigma_{GP}$ is the amplitude of GP modulation, and $T_{GP}$ presents the characteristic timescale (Ambikasaran et al. 2015);

- a quasi-periodic GP kernel (QP-GP, Foreman-Mackey et al. 2017) with the form of $\kappa_{i,j} = B_{GP}/(2 + C_{GP}) e^{(-|t_i - t_j|/L_{GP})}(\cos(2\pi|t_i - t_j|/P_{rot;GP}) + (1 + C_{GP}))$, where $B_{GP}$ amplifies the kernel, $C_{GP}$ is a constant scaling term, $L_{GP}$ is a correlation time-scale component, and finally $P_{rot;GP}$ is the rotational modulation.

In addition to the models above, we also executed three no-planet models wherein we assumed the absence of any planetary signal in the RVs (i.e., in Eq. (1): $K(t) = 0$). Table 2 provides a summary of the results obtained from testing various models. Furthermore, Table F.1 presents the priors employed in the analysis, along with detailed descriptions of all parameters. Throughout this paper, we conducted juliet runs for each model using a consistent setup, employing a configuration with 3× number of free parameters as the number of walkers, executed 10 000 steps per walker, and discarded the initial 3000 steps as burn-in.

### 5.3. Detection of the sub-Neptune HD 88986 b

We set a uniform prior to the planetary period between 135 d and 155 d. For the mid-transit time, we applied a uniform prior defined by a time window of 146 d which is consistent with the planetary period duration. For the other parameters, we used fairly broad uniform priors. The results of the no-planet model as well as the circular and eccentric fits are shown in the Table 2. We note that this table includes $T_c$, the derived center time of the inferior conjunction, as a transit of this planet is reported below in Sect. 6. The results of the circular and eccentric models are consistent and both are statistically significant compared to the no-planet model ($\Delta\ln Z \ge 97$) which confirms a clear detection of the planetary signal. The eccentricity-free model, however, exhibits bimodality on the periastron argument $\omega$. Given that the two models are statistically indistinguishable ($\Delta\ln Z \le 2$), and their posterior distributions are consistent, we continued to model the HD 88986 b planet with a circular orbit as it has fewer free parameters.

The last three lines of Table 2 present the results of the fits adopting the three different models for the stellar activity (see Sect. 5.2). In the sinusoidal model, we employed a uniform prior for the activity period ranging from 25 to 35 d, as we detected a potential stellar signal at 29 d (see Sect. 4). Regarding the QP-GP hyperparameter, we imposed some constraints. On $P_{rot;GP}$, we used a Gaussian prior centered at 29 d with a standard deviation of 3 d. On the $C_{GP}$, we initially tried a wide range of Jeffreys priors from $10^{-20}$ to 100. The posterior distribution of $C_{GP}$ reached the prior boundary at $10^{-20}$, indicating that this parameter converged to zero. This result motivates us to set the $C_{GP}$ to $10^{-20}$, which is consistent with zero. We note that the models with fixed/unfixed $C_{GP}$ are statistically indistinguishable ($\Delta\ln Z \le 2$). Therefore, we fixed $C_{GP}$ and continued using the QP-GP model with three free parameters. A broad uniform prior is taken into account for all remaining parameters of different models, as presented in Table F.1.

The results of the five different planet models reported in Table 2 are consistent. This strongly argues in favor of the detection of the planet HD 88986 b with those parameters. Notably, all planet models accounting for the potential activity signal are strongly favored statistically ($\Delta\ln Z \ge 5$) when compared to models that do not consider the stellar activity. Among these models, the planet models with simultaneous GP kernels (1Co+EXP-GP and 1Co+QP-GP) are strongly favored over the sinusoidal model (1Co+ sinusoidal). Additionally, the model incorporating the simultaneous QP-GP kernel displays moderate statistical favorability compared to the simultaneous EXP-GP kernel. We note that for the 1Co+sinusoidal model, we have a bimodality in the posterior distribution of the stellar rotation period at 29 d and a much smaller peak at about 31 d, which might be explained by differential stellar rotation.

In addition to the different planetary models and one no-planet model explored above, we extended our analysis to encompass two GP-only (EXP-GP and QP-GP) models (see Table 2). These models exhibit a higher statistical strength than those when considering the planet alone (1Co and 1eccentricity-free), suggesting that the presence of RV variabilities is primarily driven by dominant stellar activity signals rather than a planetary influence. However, a detailed examination outlined in Table 2 reveals a strong preference: models incorporating both the planet candida with the GP components are strongly favored ($\Delta\ln Z \ge 5$) over GP-only models. This robust preference, the consistency of planet parameters across various models, and the detection of a strong periodicity in our periodogram analysis (refer to Sect. 5), coupled with the consistent phase and amplitude of the signal (refer to Sect. 5), ensures the credibility of our detection.

Finally, for the sake of completeness, we incorporated all RV data presented in this paper, with the exception of ELODIE, owing to its significantly larger error bars (approximately 8 times larger) in comparison to other datasets. This additional model aimed to assess the consistency of the remaining RV data with the results obtained from SOPHIE+ RVs. To account for the long-term curvature in the data, we applied a two-Keplerian model. Our analysis yielded consistent results (Period = 146.9 ± 0.3 d, $K = 1.7 \pm 0.2$ m s⁻¹, $T_c = 58 903 \pm 3$) when compared with other models listed in Table 2. After removing the second Keplerian orbit, the RV RMS values were 5.5 m s⁻¹, 5.6 m s⁻¹, 6.2 m s⁻¹, and 3.2 m s⁻¹ for HIRES, HIRES+, SOPHIE, and SOPHIE+ respectively. These RMS values, coupled with the limited number of available data points and sparse sampling of data from instruments other than SOPHIE+, indicate the challenges faced in detecting signals of such shallow amplitude in the other datasets beyond SOPHIE+. Additionally, this reaffirms the
<table><caption>**Table 2.** Different tested models on the SOPHIE+ RV-only data along with model comparisons with juliet.</caption><thead><tr><th>Model</th><th>Δ ln Z</th><th>period (d)</th><th>K (m s<sup>-1</sup>)</th><th>T<sub>c</sub> (BJD-2400000 d)</th></tr></thead><tbody><tr><td>No planet</td><td>-97.0</td><td>-</td><td>-</td><td>-</td></tr><tr><td>No planet+EXP-GP</td><td>10.2</td><td>-</td><td>-</td><td>-</td></tr><tr><td>No planet+QP-GP</td><td>17.5</td><td>-</td><td>-</td><td>-</td></tr><tr><td>1Co</td><td>0</td><td>146.0 ± 0.5</td><td>1.7 ± 0.2</td><td>58 897 ± 3</td></tr><tr><td>1eccentricity-free</td><td>1.6</td><td>145.8<sup>+0.4</sup><sub>-0.3</sub></td><td>1.9 ± 0.3</td><td>58 884<sup>+8</sup><sub>-7</sub></td></tr><tr><td>1Co+sinusoidal</td><td>5.3</td><td>146.2 ± 0.4</td><td>1.6 ± 0.2</td><td>58 899 ± 3</td></tr><tr><td>1Co+EXP-GP</td><td>20.0</td><td>146.3 ± 0.6</td><td>1.6 ± 0.4</td><td>58 897<sup>+5</sup><sub>-6</sub></td></tr><tr><td>1Co+QP-GP</td><td>22.3</td><td>146.3 ± 0.6</td><td>1.7 ± 0.4</td><td>58 897 ± 6</td></tr></tbody></table>

**Notes.** K refers to the RV semi-amplitude and T<sub>c</sub> is the time of mid-transit from RV fit. Additionally, in the model names, “Co” represents a circular orbit. The final choice model is indicated in bold.

is tested for both scenarios of the eccentricity-free and circular orbit (Co). To explore the possible effect of stellar activity on the planet's parameters, we model the stellar activity in three different ways:

- a sinusoidal orbit;

- an exponential GP kernel (EXP-GP) with the form of $k_{i,j} = \sigma_{GP}^2 \exp(-|t_i - t_j|/T_{GP})$, where $\sigma_{GP}$ is the amplitude of GP modulation, and $T_{GP}$ presents the characteristic timescale (Ambikasaran et al. 2015);

- a quasi-periodic GP kernel (QP-GP, Foreman-Mackey et al. 2017) with the form of $\kappa_{i,j} = B_{GP}/(2 + C_{GP}) e^{(-|t_i - t_j|/L_{GP})}(\cos(2\pi|t_i - t_j|/P_{rot;GP}) + (1 + C_{GP}))$, where $B_{GP}$ amplifies the kernel, $C_{GP}$ is a constant scaling term, $L_{GP}$ is a correlation time-scale component, and finally $P_{rot;GP}$ is the rotational modulation.

In addition to the models above, we also executed three no-planet models wherein we assumed the absence of any planetary signal in the RVs (i.e., in Eq. (1): $K(t) = 0$). Table 2 provides a summary of the results obtained from testing various models. Furthermore, Table F.1 presents the priors employed in the analysis, along with detailed descriptions of all parameters. Throughout this paper, we conducted juliet runs for each model using a consistent setup, employing a configuration with 3× number of free parameters as the number of walkers, executed 10 000 steps per walker, and discarded the initial 3000 steps as burn-in.

### 5.3. Detection of the sub-Neptune HD 88986 b

We set a uniform prior to the planetary period between 135 d and 155 d. For the mid-transit time, we applied a uniform prior defined by a time window of 146 d which is consistent with the planetary period duration. For the other parameters, we used fairly broad uniform priors. The results of the no-planet model as well as the circular and eccentric fits are shown in the Table 2. We note that this table includes $T_c$, the derived center time of the inferior conjunction, as a transit of this planet is reported below in Sect. 6. The results of the circular and eccentric models are consistent and both are statistically significant compared to the no-planet model ($\Delta\ln Z \ge 97$) which confirms a clear detection of the planetary signal. The eccentricity-free model, however, exhibits bimodality on the periastron argument $\omega$. Given that the two models are statistically indistinguishable ($\Delta\ln Z \le 2$), and their posterior distributions are consistent, we continued to model the HD 88986 b planet with a circular orbit as it has fewer free parameters.

The last three lines of Table 2 present the results of the fits adopting the three different models for the stellar activity (see Sect. 5.2). In the sinusoidal model, we employed a uniform prior for the activity period ranging from 25 to 35 d, as we detected a potential stellar signal at 29 d (see Sect. 4). Regarding the QP-GP hyperparameter, we imposed some constraints. On $P_{rot;GP}$, we used a Gaussian prior centered at 29 d with a standard deviation of 3 d. On the $C_{GP}$, we initially tried a wide range of Jeffreys priors from $10^{-20}$ to 100. The posterior distribution of $C_{GP}$ reached the prior boundary at $10^{-20}$, indicating that this parameter converged to zero. This result motivates us to set the $C_{GP}$ to $10^{-20}$, which is consistent with zero. We note that the models with fixed/unfixed $C_{GP}$ are statistically indistinguishable ($\Delta\ln Z \le 2$). Therefore, we fixed $C_{GP}$ and continued using the QP-GP model with three free parameters. A broad uniform prior is taken into account for all remaining parameters of different models, as presented in Table F.1.

The results of the five different planet models reported in Table 2 are consistent. This strongly argues in favor of the detection of the planet HD 88986 b with those parameters. Notably, all planet models accounting for the potential activity signal are strongly favored statistically ($\Delta\ln Z \ge 5$) when compared to models that do not consider the stellar activity. Among these models, the planet models with simultaneous GP kernels (1Co+EXP-GP and 1Co+QP-GP) are strongly favored over the sinusoidal model (1Co+ sinusoidal). Additionally, the model incorporating the simultaneous QP-GP kernel displays moderate statistical favorability compared to the simultaneous EXP-GP kernel. We note that for the 1Co+sinusoidal model, we have a bimodality in the posterior distribution of the stellar rotation period at 29 d and a much smaller peak at about 31 d, which might be explained by differential stellar rotation.

In addition to the different planetary models and one no-planet model explored above, we extended our analysis to encompass two GP-only (EXP-GP and QP-GP) models (see Table 2). These models exhibit a higher statistical strength than those when considering the planet alone (1Co and 1eccentricity-free), suggesting that the presence of RV variabilities is primarily driven by dominant stellar activity signals rather than a planetary influence. However, a detailed examination outlined in Table 2 reveals a strong preference: models incorporating both the planet candida with the GP components are strongly favored ($\Delta\ln Z \ge 5$) over GP-only models. This robust preference, the consistency of planet parameters across various models, and the detection of a strong periodicity in our periodogram analysis (refer to Sect. 5), coupled with the consistent phase and amplitude of the signal (refer to Sect. 5), ensures the credibility of our detection.

Finally, for the sake of completeness, we incorporated all RV data presented in this paper, with the exception of ELODIE, owing to its significantly larger error bars (approximately 8 times larger) in comparison to other datasets. This additional model aimed to assess the consistency of the remaining RV data with the results obtained from SOPHIE+ RVs. To account for the long-term curvature in the data, we applied a two-Keplerian model. Our analysis yielded consistent results (Period = 146.9 ± 0.3 d, $K = 1.7 \pm 0.2$ m s⁻¹, $T_c = 58 903 \pm 3$) when compared with other models listed in Table 2. After removing the second Keplerian orbit, the RV RMS values were 5.5 m s⁻¹, 5.6 m s⁻¹, 6.2 m s⁻¹, and 3.2 m s⁻¹ for HIRES, HIRES+, SOPHIE, and SOPHIE+ respectively. These RMS values, coupled with the limited number of available data points and sparse sampling of data from instruments other than SOPHIE+, indicate the challenges faced in detecting signals of such shallow amplitude in the other datasets beyond SOPHIE+. Additionally, this reaffirms the
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** TESS observation of HD 88986 in sector 21 in 2020 February. Top: normalized TESS PDC-SAP light curve of sector 21 (black dots) along with the best-fit trend (red curve) to the data. The green vertical lines represent the times of the spacecraft's momentum dumps. Middle: normalized re-extracted light curve of sector 21 (black dots). See the text for more information. Bottom: final detrended light curve. The expected HD 88986 b transit event from SOPHIE+ RVs (Sect. 5.3), with 1 sigma uncertainties, is highlighted in blue, and a single transit event is found in the TESS photometric data within this region.</figcaption></figure>

right selection of only SOPHIE+ data for presenting the detection and characterization of this low-mass planet, underscoring the accuracy of our results and ensuring the absence of additional instrumental offsets among the instruments. It is pertinent to highlight that satisfactory convergence was not achieved for certain parameters of the second Keplerian orbit, such as eccentricity and $\omega$, as the orbit of this massive companion remains incomplete (see Sect. 9).

To conclude, the results of all investigated models agree within the error bars. Accordingly, we chose the 1Co+ QP-GP model on SOPHIE+ data as our final choice because it is strongly favored statistically over the model without GP and has moderate favorability compared to the Co+ EXP-GP model.

## 6. First transit detected in TESS sector 21 in February 2020

HD 88986 was observed in TESS sector 21⁴ with camera 1 in a 2-minute cadence from 2020 January 21 to 2020 February 18. The photometric data were produced by the Presearch Data Conditioning-Simple Aperture Photometry (PDC-SAP) pipeline (Stumpe et al. 2012, 2014; Smith et al. 2012), provided by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) at NASA Ames Research Center. The normalized raw TESS photometric data are plotted in the top panel of Fig. 6.

TESS data from sector 21 revealed a potential single transit candidate with $T_c = 2.458\,891.6$ (corresponding to 2020 February 12), a duration of about 16 h, and a depth of ~220 ppm. Remarkably, this $T_c$ is in agreement, within uncertainty, with all the $T_c$ values predicted above from the RV fits of HD 88986 b (Sect. 5.3, Table 2). We note that, neither the TESS SPOC nor Quick Look pipelines (QLP) detected this transit signature, as they require at least two transits to generate a Threshold Crossing Event (TCE) that would be vetted by the TESS Science Office.

To investigate whether that potential single, shallow transit in sector 21 is a false positive scenario, we performed a test by calculating the mean in-transit and mean out-of-transit flux, along

with the difference between them (see Fig. 6). This approach allows us to examine the offset between the different image positions and the actual position of the target star, providing valuable insights into false positive scenarios. While interpreting the difference images from saturated stars like HD 88986 is particularly challenging, we observed that most of the energy in the transit feature is associated with the upper end of the bleed of the saturated pixels in the core of the stellar image. Therefore, it is likely that the transit feature is indeed associated with the host star.

Additionally, as Fig. 6 shows, a telescope reaction wheel momentum dump occurred during the transiting event (Fausnaugh et al. 2020). The sector 21 light curve seems to have been only minimally impacted by the other momentum dumps occurring within this sector. This suggests the robustness of the applied momentum dump correction. However, to ensure that the signal is not the cause of a pipeline's imperfect momentum correction and produces an ingress-like feature, we re-extracted the light curve.

To do so, we use the TPFED/FFIED tool (hereafter TPFED) recently developed by Wilson et al. (2023) to conduct custom extractions of the TESS sector 21 data using the calibrated target pixel files (TPFs) with the default quality bitmask. In brief, we extracted target fluxes for a range of custom aperture masks created with radii of two to four pixels in steps of 0.1 pixels centered on the target. It should be noted that as the target does not fall in the exact center of a pixel increasing the aperture mask radius by 0.1 pixels can result in unique noncircular masks. All produced light curves were background-corrected after determining the sky level using custom background masks. We then detrended the data using two methods. Firstly, we conducted Principal Component Analyses (PCA) on the pixel values within our custom background masks to determine the scattered-light flux contribution to the light curves and then removed these systematics by using the five prime principal components as basis vectors in a linear model. Secondly, we corrected flux modulation due to spacecraft jitter by retrieving the co-trending basis vectors (CBVs), and two-second cadence engineering quaternion measurements for the cameras that observed HD 88986. Following the method used in Delrez et al. (2021), we computed the means and averages of the quaternions over the scientific observational cadences and subsequently used these

<footnote n="4">⁴ https://heasarc.gsfc.nasa.gov/cgi-bin/tess/webtess/wtv.py</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** TESS observation of HD 88986 in sector 21 in 2020 February. Top: normalized TESS PDC-SAP light curve of sector 21 (black dots) along with the best-fit trend (red curve) to the data. The green vertical lines represent the times of the spacecraft's momentum dumps. Middle: normalized re-extracted light curve of sector 21 (black dots). See the text for more information. Bottom: final detrended light curve. The expected HD 88986 b transit event from SOPHIE+ RVs (Sect. 5.3), with 1 sigma uncertainties, is highlighted in blue, and a single transit event is found in the TESS photometric data within this region.</figcaption></figure>

right selection of only SOPHIE+ data for presenting the detection and characterization of this low-mass planet, underscoring the accuracy of our results and ensuring the absence of additional instrumental offsets among the instruments. It is pertinent to highlight that satisfactory convergence was not achieved for certain parameters of the second Keplerian orbit, such as eccentricity and $\omega$, as the orbit of this massive companion remains incomplete (see Sect. 9).

To conclude, the results of all investigated models agree within the error bars. Accordingly, we chose the 1Co+ QP-GP model on SOPHIE+ data as our final choice because it is strongly favored statistically over the model without GP and has moderate favorability compared to the Co+ EXP-GP model.

## 6. First transit detected in TESS sector 21 in February 2020

HD 88986 was observed in TESS sector 21⁴ with camera 1 in a 2-minute cadence from 2020 January 21 to 2020 February 18. The photometric data were produced by the Presearch Data Conditioning-Simple Aperture Photometry (PDC-SAP) pipeline (Stumpe et al. 2012, 2014; Smith et al. 2012), provided by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) at NASA Ames Research Center. The normalized raw TESS photometric data are plotted in the top panel of Fig. 6.

TESS data from sector 21 revealed a potential single transit candidate with $T_c = 2.458\,891.6$ (corresponding to 2020 February 12), a duration of about 16 h, and a depth of ~220 ppm. Remarkably, this $T_c$ is in agreement, within uncertainty, with all the $T_c$ values predicted above from the RV fits of HD 88986 b (Sect. 5.3, Table 2). We note that, neither the TESS SPOC nor Quick Look pipelines (QLP) detected this transit signature, as they require at least two transits to generate a Threshold Crossing Event (TCE) that would be vetted by the TESS Science Office.

To investigate whether that potential single, shallow transit in sector 21 is a false positive scenario, we performed a test by calculating the mean in-transit and mean out-of-transit flux, along

with the difference between them (see Fig. G.1). This approach allows us to examine the offset between the different image positions and the actual position of the target star, providing valuable insights into false positive scenarios. While interpreting the difference images from saturated stars like HD 88986 is particularly challenging, we observed that most of the energy in the transit feature is associated with the upper end of the bleed of the saturated pixels in the core of the stellar image. Therefore, it is likely that the transit feature is indeed associated with the host star.

Additionally, as Fig. 6 shows, a telescope reaction wheel momentum dump occurred during the transiting event (Fausnaugh et al. 2020). The sector 21 light curve seems to have been only minimally impacted by the other momentum dumps occurring within this sector. This suggests the robustness of the applied momentum dump correction. However, to ensure that the signal is not the cause of a pipeline's imperfect momentum correction and produces an ingress-like feature, we re-extracted the light curve.

To do so, we use the TPFED/FFIED tool (hereafter TPFED) recently developed by Wilson et al. (2023) to conduct custom extractions of the TESS sector 21 data using the calibrated target pixel files (TPFs) with the default quality bitmask. In brief, we extracted target fluxes for a range of custom aperture masks created with radii of two to four pixels in steps of 0.1 pixels centered on the target. It should be noted that as the target does not fall in the exact center of a pixel increasing the aperture mask radius by 0.1 pixels can result in unique noncircular masks. All produced light curves were background-corrected after determining the sky level using custom background masks. We then detrended the data using two methods. Firstly, we conducted Principal Component Analyses (PCA) on the pixel values within our custom background masks to determine the scattered-light flux contribution to the light curves and then removed these systematics by using the five prime principal components as basis vectors in a linear model. Secondly, we corrected flux modulation due to spacecraft jitter by retrieving the co-trending basis vectors (CBVs), and two-second cadence engineering quaternion measurements for the cameras that observed HD 88986. Following the method used in Delrez et al. (2021), we computed the means and averages of the quaternions over the scientific observational cadences and subsequently used these

<footnote n="4">⁴ https://heasarc.gsfc.nasa.gov/cgi-bin/tess/webtess/wtv.py</footnote>
vectors along with the CBVs to remove any flux trends. The final light curve is presented in the middle panel of Fig. 6, and the potential single transit is clearly seen within the dataset.

In a quest to further explore possible sources of instrumental noise that could impact the detection or shape of the potential, shallow single transit, we explored alternative methodologies following Rapetti et al. (in preparation). We used the adaptation of the technique Pixel Level Decorrelation (hereafter PLD; Deming et al. 2015; Luger et al. 2016, 2018) implemented in the PLDCorrector class of the community Python package Lightkurve<fnref n="5" />. This method employs (i) a spline polynomial fit to describe stellar variability, (ii) PCA eigenmodes to model the background light, and (iii) the PLD technique to account for pointing and mechanical effects. As an additional approach, we employed a version of PDC as adapted in the CBVCorrector class of Lightkurve, utilizing the CBV technique that the PDC method of the SPOC pipeline employs (hereafter we refer to this corrector as CBV).

Before applying the PLD, we added the background flux and errors estimated by the TESS SPOC pipeline back onto the SAP light curve. Flux level, fraction, and crowding adjustments are applied to the corrected light curves. To automatically optimize the selection of parameter values for the correctors, we evaluate the resulting light curve using the Savitsky-Golay Combined Differential Photometric Precision (sgCDPP) proxy algorithm (Gilliland et al. 2011; Van Cleve et al. 2016) implemented in Lightkurve, for durations of 30, 60, 120, 160, and 200 min (see the legend of Fig. G.2). For a grid of corrector parameter values (for further details on the parameters and the grid, see Rapetti et al., in prep.), we calculated the harmonic mean (HM) of these quantities and selected the corrected light curve that minimizes the HM. In addition, the final sgCDPP metrics can be compared to those obtained for the SPOC PDC-SAP corrected light curve (see Fig. G.2). For this comparison, we also calculated the over-fitting metric implemented in Lightkurve (see Fig. G.2) to measure the broad-band power spectrum via a Lomb-Scargle periodogram and assess the level of introduced noise in the corrected light curves.

In addition to the methods detailed above for extracting data from the 2-min TPFs, we extracted the full frame image (FFI) light curve from the TESS image using a strategy similar to Vanderburg et al. (2016). In particular, we created 20 different photometric apertures, 10 circular apertures, and 10 shaped like the TESS point spread function at the star location on the detector. We then calculated light curves from each aperture and corrected for systematics by performing a linear least-squares fit modeling on the light curve with time series the mean, standard deviation, kurtosis, and skew of the spacecraft quaternion measurements within each exposure (e.g., Vanderburg et al. 2019), the SPOC PDC CBVs, the background flux time series, and a basis spline to model slow variability. We performed the least-squares linear fit iteratively, removing outliers until convergence. Once the light curves were corrected for systematics, we corrected for dilution from other nearby stars and identified the one with the best photometric precision, which we used for our FFI analysis. The final resulting light curve is plotted in Fig. G.3.

The TESS light curves of sector 21, produced using different approaches, all find a feature where the potential single, shallow transit was identified (see Figs. 6, G.2 and G.3). Additionally, its mid-time always agrees with the predicted $T_c$ reported in Sect. 5 for HD 88986 b, whatever the chosen RV model is (see Table 2). Furthermore, the feature exhibits fair similarity and

robustness across various tested approaches for correcting the instrumental variation in the TESS data. The consistent detection of this feature through various methodologies, combined with its robust nature, supports the fact that this feature is unlikely to be attributed to instrumental effects. However, the precise shape of the transit is not clearly discernible in the SAP data (as seen in the first panel of Fig. G.2). This lack of clarity is expected due to a momentum dump occurring at the time of the transit, in particular for such a shallow transit. Similar to other instances of momentum dumps observed in the SAP data, this event introduces instrumental variations. Additionally, the SAP data is not fully corrected for the scattered light that might affect the exact shape of the potential single transit. Given the robustness of the feature, as well as the fact that the observed $T_c$ of the potential single transit agrees, within uncertainties, with the $T_c$ obtained from all the RV-only models (see Sect. 5), it is likely that this feature corresponds to a single transit attributed to HD 88986 b.

We chose PDC-SAP data for our final analysis as all the newly extracted light curves are fairly consistent with the PDC-SAP light curve. We detrended the light curve using GP with an approximate Matern kernel introduced in Foreman-Mackey et al. (2017). The reason for this choice is that there is no evidence of existing quasi-periodic oscillations in the TESS light curves. Before applying this method to more accurately measure the planet's radii, we masked out the in-transit and immediately surrounding data points. The final detrended light curve is shown in Fig. 6 bottom panel, which we use for the rest of this work. We note that no additional signal was found by performing a Transit Least Squares (Hippke & Heller 2019) algorithm on this data.

Additionally, we searched for any potential light contamination caused by neighboring stars on the light curve of HD 88986. Within the aperture set by PDC-SAP, there is only one neighboring star. Because this star is one of the ~1 million new Gaia DR3 sources, its light contamination is not corrected by SPOC. This star has a magnitude of $G = 12.3$ ($\Delta G_{mag} = 6$ compared to HD 88986) and is located 1.4" west of HD 88986. There are no values for RP magnitude and renormalized unit weight error (RUWE), which is a measurement of the goodness of the star's astrometric solution. The poor behavior of this star could be due to blending with HD 88986. Because the Gaia RP bandpass is comparable to the TESS bandpass, one can assess the level of contamination using the Gaia RP fluxes of these stars. In our case, however, the lack of RP flux data for the neighbor star prevents us from estimating the contamination effect.

## 7. Joint analysis of SOPHIE+ and TESS sector 21

We performed a joint modeling of TESS photometric data of sector 21 and SOPHIE+ RVs using juliet. juliet, in addition to the packages mentioned above for RV modeling (see Sect. 5.2), employs batman (Kreidberg 2015) for transit fitting. Following the RV-only study in Sect. 5, we used the planet model with a simultaneous QP-GP model for modeling RVs in our joint modeling. Here, we tested both eccentricity free and zero models, as combining RV and transit data provides more constraints to fitting the orbital parameters including eccentricity and the argument of periastron. We employed the same priors for RV-related parameters, as detailed in Sect. 5, with the exception of $T_c$. For $T_c$, we set a Gaussian prior centered at 58 891.6 with a standard deviation of 5 d. To parameterize the limb darkening coefficient for TESS photometry, we applied a linear law through a parameter of $q$. This choice is motivated by the limited number of informative in-transit data points when modeling

<footnote n="5">⁵ https://docs.lightkurve.org</footnote>
vectors along with the CBVs to remove any flux trends. The final light curve is presented in the middle panel of Fig. 6, and the potential single transit is clearly seen within the dataset.

In a quest to further explore possible sources of instrumental noise that could impact the detection or shape of the potential, shallow single transit, we explored alternative methodologies following Rapetti et al. (in preparation). We used the adaptation of the technique Pixel Level Decorrelation (hereafter PLD; Deming et al. 2015; Luger et al. 2016, 2018) implemented in the PLDCorrector class of the community Python package Lightkurve<fnref n="5" />. This method employs (i) a spline polynomial fit to describe stellar variability, (ii) PCA eigenmodes to model the background light, and (iii) the PLD technique to account for pointing and mechanical effects. As an additional approach, we employed a version of PDC as adapted in the CBVCorrector class of Lightkurve, utilizing the CBV technique that the PDC method of the SPOC pipeline employs (hereafter we refer to this corrector as CBV).

Before applying the PLD, we added the background flux and errors estimated by the TESS SPOC pipeline back onto the SAP light curve. Flux level, fraction, and crowding adjustments are applied to the corrected light curves. To automatically optimize the selection of parameter values for the correctors, we evaluate the resulting light curve using the Savitsky-Golay Combined Differential Photometric Precision (sgCDPP) proxy algorithm (Gilliland et al. 2011; Van Cleve et al. 2016) implemented in Lightkurve, for durations of 30, 60, 120, 160, and 200 min (see the legend of Fig. G.2). For a grid of corrector parameter values (for further details on the parameters and the grid, see Rapetti et al., in prep.), we calculated the harmonic mean (HM) of these quantities and selected the corrected light curve that minimizes the HM. In addition, the final sgCDPP metrics can be compared to those obtained for the SPOC PDC-SAP corrected light curve (see Fig. G.2). For this comparison, we also calculated the over-fitting metric implemented in Lightkurve (see Fig. G.2) to measure the broad-band power spectrum via a Lomb-Scargle periodogram and assess the level of introduced noise in the corrected light curves.

In addition to the methods detailed above for extracting data from the 2-min TPFs, we extracted the full frame image (FFI) light curve from the TESS image using a strategy similar to Vanderburg et al. (2016). In particular, we created 20 different photometric apertures, 10 circular apertures, and 10 shaped like the TESS point spread function at the star location on the detector. We then calculated light curves from each aperture and corrected for systematics by performing a linear least-squares fit modeling on the light curve with time series the mean, standard deviation, kurtosis, and skew of the spacecraft quaternion measurements within each exposure (e.g., Vanderburg et al. 2019), the SPOC PDC CBVs, the background flux time series, and a basis spline to model slow variability. We performed the least-squares linear fit iteratively, removing outliers until convergence. Once the light curves were corrected for systematics, we corrected for dilution from other nearby stars and identified the one with the best photometric precision, which we used for our FFI analysis. The final resulting light curve is plotted in Fig. G.3.

The TESS light curves of sector 21, produced using different approaches, all find a feature where the potential single, shallow transit was identified (see Figs. 6, G.2 and G.3). Additionally, its mid-time always agrees with the predicted $T_c$ reported in Sect. 5 for HD 88986 b, whatever the chosen RV model is (see Table 2). Furthermore, the feature exhibits fair similarity and

robustness across various tested approaches for correcting the instrumental variation in the TESS data. The consistent detection of this feature through various methodologies, combined with its robust nature, supports the fact that this feature is unlikely to be attributed to instrumental effects. However, the precise shape of the transit is not clearly discernible in the SAP data (as seen in the first panel of Fig. G.2). This lack of clarity is expected due to a momentum dump occurring at the time of the transit, in particular for such a shallow transit. Similar to other instances of momentum dumps observed in the SAP data, this event introduces instrumental variations. Additionally, the SAP data is not fully corrected for the scattered light that might affect the exact shape of the potential single transit. Given the robustness of the feature, as well as the fact that the observed $T_c$ of the potential single transit agrees, within uncertainties, with the $T_c$ obtained from all the RV-only models (see Sect. 5), it is likely that this feature corresponds to a single transit attributed to HD 88986 b.

We chose PDC-SAP data for our final analysis as all the newly extracted light curves are fairly consistent with the PDC-SAP light curve. We detrended the light curve using GP with an approximate Matern kernel introduced in Foreman-Mackey et al. (2017). The reason for this choice is that there is no evidence of existing quasi-periodic oscillations in the TESS light curves. Before applying this method to more accurately measure the planet's radii, we masked out the in-transit and immediately surrounding data points. The final detrended light curve is shown in Fig. 6 bottom panel, which we use for the rest of this work. We note that no additional signal was found by performing a Transit Least Squares (Hippke & Heller 2019) algorithm on this data.

Additionally, we searched for any potential light contamination caused by neighboring stars on the light curve of HD 88986. Within the aperture set by PDC-SAP, there is only one neighboring star. Because this star is one of the ~1 million new Gaia DR3 sources, its light contamination is not corrected by SPOC. This star has a magnitude of $G = 12.3$ ($\Delta G_{mag} = 6$ compared to HD 88986) and is located 1.4" west of HD 88986. There are no values for RP magnitude and renormalized unit weight error (RUWE), which is a measurement of the goodness of the star's astrometric solution. The poor behavior of this star could be due to blending with HD 88986. Because the Gaia RP bandpass is comparable to the TESS bandpass, one can assess the level of contamination using the Gaia RP fluxes of these stars. In our case, however, the lack of RP flux data for the neighbor star prevents us from estimating the contamination effect.

## 7. Joint analysis of SOPHIE+ and TESS sector 21

We performed a joint modeling of TESS photometric data of sector 21 and SOPHIE+ RVs using juliet. juliet, in addition to the packages mentioned above for RV modeling (see Sect. 5.2), employs batman (Kreidberg 2015) for transit fitting. Following the RV-only study in Sect. 5, we used the planet model with a simultaneous QP-GP model for modeling RVs in our joint modeling. Here, we tested both eccentricity free and zero models, as combining RV and transit data provides more constraints to fitting the orbital parameters including eccentricity and the argument of periastron. We employed the same priors for RV-related parameters, as detailed in Sect. 5, with the exception of $T_c$. For $T_c$, we set a Gaussian prior centered at 58 891.6 with a standard deviation of 5 d. To parameterize the limb darkening coefficient for TESS photometry, we applied a linear law through a parameter of $q$. This choice is motivated by the limited number of informative in-transit data points when modeling

<footnote n="5">⁵ https://docs.lightkurve.org</footnote>
<table><caption>**Table 3.** Median values and 68% confidence interval of parameters for HD 88986 b based on the joint analysis of the photometric and RV data by juliet (see Sect. 7 and Fig. 8).</caption><thead><tr><th>Parameter (unit)</th><th>Posterior HD 88986 b</th></tr></thead><tbody><tr><td>Stellar parameters:</td><td></td></tr><tr><td>ρ<sub>*</sub> (kg m<sup>-3</sup>)</td><td>472<sup>+38</sup><sub>-35</sub></td></tr><tr><td>Planet parameters</td><td></td></tr><tr><td>P<sup>*</sup> (d)</td><td>146.05<sup>+0.43</sup><sub>-0.40</sub></td></tr><tr><td>T<sub>c</sub> (BJD-2 400 000)</td><td>58 891.690 ± 0.003</td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>1.85 ± 0.34</td></tr><tr><td>e</td><td>0.24 ± 0.05</td></tr><tr><td>ω (°)</td><td>306<sup>+10</sup><sub>-11</sub></td></tr><tr><td>R<sub>p</sub>/R<sub>*</sub></td><td>0.0148 ± 0.0004</td></tr><tr><td>b</td><td>0.21<sup>+0.17</sup><sub>-0.14</sub></td></tr><tr><td>a/R<sub>*</sub></td><td>81.1 ± 2.1</td></tr><tr><td>TESS instrumental parameters:</td><td></td></tr><tr><td>M<sub>TESS</sub> (ppm)</td><td>-0.0000001 ± 0.0000014</td></tr><tr><td>σ<sub>ω,TESS</sub> (ppm)</td><td>143.4 ± 1.3</td></tr><tr><td>q<sup>1</sup></td><td>0.27<sup>+0.17</sup><sub>-0.14</sub></td></tr><tr><td>SOPHIE+ instrumental parameters:</td><td></td></tr><tr><td>σ<sub>SOPHIE+</sub> (m s<sup>-1</sup>)</td><td>2.09<sup>+0.16</sup><sub>-0.17</sub></td></tr><tr><td>μ<sub>SOPHIE+</sub> (m s<sup>-1</sup>)</td><td>29 090.6 ± 0.4</td></tr><tr><td>Drift on SOPHIE+:</td><td></td></tr><tr><td>A (m s<sup>-1</sup>)</td><td>0.0165 ± 0.0003</td></tr><tr><td>Q (m s<sup>-1</sup>)</td><td>0.0000039 ± 0.0000002</td></tr><tr><td>QP-GP on SOPHIE+:</td><td></td></tr><tr><td>B<sub>GP</sub> (m s<sup>-1</sup>)</td><td>4.0<sup>+1.1</sup><sub>-0.8</sub></td></tr><tr><td>C<sub>GP</sub> (m s<sup>-1</sup>)</td><td>10<sup>-20</sup> (fixed)</td></tr><tr><td>P<sub>rot</sub> (d)</td><td>30.0 ± 2.0</td></tr><tr><td>L<sub>GP</sub> (d)</td><td>29<sup>+37</sup><sub>-18</sub></td></tr><tr><td>Derived planet parameters:</td><td></td></tr><tr><td>a (au)</td><td>0.58 ± 0.04</td></tr><tr><td>i (°)</td><td>89.9 ± 0.1</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>2.49 ± 0.18</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>17.2<sup>+4.0</sup><sub>-3.8</sub></td></tr><tr><td>ρ<sub>p</sub> (g cm<sup>-3</sup>)</td><td>6.1<sup>+3.3</sup><sub>-2.3</sub></td></tr><tr><td>T<sub>eq</sub> (K)</td><td>460 ± 8</td></tr></tbody></table>

**Notes.** (*) See the refined period in Table 4.

single transit events (Sandford et al. 2019). We note that applying a quadratic limb-darkening law had no effect on our final results. Additionally, the results of the spectral analysis provided in Sect. 3 were used to set a Gaussian prior to the star density ρ*. Finally, we added a jitter term of σ to TESS photometric data. Table H.1 shows all of the employed priors and the description of the parameters.

The joint modeling of TESS sector 21 and SOPHIE+ RVs with eccentricity-free and circular orbits yields consistent results. However, the eccentricity-free model exhibits strong statistical favorability (Δln Z > 9). Therefore, we present the resulting parameters of this model in Table 3 and depict its best-fit on combined photometry and RVs in Fig. 7. A corner plot of all parameters is included in Appendix J.1. We note that the

dilution factor was not taken into account in our analysis. This is due to the fact that when a wide uniform prior range of 0–1 is used, the dilution value tends to converge toward the lower edge of the prior 0. This outcome is unacceptable because there is only one faint nearby star in the TESS aperture, and it is approximately 6 mag fainter than the primary star in the G band (refer to Eq. (6) of Espinoza et al. 2019 for more details). Such behavior can be expected in transits with low S/Ns (Espinoza et al. 2019). Additionally, due to the unavailability of Gaia RP magnitude for the neighbor star, we were unable to estimate and constrain the dilution factor accurately. However, considering that the neighbor star is faint, we expect any contamination effect from it to be negligible.

Based on our final parameters, the transiting planet HD 88986 b is a sub-Neptune with a period of 146.05<sup>+0.43</sup><sub>-0.40</sub> d. It exhibits an eccentricity of 0.23 ± 0.06, alongside a radius of 2.49 ± 0.18 R<sub>⊕</sub> and a mass of 17.2<sup>+4.0</sup><sub>-3.8</sub> M<sub>⊕</sub>, corresponding to a high mean density of 6.1<sup>+3.3</sup><sub>-2.3</sub> g cm<sup>-3</sup>. Additionally, the planet has an equilibrium temperature of 460 ± 8 K (T<sub>eq</sub> = T<sub>*</sub>√(R<sub>*</sub>/2a), Méndez & Rivera-Valentín 2017), making it a relatively cool planet.

## 8. Search for a second transit with additional photometric data

The analysis presented above in Sect. 7 predicts another transit of HD 88986 b should have occurred in 2022 February. We used CHEOPS and TESS sector 48 to attempt the detection of the second transit.

### 8.1. CHEOPS photometry

The CHEOPS spacecraft is a 30 cm ESA space telescope (Benz et al. 2021) that conducts ultra-high-precision photometry to characterize planets (Bonfanti et al. 2021; Delrez et al. 2021; Lacedelli et al. 2022) and their atmospheres (Lendl et al. 2020; Hooton et al. 2022), but it has also been used to aid in the discovery of new planets (Leleu et al. 2021; Osborn et al. 2022; Serrano et al. 2022; Wilson et al. 2022).

As derived in Sect. 7, the period of HD 88986 b is 146.05<sup>+0.43</sup><sub>-0.40</sub> d. We predicted that the next transiting event would fall within the region of the TESS sector 48 which includes the gap between orbits. Therefore, to search for a second transit event, we obtained one visit of CHEOPS observation (PI: N. Heidari) spanning 167.4 h between 2022-02-08 and 2022-02-15 with an exposure time of 3.4 s. This allowed us to cover the transit period's uncertainty from joint analysis of RVs and sector 21 by ~2σ and the TESS gap to be covered.

The data were processed with the latest version of the CHEOPS Data Reduction Pipeline (DRP v13; Hoyer et al. 2020) that conducts frame calibration, instrumental and environmental correction, and aperture photometry using predefined radii (R = 22.5", 25.0", and 30.0") as well as a noise-optimized radius. The DRP-produced flux contamination was subtracted from the light curves. We retrieved the data and corresponding instrumental basis vectors and assessed the quality using the PYCHEOPS Python package (Maxted et al. 2022) and found that the DEFAULT aperture minimized the root mean square (RMS) noise. Therefore, we used this aperture for our further analysis. These are plotted in the upper panel of Fig. 8, on which the expected time of the transit, derived from the SOPHIE+ RVs and TESS sector 21 joint analysis (Sect. 7), is indicated in blue.
<table><caption>**Table 3.** Median values and 68% confidence interval of parameters for HD 88986 b based on the joint analysis of the photometric and RV data by juliet (see Sect. 7 and Fig. 8).</caption><thead><tr><th>Parameter (unit)</th><th>Posterior HD 88986 b</th></tr></thead><tbody><tr><td colspan="2"><strong>Stellar parameters:</strong></td></tr><tr><td>ρ<sub>*</sub> (kg m<sup>-3</sup>)</td><td>472<sup>+38</sup><sub>-35</sub></td></tr><tr><td colspan="2"><strong>Planet parameters</strong></td></tr><tr><td>P<sup>*</sup> (d)</td><td>146.05<sup>+0.43</sup><sub>-0.40</sub></td></tr><tr><td>T<sub>c</sub> (BJD-2 400 000)</td><td>58 891.690 ± 0.003</td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>1.85 ± 0.34</td></tr><tr><td>e</td><td>0.24 ± 0.05</td></tr><tr><td>ω (°)</td><td>306<sup>+10</sup><sub>-11</sub></td></tr><tr><td>R<sub>p</sub>/R<sub>*</sub></td><td>0.0148 ± 0.0004</td></tr><tr><td>b</td><td>0.21<sup>+0.17</sup><sub>-0.14</sub></td></tr><tr><td>a/R<sub>*</sub></td><td>81.1 ± 2.1</td></tr><tr><td colspan="2"><strong>TESS instrumental parameters:</strong></td></tr><tr><td>M<sub>TESS</sub> (ppm)</td><td>-0.0000001 ± 0.0000014</td></tr><tr><td>σ<sub>ω,TESS</sub> (ppm)</td><td>143.4 ± 1.3</td></tr><tr><td>q<sup>1</sup></td><td>0.27<sup>+0.17</sup><sub>-0.14</sub></td></tr><tr><td colspan="2"><strong>SOPHIE+ instrumental parameters:</strong></td></tr><tr><td>σ<sub>SOPHIE+</sub> (m s<sup>-1</sup>)</td><td>2.09<sup>+0.16</sup><sub>-0.17</sub></td></tr><tr><td>μ<sub>SOPHIE+</sub> (m s<sup>-1</sup>)</td><td>29 090.6 ± 0.4</td></tr><tr><td colspan="2"><strong>Drift on SOPHIE+:</strong></td></tr><tr><td>A (m s<sup>-1</sup>)</td><td>0.0165 ± 0.0003</td></tr><tr><td>Q (m s<sup>-1</sup>)</td><td>0.0000039 ± 0.0000002</td></tr><tr><td colspan="2"><strong>QP-GP on SOPHIE+:</strong></td></tr><tr><td>B<sub>GP</sub> (m s<sup>-1</sup>)</td><td>4.0<sup>+1.1</sup><sub>-0.8</sub></td></tr><tr><td>C<sub>GP</sub> (m s<sup>-1</sup>)</td><td>10<sup>-20</sup> (fixed)</td></tr><tr><td>P<sub>rot</sub> (d)</td><td>30.0 ± 2.0</td></tr><tr><td>L<sub>GP</sub> (d)</td><td>29<sup>+37</sup><sub>-18</sub></td></tr><tr><td colspan="2"><strong>Derived planet parameters:</strong></td></tr><tr><td>a (au)</td><td>0.58 ± 0.04</td></tr><tr><td>i (°)</td><td>89.9 ± 0.1</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>2.49 ± 0.18</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>17.2<sup>+4.0</sup><sub>-3.8</sub></td></tr><tr><td>ρ<sub>p</sub> (g cm<sup>-3</sup>)</td><td>6.1<sup>+3.3</sup><sub>-2.3</sub></td></tr><tr><td>T<sub>eq</sub> (K)</td><td>460 ± 8</td></tr></tbody></table>

**Notes.** (*)See the refined period in Table 4.

single transit events (Sandford et al. 2019). We note that applying a quadratic limb-darkening law had no effect on our final results. Additionally, the results of the spectral analysis provided in Sect. 3 were used to set a Gaussian prior to the star density ρ*. Finally, we added a jitter term of σ to TESS photometric data. Table H.1 shows all of the employed priors and the description of the parameters.

The joint modeling of TESS sector 21 and SOPHIE+ RVs with eccentricity-free and circular orbits yields consistent results. However, the eccentricity-free model exhibits strong statistical favorability (Δln Z > 9). Therefore, we present the resulting parameters of this model in Table 3 and depict its best-fit on combined photometry and RVs in Fig. 7. A corner plot of all parameters is included in Appendix J.1. We note that the

dilution factor was not taken into account in our analysis. This is due to the fact that when a wide uniform prior range of 0–1 is used, the dilution value tends to converge toward the lower edge of the prior 0. This outcome is unacceptable because there is only one faint nearby star in the TESS aperture, and it is approximately 6 mag fainter than the primary star in the G band (refer to Eq. (6) of Espinoza et al. 2019 for more details). Such behavior can be expected in transits with low S/Ns (Espinoza et al. 2019). Additionally, due to the unavailability of Gaia RP magnitude for the neighbor star, we were unable to estimate and constrain the dilution factor accurately. However, considering that the neighbor star is faint, we expect any contamination effect from it to be negligible.

Based on our final parameters, the transiting planet HD 88986 b is a sub-Neptune with a period of 146.05<sup>+0.43</sup><sub>-0.40</sub> d. It exhibits an eccentricity of 0.23 ± 0.06, alongside a radius of 2.49 ± 0.18 R<sub>⊕</sub> and a mass of 17.2<sup>+4.0</sup><sub>-3.8</sub> M<sub>⊕</sub>, corresponding to a high mean density of 6.1<sup>+3.3</sup><sub>-2.3</sub> g cm<sup>-3</sup>. Additionally, the planet has an equilibrium temperature of 460 ± 8 K (T<sub>eq</sub> = T<sub>*</sub>√(R<sub>*</sub>/2a), Méndez & Rivera-Valentín 2017), making it a relatively cool planet.

## 8. Search for a second transit with additional photometric data

The analysis presented above in Sect. 7 predicts another transit of HD 88986 b should have occurred in 2022 February. We used CHEOPS and TESS sector 48 to attempt the detection of the second transit.

### 8.1. CHEOPS photometry

The CHEOPS spacecraft is a 30 cm ESA space telescope (Benz et al. 2021) that conducts ultra-high-precision photometry to characterize planets (Bonfanti et al. 2021; Delrez et al. 2021; Lacedelli et al. 2022) and their atmospheres (Lendl et al. 2020; Hooton et al. 2022), but it has also been used to aid in the discovery of new planets (Leleu et al. 2021; Osborn et al. 2022; Serrano et al. 2022; Wilson et al. 2022).

As derived in Sect. 7, the period of HD 88986 b is 146.05<sup>+0.43</sup><sub>-0.40</sub> d. We predicted that the next transiting event would fall within the region of the TESS sector 48 which includes the gap between orbits. Therefore, to search for a second transit event, we obtained one visit of CHEOPS observation (PI: N. Heidari) spanning 167.4 h between 2022-02-08 and 2022-02-15 with an exposure time of 3.4 s. This allowed us to cover the transit period's uncertainty from joint analysis of RVs and sector 21 by ~2σ and the TESS gap to be covered.

The data were processed with the latest version of the CHEOPS Data Reduction Pipeline (DRP v13; Hoyer et al. 2020) that conducts frame calibration, instrumental and environmental correction, and aperture photometry using predefined radii (R = 22.5", 25.0", and 30.0") as well as a noise-optimized radius. The DRP-produced flux contamination was subtracted from the light curves. We retrieved the data and corresponding instrumental basis vectors and assessed the quality using the PYCHEOPS Python package (Maxted et al. 2022) and found that the DEFAULT aperture minimized the root mean square (RMS) noise. Therefore, we used this aperture for our further analysis. These are plotted in the upper panel of Fig. 8, on which the expected time of the transit, derived from the SOPHIE+ RVs and TESS sector 21 joint analysis (Sect. 7), is indicated in blue.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Joint analysis of SOPHIE+ and TESS sector 21 observations of HD 88986 b. Top: SOPHIE+ data overplotted by the best-fit orbit model. Bottom-left: phase-folded TESS PDC-SAP photometric data of sector 21. The data are binned (red points) in 1 hour. The black line shows the best-fit transit model. Bottom-right: phase-folded SOPHIE+ RVs of HD 88986 b at the period of 146.05 d. The red points depict the binned data, utilizing a bin size of 0.05 in orbital phase units. The black line represents the best-fit orbit model.</figcaption></figure>

In previous studies, it has been noted that environmental effects (i.e., spacecraft temperature and illumination) and the presence of nearby contaminants can induce flux modulation in light curves (Morris et al. 2021; Maxted et al. 2022; Wilson et al. 2022). In order to correct for these effects and search for the smallest transit signals in our transit search analysis, we conduct a principal component analysis on the auto-correlation function of the CHEOPS frames using the methodology detailed in Wilson et al. (2022). The process has been shown to monitor PSF shape changes, and so any effects that alter the CHEOPS PSF, such as environmental and contamination effects, are measured by this tool and can be removed by using the produced principal components as the basis vectors in a linear model detrending. Further examples of applications of this tool can be seen in Hoyer et al. (2022); Ehrenreich et al. (2023); Hawthorn et al. (2023).

To assess the existence of a transiting body with the CHEOPS observations, we conduct a statistical analysis using a newly developed tool Wilson et al. (2023). In brief, we use the PSF-based PCA components produced following (Wilson et al. 2022) above in combination with the instrumental basis vectors to construct a linear noise model that is fit simultaneously with either a 0 or 1 planet transit model that allows us to compute the True and False Inclusion Probabilities (TIP and FIP; Hara et al. 2022b) for the presence of transit in the data. These are calculated using the Bayes Evidences and posterior distributions of the 0 and 1 planet fits. For this study, we conduct this analysis twice: one with a period prior constrained by the transit model from the RV data and the second time with no period prior. For both cases and for all transit $T_0$ values within the CHEOPS dataset, we find FIP ~ 1, which statistically means that there is no transit in the lightcurve.

This non-detection by CHEOPS leads to two potential conclusions: (1) the rejection of the presumed transit event in TESS sector 21 as a spurious feature, or (2) the transit might be occurring at a time that deviates more than ~2σ away from our predicted ephemeris. We note that although the CHEOPS data did not reveal any transit features, this observation was valuable in covering the gap in TESS sector 48 data (see Sect. 8.2 and Fig. 8), substantially contributing to the refinement of the planet's period.

## 8.2. TESS photometry sector 48

TESS conducted a second observation of this star with a cadence of 2 minutes, spanning from 2022 January 28 to February 26. The PDC-SAP data provided by SPOC from this observation is depicted in Fig. 8 (second panel). Remarkably proximate to the 3σ anticipated transit region, as determined through the joint modeling of RVs and TESS sector 21 data, a second potential transit-like feature with a $T_c$ of about ~59 628.8 and a period of 147.4 d from the first transit is observed in the figure. We note that it was particularly lucky to cover two transits of that long-period planet with TESS, as only two sectors of TESS covered that star.

Similar to the first potential transit in sector 21, we tested the mean in-transit and mean out-of-transit flux, along with the difference between them (see Fig. G.1). This test confirmed that the feature in sector 48 is likely related to the host star. However, in contrast to sector 21 with a standard deviation of $\sigma$ = 204.4 ppm, TESS sector 48 displays a substantial dispersion with $\sigma$ = 635.8 ppm. This noticeable difference might be linked to the
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Joint analysis of SOPHIE+ and TESS sector 21 observations of HD 88986 b. Top: SOPHIE+ data overplotted by the best-fit orbit model. Bottom-left: phase-folded TESS PDC-SAP photometric data of sector 21. The data are binned (red points) in 1 hour. The black line shows the best-fit transit model. Bottom-right: phase-folded SOPHIE+ RVs of HD 88986 b at the period of 146.05 d. The red points depict the binned data, utilizing a bin size of 0.05 in orbital phase units. The black line represents the best-fit orbit model.</figcaption></figure>

In previous studies, it has been noted that environmental effects (i.e., spacecraft temperature and illumination) and the presence of nearby contaminants can induce flux modulation in light curves (Morris et al. 2021; Maxted et al. 2022; Wilson et al. 2022). In order to correct for these effects and search for the smallest transit signals in our transit search analysis, we conduct a principal component analysis on the auto-correlation function of the CHEOPS frames using the methodology detailed in Wilson et al. (2022). The process has been shown to monitor PSF shape changes, and so any effects that alter the CHEOPS PSF, such as environmental and contamination effects, are measured by this tool and can be removed by using the produced principal components as the basis vectors in a linear model detrending. Further examples of applications of this tool can be seen in Hoyer et al. (2022); Ehrenreich et al. (2023); Hawthorn et al. (2023).

To assess the existence of a transiting body with the CHEOPS observations, we conduct a statistical analysis using a newly developed tool Wilson et al. (2023). In brief, we use the PSF-based PCA components produced following (Wilson et al. 2022) above in combination with the instrumental basis vectors to construct a linear noise model that is fit simultaneously with either a 0 or 1 planet transit model that allows us to compute the True and False Inclusion Probabilities (TIP and FIP; Hara et al. 2022b) for the presence of transit in the data. These are calculated using the Bayes Evidences and posterior distributions of the 0 and 1 planet fits. For this study, we conduct this analysis twice: one with a period prior constrained by the transit model from the RV data and the second time with no period prior. For both cases and for all transit $T_0$ values within the CHEOPS dataset, we find FIP ~ 1, which statistically means that there is no transit in the lightcurve.

This non-detection by CHEOPS leads to two potential conclusions: (1) the rejection of the presumed transit event in TESS sector 21 as a spurious feature, or (2) the transit might be occurring at a time that deviates more than ~2σ away from our predicted ephemeris. We note that although the CHEOPS data did not reveal any transit features, this observation was valuable in covering the gap in TESS sector 48 data (see Sect. 8.2 and Fig. 8), substantially contributing to the refinement of the planet's period.

## 8.2. TESS photometry sector 48

TESS conducted a second observation of this star with a cadence of 2 minutes, spanning from 2022 January 28 to February 26. The PDC-SAP data provided by SPOC from this observation is depicted in Fig. 8 (second panel). Remarkably proximate to the 3σ anticipated transit region, as determined through the joint modeling of RVs and TESS sector 21 data, a second potential transit-like feature with a $T_c$ of about ~59 628.8 and a period of 147.4 d from the first transit is observed in the figure. We note that it was particularly lucky to cover two transits of that long-period planet with TESS, as only two sectors of TESS covered that star.

Similar to the first potential transit in sector 21, we tested the mean in-transit and mean out-of-transit flux, along with the difference between them (see Fig. G.1). This test confirmed that the feature in sector 48 is likely related to the host star. However, in contrast to sector 21 with a standard deviation of $\sigma$ = 204.4 ppm, TESS sector 48 displays a substantial dispersion with $\sigma$ = 635.8 ppm. This noticeable difference might be linked to the
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** CHEOPS and TESS observations of HD 88986 in 2022 February. The predicted time of HD 88986 b's second transit event, based on the best-fit model of combining RVs and the photometric light curve of sector 21 (see Sect. 7), is highlighted in blue. Top: CHEOPS photometric data. Second: TESS (black dots) PDC-SAP light curve of sector 48. The TESS data are binned (red points) in 1-hour increments. The potential transit-like feature is marked by a green triangle. The vertical dashed, green lines are the spacecraft momentum dumps. Third: re-extracted TESS light curve (see the text for more explanation). Bottom: zoomed on the potential transit event of HD 88986 b on PDC-SAP (left) and re-extracted (right) data. As the plot indicates, the two light curves are noticeably different. As a result, we did not include those data in our joint analysis presented in Sect. 7.</figcaption></figure>

presence of residual systematics that may have persisted even after the SPOC correction.

Therefore, following the methodology outlined in Sect. 6 using the TPFED tool, we conducted a customized extraction of TESS sector 48. The photometric results obtained from this approach are presented in Fig. 8 (third panel), with the transit-like feature zoomed in for better visualization (fourth panel). While the resulting custom extraction of the TESS light curve for sector 21 was fairly similar to the SPOC light curve, they are noticeably different for sector 48 (fourth panels).

To further explore the potential source of instrumental noises, similar to sector 21, we also extracted the light curve using PLD and CBV approaches from 2-min TPFs cadence, along with the FFI data (see Sect. 6 for more information about the methods). The resulting extracted data (see Figs. G.2 and G.3) following these methods also confirmed the noticeable difference between PDC-SAP data and the independently extracted light curves. Moreover, the variations within the transit further complicate our understanding of the TESS photometric data in sector 48, leading to its exclusion from our analysis. Additionally, we note that by the inclusion of PDC-SAP data of sector 48 in our joint modeling, the final results remain consistent compared to our model with only RVs and sector 21 in Sect. 7.

Figure G.4 shows the phase-folded TESS PDC-SAP and FFI data for sectors 21 and 48 corresponding to the 147.4 d period. We note that the FFI light curve is detrended by the

<table><caption>**Table 4.** Possible solution for HD 88986 b's period.</caption><thead><tr><td>Models</td><td>HD 88986 b's period</td></tr></thead><tbody><tr><td>RVs-only</td><td>146.3 ± 0.6</td></tr><tr><td>RVs+ sector 21</td><td>146.1 ± 0.4</td></tr><tr><td>RVs+ sector 21+ sector 48</td><td>147.4 ± 0.1</td></tr></tbody></table>

spline approach using the Wotan package (Hippke et al. 2019). While the consistency between the two potential transits in the PDC-SAP data remains uncertain, the two transits exhibit good consistency within the FFI data, particularly concerning transit depth and duration.

In Table 4, we summarize all the possible period solutions for HD 88986 b, including the periods obtained through the final choice of the RV-only model (see Sect. 5.2), combined RVs with TESS sector 21 (see Sect. 7), and combined RVs with both TESS sector 21 and 48. The orbital period derived from RV data combined with two potential transits agrees (at $3\sigma$) with the period calculated using RVs and single transit in sector 21, and also agrees (at $2\sigma$) with the RV-only period. One could expect an even better agreement; possible persistent instrumental effects not perfectly taken into account in our models might be the cause. This agreement between all period solutions, arguing here for an actual detection of a transit of HD 88986 b. Still, the noise
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** CHEOPS and TESS observations of HD 88986 in 2022 February. The predicted time of HD 88986 b's second transit event, based on the best-fit model of combining RVs and the photometric light curve of sector 21 (see Sect. 7), is highlighted in blue. Top: CHEOPS photometric data. Second: TESS (black dots) PDC-SAP light curve of sector 48. The TESS data are binned (red points) in 1-hour increments. The potential transit-like feature is marked by a green triangle. The vertical dashed, green lines are the spacecraft momentum dumps. Third: re-extracted TESS light curve (see the text for more explanation). Bottom: zoomed on the potential transit event of HD 88986 b on PDC-SAP (left) and re-extracted (right) data. As the plot indicates, the two light curves are noticeably different. As a result, we did not include those data in our joint analysis presented in Sect. 7.</figcaption></figure>

presence of residual systematics that may have persisted even after the SPOC correction.

Therefore, following the methodology outlined in Sect. 6 using the TPFED tool, we conducted a customized extraction of TESS sector 48. The photometric results obtained from this approach are presented in Fig. 8 (third panel), with the transit-like feature zoomed in for better visualization (fourth panel). While the resulting custom extraction of the TESS light curve for sector 21 was fairly similar to the SPOC light curve, they are noticeably different for sector 48 (fourth panels).

To further explore the potential source of instrumental noises, similar to sector 21, we also extracted the light curve using PLD and CBV approaches from 2-min TPFs cadence, along with the FFI data (see Sect. 6 for more information about the methods). The resulting extracted data (see Figs. G.2 and G.3) following these methods also confirmed the noticeable difference between PDC-SAP data and the independently extracted light curves. Moreover, the variations within the transit further complicate our understanding of the TESS photometric data in sector 48, leading to its exclusion from our analysis. Additionally, we note that by the inclusion of PDC-SAP data of sector 48 in our joint modeling, the final results remain consistent compared to our model with only RVs and sector 21 in Sect. 7.

Figure G.4 shows the phase-folded TESS PDC-SAP and FFI data for sectors 21 and 48 corresponding to the 147.4 d period. We note that the FFI light curve is detrended by the

<table><caption>**Table 4.** Possible solution for HD 88986 b's period.</caption><thead><tr><td>Models</td><td>HD 88986 b's period</td></tr></thead><tbody><tr><td>RVs-only</td><td>146.3 ± 0.6</td></tr><tr><td>RVs+ sector 21</td><td>146.1 ± 0.4</td></tr><tr><td>RVs+ sector 21+ sector 48</td><td>147.4 ± 0.1</td></tr></tbody></table>

spline approach using the Wotan package (Hippke et al. 2019). While the consistency between the two potential transits in the PDC-SAP data remains uncertain, the two transits exhibit good consistency within the FFI data, particularly concerning transit depth and duration.

In Table 4, we summarize all the possible period solutions for HD 88986 b, including the periods obtained through the final choice of the RV-only model (see Sect. 5.2), combined RVs with TESS sector 21 (see Sect. 7), and combined RVs with both TESS sector 21 and 48. The orbital period derived from RV data combined with two potential transits agrees (at $3\sigma$) with the period calculated using RVs and single transit in sector 21, and also agrees (at $2\sigma$) with the RV-only period. One could expect an even better agreement; possible persistent instrumental effects not perfectly taken into account in our models might be the cause. This agreement between all period solutions, arguing here for an actual detection of a transit of HD 88986 b. Still, the noise
in sector 48 light curve and the differences between reduction methods keep us prudent about the transit detection in sector 48, which we chose not to include in our final fit (Sect. 7). Conducting follow-up photometric observations for this system, with the goal of identifying HD 88986 b's second transit event, would strongly confirm that the planet is transiting while also providing a much better constraint on the planet's period.

# 9. Constraining a long-term companion

In this section, we examined different scenarios to determine the origin of long-term curvature seen in the RVs in addition to HD 88986 b. Stellar activity or a wide-orbit companion are two possibilities. To explore them, we used long-term photometric data obtained with the T8 APT, the combined RVs from ELODIE, HIRES, and SOPHIE instruments, as well as astrometric data from Gaia and HIPPARCOS.

## 9.1. APT photometric observations

To characterize the origin of long-term curvature observed in RVs (see Fig. 1), we used 1335 photometric observations of HD 88986 covering 21 observing seasons from 1995–1996 to 2019–2020, except the four observing seasons 2015–16 through 2018–19, during which the star was not observed. The observations were acquired with the T8 0.80 m APT at Fairborn Observatory in southern Arizona. The T8 APT is equipped with a two-channel photometer that uses two EMI 9124QB bi-alkali photomultiplier tubes to measure the stellar brightness simultaneously in the Strömgren *b* and *y* passbands.

The observations are made differentially with respect to three nearby comparison stars. We measured the difference in brightness between our program star HD 88986 (star d) and the comparison stars (stars a: HD 89557 (G = 7.3 mag, G8 III), b: HD 87667 (G = 7.3, F5), and c: HD 88476 (G = 6.6, G8 III)) and created differential magnitudes in the following six combinations: d-a, d-b, d-c, c-a, c-b, and b-a. Intercomparison of these six light curves shows that the comparison star a (HD 89557) is the only one that appears to be constant to the limit of our precision, so we present our results as differential magnitudes in the sense star d minus star a, which we designate as d-a.

To improve the photometric precision of the individual nightly observations, we combined the differential *b* and *y* magnitudes into a single (*b* + *y*)/2 passband. The precision of a single differential observation with T8, as measured from pairs of constant comparison stars, typically ranges between 0.001 mag and 0.0015 mag on good nights. The T8 APT is described in Henry (1999), where further details of the telescope, precision photometer, and observing and data reduction procedures can be found.

Figure 9 plots the 1335 nightly observations in the (*b* + *y*)/2 passband photometry from d-a obtained across the 21 observing seasons as small filled circles. The mean of all the nightly observations, -1.16492 mag, is plotted as the dashed line in the figure. The standard deviation of the nightly observations from their mean is 0.00118 mag, consistent with the precision of the measurements. The 21 seasonal means of these data are plotted as large filled circles. The standard deviations of the individual seasonal means are roughly the size of the plot symbols. The standard deviation of the 21 seasonal means from the mean of the seasonal means is 0.00028 mag, indicating that there is no long-term variability in HD 88986 to the limit of our photometric precision.

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Nightly Strömgren (*b* + *y*)/2 band photometry of HD 88986 from 21 observing seasons from 1995–96 to 2020–21 (small filled circles) scatter about their mean (dashed line) with a standard deviation of 0.00118 mag. Seasonal means from the 21 seasons (large filled circles) scatter about their mean with a standard deviation of 0.00028 mag. No significant variations nor periodicities are detected.</figcaption></figure>

Table L.1 summarizes observations in the (*b* + *y*)/2 passband photometry from d-a. The standard deviations of the nightly observations for each observing season indicate little or no short-term variability within each observing season. Frequency analysis of each individual observing season using the method of Vaníček (1971) confirms the lack of any periodic variability. Henry et al. (2022) show extensive examples of this method of period analysis.

In a study conducted by Lovis et al. (2011), the activity cycles of 311 FGK stars were analyzed. It shows that RV semi-amplitudes can be induced up to approximately 25 m s⁻¹ by the stellar long-period activity. However, in the case of HD 88986, the long-period RV semi-amplitude is at least about 40 m. s⁻¹ (see the SOPHIE+ RVs in Fig. 7), and also as stated in this section, there is no evidence of long-term photometric variability up to the limits of our precision. This suggests that the origin of this long-term curvature is likely unrelated to stellar activity and instead points toward the presence of a third body in the system.

## 9.2. Combining RV and HIPPARCOS/Gaia astrometry data

In order to improve the characterization of the outer massive companion, we modeled simultaneously the available RV and absolute astrometry data. The use of an MCMC algorithm enables us to explore the different solutions for each orbital parameter and for the companion mass compatible with the data. This algorithm, introduced by Philipot et al. (2023), is used to fit Keplerian orbits based on the emcee 3.0 (Foreman-Mackey et al. 2013) and HTOF (Brandt et al. 2021a) packages. The likelihood computation is similar to that of the ORVARA code (Brandt et al. 2021b).

We considered the ELODIE, SOPHIE, SOPHIE+, HIRES, and HIRES+ RV data, presented previously, coupled with the proper motion and position values calculated in the HIPPARCOS-Gaia Catalog of Accelerations (Brandt 2021) from HIPPARCOS (Perryman et al. 1997; van Leeuwen 2007) and Gaia data release 3 (DR3; Gaia Collaboration 2016, 2021) measurements. For the fit, we considered a Gaussian prior for the stellar mass and parallax, based on the values published
in sector 48 light curve and the differences between reduction methods keep us prudent about the transit detection in sector 48, which we chose not to include in our final fit (Sect. 7). Conducting follow-up photometric observations for this system, with the goal of identifying HD 88986 b's second transit event, would strongly confirm that the planet is transiting while also providing a much better constraint on the planet's period.

## 9. Constraining a long-term companion

In this section, we examined different scenarios to determine the origin of long-term curvature seen in the RVs in addition to HD 88986 b. Stellar activity or a wide-orbit companion are two possibilities. To explore them, we used long-term photometric data obtained with the T8 APT, the combined RVs from ELODIE, HIRES, and SOPHIE instruments, as well as astrometric data from Gaia and HIPPARCOS.

### 9.1. APT photometric observations

To characterize the origin of long-term curvature observed in RVs (see Fig. 1), we used 1335 photometric observations of HD 88986 covering 21 observing seasons from 1995–1996 to 2019–2020, except the four observing seasons 2015–16 through 2018–19, during which the star was not observed. The observations were acquired with the T8 0.80 m APT at Fairborn Observatory in southern Arizona. The T8 APT is equipped with a two-channel photometer that uses two EMI 9124QB bi-alkali photomultiplier tubes to measure the stellar brightness simultaneously in the Strömgren *b* and *y* passbands.

The observations are made differentially with respect to three nearby comparison stars. We measured the difference in brightness between our program star HD 88986 (star d) and the comparison stars (stars a: HD 89557 (G = 7.3 mag, G8 III), b: HD 87667 (G = 7.3, F5), and c: HD 88476 (G = 6.6, G8 III)) and created differential magnitudes in the following six combinations: d-a, d-b, d-c, c-a, c-b, and b-a. Intercomparison of these six light curves shows that the comparison star a (HD 89557) is the only one that appears to be constant to the limit of our precision, so we present our results as differential magnitudes in the sense star d minus star a, which we designate as d-a.

To improve the photometric precision of the individual nightly observations, we combined the differential *b* and *y* magnitudes into a single (*b* + *y*)/2 passband. The precision of a single differential observation with T8, as measured from pairs of constant comparison stars, typically ranges between 0.001 mag and 0.0015 mag on good nights. The T8 APT is described in Henry (1999), where further details of the telescope, precision photometer, and observing and data reduction procedures can be found.

Figure 9 plots the 1335 nightly observations in the (*b* + *y*)/2 passband photometry from d-a obtained across the 21 observing seasons as small filled circles. The mean of all the nightly observations, -1.16492 mag, is plotted as the dashed line in the figure. The standard deviation of the nightly observations from their mean is 0.00118 mag, consistent with the precision of the measurements. The 21 seasonal means of these data are plotted as large filled circles. The standard deviations of the individual seasonal means are roughly the size of the plot symbols. The standard deviation of the 21 seasonal means from the mean of the seasonal means is 0.00028 mag, indicating that there is no long-term variability in HD 88986 to the limit of our photometric precision.

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Nightly Strömgren (*b* + *y*)/2 band photometry of HD 88986 from 21 observing seasons from 1995–96 to 2020–21 (small filled circles) scatter about their mean (dashed line) with a standard deviation of 0.00118 mag. Seasonal means from the 21 seasons (large filled circles) scatter about their mean with a standard deviation of 0.00028 mag. No significant variations nor periodicities are detected.</figcaption></figure>

Table L.1 summarizes observations in the (*b* + *y*)/2 passband photometry from d-a. The standard deviations of the nightly observations for each observing season indicate little or no short-term variability within each observing season. Frequency analysis of each individual observing season using the method of Vaníček (1971) confirms the lack of any periodic variability. Henry et al. (2022) show extensive examples of this method of period analysis.

In a study conducted by Lovis et al. (2011), the activity cycles of 311 FGK stars were analyzed. It shows that RV semi-amplitudes can be induced up to approximately 25 m s⁻¹ by the stellar long-period activity. However, in the case of HD 88986, the long-period RV semi-amplitude is at least about 40 m. s⁻¹ (see the SOPHIE+ RVs in Fig. 7), and also as stated in this section, there is no evidence of long-term photometric variability up to the limits of our precision. This suggests that the origin of this long-term curvature is likely unrelated to stellar activity and instead points toward the presence of a third body in the system.

### 9.2. Combining RV and HIPPARCOS/Gaia astrometry data

In order to improve the characterization of the outer massive companion, we modeled simultaneously the available RV and absolute astrometry data. The use of an MCMC algorithm enables us to explore the different solutions for each orbital parameter and for the companion mass compatible with the data. This algorithm, introduced by Philipot et al. (2023), is used to fit Keplerian orbits based on the emcee 3.0 (Foreman-Mackey et al. 2013) and HTOF (Brandt et al. 2021a) packages. The likelihood computation is similar to that of the ORVARA code (Brandt et al. 2021b).

We considered the ELODIE, SOPHIE, SOPHIE+, HIRES, and HIRES+ RV data, presented previously, coupled with the proper motion and position values calculated in the HIPPARCOS-Gaia Catalog of Accelerations (Brandt 2021) from HIPPARCOS (Perryman et al. 1997; van Leeuwen 2007) and Gaia data release 3 (DR3; Gaia Collaboration 2016, 2021) measurements. For the fit, we considered a Gaussian prior for the stellar mass and parallax, based on the values published
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Orbital fits for HD 88986 outer massive companion. Top: fit of the HD 88986 RV data points. Bottom: fit of the HD 88986 proper motion measurements in right ascension (left) and declination (right). The black points correspond to the HIPPARCOS and Gaia EDR3 data points. In each plot, the black curve corresponds to the best fit. The color bar indicates the log-likelihood corresponding to the different fits plotted.</figcaption></figure>

by Kervella et al. (2022), and a sin(I) prior for the orbital inclination. For the semi-major axis (a), the companion mass, the eccentricity, the longitude of the ascending node, the argument of periastron, the periastron passage time, and the jitter, we set uniform priors. In addition, as we use RV data from different instruments, we added an instrument offset for each dataset, also with uniform priors.

As the mass of HD 88986 b is low and its orbital period much smaller than the HIPPARCOS and Gaia DR3 observation windows (1227 and 1038 d, respectively), the proper motion variation of HD 88986 induced by the planet HD 88986 b is negligible. We have therefore only fitted the orbit of the outer massive companion (Fig. 10). However, as the RV data covers only a small part of the RV variation due to the outer companion, the star's RV remains poorly constrained and a wide range of solutions is compatible with the data, with similar likelihood values. We thus obtain an interval, with a confidence index of $3\sigma$, between 16.7 and 38.8 au for the semi-major axis, 16 and 169° for the orbital inclination, and 68 and 284 $M_{Jup}$ for the true mass

<table><caption>**Table 5.** Median values and 68% confidence interval for parameters of the outer companion of HD 88986 based on the joint analysis of the HIPPARCOS/Gaia astrometric and RV data.</caption><thead><tr><th>Parameter (unit)</th><th>The outer companion posteriors</th></tr></thead><tbody><tr><td colspan="2"><strong>Stellar parameters</strong></td></tr><tr><td><em>M</em><sub>*</sub> (<em>M</em><sub>☉</sub>)</td><td>1.20<sup>+0.07</sup><sub>-0.06</sub></td></tr><tr><td>Parallax (mas)</td><td>30.025 ± 0.023</td></tr><tr><td><em>σ</em><sub>*</sub> (m s<sup>-1</sup>)</td><td>3.3<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td colspan="2"><strong>The outer companion parameters</strong></td></tr><tr><td><em>P</em> (yr)</td><td>116<sup>+40</sup><sub>-34</sub></td></tr><tr><td><em>T</em><sub>c</sub> (BJD-2 400 000 d)</td><td>65 000 ± 2000</td></tr><tr><td><em>a</em> (au)</td><td>26.2<sup>+6.4</sup><sub>-5.5</sub></td></tr><tr><td>√<em>e</em> cos ω</td><td>0.47<sup>+0.16</sup><sub>-0.24</sub></td></tr><tr><td>√<em>e</em> sin ω</td><td>-0.47<sup>+0.13</sup><sub>-0.10</sub></td></tr><tr><td><em>e</em></td><td>0.46 ± 0.13</td></tr><tr><td>ω (°)</td><td>314.9<sup>+14.4</sup><sub>-21.4</sub></td></tr><tr><td><em>I</em> (°)</td><td>54.5<sup>+21.3</sup><sub>-18.8</sub> or 135.5<sup>+18.4</sup><sub>-19.4</sub></td></tr><tr><td>Ω (°)</td><td>29.1<sup>+8.5</sup><sub>-3.6</sub></td></tr><tr><td><em>M</em><sub>c</sub> (<em>M</em><sub>Jup</sub>)</td><td>145<sup>+73</sup><sub>-48</sub></td></tr><tr><td colspan="2"><strong>Instrumental parameters:</strong></td></tr><tr><td><em>mu</em><sub>ELODIE</sub> (m s<sup>-1</sup>)</td><td>29 220<sup>+180</sup><sub>-110</sub></td></tr><tr><td><em>mu</em><sub>SOPHIE</sub> (m s<sup>-1</sup>)</td><td>29 320<sup>+180</sup><sub>-110</sub></td></tr><tr><td><em>mu</em><sub>SOPHIE+</sub> (m s<sup>-1</sup>)</td><td>29 320<sup>+180</sup><sub>-110</sub></td></tr><tr><td><em>mu</em><sub>HIRES</sub> (m s<sup>-1</sup>)</td><td>240<sup>+180</sup><sub>-110</sub></td></tr><tr><td><em>mu</em><sub>HIRES+</sub> (m s<sup>-1</sup>)</td><td>240<sup>+180</sup><sub>-110</sub></td></tr></tbody></table>

of the companion (Table 5). Nevertheless, these results suggest that the outer massive companion is likely a brown dwarf or a low-mass star.

As previously mentioned, there is a Gaia DR3 source situated 1.4 arcseconds west of HD 88986. Considering HD 88986's parallax value for this source, its semi-major axis deviates by approximately $3.2\sigma$ from the resulting semi-major axis of the massive companion (see Table 5). Given the compatibility of the results with a wide range of solutions, it remains uncertain whether this source is the cause of the observed acceleration. Notably, this Gaia source lacks parallax information, raising the possibility that it might be a projected neighbor, unrelated to HD 88986.

### 9.3. Other constraints from Gaia astrometric excess noise

We used the Gaia data simulator from the gaston code first developed for the Gaia DR1 (Kiefer et al. 2019, 2021; Kiefer 2019) to test whether astrometric excess noises (AEN, hereafter) from the Gaia DR3 could lead to complementary mass constraints on the outer companion of HD 88986 at the orbital period found by coupling to the proper motions of HIPPARCOS and Gaia. The AEN is a measurement of supplementary motion, beyond proper motion and parallax, in the astrometric data of a source. The AEN is obtained from the RMS of residuals after fitting out the RA–Dec position, proper motion, and parallax to the simulated astrometric Gaia measurements by the approximate formula (see also Kiefer et al. 2019)
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Orbital fits for HD 88986 outer massive companion. Top: fit of the HD 88986 RV data points. Bottom: fit of the HD 88986 proper motion measurements in right ascension (left) and declination (right). The black points correspond to the HIPPARCOS and Gaia EDR3 data points. In each plot, the black curve corresponds to the best fit. The color bar indicates the log-likelihood corresponding to the different fits plotted.</figcaption></figure>

by Kervella et al. (2022), and a sin(I) prior for the orbital inclination. For the semi-major axis (a), the companion mass, the eccentricity, the longitude of the ascending node, the argument of periastron, the periastron passage time, and the jitter, we set uniform priors. In addition, as we use RV data from different instruments, we added an instrument offset for each dataset, also with uniform priors.

As the mass of HD 88986 b is low and its orbital period much smaller than the HIPPARCOS and Gaia DR3 observation windows (1227 and 1038 d, respectively), the proper motion variation of HD 88986 induced by the planet HD 88986 b is negligible. We have therefore only fitted the orbit of the outer massive companion (Fig. 10). However, as the RV data covers only a small part of the RV variation due to the outer companion, the star's RV remains poorly constrained and a wide range of solutions is compatible with the data, with similar likelihood values. We thus obtain an interval, with a confidence index of $3\sigma$, between 16.7 and 38.8 au for the semi-major axis, 16 and 169° for the orbital inclination, and 68 and 284 $M_{Jup}$ for the true mass

<table><caption>**Table 5.** Median values and 68% confidence interval for parameters of the outer companion of HD 88986 based on the joint analysis of the HIPPARCOS/Gaia astrometric and RV data.</caption><thead><tr><th>Parameter (unit)</th><th>The outer companion posteriors</th></tr></thead><tbody><tr><td colspan="2"><strong>Stellar parameters</strong></td></tr><tr><td><em>M</em><sub>*</sub> (<em>M</em><sub>☉</sub>)</td><td>1.20<sup>+0.07</sup><sub>-0.06</sub></td></tr><tr><td>Parallax (mas)</td><td>30.025 ± 0.023</td></tr><tr><td><em>σ</em><sub>*</sub> (m s<sup>-1</sup>)</td><td>3.3<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td colspan="2"><strong>The outer companion parameters</strong></td></tr><tr><td><em>P</em> (yr)</td><td>116<sup>+40</sup><sub>-34</sub></td></tr><tr><td><em>T</em><sub>c</sub> (BJD-2 400 000 d)</td><td>65 000 ± 2000</td></tr><tr><td><em>a</em> (au)</td><td>26.2<sup>+6.4</sup><sub>-5.5</sub></td></tr><tr><td>√<em>e</em> cos ω</td><td>0.47<sup>+0.16</sup><sub>-0.24</sub></td></tr><tr><td>√<em>e</em> sin ω</td><td>-0.47<sup>+0.13</sup><sub>-0.10</sub></td></tr><tr><td><em>e</em></td><td>0.46 ± 0.13</td></tr><tr><td>ω (°)</td><td>314.9<sup>+14.4</sup><sub>-21.4</sub></td></tr><tr><td><em>I</em> (°)</td><td>54.5<sup>+21.3</sup><sub>-18.8</sub> or 135.5<sup>+18.4</sup><sub>-19.4</sub></td></tr><tr><td>Ω (°)</td><td>29.1<sup>+8.5</sup><sub>-3.6</sub></td></tr><tr><td><em>M</em><sub>c</sub> (<em>M</em><sub>Jup</sub>)</td><td>145<sup>+73</sup><sub>-48</sub></td></tr><tr><td colspan="2"><strong>Instrumental parameters:</strong></td></tr><tr><td><em>mu</em><sub>ELODIE</sub> (m s<sup>-1</sup>)</td><td>29 220<sup>+180</sup><sub>-110</sub></td></tr><tr><td><em>mu</em><sub>SOPHIE</sub> (m s<sup>-1</sup>)</td><td>29 320<sup>+180</sup><sub>-110</sub></td></tr><tr><td><em>mu</em><sub>SOPHIE+</sub> (m s<sup>-1</sup>)</td><td>29 320<sup>+180</sup><sub>-110</sub></td></tr><tr><td><em>mu</em><sub>HIRES</sub> (m s<sup>-1</sup>)</td><td>240<sup>+180</sup><sub>-110</sub></td></tr><tr><td><em>mu</em><sub>HIRES+</sub> (m s<sup>-1</sup>)</td><td>240<sup>+180</sup><sub>-110</sub></td></tr></tbody></table>

of the companion (Table 5). Nevertheless, these results suggest that the outer massive companion is likely a brown dwarf or a low-mass star.

As previously mentioned, there is a Gaia DR3 source situated 1.4 arcseconds west of HD 88986. Considering HD 88986's parallax value for this source, its semi-major axis deviates by approximately $3.2\sigma$ from the resulting semi-major axis of the massive companion (see Table 5). Given the compatibility of the results with a wide range of solutions, it remains uncertain whether this source is the cause of the observed acceleration. Notably, this Gaia source lacks parallax information, raising the possibility that it might be a projected neighbor, unrelated to HD 88986.

### 9.3. Other constraints from Gaia astrometric excess noise

We used the Gaia data simulator from the gaston code first developed for the Gaia DR1 (Kiefer et al. 2019, 2021; Kiefer 2019) to test whether astrometric excess noises (AEN, hereafter) from the Gaia DR3 could lead to complementary mass constraints on the outer companion of HD 88986 at the orbital period found by coupling to the proper motions of HIPPARCOS and Gaia. The AEN is a measurement of supplementary motion, beyond proper motion and parallax, in the astrometric data of a source. The AEN is obtained from the RMS of residuals after fitting out the RA–Dec position, proper motion, and parallax to the simulated astrometric Gaia measurements by the approximate formula (see also Kiefer et al. 2019)
<table><caption>**Table 6.** Resulting constraints on the orbital inclination and mass of companion “c”, and on the predicted photocenter semi-major axis of HD 88986, using the AEN from Gaia DR3.</caption><thead><tr><th>Parameter</th><th>1σ</th><th>3σ limits</th></tr></thead><tbody><tr><td>I<sub>c</sub> (deg)</td><td>33<sup>+33</sup><sub>-18</sub></td><td>>2.7</td></tr><tr><td>M<sub>c</sub> (M<sub>jup</sub>)</td><td>180<sup>+220</sup><sub>-80</sub></td><td>&lt;810</td></tr><tr><td>a<sub>phot</sub> (mas)</td><td>91<sup>+107</sup><sub>-47</sub></td><td>&lt;440</td></tr></tbody></table>

and references therein):

$$ \sigma_{\text{AL}}^2 + \sigma_{\text{attitude}}^2 + \text{AEN}_{\text{DR3}}^2 = \frac{\sum_j R_j^2}{N-5} \quad (2) $$

where $R_j$ are the N along-scan (AL) angle residuals of the astrometric fit; $\sigma_{\text{AL}}$ is the typical AL angle measurement noise; $\sigma_{\text{attitude}}$ is the spacecraft attitude excess noise, and $\text{AEN}_{\text{DR3}}$ is the AEN. The AL angle measurement noise has a value of $\sigma_{\text{AL}} = 0.05$ mas for targets with a G-magnitude of 6.3 (Fig. A.1 from Lindegren et al. 2021), and the typical attitude noise in the DR3 is $\sigma_{\text{attitude}} = 0.076$ mas (Lindegren et al. 2021).

HD 88986 has a magnitude of $G_{\text{mag}} \sim 6.3$ and a color Gb-Gr of $\sim 0.8$. In the Gaia-DR3 catalog, the typical AEN of single stars at that magnitude and color for sources fitted with five parameters, as HD 88986, is 0.14 mas. This nonzero AEN is present for nearly all single stars and is due to a systematic jitter, including instrumental and global modelization noises, that is accounted for in the formal errors used to calculate the $\chi^2$ (Lindegren et al. 2021). The AEN of HD 88986 is $\epsilon_{\text{DR3}} = 0.135$ mas. The Gaia DR3 astrometry of this target is thus compatible with a single star without a companion, but it also allows deriving an upper-limit constraint on the mass of the RV-detected companion, given a range of possible orbital periods.

We follow the method from Kiefer et al. (2019, 2021), using the code gaston adapted to the (E)DR3. The general principle of the method is the same as with the DR1. Fixing $P$, $m \sin i$, $e$, $\omega$ and $T_0$ within their priors derived from combined RVs in Table 5, we run several simulations of Gaia measurements of the target along a model of the orbital motion of the system due to the outer massive companion and derive simulated values of AEN that we compare with the actual $\text{AEN}_{\text{DR3}}$.

We sample orbital inclination uniformly between 0 and 90° by an MCMC routine based on the emcee code (Foreman-Mackey et al. 2013) and thoroughly explained in Kiefer et al. (2019, 2021). The orbital inclination changes the amplitude of the astrometric motion due to a different mass of the companion determined from $M = m \sin i/\sin i$ and thus changes the value of the AEN allowing us to match a range of orbital inclinations to the observed AEN.

Noises, epochs, scan angles, and the number of measurements used in the simulations are updated with respect to the new data reduction of DR3 (Gaia Collaboration 2021). An epoch is a date when the star is transiting the Gaia field of view; several measurements, typically 9, are performed during a single transit. Those epochs can be found for any target in the Gaia Observation Forecast Tool (or GOST<fnref n="6" />). We add in our simulated model a jitter of 0.16 mas, allowing us to reproduce a median AEN of 0.14 mas for single sources at $G = 6.3$ and Gb-Gr = 0.8. It is

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Illustration of HD 88986 outer companion analysis using Gaia DR3 analysis showing the companion mass posterior distribution running gaston. The dotted line shows the 3σ upper-limit, the dashed lines show the 1σ confidence interval, the solid blue line is the median mass, and the solid red line shows the RV $m \sin i$.</figcaption></figure>

modeled as a Gaussian noise changing every epoch of observation. The spacecraft attitude noise is also added to the model as a systematic Gaussian dispersion that changes every observation epoch with a standard deviation of 0.076 mas. A Gaussian measurement noise of $\sigma_{\text{AL}} = 0.05$ mas is added to each of the $N_{\text{AL}}$ astrometric measurements performed at a given epoch.

Table 6 summarises the results of AEN fitting for this star. Fig. 11 shows the relation between AEN and inclination in the simulations and plots the posterior distribution of companion mass. The posterior distribution on mass gives an upper limit on the mass of the companion below 810 $M_{\text{jup}}$ at 3σ. This result agrees with the one presented in Sect. 9.2. Finally, we adopted the results from Sect. 9.2 as it provides a higher level of precision in the mass and orbital inclination of the outer massive companion.

## 10. Discussion and conclusion

We discovered and have characterized HD 88986 b, a sub-Neptune in orbit around a subgiant star, which stands as one of the nearest and brightest ($V = 6.47$ mag) exoplanet host stars (see Fig. 12 top right). Our analysis indicates that this planet is transiting, based on two potential single transit detections in TESS sectors 21 and 48, both of which are consistent with the anticipated transit time from the RV model. By combining data from SOPHIE+ RV measurements and TESS sector 21 photometric data, we determined the following parameters for HD 88986 b: a period of $P_p = 146.05_{-0.40}^{+0.43}$ d, a mass of $M_p = 17.2_{-3.8}^{+4.0} M_{\oplus}$, and a radius of $R_p = 2.49 \pm 0.18 R_{\oplus}$, resulting in a high mean density of $\rho_p = 6.1_{-2.3}^{+3.3} \text{ g cm}^{-3}$. The two-layer theoretical composition model developed by Zeng et al. (2016) indicates that the planet is composed predominantly of rock, accounting for approximately 75% of its mass, while water makes up the remaining 25% (Fig. 12 bottom left). Additional photometric observations of the system targeting another transit event of HD 88986 b are needed. Such observations would provide a strong confirmation of the planet's transiting nature and yield better estimates for its period and radius.

Additionally, we identified a clear long-term curvature in the RV caused by the presence of a massive companion in the system. The nature of this companion has yet to be confirmed. A

<footnote n="6">⁶ https://gaia.esac.esa.int/gost/index.jsp</footnote>
<table><caption>**Table 6.** Resulting constraints on the orbital inclination and mass of companion “c”, and on the predicted photocenter semi-major axis of HD 88986, using the AEN from Gaia DR3.</caption><thead><tr><th>Parameter</th><th>1σ</th><th>3σ limits</th></tr></thead><tbody><tr><td>I<sub>c</sub> (deg)</td><td>33<sup>+33</sup><sub>-18</sub></td><td>>2.7</td></tr><tr><td>M<sub>c</sub> (M<sub>jup</sub>)</td><td>180<sup>+220</sup><sub>-80</sub></td><td>&lt;810</td></tr><tr><td>a<sub>phot</sub> (mas)</td><td>91<sup>+107</sup><sub>-47</sub></td><td>&lt;440</td></tr></tbody></table>

and references therein):

$$ \sigma_{\text{AL}}^2 + \sigma_{\text{attitude}}^2 + \text{AEN}_{\text{DR3}}^2 = \frac{\sum_j R_j^2}{N-5} \qquad (2) $$

where $R_j$ are the N along-scan (AL) angle residuals of the astrometric fit; $\sigma_{\text{AL}}$ is the typical AL angle measurement noise; $\sigma_{\text{attitude}}$ is the spacecraft attitude excess noise, and $\text{AEN}_{\text{DR3}}$ is the AEN. The AL angle measurement noise has a value of $\sigma_{\text{AL}} = 0.05$ mas for targets with a G-magnitude of 6.3 (Fig. A.1 from Lindegren et al. 2021), and the typical attitude noise in the DR3 is $\sigma_{\text{attitude}} = 0.076$ mas (Lindegren et al. 2021).

HD 88986 has a magnitude of $G_{\text{mag}} \sim 6.3$ and a color Gb-Gr of $\sim 0.8$. In the Gaia-DR3 catalog, the typical AEN of single stars at that magnitude and color for sources fitted with five parameters, as HD 88986, is 0.14 mas. This nonzero AEN is present for nearly all single stars and is due to a systematic jitter, including instrumental and global modelization noises, that is accounted for in the formal errors used to calculate the $\chi^2$ (Lindegren et al. 2021). The AEN of HD 88986 is $\epsilon_{\text{DR3}} = 0.135$ mas. The Gaia DR3 astrometry of this target is thus compatible with a single star without a companion, but it also allows deriving an upper-limit constraint on the mass of the RV-detected companion, given a range of possible orbital periods.

We follow the method from Kiefer et al. (2019, 2021), using the code gaston adapted to the (E)DR3. The general principle of the method is the same as with the DR1. Fixing $P$, $m \sin i$, $e$, $\omega$ and $T_0$ within their priors derived from combined RVs in Table 5, we run several simulations of Gaia measurements of the target along a model of the orbital motion of the system due to the outer massive companion and derive simulated values of AEN that we compare with the actual $\text{AEN}_{\text{DR3}}$.

We sample orbital inclination uniformly between 0 and 90° by an MCMC routine based on the emcee code (Foreman-Mackey et al. 2013) and thoroughly explained in Kiefer et al. (2019, 2021). The orbital inclination changes the amplitude of the astrometric motion due to a different mass of the companion determined from $M = m \sin i/\sin i$ and thus changes the value of the AEN allowing us to match a range of orbital inclinations to the observed AEN.

Noises, epochs, scan angles, and the number of measurements used in the simulations are updated with respect to the new data reduction of DR3 (Gaia Collaboration 2021). An epoch is a date when the star is transiting the Gaia field of view; several measurements, typically 9, are performed during a single transit. Those epochs can be found for any target in the Gaia Observation Forecast Tool (or GOST<fnref n="6" />). We add in our simulated model a jitter of 0.16 mas, allowing us to reproduce a median AEN of 0.14 mas for single sources at $G = 6.3$ and Gb-Gr = 0.8. It is

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Illustration of HD 88986 outer companion analysis using Gaia DR3 analysis showing the companion mass posterior distribution running gaston. The dotted line shows the $3\sigma$ upper-limit, the dashed lines show the $1\sigma$ confidence interval, the solid blue line is the median mass, and the solid red line shows the RV $m \sin i$.</figcaption></figure>

modeled as a Gaussian noise changing every epoch of observation. The spacecraft attitude noise is also added to the model as a systematic Gaussian dispersion that changes every observation epoch with a standard deviation of 0.076 mas. A Gaussian measurement noise of $\sigma_{\text{AL}} = 0.05$ mas is added to each of the $N_{\text{AL}}$ astrometric measurements performed at a given epoch.

Table 6 summarises the results of AEN fitting for this star. Fig. 11 shows the relation between AEN and inclination in the simulations and plots the posterior distribution of companion mass. The posterior distribution on mass gives an upper limit on the mass of the companion below 810 $M_{\text{jup}}$ at $3\sigma$. This result agrees with the one presented in Sect. 9.2. Finally, we adopted the results from Sect. 9.2 as it provides a higher level of precision in the mass and orbital inclination of the outer massive companion.

# 10. Discussion and conclusion

We discovered and have characterized HD 88986 b, a sub-Neptune in orbit around a subgiant star, which stands as one of the nearest and brightest ($V = 6.47$ mag) exoplanet host stars (see Fig. 12 top right). Our analysis indicates that this planet is transiting, based on two potential single transit detections in TESS sectors 21 and 48, both of which are consistent with the anticipated transit time from the RV model. By combining data from SOPHIE+ RV measurements and TESS sector 21 photometric data, we determined the following parameters for HD 88986 b: a period of $P_p = 146.05^{+0.43}_{-0.40}$ d, a mass of $M_p = 17.2^{+4.0}_{-3.8} M_{\oplus}$, and a radius of $R_p = 2.49 \pm 0.18 R_{\oplus}$, resulting in a high mean density of $\rho_p = 6.1^{+3.3}_{-2.3} \text{ g cm}^{-3}$. The two-layer theoretical composition model developed by Zeng et al. (2016) indicates that the planet is composed predominantly of rock, accounting for approximately 75% of its mass, while water makes up the remaining 25% (Fig. 12 bottom left). Additional photometric observations of the system targeting another transit event of HD 88986 b are needed. Such observations would provide a strong confirmation of the planet's transiting nature and yield better estimates for its period and radius.

Additionally, we identified a clear long-term curvature in the RV caused by the presence of a massive companion in the system. The nature of this companion has yet to be confirmed. A

<footnote n="6">⁶ https://gaia.esac.esa.int/gost/index.jsp</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Position of HD 88986 b and its host star among known transiting exoplanetary systems, with precise planet mass and radius measurements ($\sigma M/M = 25\%$ and $\sigma R/R = 8\%$, Otegi et al. 2020) from the NASA Exoplanet Data Archive (June 7, 2023). Top left: radius-period diagram of exoplanets. Top right: brightness in the K band versus distance of the exoplanet host star for the same planets. Bottom left: mass-radius diagrams of small planets ($R_p < 4 R_⊕$). The colored curves are the two-layer theoretical composition models of Zeng et al. (2016). Bottom right: equilibrium temperature-radius diagram of planets within the radius range of $1.5 R_⊕$ to $2.75 R_⊕$ with each planet's size scaled according to the propositional TSM introduced by Kempton et al. (2018). The gray area is the proposed hazy atmosphere zone by Yu et al. (2021). The purple mark indicates HD 88986 b. These planets possess either a known equilibrium temperature or information about the radius and temperature of their host star to estimate their temperature ($T_{eq} \sim T_* \sqrt{(R_*/2a)}$, Méndez & Rivera-Valentín 2017). These figures highlight the unique position of the HD 88986 system.</figcaption></figure>

joint analysis of RV, HIPPARCOS, and Gaia astrometric data shows that with a $3\sigma$ confidence interval, its semi-major axis is between 16.7 and 38.8 au and its mass is between 68 and 284 $M_{Jup}$. SOPHIE+ observations are being conducted to disclose the nature of this massive companion, and in particular better constrain its period and eccentricity. Furthermore, given its extensive semi-major axis, this outer massive companion presents an opportunity for being directly imaged, aiming to provide a more precise characterization of its orbit and mass. This study can be facilitated using the current generation of high-contrast imaging instruments, such as SPHERE.

The top-left panel of Fig. 12 highlights the unique position of HD 88986 b in the radius-period diagram among other known planets from the NASA Exoplanet Data Archive<fnref n="7" /> (as of June 7, 2023) with precise mass and radius measurements ($\sigma M/M = 25\%$ and $\sigma R/R = 8\%$, Otegi et al. 2020). Notably, HD 88986 b

has the longest orbital period among the discovered small transiting planets ($R_p < 4 R_⊕$). This wide orbit suggests that the planet did not undergo significant mass loss due to extreme-ultraviolet radiation, and hence it probably retains its primordial composition (Kubyshkina & Fossati 2022). Consequently, HD 88986 b is an excellent candidate for investigating the planet's internal structure and formation conditions.

The bottom-right panel of Fig. 12 compares the equilibrium temperatures of HD 88986 b with those of other known small planets possessing a precise mass and radius (Otegi et al. 2020). HD 88986 b, thanks to its long orbital period, is a cold planet ($T_{eq} = 460 \pm 8$ K). The study of atmospheric characteristics of cold planets ($\le 500$ K) transiting bright hosts is extremely limited by the lack of such planets. In terms of atmospheric chemistry, colder atmospheres may contain disequilibrium chemistry that is useful for understanding atmospheric physics (Fortney et al. 2020). Remarkably, HD 88986 b's equilibrium temperature places it within the proposed hazy atmosphere zone, ranging between 270 and 600 K, as suggested by Yu et al. (2021).

<footnote n="7">7 https://exoplanetarchive.ipac.caltech.edu/</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Position of HD 88986 b and its host star among known transiting exoplanetary systems, with precise planet mass and radius measurements ($\sigma M/M = 25\%$ and $\sigma R/R = 8\%$, Otegi et al. 2020) from the NASA Exoplanet Data Archive (June 7, 2023). Top left: radius-period diagram of exoplanets. Top right: brightness in the K band versus distance of the exoplanet host star for the same planets. Bottom left: mass-radius diagrams of small planets ($R_p < 4 R_⊕$). The colored curves are the two-layer theoretical composition models of Zeng et al. (2016). Bottom right: equilibrium temperature-radius diagram of planets within the radius range of $1.5 R_⊕$ to $2.75 R_⊕$ with each planet's size scaled according to the propositional TSM introduced by Kempton et al. (2018). The gray area is the proposed hazy atmosphere zone by Yu et al. (2021). The purple mark indicates HD 88986 b. These planets possess either a known equilibrium temperature or information about the radius and temperature of their host star to estimate their temperature ($T_{eq} \sim T_* \sqrt{(R_*/2a)}$, Méndez & Rivera-Valentín 2017). These figures highlight the unique position of the HD 88986 system.</figcaption></figure>

joint analysis of RV, HIPPARCOS, and Gaia astrometric data shows that with a $3\sigma$ confidence interval, its semi-major axis is between 16.7 and 38.8 au and its mass is between 68 and 284 $M_{Jup}$. SOPHIE+ observations are being conducted to disclose the nature of this massive companion, and in particular better constrain its period and eccentricity. Furthermore, given its extensive semi-major axis, this outer massive companion presents an opportunity for being directly imaged, aiming to provide a more precise characterization of its orbit and mass. This study can be facilitated using the current generation of high-contrast imaging instruments, such as SPHERE.

The top-left panel of Fig. 12 highlights the unique position of HD 88986 b in the radius-period diagram among other known planets from the NASA Exoplanet Data Archive<fnref n="7" /> (as of June 7, 2023) with precise mass and radius measurements ($\sigma M/M = 25\%$ and $\sigma R/R = 8\%$, Otegi et al. 2020). Notably, HD 88986 b

has the longest orbital period among the discovered small transiting planets ($R_p < 4 R_⊕$). This wide orbit suggests that the planet did not undergo significant mass loss due to extreme-ultraviolet radiation, and hence it probably retains its primordial composition (Kubyshkina & Fossati 2022). Consequently, HD 88986 b is an excellent candidate for investigating the planet's internal structure and formation conditions.

The bottom-right panel of Fig. 12 compares the equilibrium temperatures of HD 88986 b with those of other known small planets possessing a precise mass and radius (Otegi et al. 2020). HD 88986 b, thanks to its long orbital period, is a cold planet ($T_{eq} = 460 \pm 8$ K). The study of atmospheric characteristics of cold planets ($\le 500$ K) transiting bright hosts is extremely limited by the lack of such planets. In terms of atmospheric chemistry, colder atmospheres may contain disequilibrium chemistry that is useful for understanding atmospheric physics (Fortney et al. 2020). Remarkably, HD 88986 b's equilibrium temperature places it within the proposed hazy atmosphere zone, ranging between 270 and 600 K, as suggested by Yu et al. (2021).

<footnote n="7">7 https://exoplanetarchive.ipac.caltech.edu/</footnote>
This intriguing positioning opens up exciting opportunities for studying the haze layer and atmospheric composition above it (Kawashima et al. 2019). Hence, for the bottom-right panel of Fig. 12, we examined the proportional transmission spectroscopy metric (TSM) $\propto (T_{eq} \times R_j^2 / M_p \times R_s^2) \times 10^{-mJ/5}$ introduced by Kempton et al. (2018) and scaled the size of each planet accordingly. This metric considers the planet's equilibrium temperature $T_{eq}$, radius $R_p$, mass $M_p$, host star radius $R_s$, and host star magnitude in the J band $m_J$. We selected a radius bin of 1.5 $R_⊕$ to 2.75 $R_⊕$ based on Kempton et al. (2018)'s assumption that planets within the same size bin share similar atmospheric compositions. We exclusively employed the proportional TSM because the constant scale factor in Kempton et al. (2018)'s TSM is intended for stars with $m_J > 9$, while HD 88986 has a magnitude of $m_J = 5.2$ mag. HD 88986 b ranked 19th in comparison with proportional TSM S/N of only 24 planets detected in the hazy atmosphere zone (see Fig. 12 bottom-right panel). This relatively low rank is to some extent compensated for by the exceptionally long duration of the transit (16 h). Therefore, HD 88986 b's unique characteristics, such as being a cold exoplanet orbiting a bright star at a distance of 33 pc, make it a good target for atmospheric characterization studies of cold planets residing in the hazy zone.

Furthermore, with a mass of $M_p = 17.2_{-3.8}^{+4.0} M_⊕$, HD 88986 b surpasses the critical mass threshold of ~10 $M_⊕$ required for envelope accretion (Johnson et al. 2010). This indicates that it likely formed similarly to the cores of giant planets in our Solar System. However, HD 88986 b failed to accumulate much gas during its formation process. One possible scenario is that HD 88986 b formed at a late stage in the protoplanetary disk when there was little gas present during core assembly (Lee & Chiang 2016). Moreover, according to the minimum mass solar nebular model, it is unlikely for such a massive planet to form in situ at its current location, situated 0.6 au from its host star (Schlichting 2014). Instead, it likely formed farther away and subsequently migrated inward over time, potentially influenced by interactions with the detected massive companion in the system. However, to gain a comprehensive understanding of the HD 88986 planetary system, additional photometric and spectroscopic observations are required.

**Acknowledgements.** We warmly thank the OHP staff for their support on the observations. We received funding from the French Programme National de Physique Stellaire (PNPS) and the Programme National de Planétologie (PNP) of CNRS (INSU). N.H. acknowledges CNES postdoctoral funding fellowship. N.H. also acknowledges the financial support of the French embassy in Tehran as well as the Iran Ministry of Science Research and Technology. J.S.J. acknowledges support by FONDECYT grant 1201371 and from the ANID BASAL project FB210003. This paper made use of data collected by the TESS mission which is publicly available from the Mikulski Archive for Space Telescopes (MAST) operated by the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by NASA's Science Mission Directorate. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. NASA supported DR under award number NNA16BD14C for NASA Academic Mission Services. A.C., X.D., and T.F. acknowledge support by the French National Research Agency in the framework of the Investissement d'Avenir program (ANR-15-IDEX-02), through the funding of the « Origin of Life » project of the Grenoble-Alpes University. We acknowledge funding from the French ANR under contract number ANR18CE310019 (SPlaSH). This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UID/FIS/04434/2019, UIDB/04434/2020, UIDP/04434/2020, PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113, PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953, PTDC/ FIS-AST/28987/2017 & POCI-01-0145-FEDER-028987. N.C.S. further

acknowledges funding by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. O.V. acknowledges funding from the ANR project "EXACT" (ANR-21-CE49-0008-01), from the Centre National d'Études Spatiales (CNES), and from the CNRS/INSU Programme National de Planétologie (PNP). M.H. acknowledges support from ANID-Millennium Science Initiative-ICN12_009. S.D. is funded by the UK Science and Technology Facilities Council (grant number ST/V004735/1). 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 project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (project SPICE DUNE, grant agreement No 947634).

## References

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

Baluev, R. V. 2008, MNRAS, 385, 1279

Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373

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

Berger, T. A., Huber, D., Gaidos, E., & van Saders, J. L. 2018, ApJ, 866, 99

Bijaoui, A. 1980, A&A, 84, 81

Boisse, I., Eggenberger, A., Santos, N., et al. 2010, A&A, 523, A88

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

Borgniet, S., Meunier, N., & Lagrange, A. M. 2015, A&A, 581, A133

Bouchy, F., Hébrard, G., Udry, S., et al. 2009, A&A, 505, 853

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

Brandt, T. D. 2021, ApJS, 254, 42

Brandt, G. M., Michalik, D., Brandt, T. D., et al. 2021a, AJ, 162, 230

Brandt, T. D., Dupuy, T. J., Li, Y., et al. 2021b, AJ, 162, 186

Butler, R. P., Vogt, S. S., Laughlin, G., et al. 2017, AJ, 153, 208

Courcol, B., Bouchy, F., Pepe, F., et al. 2015, A&A, 581, A38

Cutri, R., Skrutskie, M., Van Dyk, S., et al. 2003, The IRSA 2MASS All-Sky Point Source Catalog

Da Silva, J. G., Santos, N., Bonfils, X., et al. 2011, A&A, 534, A30

Delisle, J.-B., Ségransan, D., Buchschacher, N., & Alesina, F. 2016, A&A, 590, A134

Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nat. Astron., 5, 775

Deming, D., Knutson, H., Kammer, J., et al. 2015, ApJ, 805, 132

Díaz, R., Santerne, A., Sahlmann, J., et al. 2012, A&A, 538, A113

Díaz, R. F., Ségransan, D., Udry, S., et al. 2016, A&A, 585, A134

Dobos, V., Charnoz, S., Pál, A., Roque-Bernard, A., & Szabó, G. M. 2021, PASP, 133, 094401

Ehrenreich, D., Delrez, L., Akinsanmi, B., et al. 2023, A&A, 671, A154

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

Fausnaugh, M. M., Burke, C. J., Caldwell, D. A., et al. 2020, TESS Data Release Notes: Sector 21, DR29, Technical Report

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

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220

Fortney, J. J., Visscher, C., Marley, M. S., et al. 2020, AJ, 160, 288

Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109

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

Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1

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

Gilliland, R. L., Chaplin, W. J., Dunham, E. W., et al. 2011, ApJS, 197, 6

Hall, J. C., Lockwood, G., & Skiff, B. A. 2007, AJ, 133, 862

Hara, N. C., & Delisle, J.-B. 2023, A&A, submitted [arXiv:2304.08489]

Hara, N. C., Bouchy, F., Stalport, M., et al. 2020, A&A, 636, L6

Hara, N. C., Delisle, J.-B., Unger, N., & Dumusque, X. 2022a, A&A, 658, A177

Hara, N. C., Unger, N., Delisle, J.-B., Díaz, R. F., & Ségransan, D. 2022b, A&A, 663, A14

Hawthorn, F., Bayliss, D., Wilson, T. G., et al. 2023, MNRAS, 520, 3649

Haywood, R., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517

Hébrard, G., Bouchy, F., Pont, F., et al. 2008, A&A, 488, 763

Heidari, N. 2022, PhD thesis, Université Côte d'Azur; Shahid Beheshti University (Tehran), Iran

Heidari, N., Boisse, I., Orell-Mique, J., et al. 2022, A&A 658, A176

Henry, G. W. 1999, PASP, 111, 845

Henry, G. W., Fekel, F. C., & Williamson, M. H. 2022, AJ, 163, 180

Hippke, M., & Heller, R. 2019, A&A, 623, A39

Hippke, M., David, T. J., Mulders, G. D., & Heller, R. 2019, AJ, 158, 143

Hobson, M. J. 2019, Ph.D. thesis, Aix-Marseille, France