Press Release
Open Access
Issue
A&A
Volume 681, January 2024
Article Number A55
Number of page(s) 37
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202347897
Published online 11 January 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 exoplan-ets 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 (Teff > 2800 K), is 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 astro-metric 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 Planetes Extrasolaires (RPE) 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 λλ ≈ 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σ 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−1, with a root mean square (RMS) of 15.30 m s−1. The final SOPHIE dataset is provided in Table D.1.

thumbnail Fig. 1

Radial velocity measurements of HD 88986 from ELODIE, HIRES, HIRES+, SOPHIE and SOPHIE+.

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 ms−1 resulting from a CCD change in August 2004, as well as a long-term drift (≲1 m s−1) 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 tool2 following the methodology outlined in Tal-Or et al. (2019). The mean uncertainty of the HIRES data is 1.2 m s−1, with a RMS of 11.1 m s−1. 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−1 and RMS of 13.0 m s−1, 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 (Teff), metallicity ([Fe/H]), and surface gravity (log ɡ), using the procedure described in Santos et al. (2013) and Sousa et al. (2018). The resulting Teff, [Fe/H], and log ɡ 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 Stassun & Torres (2016); Stassun et al. (2017); Stassun & Torres (2018). We obtained the BT VT magnitudes from Tycho-2 (Høg 2001), the JHKS 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 GBP GRP magnitudes from Gaia (Gaia Collaboration 2016, 2021). We also used the UV measurement at 274 nm from the TB1 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 Teff, log ɡ, and [Fe/H] adopted from the spectroscopic analysis. The remaining free parameter is the extinction AV, which we fixed at zero due to the proximity of the system to Earth.

The resulting fit (Fig. 2) has a reduced χ2 of 1.4. Integrating the model SED gives the bolometric flux at Earth, Fbol = 7.28 ± 0.17 × 10−8 erg s−1 cm−2. Taking the Fbol and Teff together with the Gaia parallax, gives the stellar radius, R = 1.543 ± 0.010 R, placing the star within the subgiant range (1.5–3 R; 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 = 1.19 ± 0.07 M, which is consistent with the value of 1.25 ± 0.05 M determined empirically via R and log ɡ. 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 RHK (see Sect. 4) to estimate the stellar age via empirical activity-age relations. We obtain an age of τ = 7.9 ± 1.3 Gyr via the empirical relations of Mamajek & Hillenbrand (2008).

Table 1

Stellar properties of HD 88986.

thumbnail 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).

4 Stellar rotation and activity

To study the host star’s activity and rotational period, we used log RHK, 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 RHK 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 RHK 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, λ ~ 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 RHK = −5.04 ± 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 RHK parameter, with values transitioning from −5.12 ± 0.10 during the initial subset to −5.05 ± 0.10 in the second subset, and ultimately stabilizing at −5.00 ± 0.10 in the final subset. All these values are in good agreement with the value of log RHK = −5.22 and log RHK = −5.07 reported by Radick et al. (2018) and Hall et al. (2007), respectively. Finally, we estimated a rotational period of d following Noyes et al. (1984), a value consistent with the Mamajek & Hillenbrand (2008) recipe, which yields 26.3 ± 3.1 d.

We searched for rotational modulation in the Simple Aperture 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, λ ~ 3955 Å) and significant contamination caused by the calibration lamp. Additionally, one data point was excluded due to its identification as a 5σ 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, λ ~ 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 d. Furthermore, there is an activity signal at 141.2 d in the periodogram of the HIRES and HIRES+ S-index, which will be further discussed in the following section.

thumbnail 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).

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−1 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 ms−1, τ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 ms−1, σR = 1.5 ms−1, τ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 ɡ 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; 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 Tobs/3 = 1211 days and Tobs/9 = 411 days where Tobs 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 ( 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 , 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.

thumbnail Fig. 4

1 periodogram of SOPHIE+ data. The identified periods are shown in red.

thumbnail 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 Tobs the total time span of observations, a) and b) are obtained with windows of size Tobs/3 = 1234.2 d and Tobs/9 = 411.42 d, respectively.

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 (Aln 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 Aln Z ≤ 2. In this case, the model with the fewest free parameters would be chosen.

We defined our RV model as follows: (1)

where K(t) is the Keplerian model and the , is a normal distribution (𝒩) of white-Gaussian noise where σ (t) is the uncertainty of each RV point at time t, and σw is a jitter term. Additionally, µ 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 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 , where σGP is the amplitude of GP modulation, and TGP presents the characteristic timescale (Ambikasaran et al. 2015);

  • a quasi-periodic GP kernel (QP-GP, Foreman-Mackey et al. 2017) with the form of , where BGP amplifies the kernel, CGP is a constant scaling term, LGP is a correlation time-scale component, and finally Prot;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.

Table 2

Different tested models on the SOPHIE+ RV-only data along with model comparisons with juliet.

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 Tc, 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 (Aln Z ≥ 97) which confirms a clear detection of the planetary signal. The eccentricity-free model, however, exhibits bimodality on the periastron argument ω. Given that the two models are statistically indistinguishable (Aln Z ≤ 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 Prot;Gp, we used a Gaussian prior centered at 29 d with a standard deviation of 3 d. On the CGP, we initially tried a wide range of Jeffreys priors from 10−20 to 100. The posterior distribution of CGP reached the prior boundary at 10−20, indicating that this parameter converged to zero. This result motivates us to set the CGP to 10−20, which is consistent with zero. We note that the models with fixed/unfixed CGP are statistically indistinguishable (∆ln Z ≤ 2). Therefore, we fixed CGP 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 (∆ ln Z ≥ 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 1eccentricityfree), 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 (∆ln Z ≥ ~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± 0.2 m s−1, Tc = 58903 ± 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−1, 5.6 m s−1, 6.2 m s−1, and 3.2 m s−1 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 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 ω, 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.

thumbnail 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.

6 First transit detected in TESS sector 21 in February 2020

HD 88986 was observed in TESS sector 214 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 Tc = 2 458 891.6 (corresponding to 2020 February 12), a duration of about 16 h, and a depth of ~ 220 ppm. Remarkably, this Tc is in agreement, within uncertainty, with all the Tc 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 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 Lightkurve5. 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 Tc 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 Tc of the potential single transit agrees, within uncertainties, with the Tc 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 (∆Gmag = 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 Tc. For Tc, 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 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 d. It exhibits an eccentricity of 0.23 ± 0.06, alongside a radius of 2.49 ± 0.18 R and a mass of M, corresponding to a high mean density of . Additionally, the planet has an equilibrium temperature of , Méndez & Rivera-Valentín 2017), making it a relatively cool planet.

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).

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 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.4s. 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 v 13; 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.

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 T0 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.

thumbnail 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.

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 Tc 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 σ = 204.4 ppm, TESS sector 48 displays a substantial dispersion with σ = 635.8 ppm. This noticeable difference might be linked to the 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 transitlike 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 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σ) with the period calculated using RVs and single transit in sector 21, and also agrees (at 2σ) 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.

thumbnail 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.

Table 4

Possible solution for HD 88986 b’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 bialkali 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.

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−1 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−1 (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.

thumbnail 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.

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 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 88986b 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σ, between 16.7 and 38.8 au for the semi-major axis, 16 and 169° for the orbital inclination, and 68 and 284 MJup for the true mass 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σ 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.

thumbnail 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.

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.

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 and references therein): (2)

where Rj are the N along-scan (AL) angle residuals of the astrometric fit; σAL is the typical AL angle measurement noise; σattitude is the spacecraft attitude excess noise, and AENDR3 is the AEN. The AL angle measurement noise has a value of σ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 σattitude = 0.076 mas (Lindegren et al. 2021).

HD 88986 has a magnitude of Gmag ~ 6.3 and a color Gb−Gr of ~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 χ2 (Lindegren et al. 2021). The AEN of HD 88986 is 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, ω and T0 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 AENDR3.

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 GOST6). 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 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 σAL = 0.05 mas is added to each of the NAL 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 Mjup 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.

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.

thumbnail 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.

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 , a mass of , and a radius of Rp = 2.49 ± 0.18 R, resulting in a high mean density of . 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 joint analysis of RV, HIPPARCOS, and Gaia astrometric data shows that with a 3σ confidence interval, its semi-major axis is between 16.7 and 38.8 au and its mass is between 68 and 284 MJup. 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 Archive7 (as of June 7, 2023) with precise mass and radius measurements (σM/M = 25% and σR/R = 8%, Otegi et al. 2020). Notably, HD 88986b has the longest orbital period among the discovered small transiting planets (Rp< 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 88986b 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 (Teq = 460 ± 8 K). The study of atmospheric characteristics of cold planets (≤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).

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 introduced by Kempton et al. (2018) and scaled the size of each planet accordingly. This metric considers the planet's equilibrium temperature Teq, radius Rp, mass Mp, host star radius Rs, and host star magnitude in the J band mJ. 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 mJ > 9, while HD 88986 has a magnitude of mJ = 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 88986b’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 , 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.

thumbnail Fig. 12

Position of HD 88986 b and its host star among known transiting exoplanetary systems, with precise planet mass and radius measurements (σ-M/M = 25% and σ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 (Rp < 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 (, Méndez & Rivera-Valentín 2017). These figures highlight the unique position of the HD 88986 system.

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 Internacionalizacã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).

Appendix A Atmospheric dispersion correction

Atmospheric dispersion can introduce a slope in the spectral continuum, leading to a shift in the mean RV values of the observed targets (Pepe & Lovis 2008; Wehbe et al. 2020). In order to achieve a higher RV precision necessary for detecting low-mass planets, it is imperative to consider the impact of atmospheric dispersion. To address this effect in SOPHIE, we implemented a correction method based on the HARPS8 and ESPRESSO DRS (Modigliani et al. 2019). This correction involves scaling the target spectrum to match its flux distribution with that of a template. Specifically, we multiplied the target spectrum by the flux ratio between the target spectrum and the template. The template is constructed from a high S/N spectrum of a standard star with the same spectral type, acquired at low air mass conditions. Through the application of this method, the flux distributions of star spectra always have the same distribution, thereby minimizing the influence of atmospheric conditions on the computation of the CCF.

The method described herein yielded an enhancement in the precision of SOPHIE RV measurements by 8 cm s−1 when applied to the high-precision RV measurements from the SOPHIE SP1 star catalog (which comprises 96 stars that have more than 10 observations). This improvement equates to 7% of the mean error bar of 1.2 m s−1 associated with these stars. With the application of this method to the full width at half maximum (FWHM) of the same set of stars, we observed a substantial precision enhancement of 15 m s−1 on the FWHM. Notably, this enhancement is nearly five times greater than the average error bars of 2.8 m s−1 for the FWHM of these stars. This correction has been incorporated into the SOPHIE DRS and will be used in future planet detections conducted by the SOPHIE instrument. Moreover, the correction method is applied to all SOPHIE RV constant stars (4 super constant stars in addition to ~ 26 other stars with low RV dispersion), resulting in the creation of a more robust RV master constant time series (see Appendix C).

Appendix B Moon contaminated spectra

To effectively detect and characterize planets through the RV method, it is imperative to meticulously manage systematic noise sources and eliminate any outliers. One such source of outliers in RV data sets is the spectra contaminated by Moonlight (e.g., Hébrard et al. 2008; Santerne et al. 2011). The Moon’s reflected light can cause spectral contamination, leading to potential masking of the planet’s signals or presenting systematic errors in the properties of the detected planets. Therefore, it is crucial to carefully identify and exclude any contamination by the lunar light from our data sets.

To achieve this, we developed a recipe that can identify Moon-polluted spectra for star observations made through the simultaneous calibration lamp, where no sky observation is available. This is achieved by taking into account the phase and position of the Moon at the time of observation, using two empirical criteria.

The first criterion considers whether moonlight contributes significantly to the target spectrum. This can be assumed when either:

  • 1.

    the Moon phase is more than 68% at the time of observation and the sky-level is above the mean of the sky-level of all observations,

  • 2.

    or when the separation between the target and Moon is less than 30°. We note that the sky-level is a criterion for estimating sky background light (see Bijaoui (1980); Ji et al. (2018)). In SOPHIE, it is calculated using SOPHIE DRS and is available in each FITS file spectrum header.

The second criterion we utilized in our study takes into account the proximity of the targets’ RV to the Moon’s RV. The Barycentric Earth radial velocity (BERV) at the time of observation and in the direction of the target is within approximately 1 kms−1 of Moon RVs (Díaz et al. 2012). Therefore, it is reasonable to consider it as a Moon RV. If the BERV is close to the target radial velocity, with |RVtargetBERV| < 2*FWHM, then a spectrum can be considered moonlight polluted.

To test these empirical criteria, we used observations of three stars with simultaneously recorded sky observations where recorded sky spectra were available. Over a total of 59 spectra, 13 data were Moon contaminated. Our criteria allowed us to successfully detect 8 of these contaminated spectra but also flagged three uncontaminated data. So we conclude our criteria are conservative.

Overall, our recipe provides a useful tool for identifying suspected Moon-polluted spectra in RV observations made through simultaneous calibration lamps. This can help improve data quality and ensure accurate scientific results. In the case of HD 88986, following these criteria, we identified 31 spectra that were suspected to be contaminated by moonlight. These data conservatively were discarded from our analysis. We note that including or excluding these data had no significant effects on our final results within the uncertainties, showing our criteria indeed are conservative.

Appendix C Update on constructing RV master constant timeseries

To monitor long-term instrumental variations, we have been conducting nightly observations of a few constant stars using SOPHIE since 2012. These constant stars include four super constant stars and approximately 25 additional stars with low RV dispersion (σ < 3.5 m s−1). Then, following the method outlined in Courcol et al. (2015), we create a master constant time series and subtract it from the RVs of each star. However, our latest analysis suggests that one of our constant stars, HD185144, exhibits activity over an extended period (see Fig. C.1). This activity could potentially affect our master constant time series and consequently impact the RVs of other stars.

We note that this star’s activity cycle was already known thanks to HIRES data (Isaacson & Fischer 2010). SOPHIE later confirmed this when we built the first master constant correction in 2015 (Courcol et al. 2015). At that time, we estimated the semi-amplitude of the signal to be less than 1.5 m s−1, which was negligible given the spectrograph’s precision. With more data points in 2023, the effect is estimated to have a semi-amplitude of 2.7 ± 0.1 m s−1. Given SOPHIE’s improved RV precision, it has become necessary to correct the impact of HD185144’s stellar activity on its RVs. Since this star is one of the most frequently observed constant stars by SOPHIE with more than 1100 data points, removing the RVs of this star from our master constant time series is not an option.

thumbnail Fig. C.1

log RHK (top) and RVs (bottom) of constant star HD185144. The best fit of a third-order polynomial model (black line bottom) is overplotted to the stellar long-period activity.

To correct this stellar activity, we utilized a four-step approach. Firstly, we corrected the HD185144 RVs for SOPHIE instrumental variations. To accomplish this, we subtracted the master constant time series, derived from SOPHIE constant stars excluding HD185144, from the RV of HD185144. Next, we fit a third-order polynomial model on the HD185144 RV time series to determine its activity phase (see Fig. C.1 bottom). Subsequently, we subtracted the same polynomial model from the raw RVs of HD185144 (prior to master constant correction). Finally, we used these corrected RVs to construct the final master constant time series along with other stars.

Appendix D RVs

Table D.1

SOPHIE RVs for HD 88986. We note that a value of 999.0 indicates invalid data on the corresponding date.

Table D.2

HIRES RVs for HD88986

Table D.3

ELODIE RVs for HD 88986

Appendix E Background contamination of the calibration lamp from the SOPHIE spectra

The combined analysis of RVs and activity indicators plays a crucial role in determining the origin of a signal. To do so, having accurate activity indicators is essential for effectively interpreting the data. In CCD spectrograph images, the recorded light from fibers A and B in a spectral order are located close to each other (e.g., within approximately 17 pixels for SOPHIE). This proximity introduces a small but non-negligible amount of light diffusion, primarily caused by the calibration lamp’s light from fiber B to fiber A (Lovis et al. 2011). Consequently, before deriving the activity indexes such as log RHK and Hα, it becomes imperative to subtract the background light originating from the diffuse light emitted by either the Th-Ar or FP calibration lamp present in the star spectrum.

To address this contamination issue, the SOPHIE DRS employs various methods, depending on the calibration lamps used. In the case of the Th-Ar lamp, a background is estimated from the flux of fiber B in the same spectral order by fitting a polynomial function on local minima (Boisse et al. 2010). On the other hand, for the FP lamp, which has been installed since semester 2017B, the background is directly measured using a Dark-FP frame, that is, illumination of fiber B with the FP calibration lamp while keeping fiber A completely dark (Hobson 2019; Lovis et al. 2011). However, as more years of observations were conducted, it became evident that there was a noticeable discrepancy in the data obtained from the two calibration lamps. This discrepancy can be attributed to the utilization of different background correction methods and likely over-estimation of background contamination in the method used in Th-Ar data. In order to rectify this issue, we employed a direct measuring method (Hobson 2019; Lovis et al. 2011) for both calibration lamps. By implementing this approach, we successfully corrected this discrepancy and significantly improved the consistency between data sets obtained using different calibration lamps. This method has been implemented into the SOPHIE DRS and will be utilized in forthcoming SOPHIE planet publications. We applied this method to the SOPHIE spectra prior to deriving the Hα and log RHK (or S-index) activity indices for HD 88986. We note that we did not use the Hα activity indexes in our final analysis due to its high contamination by the telluric.

Appendix F Priors on HD88986 b for RV-only model

Table F.1

Priors and description of parameters used within juliet to model RVs of HD88986 in Sect. 5.2.

Appendix G False positive tests

thumbnail Fig. G.1

Centroid analysis for HD88986. Difference images (panels a) for the potential transits in sector 21 (left) and sector 48 (right), along with the mean out-of-transit image (panel b) and the mean in-transit image (panel c). The difference images are obtained by subtracting the mean in-transit image from the mean out-of-transit image and ideally appear as an isolated stellar image of the host star. We note that in sector 48, the host star is located just beside a bad column. While interpreting the difference images from saturated stars such as HD88986 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. Thus, the transit features are likely associated with the host star in both sectors.

thumbnail Fig. G.2

TESS light curves reproduced using PLD and CBV approaches (left panels) and zoomed in the potential single transits (right panels) for sector 21 (top) and sector 48 (bottom). The SAP and PDC-SAP data are also plotted to provide a reference for comparison. The data are binned (black points) in 1 hour. The legend includes an overfitting score and the sgCDPP metric to facilitate an assessment of the different light curves. The brown vertical lines are the telescope momentum dumps.

thumbnail Fig. G.3

TESS data extracted from FFI for sector 21 (top) and sector 48 (bottom). The data are zoomed in on the potential single transits in the right panels.

thumbnail Fig. G.4

Raw PDC-SAP (left) and detrended FFI (right) phase-folded data from sector 21 and sector 48 with a period of 147.4 d. The PDC-SAP data were binned in 30 minutes to ensure compatibility with the FFI data.

Appendix H Priors on joint modeling of RV and photometric data

Table H.1

Priors for the joint modeling of RV and photometric data with juliet in Sect. 7

Appendix I Priors on joint modeling of RV and Hipparcos/Gaia astrometric data

Table I.1

Priors for the joint modeling of RV and Hipparcos/Gaia astrometric data.

Appendix J Corner plot of the joint modeling SOPHIE+ RVs and TESS sector 21 photometry

thumbnail Fig. J.1

Nested samples distribution of fitted parameters of HD88986 b on joint modeling of SOPHIE+ RVs and TESS sector 21 light curve. The 1, 2, and 3σ confidence levels of the posterior samples are shown by the contours.

Appendix K Corner plot of the joint modeling of combined HD88986 RVs, Hipparcos and Gaia astrometry

thumbnail Fig. K.1

MCMC samples distribution of fitted parameters of HD88986 outer massive companion on joint modeling of HD88986 RVs and Hipparcos/Gaia astrometry. V0, V1, V2, V3, and V4 correspond to the ELODIE, SOPHIE, SOPHIE+, HIRES, and HIRES+ RV dataset, respectively.

Appendix L APT photometric data

Table L.1

Summary of APT photometric observation for HD88986

References

  1. Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 252 [Google Scholar]
  2. Baluev, R. V. 2008, MNRAS, 385, 1279 [Google Scholar]
  3. Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron., 51, 109 [Google Scholar]
  5. Berger, T. A., Huber, D., Gaidos, E., & van Saders, J. L. 2018, ApJ, 866, 99 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bijaoui, A. 1980, A&A, 84, 81 [NASA ADS] [Google Scholar]
  7. Boisse, I., Eggenberger, A., Santos, N., et al. 2010, A&A, 523, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Bonfanti, A., Delrez, L., Hooton, M. J., et al. 2021, A&A, 646, A157 [EDP Sciences] [Google Scholar]
  9. Borgniet, S., Meunier, N., & Lagrange, A. M. 2015, A&A, 581, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bouchy, F., Hébrard, G., Udry, S., et al. 2009, A&A, 505, 853 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bouchy, F., Díaz, R., Hébrard, G., et al. 2013, A&A, 549, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Brandt, T. D. 2021, ApJS, 254, 42 [NASA ADS] [CrossRef] [Google Scholar]
  13. Brandt, G. M., Michalik, D., Brandt, T. D., et al. 2021a, AJ, 162, 230 [NASA ADS] [CrossRef] [Google Scholar]
  14. Brandt, T. D., Dupuy, T. J., Li, Y., et al. 2021b, AJ, 162, 186 [NASA ADS] [CrossRef] [Google Scholar]
  15. Butler, R. P., Vogt, S. S., Laughlin, G., et al. 2017, AJ, 153, 208 [Google Scholar]
  16. Courcol, B., Bouchy, F., Pepe, F., et al. 2015, A&A, 581, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Cutri, R., Skrutskie, M., Van Dyk, S., et al. 2003, The IRSA 2MASS All-Sky Point Source Catalog [Google Scholar]
  18. Da Silva, J. G., Santos, N., Bonfils, X., et al. 2011, A&A, 534, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Delisle, J.-B., Ségransan, D., Buchschacher, N., & Alesina, F. 2016, A&A, 590, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nat. Astron., 5, 775 [NASA ADS] [CrossRef] [Google Scholar]
  21. Deming, D., Knutson, H., Kammer, J., et al. 2015, ApJ, 805, 132 [Google Scholar]
  22. Díaz, R., Santerne, A., Sahlmann, J., et al. 2012, A&A, 538, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Díaz, R. F., Ségransan, D., Udry, S., et al. 2016, A&A, 585, A134 [Google Scholar]
  24. Dobos, V., Charnoz, S., Pál, A., Roque-Bernard, A., & Szabó, G. M. 2021, PASP, 133, 094401 [NASA ADS] [CrossRef] [Google Scholar]
  25. Ehrenreich, D., Delrez, L., Akinsanmi, B., et al. 2023, A&A, 671, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262 [Google Scholar]
  27. Fausnaugh, M. M., Burke, C. J., Caldwell, D. A., et al. 2020, TESS Data Release Notes: Sector 21, DR29, Technical Report [Google Scholar]
  28. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  29. Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220 [Google Scholar]
  30. Fortney, J. J., Visscher, C., Marley, M. S., et al. 2020, AJ, 160, 288 [NASA ADS] [CrossRef] [Google Scholar]
  31. Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [Google Scholar]
  32. Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504 [Google Scholar]
  33. Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Gilliland, R. L., Chaplin, W. J., Dunham, E. W., et al. 2011, ApJS, 197, 6 [NASA ADS] [CrossRef] [Google Scholar]
  36. Hall, J. C., Lockwood, G., & Skiff, B. A. 2007, AJ, 133, 862 [NASA ADS] [CrossRef] [Google Scholar]
  37. Hara, N. C., & Delisle, J.-B. 2023, A&A, submitted [arXiv:2304.08489] [Google Scholar]
  38. Hara, N. C., Bouchy, F., Stalport, M., et al. 2020, A&A, 636, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Hara, N. C., Delisle, J.-B., Unger, N., & Dumusque, X. 2022a, A&A, 658, A177 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Hara, N. C., Unger, N., Delisle, J.-B., Díaz, R. F., & Ségransan, D. 2022b, A&A, 663, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Hawthorn, F., Bayliss, D., Wilson, T. G., et al. 2023, MNRAS, 520, 3649 [NASA ADS] [CrossRef] [Google Scholar]
  42. Haywood, R., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517 [NASA ADS] [CrossRef] [Google Scholar]
  43. Hébrard, G., Bouchy, F., Pont, F., et al. 2008, A&A, 488, 763 [Google Scholar]
  44. Heidari, N. 2022, PhD thesis, Université Côte d’Azur; Shahid Beheshti University (Tehran), Iran [Google Scholar]
  45. Heidari, N., Boisse, I., Orell-Mique, J., et al. 2022, A&A 658, A176 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Henry, G. W. 1999, PASP, 111, 845 [NASA ADS] [CrossRef] [Google Scholar]
  47. Henry, G. W., Fekel, F. C., & Williamson, M. H. 2022, AJ, 163, 180 [NASA ADS] [CrossRef] [Google Scholar]
  48. Hippke, M., & Heller, R. 2019, A&A, 623, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Hippke, M., David, T. J., Mulders, G. D., & Heller, R. 2019, AJ, 158, 143 [Google Scholar]
  50. Hobson, M. J. 2019, Ph.D. thesis, Aix-Marseille, France [Google Scholar]
  51. Høg, E. 2001, in Encyclopedia of Astronomy & Astrophysics (Boca Raton: CRC Press), 1 [Google Scholar]
  52. Hooton, M. J., Hoyer, S., Kitzmann, D., et al. 2022, A&A, 658, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Hoyer, S., Guterman, P., Demangeon, O., et al. 2020, A&A, 635, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Hoyer, S., Bonfanti, A., Leleu, A., et al. 2022, A&A, 668, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Huber, D., Zinn, J., Bojsen-Hansen, M., et al. 2017, ApJ, 844, 102 [Google Scholar]
  56. Isaacson, H., & Fischer, D. 2010, ApJ, 725, 875 [Google Scholar]
  57. Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, Proc. SPIE, 9913, 99133E [Google Scholar]
  58. Ji, I., Hasan, I., Schmidt, S. J., & Tyson, J. A. 2018, PASP, 130, 084504 [NASA ADS] [Google Scholar]
  59. Johnson, J. A., Aller, K. M., Howard, A. W., & Crepp, J. R. 2010, PASP, 122, 905 [Google Scholar]
  60. Kawashima, Y., Hu, R., & Ikoma, M. 2019, ApJ, 876, L5 [Google Scholar]
  61. Kempton, E. M.-R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401 [CrossRef] [Google Scholar]
  62. Kervella, P., Arenou, F., & Thevenin, F. 2022, A&A, 657, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Kiefer, F. 2019, A&A, 632, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Kiefer, F., Hébrard, G., Sahlmann, J., et al. 2019, A&A, 631, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Kiefer, F., Hébrard, G., Lecavelier des Étangs, A., et al. 2021, A&A, 645, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
  67. Kreidberg, L. 2015, PASP, 127, 1161 [Google Scholar]
  68. Kubyshkina, D. & Fossati, L. 2022, A&A, 668, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Lacedelli, G., Wilson, T. G., Malavolta, L., et al. 2022, MNRAS, 511, 4551 [NASA ADS] [CrossRef] [Google Scholar]
  70. Lee, E. J., & Chiang, E. 2016, ApJ, 817, 90 [Google Scholar]
  71. Leleu, A., Alibert, Y., Hara, N. C., et al. 2021, A&A, 649, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Lendl, M., Csizmadia, S., Deline, A., et al. 2020, A&A, 643, A94 [EDP Sciences] [Google Scholar]
  73. Li, G., & Batygin, K. 2014, ApJ, 790, 69 [NASA ADS] [CrossRef] [Google Scholar]
  74. Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2 [EDP Sciences] [Google Scholar]
  75. Lomb, N. R. 1976, Astrophys. Space Sci., 39, 447 [Google Scholar]
  76. Lovis, C., Dumusque, X., Santos, N., et al. 2011, arXiv e-prints [arXiv:1107.5325] [Google Scholar]
  77. Luger, R., Agol, E., Kruse, E., et al. 2016, AJ, 152, 100 [Google Scholar]
  78. Luger, R., Kruse, E., Foreman-Mackey, D., Agol, E., & Saunders, N. 2018, AJ, 156, 99 [Google Scholar]
  79. Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [Google Scholar]
  80. Maxted, P. F. L., Ehrenreich, D., Wilson, T. G., et al. 2022, MNRAS, 514, 77 [NASA ADS] [CrossRef] [Google Scholar]
  81. Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [Google Scholar]
  82. Méndez, A., & Rivera-Valentín, E. G. 2017, ApJ, 837, L1 [Google Scholar]
  83. Meunier, N. 2021, arXiv e-prints [arXiv:2104.06072] [Google Scholar]
  84. Modigliani, A., Sownsowska, D., & Lovis, C. 2019, ESPRESSO Pipeline User Manual [Google Scholar]
  85. Morris, R. L., Twicken, J. D., Smith, J. C., et al. 2020, Kepler Data Processing Handbook: Photometric Analysis, Kepler Science Document KSCI-19081-003, 6 [Google Scholar]
  86. Morris, B. M., Delrez, L., Brandeker, A., et al. 2021, A&A, 653, A173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Noyes, R., Hartmann, L., Baliunas, S., Duncan, D., & Vaughan, A. 1984, ApJ, 279, 763 [NASA ADS] [CrossRef] [Google Scholar]
  88. Osborn, H. P., Bonfanti, A., Gandolfi, D., et al. 2022, A&A, 664, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Otegi, J., Bouchy, F., & Helled, R. 2020, A&A, 634, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Paunzen, E. 2015, A&A, 580, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Pepe, F., & Lovis, C. 2008, Phys. Scr., 2008, 014007 [NASA ADS] [CrossRef] [Google Scholar]
  92. Pepe, F., Mayor, M., Galland, F., et al. 2002, A&A, 388, 632 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49 [Google Scholar]
  94. Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proc. Natl. Acad. Sci., 110, 19273 [Google Scholar]
  95. Philipot, F., Lagrange, A. M., Rubini, P., Kiefer, F., & Chomez, A. 2023, A&A, 670, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Queloz, D., Henry, G. W., Sivan, J.-P., et al. 2001, A&A, 379, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Radick, R. R., Lockwood, G. W., Henry, G. W., Hall, J. C., & Pevtsov, A. A. 2018, ApJ, 855, 75 [Google Scholar]
  98. Reynolds, R. T., Squyres, S. W., Colburn, D. S., & McKay, C. P. 1983, Icarus, 56, 246 [NASA ADS] [CrossRef] [Google Scholar]
  99. Sandford, E., Espinoza, N., Brahm, R., & Jordán, A. 2019, MNRAS, 489, 3149 [NASA ADS] [CrossRef] [Google Scholar]
  100. Santerne, A., Bonomo, A., Hébrard, G., et al. 2011, A&A, 536, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Santos, N., Sousa, S., Mortier, A., et al. 2013, A&A, 556, A150 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
  103. Schlichting, H. E. 2014, ApJ, 795, L15 [Google Scholar]
  104. Serrano, L. M., Gandolfi, D., Hoyer, S., et al. 2022, A&A, 667, A1 [Google Scholar]
  105. Silburt, A., Gaidos, E., & Wu, Y. 2015, ApJ, 799, 180 [Google Scholar]
  106. Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000 [Google Scholar]
  107. Sousa, S., Adibekyan, V., Delgado-Mena, E., et al. 2018, A&A, 620, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Stassun, K. G., & Torres, G. 2016, AJ, 152, 180 [Google Scholar]
  109. Stassun, K. G., & Torres, G. 2018, ApJ, 862, 61 [Google Scholar]
  110. Stassun, K. G., & Torres, G. 2021, ApJ, 907, L33 [NASA ADS] [CrossRef] [Google Scholar]
  111. Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136 [Google Scholar]
  112. Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985 [Google Scholar]
  113. Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100 [Google Scholar]
  114. Tabernero, H., Allende Prieto, C., Zapatero Osorio, M. R., et al. 2020, MNRAS, 498, 4222 [NASA ADS] [CrossRef] [Google Scholar]
  115. Tal-Or, L., Trifonov, T., Zucker, S., Mazeh, T., &Zechmeister, M. 2019, MNRAS, 484, L8 [NASA ADS] [CrossRef] [Google Scholar]
  116. Tayar, J., Claytor, Z. R., Huber, D., & van Saders, J. 2022, ApJ, 927, 31 [NASA ADS] [CrossRef] [Google Scholar]
  117. Thompson, G. I., Nandy, K., Jamar, C., et al. 1978, Catalogue of Stellar Ultraviolet Fluxes: A Compilation of Absolute Stellar Fluxes Measured by the Sky Survey Telescope (S2/68) aboard the ESRO Satellite TD-1 (London: The Science Research Council) [Google Scholar]
  118. Torres, G., Andersen, J., & Giménez, A. 2010, A&ARv, Rev, 18, 67 [NASA ADS] [CrossRef] [Google Scholar]
  119. Trotta, R. 2008, Contemp. Phys., 49, 71 [Google Scholar]
  120. Twicken, J. D., Clarke, B. D., Bryson, S. T., et al. 2010, Proc. SPIE, 7740, 774023 [NASA ADS] [CrossRef] [Google Scholar]
  121. Van Cleve, J. E., Howell, S. B., Smith, J. C., et al. 2016, PASP, 128, 075002 [CrossRef] [Google Scholar]
  122. van Leeuwen, F. 2007, A&A, 474, 653 [CrossRef] [EDP Sciences] [Google Scholar]
  123. Vanderburg, A., Latham, D. W., Buchhave, L. A., et al. 2016, ApJS, 222, 14 [Google Scholar]
  124. Vanderburg, A., Huang, C. X., Rodriguez, J. E., et al. 2019, ApJ, 881, L19 [Google Scholar]
  125. VanderPlas, J. T. 2018, ApJS, 236, 16 [Google Scholar]
  126. Vaníček, P. 1971, Astrophys. Space Sci., 12, 10 [CrossRef] [Google Scholar]
  127. Vogt, S. S., & Penrod, G. D. 1988, in Instrumentation for Ground-Based Optical Astronomy: Present and Future The Ninth Santa Cruz Summer Workshop in Astronomy and Astrophysics, July 13–July 24, 1987, Lick Observatory (Berlin: Springer), 68 [Google Scholar]
  128. Wehbe, B., Cabral, A., Martins, J., et al. 2020, MNRAS, 491, 3515 [CrossRef] [Google Scholar]
  129. Wilson, T. G., Goffo, E., Alibert, Y., et al. 2022, MNRAS, 511, 1043 [NASA ADS] [CrossRef] [Google Scholar]
  130. Wilson, T. G., Simpson, A. M., & Collier Cameron, A. 2023, Nature, submitted [Google Scholar]
  131. Wright, E. L., Eisenhardt, P. R., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [Google Scholar]
  132. Yu, X., He, C., Zhang, X., et al. 2021, Nat. Astron., 5, 822 [NASA ADS] [CrossRef] [Google Scholar]
  133. Zeng, L., Sasselov, D. D., & Jacobsen, S. B. 2016, ApJ, 819, 127 [Google Scholar]

3

Available at https://dace.unige.ch

All Tables

Table 1

Stellar properties of HD 88986.

Table 2

Different tested models on the SOPHIE+ RV-only data along with model comparisons with juliet.

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).

Table 4

Possible solution for HD 88986 b’s period.

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.

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.

Table D.1

SOPHIE RVs for HD 88986. We note that a value of 999.0 indicates invalid data on the corresponding date.

Table D.2

HIRES RVs for HD88986

Table D.3

ELODIE RVs for HD 88986

Table F.1

Priors and description of parameters used within juliet to model RVs of HD88986 in Sect. 5.2.

Table H.1

Priors for the joint modeling of RV and photometric data with juliet in Sect. 7

Table I.1

Priors for the joint modeling of RV and Hipparcos/Gaia astrometric data.

Table L.1

Summary of APT photometric observation for HD88986

All Figures

thumbnail Fig. 1

Radial velocity measurements of HD 88986 from ELODIE, HIRES, HIRES+, SOPHIE and SOPHIE+.

In the text
thumbnail 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).

In the text
thumbnail 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).

In the text
thumbnail Fig. 4

1 periodogram of SOPHIE+ data. The identified periods are shown in red.

In the text
thumbnail 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 Tobs the total time span of observations, a) and b) are obtained with windows of size Tobs/3 = 1234.2 d and Tobs/9 = 411.42 d, respectively.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 12

Position of HD 88986 b and its host star among known transiting exoplanetary systems, with precise planet mass and radius measurements (σ-M/M = 25% and σ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 (Rp < 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 (, Méndez & Rivera-Valentín 2017). These figures highlight the unique position of the HD 88986 system.

In the text
thumbnail Fig. C.1

log RHK (top) and RVs (bottom) of constant star HD185144. The best fit of a third-order polynomial model (black line bottom) is overplotted to the stellar long-period activity.

In the text
thumbnail Fig. G.1

Centroid analysis for HD88986. Difference images (panels a) for the potential transits in sector 21 (left) and sector 48 (right), along with the mean out-of-transit image (panel b) and the mean in-transit image (panel c). The difference images are obtained by subtracting the mean in-transit image from the mean out-of-transit image and ideally appear as an isolated stellar image of the host star. We note that in sector 48, the host star is located just beside a bad column. While interpreting the difference images from saturated stars such as HD88986 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. Thus, the transit features are likely associated with the host star in both sectors.

In the text
thumbnail Fig. G.2

TESS light curves reproduced using PLD and CBV approaches (left panels) and zoomed in the potential single transits (right panels) for sector 21 (top) and sector 48 (bottom). The SAP and PDC-SAP data are also plotted to provide a reference for comparison. The data are binned (black points) in 1 hour. The legend includes an overfitting score and the sgCDPP metric to facilitate an assessment of the different light curves. The brown vertical lines are the telescope momentum dumps.

In the text
thumbnail Fig. G.3

TESS data extracted from FFI for sector 21 (top) and sector 48 (bottom). The data are zoomed in on the potential single transits in the right panels.

In the text
thumbnail Fig. G.4

Raw PDC-SAP (left) and detrended FFI (right) phase-folded data from sector 21 and sector 48 with a period of 147.4 d. The PDC-SAP data were binned in 30 minutes to ensure compatibility with the FFI data.

In the text
thumbnail Fig. J.1

Nested samples distribution of fitted parameters of HD88986 b on joint modeling of SOPHIE+ RVs and TESS sector 21 light curve. The 1, 2, and 3σ confidence levels of the posterior samples are shown by the contours.

In the text
thumbnail Fig. K.1

MCMC samples distribution of fitted parameters of HD88986 outer massive companion on joint modeling of HD88986 RVs and Hipparcos/Gaia astrometry. V0, V1, V2, V3, and V4 correspond to the ELODIE, SOPHIE, SOPHIE+, HIRES, and HIRES+ RV dataset, respectively.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.