© The Authors 2025

1 Introduction

In 1995, Mayor and Queloz discovered 51 Peg b, the first confirmed exoplanet around a main sequence star, using the radial velocity (RV) method (Mayor & Queloz 1995). This planetary detection represented a paradigm shift and launched the exoplanet field on the path toward its current sophisticated state. Specifically, their discovery proved giant planets could be detected with the technology at the time, resulting in a huge push to find more new worlds, such as 70 Vir b (Marcy & Butler 1996) and 47 UMa b (Butler & Marcy 1996). Subsequently, more advanced telescopes with much higher resolution and precision were developed, allowing the detection of less massive planets, for example the super-Earths GJ667C b, c, and d (Anglada-Escudé et al. 2012, 2013). Additionally, astronomers searched for planets in the habitable zone (HZ), where liquid water could exist on a planet’s surface. Santos et al. (2001) discovered HD 28185, which Jones et al. (2006) confirmed to be the first planet entirely within the HZ. Eventually, Quintana et al. (2014) discovered Kepler 186f, the first of many Earth-like planets orbiting its star in the habitable zone. As more and more discoveries accumulated, multi-planet solar systems were discovered, including Upsilon Andromodae b, c, and d (Lissauer 1999). Further, there have now been many detections of multi-planet systems in which two or more planets are locked in a mean-motion orbital resonance (MMR) (Popova & Shevchenko 2014; Beaugé et al. 2006; Ramos et al. 2017), such as TOI-216b and c, two warm large exoplanets in a 2:1 MMR (Dawson et al. 2019).

As the body of planetary systems with confirmed or candidate planets increases, it remains extremely important to revisit these systems to ensure that the most accurate planetary parameters are recorded. It is not uncommon for a planet detection to be contested at a later date, due to analyses with more data or differing methods producing inconsistent results (Jones & Jenkins 2014; Wittenmyer et al. 2013, 2019a,b; Jenkins et al. 2013; Jenkins & Tuomi 2014; Santos et al. 2014). For example, it is well known that stellar activity can create periodic RV signals (Oshagh et al. 2017; Costes et al. 2021). Magnetic activity, photometric variations, stellar rotation cycles, or long-term magnetic cycles can imitate planetary RV signals (Dumusque et al. 2011; Lovis et al. 2011). Therefore, examining available stellar indices, photometric data, and any correlations with the RV data remain important steps in verifying true planetary RV signals (see, e.g., Díaz et al. 2018). On the other hand, re-examinations of previously studied star systems frequently lead to new planetary detections (Mayor et al. 2009; Vogt et al. 2010; Tuomi et al. 2013; Orell-Miquel et al. 2023).

The scientific community places a special emphasis on the search for exoplanets in nearby solar systems, a case in point being the discovery of Proxima b, a terrestrial world orbiting within the HZ of our nearest neighbor star, Proxima Centauri (Anglada-Escudé et al. 2016); these exoplanets offer immense hope to advance our understanding of planet formation and evolution. HD 208487 serves as a particularly interesting system to re-examine since it is one of the nearest (~45pc) and brightest (V = 7.47) stars in the southern sky, and it is already known to be orbited by one planet. Additionally, a second longer-period signal was detected (Gregory 2005, 2007), though later disregarded as likely activity-induced (GoZdziewski & Migaszewski 2006; Wright et al. 2007). We now have 184 new RV observations spanning more than 15 years since Tinney et al. (2005) published this original discovery of HD208487b.

Of the nearly 900 multi-planet systems¹, only a handful exhibit similarities to our own Solar System. Weiss et al. (2018) analysed the architectures of the Kepler sample of multi-planet systems, revealing a phenomenon denoted as peas in a pod. Within a multi-planet system, each planet’s size is strongly correlated with that of its neighboring planets. In systems with three or more transiting planets, an additional pattern emerges: a systematic regularity in the spacing of planets. Frequently, these systems are densely populated with planets near their stability thresholds. Extensive research on the dynamic evolution and stability indicates that gravitational scattering plays a crucial role in shaping the final architecture of such multi-planet systems (e.g., Chatterjee et al. 2008; Johansen et al. 2012; Izidoro et al. 2017; Wu & He 2023; Ghosh & Chatterjee 2024). This process triggers perturbations in planetary orbits, giving rise to close encounters between planets. Ultimately, this can lead to either collisions (whether between two planets or between a planet and its star) or the ejection of planets from the system. The former is more prevalent at smaller semimajor axes (see, e.g., Petrovich et al. 2014). This observed correlation between planet size and spacing underscores the pivotal role of dynamics in shaping the final sizes and orbital separations of planets within these systems.

In this paper, we discuss how our new analysis has led to the statistical confirmation and discovery of two new planets in this system, HD 208487c and HD 208487d. The rest of the paper is organized as follows. In Sect. 2 we discuss the RV observations in detail, including a breakdown of the measurements by telescope. In Sect. 3, we present updated stellar parameters for HD 208487 using the SPECIES and ARIADNE codes. We begin the analysis of the RV data set by observing the GLS periodograms in Sect. 4, while Sect. 5 details the detection of the three planets with the more rigorous Bayesian analysis. Next, we highlight the search for correlations between stellar activity indices and the RVs in Sect. 6. In Sect. 7 we describe the analysis of ASAS photometric data and the fitting of stellar cycles. Finally, we discuss the results and their implications in Sect. 8, followed by our conclusions in Sect. 9.

2 RV observations

We used a total of 219 high-precision Doppler spectroscopic RV measurements of HD 208487 spanning from JD 2451034.18, August 8, 1998, to JD 2459508.64, October 21, 2021 for a total baseline of 8474 days (23.2 years), as seen in Figure 1. The data comes from the University College London Échelle Spectrograph (UCLES) instrument installed on the 3.9m Anglo-Australian Telescope at the Siding Spring Observatory (SSO), the High Accuracy Radial Velocity Planet Searcher (HARPS) installed on the 3.6 m ESO Telescope (before and after the fibre upgrade in 2015) at La Silla Observatory, and the Carnegie Planet Finder Spectograph (PFS) mounted on the 6.5 m Magellan II Telescope at Las Campanas Observatory.

2.1 3.9 m University College London Échelle Spectrograph (UCLES)

A total of 48 observations of HD 208487 were collected using UCLES at the 3.9 m Anglo-Australian Telescope (AAT) between JD 2451034.18 (August 8, 1998) to JD 2456969.96 (November 8, 2014) as part of the Anglo-Australian Planet Search (AAPS). UCLES is a cross-dispersed echelle spectrograph, which uses an I₂ absorption cell to derive RVs at the ~3 m/s precision level (Butler et al. 2001). In addition to bias subtraction and flatfielding, the I₂ method utilizes a point spread function (PSF) to properly wavelength calibrate the spectra. The known I₂ spectral features are used to calibrate the wavelengths of the stellar spectra and to model the PSF. Next, the stellar spectrum is deconvolved from the densely packed I₂ features to successfully extract the correct Doppler shifts found in the observed stellar spectra, which enables the calculation of long-term stable, precise RVs (Butler et al. 1996). UCLES offers a resolution of 40 000120 000 depending on the slit width and CCD binning (Diego et al. 1990; Tinney et al. 2005). The precision of UCLES data is between 2-3 m s⁻¹ (Butler et al. 2006). These data were obtained from the NASA Exoplanet Archive² from (Butler et al. 2006) and from newly measured RV data calculated for this work (Butler private communication). These data will henceforth be referred to as UCLES.

Fig. 1 All 219 mean-subtracted radial velocity measurements of HD 208487, including 48 UCLES (gray), 20 HARPS-TERRA (red), 106 HARPS-TERRA2 (olive), and 45 PFS (teal) observations.

2.2 3.6 m High Accuracy Radial Velocity Planet Searcher (HARPS)

HARPS is a stabilized, high-precision fiber-fed echelle spectrograph, contained in a vacuum cell (variations maintained ≤0.01 mbar) surrounded by a thermally stablized enclosure in order to avoid temperature fluctuations (≤0.01K) affecting the RV measurements (see Lovis et al. 2006). This instrument utilizes simultaneous observations through two fibers, one typically used to image the target and the other a Th-Ar reference spectrum to maintain high wavelength stability and monitoring over the long term. With a fibre aperture on the sky of 1″, HARPS maintains a resolving power of around 115 000. HARPS achieves an overall radial velocity accuracy of ~1 m/s (Mayor et al. 2003; Figueira et al. 2010).

2.2.1 TERRA

Using publicly available data from HARPS obtained from the ESO HARPS archive (PIs: Mayor, M; Kuerster, M; Trifonov, T), we found 20 observations between JD 2453265.69 (September 17, 2004) to JD 2454258.94 (June 7, 2007). These observations were all processed with the built-in HARPS Data Reduction Software (DRS), which performs the standard bias, flat fielding, and wavelength calibration of the high-resolution spectra. The DRS also yields the extracted spectra, the cross-correlation functions (CCF), a measurement of the full-width at half maximum (FWHM), and the bisector inverse slope (BIS) (Mayor et al. 2003; Díaz et al. 2018). After the data passes through the DRS, radial velocity measurements are then re-computed using the HARPS-TERRA (Template-Enhanced Radial velocity Reanalysis Application) package (Anglada-Escude & Butler 2012). Using all available orders of spectra of a star, HARPS-TERRA shifts them to zero, then co-adds the spectra, yielding one high S/N spectrum of the target. Next, the individual spectra are crosscorrelated order by order against the high S/N stellar template to re-calculate the RVs. We refer to these data as TERRA hereafter.

2.2.2 TERRA2

On June 3, 2015, a new fiber link was installed in HARPS in order to increase the stability of the instrument and the radial velocity precision (Lo Curto et al. 2015). As a result, we treat the observations from HARPS post-upgrade as a distinct instrument and refer to these radial velocity measurements as TERRA2, since they were processed in a similar manner to the preupgrade HARPS data. After the upgrade, we have a total of 106 observations spanning from JD 2457559.81 (June 20, 2016) to JD 2458755.48 (September 28, 2019).

2.3 6.5 m Carnegie Planet Finder Spectrograph (PFS)

PFS is a high-resolution optical echelle spectrograph, which is outfitted with an I₂ cell. This Iodine absorption cell aids with wavelength calibrations to ensure precise radial velocity measurements are derived from the Doppler shift of certain spectral lines (Butler et al. 1996). PFS is temperature controlled at 27 ± 0.01 °C and operated at ambient atmospheric pressure. The data taken prior to 2018 has a resolution of ~80000 (0.5” slit). In early 2018, the old 4K x 4K CCD with 15 micron pixels was replaced with a modernized 10K × 10K chip with 9 micron pixels. Since the upgrade³, observations have been taken with a 0.3” slit, yielding a resolution of ~ 130 000 with a radial velocity uncertainty on the order of ~1 m/s (Crane et al. 2006, 2008; Crane et al. 2010), for targets brighter than V=10. In this study, we have 45 observations from JD 2455847.60 (October 13, 2011) until JD 2459508.64 (October 21, 2021).

3 Stellar properties

We derived the host star properties using two independent methods based on equivalent width measurements implemented in SPECIES⁴ (Soto & Jenkins 2018) and spectral energy distribution (SED) fitting with ARIADNE⁵ (Vines & Jenkins 2022).

3.1 SPECIES - Spectroscopic Parameters and atmosphEric ChemIstriEs of Stars

SPECIES computes atmospheric parameters (T_eff, [Fe/H], log g, ξ_t) as well as 11 elemental abundances from high-resolution spectra (Soto & Jenkins 2018; Soto et al. 2021). Initially, the radial velocity is measured by cross-correlation with a G2 binary mask, and the spectrum is then shifted to the laboratory reference frame, where the Fe I and Fe II lines are identified and equivalent widths measured by fitting Gaussian-like profiles to the absorption lines. A grid of models taken from interpolating ATLAS9 (Castelli & Kurucz 2003) model atmospheres, as well as the equivalent widths, are provided to MOOG (Sneden 1973) in order to solve the radiative transfer equation (RTE) and estimate the correlation between Fe line abundances as a function of excitation potential and equivalent widths under the Local Thermodynamic Equilibrium (LTE) assumption. The above process is iterative, and is executed with a new set of parameters (T_eff, [Fe/H], log g, ξ_t) until either the number of maximum iterations is reached or when no correlation is detected above an established threshold e, where convergence is assumed for the corresponding stellar parameters. SPECIES mass, radius and age are estimated using the Isochrones package (Morton 2015) by interpolating through a grid of MIST (Dotter 2016) evolutionary tracks using T_eff, [Fe/H], log g as well as astrometric and photometric information. Finally, the vsin i and macro turbulent velocities are obtained from temperature calibrators and fitting Fe absorption lines of the observed HD 208487 spectra with synthetic line profiles.

HD 208487 stellar parameters were derived from the 10 highest signal-to-noise (SN) spectra (SN~200) fetched from the ESO database⁸, where SPECIES outputs from every run being within the 1σ statistical agreement for the stellar parameters were recorded. The typical median and uncertainties for the main stellar parameters from SPECIES were T_eff = 6163 ± 50 K, [Fe/H] = 0.09 ± 0.05 dex, log g = 4.39 ± 0.07, vsini = 5.63 ± 0.28 km s⁻¹, Age = ${1.22}_{-0.55}^{0.61}$ Gyr, M_⋆ = 1.17 ± 0.02 M_⊙, and >R_⋆ = 1.15 ± 0.01 R_⊙.

3.2 ARIADNE - spectrAl eneRgy distribution bAyesian moDel averagiNg fittEr

Adjusting the SED to a stellar atmosphere model is a common method to derive stellar properties. However, the physical assumptions of the underlying atmosphere model may contribute to bias the stellar parameters. Here we used a Bayesian Model Averaging approach implemented in ARIADNE (Vines & Jenkins 2022), which fits archival photometry to several stellar atmosphere models, in order to mitigate the bias from a single atmosphere model. The stellar parameters were the results of an averaged posterior distributions from the Phoenix V2 (Husser et al. 2013), BT-Settl (Hauschildt et al. 1999; Allard et al. 2012; Kurucz 1993; Castelli & Kurucz 2003) SED models weighted by their respective Bayesian evidence estimates.

To derive the parameters from Table 1, SPECIES T_eff, log g, [Fe/H] were used as priors along with archival photometric data, which was automatically fetched by ARIADNE for the SED fitting. ARIADNE outputs are used to derive the stellar age, mass and the equal evolutionary points from the isochrone package.

4 Periodogram analysis

To begin the analysis of the RV data, we utilized the EXOSTRIKER package, which is a “transit and radial velocity interactive fitting tool for orbital analysis and N-body simulations” built as an accessible GUI (Trifonov 2019). We searched for embedded periodic signals in the data using the Generalized Lomb-Scargle (GLS) periodogram with a minimum period of 1 day and a maximum period of 10 000 days (Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009).

4.1 Full RV data set

After inserting the combined RV data from UCLES, HARPS, and PFS into EXOSTRIKER, we mean-subtracted the RVs and then fit a linear trend to the data. The top panel in Figure 2 shows the results of the original GLS periodogram, after the initial mean-subtraction and linear trend fit (with the Window function in green). The next panel displays the GLS periodogram of the residuals after fitting the first Keplerian (k1) signal, where the result shows the remaining embedded signals. Similarly, the third panel illustrates the GLS periodogram on the residuals of the combined fit of the k1 and k2 signals. Additionally, in that panel, we see flanking peaks that surround the ~930 d signal. Although they appear to be statistically significant peaks, they are artifacts of the main signal due the discrete sampling of the data (window function) beating with the main signal. The bottom panel reveals that there are no remaining significant peaks after fully fitting a three-planet model.

After fitting all the signals that surpass the 0.1% FAP threshold, we find periods at P₁ = 129.3 d, P₂ = 493.2 d, and P₃ = 939.7 d. This preliminary analysis helps inform what signals may appear during the more rigorous Bayesian analysis.

4.2 Split RV data set

As a first check to test whether stellar activity is adding noise to the RV data, or indeed causing the signals we are seeing, we repeated the periodogram analysis with split data sets, where the split is made after calculating the median S-index value in each data set (excluding UCLES data, which do not have S-index values recorded). We then split the RVs in half: one data set with all RVs corresponding to S-index values less than the median and one data set with all RVs corresponding to S-index values greater than or equal to the median. We found that the low Sindex data find all three signals, while the high S-index data only find the inner gas giant. Thus, the low S-index data provide a significantly better fit to the data when compared to the high Sindex data, which would be expected if the RV signals are all independent of the stellar activity.

Fig. 2 Top panel: original GLS periodogram of the combined and mean-subtracted RV data (black). The periodogram of the window function (sampling) for the combined data (green) is also shown. Second panel: GLS periodogram after fitting and removing the first detected signal. Third panel: GLS periodogram after fitting and removing both the first and second signals. Bottom panel: GLS periodogram after fitting and removing all three previously detected signals. The blue dotted horizontal lines represents the 0.1% false alarm probability level estimated from 5000 bootstrap resamplings.

4.2.1 Low S-index

The low S-index data set includes 10 TERRA, 22 PFS, and 53 TERRA2 observations, for a total of 85 RVs. Duplicating the analysis from the full data set, the GLS periodogram is used to search for any periodicity in this data set of lower stellar activity. To evaluate the models that we fit, we check the Bayesian Information Criteria (BIC). In general, a ΔBIC ≥ 5 (~150x more probable) is considered statistically significant (Kass & Raftery 1995; Feng et al. 2016). After mean-subtracting the RVs and fitting a linear trend, the initial k0 model yields a BIC of 5940.85. The GLS periodogram clearly shows the first two signals above the 0.1% FAP threshold, with the third slightly below, as seen in Figure B.1. Fitting the first signal, with a period of P₁ = 129.3 d, yields a BIC of 620.87 and a χ²_red = 4.60. Next, we add the second signal of P₂ = 536.6 d to the model, which after fully fitting ends with a period of P₂ = 526.4 d, with a BIC of 452.98 and a χ²_red = 2.25. Initially, the P₃ signal falls below the 0.1% FAP threshold. However, fitting the third period found from the highest peak and refitting all of the signals ends up reducing the BIC to 444.66, where P₃ = 937.6 d. Adding the third signal reduced the BIC by 8.32 (~4100x improvement), enabling us to conclude that the third signal still exists in this data set. After fully fitting the three-planet model, we found that rms = 2.18 m/s, χ² = 137.94, and χ²_red = 2.03, showing a good fit to the data.

4.2.2 High S-index

The high S-index data set includes 10 TERRA, 23 PFS, and 53 TERRA2 observations, for a total of 86 RVs. The RVs corresponding to the data with higher stellar activity were examined with the GLS periodogram, finding that only the first of the previously detected peaks shows up at high statistical significance (Figure B.2). A BIC of 3719.79 was found after mean-subtracting and fitting a linear trend to the RVs. Once again, we attempted to manually fit the peaks from the periods found with the full data set. After adopting the first period of P₁ = 130.3 d as the k1 signal, the best fit is found at P₁ = 130.6 d with a BIC of 1020.59 and a χ²_red = 9.67, showing the signal is significant within the BIC test. Neither of the remaining two signals are detected from the high S-index data. Thus, we concluded that the stellar activity is actually obscuring the signals, rather than contributing to them.

5 Bayesian analysis

After examining the results from the periodogram analyses, we found significant evidence for two additional planetary signals in the data, along with the previously confirmed planet, HD 208487b. In order to corroborate the initial results, we followed similar procedures to those employed in Tuomi et al. (2014); Jenkins & Tuomi (2014); Díaz et al. (2020) to model the RVs. The Markov chain Monte Carlo (MCMC) code, Exoplanet Mcmc Parallel tEmpering Radial velOcity fitteR⁹ (Peña Rojas & Jenkins 2025), searches for signals in radial velocity time series, using adaptive parallel tempering MCMC, convergence tests, and Bayesian statistics, along with various noise models. Further, the Python code provides Bayesian posterior distribution models and the parameterization for the planetary signals to describe the best-fit orbital solution to the star system. An early version of EMPEROR has already been used in a number of works to detect small signals in RV data (Wittenmyer et al. 2017; Jenkins et al. 2019a,b; Barnes et al. 2020; Jenkins et al. 2020), where the updated version employs a custom-written Adaptive Parallel Tempering MCMC engine, similar to the EMCEE package (Foreman-Mackey et al. 2013), as used in Gregory (2005, 2007).

We successively ran the EMPEROR code with different properties to find the optimal setup. The statistical models included a function describing a k-Keplerian planet model, a linear trend term (to account for any existing accelerations within the system), a red-noise model (consisting of a moving average model with exponential weighted smoothing), and linear correlations with the stellar activity indices. The full mathematical description of the model with all parts included is shown below: $$\begin{array}{r}{y}_{i,ins}={\gamma }_{ins}+\dot{\gamma }{t}_i+f_k(t_i)+\sum _{n=1}^{q_{ins}}{c}_{n,ins}{\xi }_{n,i,ins}\\ +\sum _{l=1}^p{\varphi }_{ins,l}\exp \left\{\frac{t_{i-l}-t_i}{{\tau }_{ins}}\right\}{r}_{i-l,ins}+{ϵ}_{i,ins},\end{array}$$(1)

where y_i,ins corresponds to the observation at time t_i for each of the distinct instruments, γ_ins is the velocity offset for each of the instruments, γ˙ is a linear trend term (for any existing acceleration), and c_n,ins describes the linear correlations with q stellar activity indices ξ_n,i,ins for each instrument. r_i,ins denotes the residuals once the model has been subtracted from the measurement. The function f_k is a superposition of k-Keplerian signals denoted by $$f_k\left(t_i\right)=\sum _{m=1}^k{K}_m[\cos ({\omega }_m+{\nu }_m\left(t_i\right))+e_m\cos \left({\omega }_m\right)],$$(2)

where K_m is the velocity semi-amplitude, ω_m is the longitude of periastron, ν_m is the true anomaly and e_m is the eccentricity. Since ν_m is a function of the orbital period and the mean anomaly M_0,m at epoch, we can fully describe f_k by Km, ω_m, e_m, M_0,m, P_m; m ∈ {1,..., k}.

The white noise term is denoted by the additive random variable ε_i,ins, shown by $${ϵ}_{i,ins}\sim \mathcal{N}(0,{\sigma }_i^2+{\sigma }_{ins}^2),$$(3)

where σ_i is the uncertainty reported with the measurement y_i,ins and σ_ins is the excess white noise or jitter modelled after each instrument. This nuisance parameter, which is required to model the mechanisms that create the data, has a Normal prior bound (Table 2), in the model. The remaining terms define the rest of the noise model, as described in Díaz et al. (2018).

5.1 Posterior samplings and signal detection

The adaptive parallel tempering provides a ladder of replicas of the system, where each step has the likelihood at a different temperature. The hotter the chain, the flatter the posterior, allowing the sampler to jump from one local maxima to another, whilst the cold chains explore the phase space in detail. Swapping states for random members of different systems makes the cold chain fully characterize each of the maxima. The adaptive scheme utilizes varying temperatures to serve as a ‘good ladder.’ This ladder allows for the calculation of the Bayesian evidence via thermodynamic integration as well (more details in Vousden et al. (2016). The prior distributions for the orbital and instrumental parameters for our model are listed in Table 2.

5.2 Model selection

The final piece of the EMPEROR code deals with choosing the proper model. Ideally, the final model represents the most accurate description of the orbital and planetary parameters. As a result, EMPEROR performs model selection at each stage to ensure the current model actually produces an improvement over the previous model. In order for the new k + 1 model to be preferred over the existing k model, we require an improvement over a statistic of choice. In this scenario, for the model comparison, the log of the evidence (Z) was chosen, and a Bayes factor of at least 5 was required, which represents a solution that is 150× more probable than the previous model. The model evidence, or marginalized likelihood, corresponds to the likelihood of generating the data from the prior, given a correct model. In other words, it could be understood as the probability of the model itself, hence its use for model comparison (Robert 2007). The ratio of evidences for two competing models is called the Bayes factor and is used for model comparison.

The EMPEROR MCMC code employs a burn-in stage, in which the chains search for the high-probability zones, whilst adapting their temperatures to have an efficient swap rate. When the burn-in stage is over, the samples are discarded, and temperatures stop adapting. This ensures that the chain remains ergodic whilst satisfying detailed balance. This process lasts for 80% of the total iterations. In our analysis, the Markov chains comprised 12 temperatures, 1200 walkers, and 20 000 steps, combining for 288 million iterations per model (with a 230 million iteration burn-in). Then, the chain must be longer than 50 times the autocorrelation time, or it is extended. The models we ran, beyond the straight Keplerian models include:

1. a white-noise model (WNO, as shown in Eq. (3)).

2. correlated noise, modelled by an exponentially weighted moving average (MAO, see last term of Eq. (1)).

3. white noise with linear stellar activity correlations (CorrO, see the first sum of Eq. (1)).

4. correlated noise plus the linear stellar activity correlations (CorrMA).

Initially, the Corr models were run with every available stellar index; however, only the indices with correlations with the RVs that were inconsistent with zero (no matter how low the rank) were included in the final set of runs. As a result, the Corr models are based on the S-index from PFS and the BIS from TERRA1 and TERRA2.

Additionally, we ran each of these models with a fixed eccentricity of zero, a tightly constrained eccentricity described by a normal distribution centered at zero with a standard deviation of 0.1, and a constrained eccentricity described by a normal distribution centered at zero with a standard deviation of 0.3 (as explained in Jenkins et al. 2025), which enables EMPEROR to explore a wide variety of solutions based on the potential planetary eccentricities. In general, we found that fixing the eccentricity prior to zero would cause additional signals to appear that were harmonics of the real signals, primarily a consequence of the fact that the largest signal in the data corresponds to a massive planet on a non-circular orbit, and hence is poorly explained by a strict cosine-like function. All circular models that caused the 1st harmonic of this primary signal (HD 208487b) to appear as a secondary signal in the data performed statistically much worse, enabling us to securely disregard these harmonic models. On the other hand, most of the models with the constrained eccentricity prior or the tightly constrained eccentricity prior yield very similar results. The main difference between the constrained and tightly constrained models are the predicted eccentricities of the outer two planets. The periods and semi-amplitudes are similar, but the constrained models tend to terminate with higher eccentricities, which results in both worse final statistics and long-term unstable planetary dynamics. In the theoretical limit of infinite walkers and steps, we would expect all of the models to converge to the same parameter configuration, even though we extended our analyses as far as reasonably possible to 288 million iterations.

We found that the CorrO with the tightly constrained eccentricity model performed the best, yielding a dynamically stable configuration. An independent verification method was performed to validate our results, which was in agreement with our best-fit model (see Sect. 5.3). The final orbital and instrumental parameters for the best-fitting model (CorrO-tightly constrained eccentricity) are shown in Table 3. The best-fit Keplerians are shown in Figure 3 and the posterior distributions for the periods and semi-amplitudes are shown in Figure 4.

5.2.1 Instrument contributions to the model

Due to the varying number of RVs in each data set as well as the differences in quality of each instrument, it is worthwhile to determine each data set’s overall contribution to the final model based on each instrument’s RMS. Each data set’s contribution is approximated by dividing its own number of observations by the square of its RMS, ($\frac{N_i}{RM{S}_i^2}$) and then finding the weighted sum. As shown in Table 4, although AAT and PFS comprise around one-fifth of the overall number of RVs, they both contribute a disproportionately low amount to the overall best-fit model. In fact, the AAT data contributes seven times less value than expected based on the number of RVs. On the other side, the TERRA2 data contributes most to the model due to its high precision, and large data volume. However on the other hand, including the AAT data increases the duration of the baseline from about three to five full periods of the outer planet, which lowers the uncertainties of the orbital parameters, especially for the longer-period planets. Therefore, we find an interesting balance between modeling long-baseline-lower-precision RV data sets against short-baseline-high-precision measurements; even though some data sets can be an order of magnitude lower in precision, their overall observing baselines can overcome these shortcomings when targeting longer-period worlds.

5.2.2 Comparison of earlier data and later data

Checking the stability of the RV signals over time provides another important determination of the veracity of the signals. We split the complete RV data set into first half and second half data, separated by instrument. Thus, the first half data includes 68 RVs, 48 data points from AAT and 20 data points from TERRA. The second half data includes 151 data points, 45 RVs from PFS and 106 RVs from TERRA2. We executed the EMPEROR code on the split data sets, testing the various types of models previously discussed. In the first half data, the best-fit model finds three statistically significant signals, $P_1={129.39}_{-0.22}^{+0.15}$ d, $P_2={934.56}_{-42.29}^{+17.97}$ d, and $P_3={1294.65}_{-38.62}^{+1191.79}$ d (see Table 3).

The second half data, reveals similarities to the orbital parameters derived from the first half data set. The periods and semi-amplitudes of the fit to the second-half data are as follows: $P_1={129.34}_{-0.03}^{+0.01}$ d, $P_2={485.60}_{-6.68}^{+4.71}d$ d, and $P_3={979.74}_{-36.32}^{+33.79}$ d (see Table 3). The period and semi-amplitude of HD 208487b remain almost identical between the early and late data. The period and semi-amplitude of HD 208487c are also consistent within uncertainties. The main difference between the two split data sets is the appearance of HD 208487d as an intermediate-period signal rather than a long-period signal. Interestingly, the later data set, with a baseline of only 3661d, is in close agreement with the initial EMPEROR model, which performs worse statistically than the chosen model. However, the model based on the earlier data, with a longer baseline of 5936d, more closely resembles the chosen CorrO model, as evidenced by the similar periods for all three signals. We note that since this intermediate signal repeatedly appears, just with lower statistical significance, this may indicate the reality of a further Doppler signal that requires more high-precision RVs to confirm.

5.2.3 AAT-only run

Ensuring consistency between earlier research and the current analysis remains an essential aspect of follow-up studies. Initially, Tinney et al. (2005) detected HD 208487b utilizing a Periodogram analysis. Later, Gregory (2005, 2007) claimed the detection of HD 208487c with a more statistically robust MCMC method compared to Tinney’s Periodogram analysis. To maintain consistency, we re-analyzed the AAT data alone with EMPEROR. The results from these MCMC runs, summarized in Table 3, indicate very similar findings, adding confidence to the past research.

5.3 Independent verification of results

In order to increase the robustness and reliability of the results, we analyzed the data with independent sampling methods and codes. We applied the delayed-rejection adaptive Metropolis algorithm (Haario et al. 2006) to detect maxima in the posterior probability density and the simple adaptive-Metropolis algorithm (Haario et al. 2001) to estimate parameters and obtain solutions (see, e.g., Butler et al. 2017; Tuomi et al. 2018).

We also increased the set of available activity-indices by calculating the “differential velocities” of Feng et al. (2017) for the TERRA reduction of HARPS. This is to say that we divided the 72 HARPS orders into three equal subsets, corresponding to the red-most, middle, and blue-most parts of the corresponding spectra. These were then used independently to obtain three sets of radial velocities denoted as y_i,1, y_i,2, and y_i,3 for i = 1,...,N_HARPS. The differential velocities are thus ξ_i,1 = y_i,1 - y_i,2 and ξ_i,2 = y_i,2 - y_i,3 that can be used as additional activity indices. Because any and all Doppler signals of planets have been subtracted when producing these indices, the remaining variability simply corresponds to noise and activity-induced variations. In particular, wavelength-dependent variations that are present in these indices enable one to filter them out from the data by accounting for the respective correlations with these indices in accordance with Eq. (1). The models are thus equivalent except for the fact that the additional activity indices enable removing some additional activity-induced variability. Such variability indeed was found to be only weakly present in the HARPS data. However, since the results were consistent, the postulated planetary signals presented in the manuscript show no evidence for variations as a function of spectral wavelength.

The results from these analyses are consistent with those presented in Table 3. The obtained solution for the three-Keplerian model contains periods of 129.40±0.04, 917.0±6.4, and 1396±27 days, respectively. These estimates are consistent with the results obtained with the EMPEROR code, albeit slightly more uncertain. This excess uncertainty may be caused by the fact that the model contains additional activity indices and is thus more flexible. It is also clear that wavelength-dependent variations in the HARPS radial velocities, as estimated by the differential velocities, are not responsible for the observed Keplerian periodicities.

Fig. 3 Top to bottom, left to right: Phase-folded best-fit model for the first, second, third, and complete Keplerian model.

Fig. 4 Top, left to right: final posterior distribution of the period for the first, second, and third Keplerians from the best-fitting parallel-tempered MCMC run. Bottom, left to right: final posterior distribution of the semi-amplitude for the first, second, and third Keplerians from the best-fitting parallel-tempered MCMC run. The numbers at the top of each figure correspond to the mean, variance, skewness, kurtosis, median, and mode, respectively. The solid line represents a Gaussian curve with the same mean and variance.

Fig. 5 GLS periodogram of the available stellar activity indices. From top to bottom: PFS S-index, TERRA S-index, TERRA FWHM, TERRA BIS, TERRA2 S-index, TERRA2 FWHM, TERRA2 BIS, Combined Sindex, Combined BIS, Combined FWHM. The green solid vertical lines show the best-fit periods of the three planetary signals. The blue dotted horizontal lines represents the 0.1% false alarm probability (FAP) level estimated from 5000 bootstrap resamplings.

6 Stellar activity and correlations

6.1 GLS periodograms and correlations of stellar activity

Following a similar strategy as in Santos et al. (2014) and Díaz et al. (2018), we investigated how stellar indices might be affecting or contributing to the periodic signals found within the RVs. First, we checked the GLS periodograms of every available stellar index to see if any stellar activity peak is associated with a previously found periodic signal, or its harmonics, as shown in Figure 5. No signals are significant above the 0.1% FAP threshold for any index in TERRA or PFS, as well as the Bisector Inverse Slope (BIS) of TERRA2. The S-index in TERRA2 has significant peaks in the mid 20 days and also around 800 days. Also, the Full Width at Half Maximum (FWHM) from TERRA2 has a significant peak at ~22 d. However, none of these peaks approach the previously found periodic signals. The Combined S-index and Combined FWHM both had significant peaks at long periods around ~4000 d. Additionally, the RV vs stellar index correlations were plotted in Figure 6. Table 5 shows the calculated Pearson rank and p-value of the RV/index correlations. Generally, in order for a correlation to be statistically significant, the p-value must be p ≤ 0.05. We note that the TERRA S-index, TERRA2 S-index, and TERRA2 FWHM have p-values approaching significance; however, none have a rank corresponding to even a moderate correlation. A rank of |r| ≥ 0.5 is required to be considered a moderate or stronger correlation. Since none of the stellar index correlations have p ≤ 0.05 nor |r| ≥ 0.5, there are no significant nor strong/moderate correlations between the RVs and the stellar indices.

6.2 Fitting stellar activity with EMPEROR

To expand our search for periodic RV modulations induced by stellar activity, we incorporated Bayesian statistical analysis. As a novel approach, we executed the EMPEROR code on the available stellar activity indices (the authors are unaware of any previous procedures involving Bayesian statistical analysis on stellar activity indices). Utilizing the EMPEROR code enables a more rigorous analysis of the stellar activity signals compared with the standard GLS and correlation analysis. Moreover, the posterior distributions derived from the Markov chains afford a comprehensive assessment of certainty across the entire period space, facilitating a thorough evaluation of the distinctness of the signals.

We completed several EMPEROR runs for each stellar activity index. Each index was initially analyzed individually, before joining the data sets of similar indices from different instruments into one run to test the stability and significance of the signals over the full baseline of the activity data. Further, we tested the differences between how EMPEROR fits the stellar activity indices with the fixed zero eccentricity and the constrained (N ~ (0,0.3,0,1)) eccentricity prior. The results of the noteworthy EMPEROR run is summarized in Table 6.

The result needing further inspection arose from the highest-probability EMPEROR run (constrained eccentricity prior) of the combined FWHM (TERRA and TERRA2). The FWHM period found at $P_{FWHM}={860.7}_{-9.3}^{+8.4}$ d somewhat near to the best-fit period for HD 208487c at $P_{RV}={923.06}_{-2.76}^{+2.02}$ d. Propagating the errors yields a sigma separation of ~7.1σ away from HD 208487c. Thus, based on the current data, the combined FWHM signal is again statistically distinct from the long-period RV signal. We can therefore say to a high level of statistical confidence that none of the three significant signals we find in the RV data have periods that overlap with signal frequencies in any of the activity indices we analyzed.

7 Photometric analysis

7.1 ASAS photometry

In order to provide additional analysis of this star/planetary system, All-Sky Automated Survey (ASAS) photometric data were analysed (Pojmanski 1997). Searching through the ASAS archive, we found 608 V-band observations spanning just over nine years, between JD 2451870.55 (November 22, 2000) and JD 2455168.58 (December 3, 2009). Following the ASAS guidelines, only data marked with an “A” or “B” quality flagged are kept (henceforth all data), leaving 509 usable data points. However, when examining the photometric time series, we see evidence of instrumental issues in the first ~300 d of observations, see Figure 7 top. In order to clean the data, we calculated the mean and standard deviation of all data after JD 2452300. Then all data falling within three standard deviations of the calculated mean for the entire baseline were kept, as seen in Figure 7 middle. Finally, the GLS periodogram was calculated from the cleaned data, where we sampled the period space from 1 day to 10 000 days, performing 80 000 period samples.

Table 7 shows the top five peaks and their corresponding powers; however, all of them fall significantly below the 0.1% FAP line. It is interesting to note that Wright et al. (2007) posit that a ~28-29 d signal they found in earlier RV data could be related to a planet; however, given the relatively strong peak around that period in the ASAS data, we believe the planetary hypothesis is unlikely. Although the top two peaks are under the 0.1% FAP line, they fall in a period range that could be related to either the stellar rotation cycle or the lunar cycle. Additionally, using the same setup, we ran the GLS Periodogram for the entire photometric data set, also showing no significant periods.

Fig. 5 Correlations of the available stellar activity indices. From top to bottom, left to right: PFS S-index, TERRA S-index, TERRA FWHM, TERRA BIS, TERRA2 S-index, TERRA2 FWHM, TERRA2 BIS. The lines show the best-fit correlation between the activity indicator and the RVs.

7.2 Fitting ASAS stellar cycles

Extending our analysis of the ASAS photometry, we searched for periodic stellar cycles in the cleaned data. After noting that none of the GLS periodogram peaks align with the planetary signals, we fit the strongest periods, as shown in Figure 8. The segment of the periodogram with the highest power is continuously increasing, with no defined peak, even as it approaches a period of 10 000 days. Thus, we conclude that the signal is likely a baseline effect, rather than a real cycle. The next two peaks constitute very similar periods of 29.51 days and 28.73 days. Based on the strong peak in ASAS and a very similar period in the GLS periodogram of the S-index of TERRA2, we conclude that the spikes in this period range are very likely to be due to the stellar rotation period, possibly even being affected by differential rotation. With continued spectral monitoring of this star we can further test this hypothesis in the future.

Fig. 7 Left (right): ASAS V-band photometry (Hipparcos photometry). Top: A and B flagged V-band ASAS photometry (all Hipparcos data) for HD 208487. Middle: cleaned V-band ASAS (Hipparcos) photometry for HD 208487. Bottom: GLS periodogram of the photometry. The horizontal lines correspond to the 0.1% (dotted blue) FAP significance level computed via 5000 bootstrap iterations on the photometry time series. The vertical green lines represent the periods of the three planets.

7.3 Hipparcos photometry

To further explore the photometric properties of HD 208487, we investigated the publicly available HIPPARCOS data (Hipparcos 1997). These data were obtained from the NASA Exoplanet Archive¹⁰, which contained 90 data points, spanning three years from JD 2447881.99 (December 21, 1989) to JD 2448977.85 (December 21, 1992). Following the data quality flag information, 1 point was identified as low-quality and unusable data. After removing the poor data point, the GLS Periodogram was calculated, in which we sampled the period space from 1 day to 10 000 days, performing 80 000 period samples. There were no peaks coming close to statistical significance, demonstrating a lack of any periodic photometric signals in the Hipparcos data.

8 Discussion

With 184 new RV observations spanning 15+ years since the previous analyses (Tinney et al. 2005; Gregory 2005; Gozdziewski & Migaszewski 2006; Gregory 2007; Wright et al. 2007; Babu et al. 2010), this system is ripe to study and to re-evaluate the reality of the planetary system. Using the SPECIES and ARIADNE codes, we provide updated stellar parameters to HD 208487. In order to fully understand the RV data set and best characterize the system, we utilized GLS periodograms and EMPEROR parallel-tempering MCMC samplings, which revealed new periodic signals embedded in the RVs. An independent, analysis accounting for wavelength-dependent variations in the HARPS-TERRA velocities yielded consistent results providing confidence in the final parameters. The best-fitting model to the current RVs reveals a candidate three-planet system.

Our best-fit model, which confirms the existence of HD 208487c, alters some of the parameters. Although the period, ${923.06}_{-2.76}^{+2.02}$ d, falls within the uncertainty published by Gregory (2005, 2007), the semi-amplitude, ${6.18}_{-0.07}^{+0.09}$ m/s, and mass, 0.32 ± 0.01 m_j, are now slightly lower leading to a classification as a Saturn-mass planet. Further, our analysis reveals a third signal, relating to a candidate outer planet, HD 208487d, with a period of ${1380.13}_{-8.25}^{+19.20}$ d, a semi-amplitude of ${2.46}_{-0.09}^{+0.30}$ m/s, and a minimum mass of ${0.15}_{-0.02}^{+0.01}{m}_j$. The two outer signals are found very close to the 3:2 period ratio (P_d/P_c = 1.495), further increasing the reliability of these as being due to Doppler signals, since it is difficult to produce 3:2 frequency ratio signals through aliasing (e.g., Hara et al. 2020), particularly with a possible 2:1 signal also present in the data, and also transit observations have shown such period ratios to be fairly common (Lissauer et al. 2011; Fabrycky et al. 2014; Lissauer et al. 2024). These tests therefore leave us with a system that presents an interesting and unique system architecture, as there are currently no other known planetary systems that contain two planets on 500-1500 d orbits that are between 0.1-0.4m_j. HD 208487 contains an inner gas giant along with a resonant pair of outer, Saturn/sub-Saturn planets that are close to a second order period resonance.

As part of our analysis, we examined all of the available stellar activity indices, including, to our knowledge, the first use of Bayesian analysis on stellar activity indices to check for activity-induced RV signals. We find that none of the stellar indices have peaks matching the RV periods of the discovered planets, nor their harmonics. Although the Bayesian analysis of the stellar activity indices revealed one signal (Combined FWHM) in proximity to the RV signals, it is statistically distinct based on the current data. Additionally, none of the stellar activity indices are statistically significantly correlated with the observed RVs, nor reach the level of even a moderate correlation, which was not the case for analyses done on stars like HD 26965, for example Díaz et al. (2018). Furthermore, during the GLS periodogram analysis with EXOSTRIKER, we split the RV data sets into two separate sets, based on their S-index values. Although only the inner gas giant was detected in the active data set, three signals were found in the inactive data set, the inner gas giant, the ~500 d signal, and the Saturn-mass planet, HD 208487c. Since the model fit of the inactive set was significantly better, the stellar activity may actually be masking the Doppler signals rather than contributing to them. Again, we reiterate that more data are needed to confirm the nature of the intermediate ~500 d signal.

In order to evaluate the stability of the signals over time, we split the data into early and later subsets. HD 208487b remains entirely consistent in both data sets. HD 208487c is also found in both data sets, with the period and semi-amplitude remaining consistent within uncertainties. The major difference between the earlier and later data sets relates to the period of the third significant signal. The earlier data yields the longest-period signal, while an intermediate-period signal arises in the later data. Both the overall best-fit EMPEROR model and an independent analysis favor the longest-period signal, leading to the adoption of HD 208487d as an outer planet. However, we note that often times our Bayesian analyses arrived at a model favoring the solution with the shorter-period planet; however, even in the best case scenarios the model probabilities were far below the models containing the outer signal. This might indicate that there is indeed another Doppler signal close to a 2:1 period ratio with the HD 208487c signal, and from a dynamical standpoint this also may be favored, since it would be another example of a dynamically packed planetary system.

Extending our investigation of this system, we checked ASAS photometry data to search for any periodic variations in the star’s visible magnitude. Examining the data showed no significant peaks within the photometry. However, the highest power peaks, paired with stellar activity indices from TERRA2, suggest a stellar rotational period around P_rot ~ 29 d.

Fig. 8 GLS periodogram of cleaned ASAS photometry. Top: original GLS periodogram of combined mean-subtracted ASAS photometry data. Second: GLS periodogram after fully fitting the first signal. Third: GLS periodogram after fully fitting the second signal. The blue dotted horizontal lines represents the 0.1% FAP level estimated from 5000 bootstrap resamplings.

8.1 Previous signal conclusions

After the discovery of HD 208487b by Tinney et al. (2005) and the subsequent detection of HD 208487c by Gregory (2005), some disagreements about the planetary configuration arose. Gregory (2007), Goździewski & Migaszewski (2006), and Wright et al. (2007) also detected HD 208487c as a moderately eccentric (e ≃ 0.4) sub-Jovian world; however, the latter two papers find equally probable solutions at P_c ~ 14.5 d, P_c ~ 28 d, or as an activity-induced RV signal, leading them to disregard the existence of HD 208487c as an outer planet. Babu et al. (2010) also claim a potential planet around P_c ~ 14.5 d.

Conclusions drawn from our analysis, applying updated analysis methods and much more data over a longer baseline and at high RV precision, posits that the stellar rotation period is around P_rot ~ 29 d, possibly with some evidence for differential rotation, and therefore this leads to the conclusion that Goździewski & Migaszewski (2006); Wright et al. (2007); Babu et al. (2010) actually found the rotation period and a likely alias at half this rotation period. Further, we believe that the detection of HD 208487c as a moderately eccentric, sub-Jovian world fits the original data, as seen in Section 5.2.3. The new, higher-precision RV data collected over the last 15+ years now enables the detection of candidates HD 208487d, which would be an outer super-Neptune-sub-Saturn planet, and HD 208487c, which would be a Saturn-mass planet in a near-circular orbit, explaining the different planetary models discussed in these works.

8.2 Habitable zone planets

As we have found new planets in this work, we can try to more deeply understand and characterize them. Without observing a transit, we do not have a way to constrain their radii, meaning we can’t determine their densities or physical parameters. However, since their orbits are well-constrained, we can gain a better understanding of this multi-planet system. Based on the definition described by Kasting et al. (1993) and Kopparapu et al. (2013), we can measure the limits of the habitable zone (HZ) of HD 208487, using the following equations: $$S_{eff}=S_{eff\odot }+a{T}_{\star }+b{T}_{\star }^2+c{T}_{\star }^3+d{T}_{\star }^4,$$(4)

where S_eff is the effective solar flux; T_⋆ = T_eff - 5780 K (T_eff is the stellar effective temperature); and a, b, c, d, and S_eff⊙ are constants whose values are given in Table 3 of Kopparapu et al. (2013); $$d=\sqrt{\frac{L_{\star }/L_{\odot }}{S_{eff}}},$$(5)

where L_⋆ and L_⊙ is the luminosity of HD 208487, in this case, and the Sun, respectively, and the radius, d, is measured in astronomical units (AU). We used the coefficient values from the recent Venus and early Mars models (Kopparapu et al. 2013). From the results of the ARIADNE code, the luminosity of HD 208487 is L_⋆ = 1.71 L_⊙. As a result we found that the HZ ranges from 0.97 AU to 2.25 AU. We find that HD 208487c would orbit the star at around 1.94 AU and HD208487d would orbit at around 2.54AU, so HD208487c would lie within the HZ, while HD208487d would orbit just outside the proposed range of the HZ (Jones et al. 2006).

8.3 Exomoons

Long-period exoplanets are also interesting as they could potentially host natural satellites (also known as exomoons). This is because, for the case of planets in close-in orbits, the extreme gravitational interactions with the host stars could destabilize potential exomoon orbits (e.g., Namouni 2010; Spalding et al. 2016; Dobos et al. 2021). Moons are thought to be formed through one of two processes, i.e., coaccretion in the circumplanetary discs around the young planets (e.g., Heller et al. 2014; Barr 2016; Nakajima et al. 2022) and planetary collisions and/or captures (Heller et al. 2014; Barr 2016; Barr & Bruck Syal 2017). Since coaccretion in the circumplanetary discs might require larger amounts of material to form moons, this process might be possible in the case of larger Jupiter-sized planets, as could be the case with HD 208487c. However, for HD 208487d, which appears to be a Neptunian-sized, planetary collisions/capture could potentially lead to a moon formation.

The critical distance from the host planet, within which a moon could have a stable orbit, depends upon the Hill radius of the planet, the eccentricity of the planet’s orbit, and the relative inclination of the moon’s orbit (i.e., if the moon’s orbit is prograde or retrograde with respect to the planet). This critical orbital distance can be expressed as 0.4895R_H(1 - 1.0305e) for a prograde orbit of the moon, and 0.9309R_H(1.0000 - 1.0764e) for a retrograde orbit of the moon (Domingos et al. 2006; Namouni 2010), where R_H is the Hill radius, and e is the orbital eccentricity of the host planet.

Since both HD 208487c and HD 208487d signals have been detected only through the RV method, we can only estimate the minimum mass, and thus, the minimum radius of the Hill sphere for them, which we have estimated to be 12.935 MKm and 13.156 MKm respectively. The minimum distance from the hoststar within which a moon can survive in a prograde orbit around HD 208487c and HD 208487d are 4.993 MKm and 5.599 MKm respectively, and those for a retrograde orbit are 9.57 MKm and 10.786 MKm respectively. For comparison, the distance of the outermost Galilean moon, Callisto, from Jupiter is 1.883 MKm, and the largest retrograde satellite in our Solar system, Triton, orbits around Neptune only at a distance of 0.354 MKm. Thus, both HD 208487c and HD 208487d are capable of hosting potential exomoons, and since they are both orbiting in the HZ around their host star, potential presence of HZ exomoons make them even more interesting. However, detecting such exomoons would require extremely high-precision photometric and spectroscopic follow-ups, which could only be possible using the next generation large telescopes (e.g., Simon et al. 2009; Kipping 2009; Timmermann et al. 2020; Saha & Sengupta 2022; Saha 2022, 2024).

8.4 Dynamics

Motivated by the apparent architecture of the HD 208487 system, and with the aim of testing if such a system could be possible and understandable, we explore the possibility that the system’s origin could be attributed to a dynamical instability of an initially more ordered system. The star hosts two moderately eccentric inner planets (e₁ ~ 0.37 and e₂ ~ 0.2), each three and two times more massive than the outer Neptune-like planet. The two outer planets are near resonance, while the innermost giant appears dynamically detached from the outer two.

We posit that this more ordered system consisted of a coplanar six-planet system orbiting a M_⋆ = 1.15 M_⊙ central star in a near-resonant chain configuration (with consecutive period ratios near 3:2). We envision a scenario in which a six-planet system was formed in a gaseous disk and got captured in resonance due to planet-disk interactions. After dispersal of the disk, the resonance chain is broken and the system becomes unstable - a scenario reminiscent to breaking of the resonance chains applied to super-Earths in the Kepler sample (Izidoro et al. 2017). However, unlike the Kepler planets where instabilities mainly lead to mergers between planets, the cold super-Neptunes in HD 208487 may lead to either collisions between planets or ejections from the system. We can quantify their relative probabilities using the so-called Safronov number for parabolic orbits as $$\begin{array}{r}{\theta }^2\equiv \left(\frac{2G{M}_p}{R_p}\right)\left(\frac{a_p}{G{M}_{\star }}\right)\\ =\left(\frac{M_p}{M_{\mathrm{J}}}\frac{R_{\mathrm{J}\mathrm{u}\mathrm{p}}}{R_p}\frac{a_p}{0.25\mathrm{a}\mathrm{u}}\right)\left(\frac{M_{\odot }}{M_{\star }}\right),\end{array}$$(6)

where for θ² ≫ 1 the planetary radius is so small that close encounters mostly lead to ejections relative to collisions, while for θ² < 1 collisions are more frequent (Ford et al. 2001; Petrovich et al. 2014). A typical super-Earth at 0.5 au has θ² ~ 0.25 leading to collisions and small eccentricity excitation (Ford et al. 2001; Petrovich et al. 2014), but a cold Neptune at ~2 au has θ² ~ 1 leading to mixed outcomes and potential mergers with significant eccentricity excitation. We numerically evaluate this possibility next.

We perform 250 simulations using the REBOUND code (Rein & Liu 2012) with the primary objective of investigating the possibility of achieving an orbital configuration resembling the observed data, thereby deviating from the conventional peas in a pod scenario. The planetary masses are m_i = 0.15 m_Jup with i = 1... 6, and we place the planets 5% away from the exact 3:2-3:2-3:2-3:2-3:2 resonant chain. We chose 5 values for the initial eccentricities for the two inner planets: e₁ = e₂ = 0.05,0.1,0.15,0.2,0.25 to speed up the phase of dynamical instability, while the outer four bodies start with a fixed value of 0.05. For each of them we randomly selected 50 values for the initial argument of periastron ω_i ∈ [0,2.π] and inclinations randomly distributed between 0° and 2°. Collisions happen whenever the distance between two planets is less than the sum of their radii: d < r_i + r_j, where the radius of each planet follows the $r_i\propto m_i^{-2.06}$ law by Lissauer et al. (2011).

Over a simulation period of 5 × 10⁵ years, none of the systems remained stable and maintained their six-planet configuration. In this compact configuration, all the simulations underwent an instability leading to planet collisions generally leaving behind 3, 4, and 5 planets. We note that the rate of stable (unstable) systems is highly sensitive to the initial conditions and will likely decrease (increase) as our simulations continue to evolve for longer timescales. Around 50% of the systems that underwent an instability, led to the formation of a three-planet system. Of these, in ~20% of cases, the three inner planets merged as a massive innermost body, and the fourth and fifth planets merged as an intermediate-mass middle planet, resulting in a three-planet system with an architecture similar to the one observed in our data, with inclinations reaching a final root mean square (RMS) of 1°. Figure 9 illustrates the temporal evolution of a representative simulation where, at ~5 × 10⁴ years, a collision occurs leading to the merging of the inner 2 bodies into a single planet, after which a second collision occurs at ~10⁵ years, finally forming the innermost body, m₁. Around ~5 × 10⁴ years, a collision of the 4^th and 5^th forms the middle planet, m₂, while the outermost body never collides with other bodies. The top panel displays the semimajor axis, peri-and apocentric distances of each planet (a_i(1 + e_i), respectively), while the bottom panel showcases the eccentricities as a function of time. The post-merger planet attains high eccentricities, providing a potential explanation for the observed eccentricities of the inner planets.

Although the approximation of treating collisions as perfectly inelastic mergers might not capture the complete realism of the scenario, it provides valuable insights as a proof of concept for understanding how such a system can form. As already stated by Izidoro et al. (2017), once the gas dissipates resonant chains may become dynamically unstable, undergoing a phase of gravitational scattering. Our simulations begin after the dissipation of the gas, exploiting the phase of dynamical instability, allowing the formation of a system featuring a fairly eccentric massive inner planet, largely separated from the outer pair.

Fig. 9 Temporal evolution of a six-planet simulation showing a merger between the first three and the fourth and fifth planets, leaving behind a three-planet system with similar masses to the observations. The upper panel displays the semimajor axis ai, the variation of the apocenter Q (defined as ai[1 + ei]), and the pericenter (q = ai[1 - ei]). In the bottom panel we show the variation of the eccentricities.

9 Conclusions

This paper provides a thorough re-investigation into the HD 208487 planetary system. First, we derived an updated set of stellar parameters utilizing the stellar modeling codes SPECIES and ARIADNE. Next, we comprehensively analyzed the currently available RV data of HD 208487 to search for signals not previously found in older investigations. Further, we developed a more rigorous methodology to analyze stellar activity indices. Finally, we evaluated ASAS and Hipparcos data to search for any periodic photometric signals.

As a result of our RV analysis, we determined that HD 208487 possibly contains at least three planets in the system (see Table 3). We provided updated parameters to the sub-Jovian HD 208487b with a period of 129.36 d. Formerly believed to be a possible sub-Jovian and moderately eccentric planet (Gregory 2005), we now show that HD 208487c is actually a candidate Saturn with a period of 923.06 d and a near-circular orbit. Additionally, we discovered a new super-Neptune-sub-Saturn candidate, HD 208487d, a planet that, if confirmed, would be orbiting with a period of 1380.13 d. HD 208487c and HD 208487d have an orbital period ratio of P_d/P_c = 1.495, indicating they are very close to a 3:2 MMR. There is also some evidence that leads us to believe that there may exist a fourth planet in this system at close to the 2:1 period ratio with HD 208487c, which would further dynamically pack the planetary system; however, more long-baseline RV data will be required to confirm this as a reality.

Our methodology to analyze the available stellar activity indices with a Bayesian statistical analysis enabled us to distinguish the RV signals from the statistically significant activity signals. We found no correlations between the RVs and stellar activity signals and no peaks in the GLS at the periods of the three planets. Using EMPEROR to provide a more accurate description, we found that the nearest detected stellar activity signal is separated by 7.1σ from P_c. Thus, we determined that there is a high probability that the RVs and stellar activity signals are distinct, based on the current data. An independent analysis of the RVs using color-dependent terms to detrend sources of color-dependent noise also did not yield any correlations, arriving at the same system configuration as that of the EMPEROR analysis.

Concluding our observational analysis, we studied the ASAS and Hipparcos photometric data. After fitting the long-period baseline-effect signal, we found a statistically significant period at ~29 d. Combining this result with our GLS and EMPEROR analysis of the stellar activity indices, we believe this signal to indicate the rotation period of the star, with potential differential rotation. This result explains the false RV signals, and a likely half-period alias, posited by Gozdziewski & Migaszewski (2006), Wright et al. (2007), and Babu et al. (2010).

Finally, we performed dynamical simulations in an effort to better understand the formation history of this complex system. By assuming the common peas in a pod configuration to begin with, including additional planets in the system close to a resonant chain, we were able to show that the inner moderately eccentric gas giant could be the remnant of a previous collision of lower-mass planets. Of the systems that underwent an instability and subsequent planetary collisions, such a scenario occurred ~20% of the time in our simulation set. This scenario originates from systems with Safronov numbers close to unity, which allows for significant scattering and eccentricity excitation before the merger occurs, possibly representing the breaking point of more ordered configurations and providing a chaotic look at the final architectures of planetary systems. Systems like HD 208487 in our simulations end up with relatively low inclinations (RMS ~ 1°), similar to the inclinations of Kepler planets within three-planet systems (Fabrycky et al. 2014; Zhu et al. 2018). However, unlike Kepler planets, HD 208487 planets orbit with semimajor axes >1 au, unreachable by Kepler. Future missions such as PLATO may spot systems with these architectures.

Data availability

Full Tables A.1–A.4 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/702/A139.

Acknowledgements

Part of this work was supported by the United States Fulbright Fellowship in partnership with the Fulbright Chile Commission. JSJ gratefully acknowledges support by FONDECYT grant 1201371, from the ANID BASAL project FB210003, and from CASSACA grant CCJRF2205. CC was supported by FONDECYT grant n° 3230283. CP acknowledges support from CATA-Basal AFB-170002, ANID BASAL project FB210003, FONDECYT Regular grant 1210425, CASSACA grant CCJRF2105, and ANID+REC Convocatoria Nacional subvencion a la instalacion en la Academia convocatoria 2020 PAI77200076. MT acknowledges support from the Jenny and Antti Wihuri Foundation. This work has also made use of the University of Hertfordshire’s high-performance computing facility.

Appendix A RV and stellar activity index data

The following tables list the first five RVs and stellar activity index data for each instrument utilized in this paper. The full data can be found at the CDS.

Appendix B Split RV GLS periodograms

Figure B.1 shows the GLS Periodogram of the Low S-index data set, while Figure B.2 shows the GLS Periodogram of the High S-index data set.

Fig. B.1 Top panel: Original GLS periodogram of the combined and mean-subtracted low S-index RV data. Second panel: GLS periodogram after fitting and removing the first signal. Third panel: GLS periodogram after fitting and removing both the first and second signals. Bottom panel: GLS periodogram after manually fitting the third signal and removing all three previously detected signals. The blue dotted horizontal lines represent the 0.1% FAP level estimated from 5000 bootstrap resamplings.

Fig. B.2 Top panel: Original GLS periodogram of the combined and mean-subtracted high S-index RV data. Second panel: GLS periodogram after fitting and removing the first signal. The blue dotted horizontal lines represent the 0.1% FAP level estimated from 5000 bootstrap resamplings.

References

All Tables

All Figures

	Fig. 1 All 219 mean-subtracted radial velocity measurements of HD 208487, including 48 UCLES (gray), 20 HARPS-TERRA (red), 106 HARPS-TERRA2 (olive), and 45 PFS (teal) observations.
In the text

Fig. 2 Top panel: original GLS periodogram of the combined and mean-subtracted RV data (black). The periodogram of the window function (sampling) for the combined data (green) is also shown. Second panel: GLS periodogram after fitting and removing the first detected signal. Third panel: GLS periodogram after fitting and removing both the first and second signals. Bottom panel: GLS periodogram after fitting and removing all three previously detected signals. The blue dotted horizontal lines represents the 0.1% false alarm probability level estimated from 5000 bootstrap resamplings.

In the text

	Fig. 3 Top to bottom, left to right: Phase-folded best-fit model for the first, second, third, and complete Keplerian model.
In the text

Fig. 4 Top, left to right: final posterior distribution of the period for the first, second, and third Keplerians from the best-fitting parallel-tempered MCMC run. Bottom, left to right: final posterior distribution of the semi-amplitude for the first, second, and third Keplerians from the best-fitting parallel-tempered MCMC run. The numbers at the top of each figure correspond to the mean, variance, skewness, kurtosis, median, and mode, respectively. The solid line represents a Gaussian curve with the same mean and variance.

In the text

Fig. 5 GLS periodogram of the available stellar activity indices. From top to bottom: PFS S-index, TERRA S-index, TERRA FWHM, TERRA BIS, TERRA2 S-index, TERRA2 FWHM, TERRA2 BIS, Combined Sindex, Combined BIS, Combined FWHM. The green solid vertical lines show the best-fit periods of the three planetary signals. The blue dotted horizontal lines represents the 0.1% false alarm probability (FAP) level estimated from 5000 bootstrap resamplings.

In the text

	Fig. 5 Correlations of the available stellar activity indices. From top to bottom, left to right: PFS S-index, TERRA S-index, TERRA FWHM, TERRA BIS, TERRA2 S-index, TERRA2 FWHM, TERRA2 BIS. The lines show the best-fit correlation between the activity indicator and the RVs.
In the text

Fig. 7 Left (right): ASAS V-band photometry (Hipparcos photometry). Top: A and B flagged V-band ASAS photometry (all Hipparcos data) for HD 208487. Middle: cleaned V-band ASAS (Hipparcos) photometry for HD 208487. Bottom: GLS periodogram of the photometry. The horizontal lines correspond to the 0.1% (dotted blue) FAP significance level computed via 5000 bootstrap iterations on the photometry time series. The vertical green lines represent the periods of the three planets.

In the text

	Fig. 8 GLS periodogram of cleaned ASAS photometry. Top: original GLS periodogram of combined mean-subtracted ASAS photometry data. Second: GLS periodogram after fully fitting the first signal. Third: GLS periodogram after fully fitting the second signal. The blue dotted horizontal lines represents the 0.1% FAP level estimated from 5000 bootstrap resamplings.
In the text

Fig. 9 Temporal evolution of a six-planet simulation showing a merger between the first three and the fourth and fifth planets, leaving behind a three-planet system with similar masses to the observations. The upper panel displays the semimajor axis ai, the variation of the apocenter Q (defined as ai[1 + ei]), and the pericenter (q = ai[1 - ei]). In the bottom panel we show the variation of the eccentricities.

In the text

Fig. B.1 Top panel: Original GLS periodogram of the combined and mean-subtracted low S-index RV data. Second panel: GLS periodogram after fitting and removing the first signal. Third panel: GLS periodogram after fitting and removing both the first and second signals. Bottom panel: GLS periodogram after manually fitting the third signal and removing all three previously detected signals. The blue dotted horizontal lines represent the 0.1% FAP level estimated from 5000 bootstrap resamplings.

In the text

	Fig. B.2 Top panel: Original GLS periodogram of the combined and mean-subtracted high S-index RV data. Second panel: GLS periodogram after fitting and removing the first signal. The blue dotted horizontal lines represent the 0.1% FAP level estimated from 5000 bootstrap resamplings.
In the text