Revised orbits of the two nearest Jupiters

ABSTRACT

With its near-to-mid-infrared high-contrast imaging capabilities, JWST is ushering us into a golden age of directly imaging Jupiter-like planets. As the two closest cold Jupiters, ε Ind A b and ε Eridani b have sufficiently wide orbits and adequate infrared emissions to be detected by JWST. To detect more Jupiter-like planets for direct imaging, we develop a gost-based method to analyse radial velocity data and multiple Gaia data releases simultaneously. Without approximating instantaneous astrometry by catalogue astrometry, this approach enables the use of multiple Gaia data releases for detection of both short-period and long-period planets. We determine a mass of |$2.96_{-0.38}^{+0.41}$| M_Jup and a period of |$42.92_{-4.09}^{+6.38}$| yr for ε Ind A b. We also find a mass of |$0.76_{-0.11}^{+0.14}$| M_Jup , a period of |$7.36_{-0.05}^{+0.04}$| yr, and an eccentricity of 0.26|$_{-0.04}^{+0.04}$| M_Jup, for ε Eridani b. The eccentricity differs from that given by some previous solutions, probably due to the sensitivity of orbital eccentricity to noise modelling. Our work refines the constraints on orbits and masses of the two nearest Jupiters and demonstrate the feasibility of using multiple Gaia data releases to constrain Jupiter-like planets.

methods: data analysis, techniques: radial velocities, astrometry, exoplanets, stars: individual: ε Eridani, stars: individual: ε Ind A

1 INTRODUCTION

The detection and characterization of Jupiter-like planets in extrasolar systems are crucial for our understanding of their formation and evolution, as well as their role in shaping the system architecture and habitability (Stevenson & Lunine 1988; Lunine 2001; Tsiganis et al. 2005; Horner et al. 2020). While the current ground-based facilities are mainly sensitive to young Jupiters, the Mid-Infrared Instrument (MIRI; Rieke et al. 2015) mounted on the JWST is optimal for imaging cold Jupiters, which are far more abundant than young Jupiters. Equipped with a coronagraph, MIRI has an imaging observing mode from 4.7 to 27.9 μm , sensitive to the mid-IR emission from cold Jupiters, with equilibrium temperature of less than 200 K (Bouchet et al. 2015). On the other hand, the JWST Near-IR Spectrograph (Jakobsen et al. 2022) conducts medium-resolution spectroscopic observation from 0.6 to 5.3 μm . The differential Doppler shift of planet relative to its host star could be used to disentangle planetary, and stellar spectrums for planet imaging (Llop-Sayson et al. 2021).

Considering that high-contrast imaging typically requires large planet-star separation and high planetary-mass, nearby cold Jupiters are optimal targets for JWST imaging. However, the orbital periods of cold Jupiters are as long as decades, beyond most of the individual radial velocity (RV), astrometric, and transit surveys. Hence, it is important to form extremely long observational time span (or baseline) by combining different types of data. One successful application of such synergistic approach is combined analyses of legacy RV data and the astrometric data from Hipparcos (Perryman et al. 1997; van Leeuwen 2007) and Gaia data releases (DRs) (Gaia Collaboration 2016, 2018, 2021, 2023). This approach not only provides decade-long baseline to constrain the orbit of a cold Jupiter, but also fully constrain planetary mass and orbit due to the complementarity between RV and astrometry. Assuming that the catalogue astrometry approximates the instantaneous astrometry at the reference epoch, one can use the difference between Hipparcos and Gaia astrometry to constrain the non-linear motion of a star (or reflex motion) induced by its companion (Brandt 2018; Snellen & Brown 2018; Kervella et al. 2019; Feng et al. 2019c). However, under this assumption, Gaia DR2 and DR3 are not independent and thus cannot be used simultaneously to constrain planetary orbit. Moreover, this assumption is not appropriate for constraining orbits with periods comparable or shorter than the observational time span of the epoch data used to produce an astrometric catalogue.

To use multiple Gaia DRs to constrain both short and long-period planets, we simulate the Gaia epoch data using the Gaia Observation Forecast Tool (gost¹) and fit the synthetic data by a linear astrometric model, corresponding to the Gaia five-parameter solution. The difference between the fitted and the catalogue astrometry constrains the non-linear reflex motion induced by planets. In this work, we apply this approach to constrain the orbits and masses of the two nearest cold Jupiters, ε Ind A b and ε Eridani b.

This paper is structured as follows. In Section 2, we introduce the RV and astrometry data measured for ε Ind A and ε Eridani. We then describe the modelling and statistical methods used in our analyses in Section 3. We give the orbital solutions for the two planets and compare our solutions with previous ones in Section 4. Finally, we discuss and conclude in Section 5.

2 RV AND ASTROMETRY DATA

ε Ind A is a K-type star with a mass of 0.754 ± 0.038 M_⊙ (Demory et al. 2009), located at a heliocentric distance of 3.7 pc. It hosts two brown dwarfs (Smith et al. 2003) and a Jupiter-like planet (Feng et al. 2019c). It has been observed by high-precision spectrographs including the High Accuracy RV Planet Searcher (HARPS; Pepe et al. 2000) mounted on the ESO La Silla 3.6 m telescope, the ESO UV-visual echelle spectrograph (UVES) on Unit 2 of the Very Large Telescope array (Dekker et al. 2000), and the Coudé Echelle Spectrograph (CES) at the 1.4 m telescope in La Silla, Chile. We use ‘HARPSpre’ and ‘HARPSpost’ to denote HARPS data obtained before and after the fibre change in 2015 (Lo Curto et al. 2015). The HARPS data are reduced by Trifonov et al. (2020) using the SpEctrum Radial Velocity AnaLyser (SERVAL) pipeline (Zechmeister et al. 2018). The RV data obtained by the Long Camera (LC) and the Very LC (VLC) were released by Zechmeister et al. (2013).

As a K-type star with a mass of 0.82 ± 0.02 M_⊙ (Gonzalez, Carlson & Tobin 2010), ε Eridani is the third closest star system to the Earth. It is likely that a Jupiter analog orbits around this star on a wide orbit (Hatzes et al. 2000; Butler et al. 2006; Mawet et al. 2019). It has been observed by the CES, LC, and VLC (Zechmeister et al. 2013), the Lick Observatory Hamilton echelle spectrometer (Vogt 1987), the Automated Planet Finder (APF), the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 2014, 1994) at the Keck observatory, the HARPS for the Northern hemisphere (HARPS-N or HARPN; Cosentino et al. 2012) installed at the Italian Telescopio Nazionale Galileo, the EXtreme PREcision Spectrograph (EXPRES; Jurgenson et al. 2016) installed at the 4.3 m Lowell Discovery Telescope (Levine et al. 2012). By accounting for the RV offets caused by the updates of Lick Hamilton spectrograph, we use ‘Lick6’, ‘Lick8’, and ‘Lick13’ to denote multiple data sets, following the convention given by Fischer, Marcy & Spronck (2013). We use the APF and HIRES data reduced by the standard CPS pipeline (Howard et al. 2010) and released by Mawet et al. (2019) as well as the other APF data reduced using the pipeline developed by Butler et al. (1996) and the HIRES data released by Butler et al. (2017). We use ‘APFh’ and ‘HIRESh’ to denote the former sets and use ‘APFp’ and ‘HIRESp’ to denote the latter. The APFp data set is presented in Table C1. The EXPRES data is from Roettenbacher et al. (2022), and the HARPS data is obtained from Trifonov et al. (2020).

We obtain the Hipparcos epoch data for ε Ind A and ε Eridani from the new Hipparcos reduction (van Leeuwen 2007). We use the Gaia second and third data releases (DR2 and DR3; Gaia Collaboration 2021, 2023) as well as the epoch data generated by gost. The Gaia first data release (DR1) is not used because it is derived partly from the Hipparcos data (Michalik, Lindegren & Hobbs 2015), and is thus not treated as independent of the Hipparcos data.

3 METHOD

Considering that the Gaia intermediate astrometric data (or epoch data) are not available in the current Gaia DRs, techniques have been developed to use the difference between Gaia and Hipparcos catalogue data to constrain the orbits of substellar companions (Brandt, Dupuy & Bowler 2019; Kervella et al. 2019; Feng et al. 2022). Without using epoch data, previous studies are limited by approximating the simultaneous astrometry at the reference epoch with the catalogue astrometry (Feng et al. 2022). In other words, a linear function is not appropriate to model the centre (position) and tangential line (proper motion) of the orbital arc of the stellar reflex motion, induced by short-period companions.

To avoid the above assumptions, we use gost to predict the Gaia observation epochs for a given star. Considering that the models of RV and reflex motion are already introduced in the previous papers authored by some of us (Feng et al. 2019c, 2021), we introduce the newly developed technique of using Hipparcos intermediate data and Gaiagost data as follows.

• Obtain data. For a given target, we obtain the revised Hipparcos intermediate data IAD from van Leeuwen (2007) and Gaiagost data, including the scan angle ψ, the along-scan (AL) parallax factor f^AL, and the observation time at barycentre. In previous studies, different versions of the Hipparcos catalogue data have been recommended (ESA 1997; van Leeuwen 2007). However, in our research, we find that the choice of the Hipparcos version has minimal impact on our orbital solutions. We attribute this to several reasons:

  • Our approach involves modelling the systematics in Hipparcos IAD using offsets and jitters instead of calibrating them a priori, as described in Brandt, Michalik & Brandt (2023). By incorporating these offsets and jitters, we account for the systematic effects in the data, making our solutions less sensitive to the specific Hipparcos version used.

  • The astrometric precision of Hipparcos is considerably inferior to that of Gaia. Additionally, the time difference between Gaia and Hipparcos is much greater than the duration covered by Gaia DRs. Consequently, when it comes to constraining long-period orbits, the crucial factor is the temporal baseline between Hipparcos and Gaia, rather than the particular version of the Hipparcos catalogue.

  • For short-period orbits, it is the curvature of the Hipparcos IAD that primarily constrains the orbit, rather than the absolute offset. Hence a calibration of the offset of Hipparcos IAD becomes less critical in determining short-period orbits.

• Model astrometry of target system barycentre (TSB) at Gaia DR3 reference epoch. We model the astrometry of the TSB at the Gaia DR3 epoch J2016.0 (t_DR3) as follows:

  $$\begin{eqnarray} \alpha ^b_{\rm DR3}=\alpha _{\rm DR3}-\frac{\Delta \alpha _*}{\cos \delta _{\rm DR3}}~, \end{eqnarray}$$

  (1)

  $$\begin{eqnarray} \delta ^b_{\rm DR3}=\delta _{\rm DR3}-\Delta \delta ~, \end{eqnarray}$$

  (2)

  $$\begin{eqnarray} \varpi ^b_{\rm DR3}=\varpi _{\rm DR3}-\Delta \varpi ~, \end{eqnarray}$$

  (3)

  $$\begin{eqnarray} \mu ^b_{\alpha {\rm DR3}}=\mu _{\alpha {\rm DR3}}-\Delta \mu _\alpha ~, \end{eqnarray}$$

  (4)

  $$\begin{eqnarray} \mu ^b_{\delta {\rm DR3}}=\mu _{\delta {\rm DR3}}-\Delta \mu _\delta ~, \end{eqnarray}$$

  (5)

  where α, δ, ϖ, μ_α, and μ_δ are R.A., decl., parallax, and proper motion in R.A and decl,² subscript _DR3 represents quantities at epoch t_DR3, superscript ^b represents TSB, and Δ means offset. Considering that the Gaia measurements of systematic RVs are quite uncertain, we use Gaia DR3 RVs to propagate astrometry instead of using them to constrain reflex motion.
• Model astrometry of TSB at Hipparcos and Gaia DR3 reference epochs. We model the TSB astrometry at the Hipparcos reference epoch t_HIP through linear propagation of state vectors in the Cartesian coordinate system as follows:

  $$\begin{eqnarray} \vec {r}_{\rm HIP}^b=\vec {r}^b_{\rm DR3}+\vec {v}^b_{\rm DR3}(t_{\rm HIP}-t_{\rm DR3})~, \end{eqnarray}$$

  (6)

  $$\begin{eqnarray} \vec {v}_{\rm HIP}^b=\vec {v}^b_{\rm DR3}~, \end{eqnarray}$$

  (7)

  where |$(\vec {r}_{\rm HIP}^b,\vec {v}_{\rm HIP}^b)$| and |$(\vec {r}_{\rm DR3}^b,\vec {v}_{\rm DR3}^b)$| are, respectively the state vectors, including location and velocity, of the TSB at the Hipparcos and Gaia DR3 epochs. We first transform TSB astrometry from equatorial coordinate system, |$\vec {\iota }^b_{\rm DR3}=(\alpha ^b_{\rm DR3},\delta ^b_{\rm DR3},\varpi ^b_{\rm DR3},\mu ^b_{\alpha {\rm DR3}},\mu ^b_{\delta {\rm DR3}},{\rm RV}_{\rm DR3})$|⁠, to Cartesian coordinate system, to get state vector at Gaia DR3 epoch, and then propagate the vector to the Hipparcos epoch, and then transform the new vector back to the astrometry at the Hipparcos epoch, |$\vec {\iota }^b_{\rm HIP}$|⁠. The whole process is: equatorial state vector at t_DR3→ Cartesian state vector at t_DR3 → linear propagation of Cartesian state vector to t_HIP→ equatorial state vector at t_HIP. By propagating state vectors in Cartesian coordinate system instead of spherical coordinate system, this approach completely solve the problem of perspective acceleration. The transformation between different coordinate systems is described in Lindegren et al. (2012) and Feng et al. (2019b).

• Simulate Gaia abscissae using GOST. Instead of simulating the Gaia epoch data precisely by considering relativistic effects, perspective acceleration, and instrumental effects (Lindegren et al. 2012), we simulate the Gaia abscissae (or along-scan coordinates) by only considering linear motion of the TSB in the equatorial coordinate system as well as the target reflex motion.³ This is justified because the various effects are independent of reflex motion, and can be estimated and subtracted from the data a priori.

  We simulate the position of the target at gost epoch t_j relative to the Gaia DR3 reference position by adding the stellar reflex motion (denoted by superscript ^r) on to the TSB motion,

  $$\begin{eqnarray} \Delta \alpha _{*i}=\Delta \alpha ^b_{*\rm DR3}+\mu _{\alpha {\rm DR3}}^b(t_i-t_{\rm DR3})+ \Delta \alpha ^r_{*i}~, \end{eqnarray}$$

  (8)

  $$\begin{eqnarray} \Delta \delta _{i}=\Delta \delta ^b_{\rm DR3} + \mu _{\delta {\rm DR3}}^b(t_i-t_{\rm DR3}) + \Delta \delta ^r_i~, \end{eqnarray}$$

  (9)

  where |$\Delta \alpha ^b_{*\rm DR3}=(\alpha ^b_{\rm DR3}-\alpha _{\rm DR3})\cos \delta ^b_{\rm DR3}$|⁠, and |$\Delta \delta ^b_{\rm DR3}=\delta ^b_{\rm DR3}-\delta _{\rm DR3}$|⁠. Because the reflex motion caused by cold Jupiters is insignificant compared with barycentric motion, the parallax at epoch t_i is approximately |$\varpi _i=\varpi ^b_i+\Delta \varpi ^r_i=\varpi ^b_{\rm DR3}+\Delta \varpi ^r_i\approx \varpi ^b_{\rm DR3}$|⁠, and the systematic RV is approximately |${\rm RV}_i={\rm RV}^b_i+\Delta {\rm RV}^r={\rm RV}^b_{\rm DR3}+\Delta {\rm RV}^r\approx {\rm RV}_{\rm DR3}$|⁠.
  Considering that Gaia pixels in along-scan (AL) direction are much smaller than that in the cross-scan direction, we only model the along-scan position of the target (or abscissa; η_i). Considering the parallax caused by Gaia’s heliocentric motion, abscissa is modelled by projecting the motion of the target on to the AL direction using

  $$\begin{eqnarray} \eta _i= \Delta \alpha _{*i}\sin \psi _i+\Delta \delta _i\cos \psi _i+\varpi ^b_{\rm DR3}f^{\rm AL}_i~. \end{eqnarray}$$

  (10)
• Fit a five-parameter model to synthetic Gaia abscissae. Considering that binary solution is only applied to a small fraction of DR3 targets and most planet-induced reflex motion is not yet available in the Gaia non-single star catalogue (Gaia Collaboration 2023), we model the simulated abscissae using a five-parameter model as follows:

  $$\begin{eqnarray} \hat{\eta }_i=\Delta \alpha _{*{\rm DR3}}^l\sin \psi _i+\Delta \delta _{\rm DR3}^l\cos \psi _i+\hat{\varpi }_{\rm DR3} f^{\rm AL}_i~, \end{eqnarray}$$

  (11)

  $$\begin{eqnarray} \Delta \alpha _{*{\rm DR3}}^l=(\hat{\alpha }_{\rm DR3}-\alpha _{\rm DR3})\cos \hat{\delta }_{\rm DR3}+\hat{\mu }_{\alpha {\rm DR3}}(t_i-t_{\rm DR3})~, \end{eqnarray}$$

  (12)

  $$\begin{eqnarray} \Delta \delta _{\rm DR3}^l=(\hat{\delta }_{\rm DR3}-\delta _{\rm DR3})+\hat{\mu }_{\delta {\rm DR3}}(t_i-t_{\rm DR3})~, \end{eqnarray}$$

  (13)

  where |$\hat{\vec \iota} _{\rm DR3}=(\hat{\alpha }_{\rm DR3}, \hat{\delta }_{\rm DR3}, \hat{\varpi }_{\rm DR3}, \hat{\mu }_{\alpha {\rm DR3}},\hat{\mu }_{\delta {\rm DR3}})$| is the set of model parameters at t_DR3, |$f_i^{\rm AL}$| is the along-scan parallax factor, and ψ_i is the scan angle at epoch t_i. This scan angle is the complementary angle of ψ in the new Hipparcos IAD (van Leeuwen 2007; Brandt et al. 2021b; Holl et al. 2022), and thus ψ will be replaced by π/2 – ψ when modelling Hipparcos IAD. The above modelling of Gaia DR3 can be applied to Gaia DR2 by changing the subscript _DR3 into _DR2. Given the limited information available through gost, we are unable to reconstruct the uncertainties of individual observations, and which epochs are actually used in producing the catalogues. As long as the astrometric uncertainties and rejected epochs are not significantly time-dependent, it is reasonable to assume that all gost epochs are used in the astrometric solution of Gaia and all abscissae have the same uncertainty. Under this assumption, we fit the five-parameter model shown in equation (11) to the simulated abscissae (η_i) for Gaia DR2 and DR3, through linear regression.
• Calculate the likelihood for Gaia DR2 and DR3. To avoid numerical errors, the catalogue astrometry at t_i relative to the Gaia DR3 epoch is defined as |$\Delta \vec {\iota }_i\equiv (\Delta \alpha _{*i},\Delta \delta _i,\Delta \varpi _i,\Delta \mu _{\alpha i},\Delta \mu _{\delta i})=((\alpha _i-\alpha _{\rm DR3})\cos \delta _i,\delta _i-\delta _{\rm DR3},\varpi _i-\varpi _{\rm DR3},\mu _{\alpha i}-\mu _{\alpha {\rm DR3}},\mu _{\delta i}-\mu _{\delta {\rm DR3}})$|⁠. The fitted astrometry for epoch t_i is |$\Delta \hat{\vec \iota }_i$|⁠. The likelihood for Gaia DR2 and DR3 is

  $$\begin{eqnarray} \mathcal {L}_{\rm gaia}=\prod _{i=1}^{N_{\rm DR}}(2\pi)^{-5/2}({\rm det}\Sigma ^{\prime }_i)^{-\frac{1}{2}}e^{-\frac{1}{2}(\Delta \hat{\vec \iota }_i -\Delta {\vec \iota }_i)^T[\Sigma _i(1+J_i)]^{-1}(\Delta \hat{\vec \iota }_i-\Delta {\vec \iota }_i)}~, \\ \end{eqnarray}$$

  (14)

  where N_DR is the number of data releases used in the analyses, and Σ_i (1 + J_i) is the jitter-corrected covariance for the five-parameter solutions of Gaia DRs, i = 1 and 2 respresent Gaia DR2 and DR3, respectively. In this work, we only use DR2 and DR3, and thus N_DR = 2.
• Calculate the likelihood for Hipparcos intermediate data. We model the abscissae of Hipparcos (ξ) by adding the reflex motion on to the linear model shown in equation (11), and calculate the likelihood as follows:

  $$\begin{eqnarray} \mathcal {L}_{\rm hip}=\prod \limits _{i=1}^{N_{\rm epoch}}\frac{1}{\sqrt{2\pi (\sigma _i^2+J_{\rm hip}^2)}}e^{-\frac{(\hat{\xi }_i -\xi _i)^2}{2(\sigma _i^2+J_{\rm hip}^2)}}, \end{eqnarray}$$

  (15)

  where N_epoch is the total number of epochs of Hipparcos IAD.

To obtain the total likelihood (⁠|$\mathcal {L}=\mathcal {L}_{\rm RV}\cdot \mathcal {L}_{\rm hip}\cdot \mathcal {L}_{\rm gaia}$|⁠), we derive the likelihoods for the Gaia and Hipparcos data through the above steps, and calculate the likelihood for the RV data (⁠|$\mathcal {L}_{\rm RV}$|⁠) following Feng, Tuomi & Jones (2017). We adopt log uniform priors for time-scale parameters such as period and correlation time-scale in the moving average (MA) model (Feng et al. 2017), and uniform priors for other parameters. Finally, we infer the orbital parameters by sampling the posterior through adaptive and parallel MCMC, developed by Haario, Saksman & Tamminen (2001) and Feng et al. (2019a).

While packages like orbitize (Blunt et al. 2020), FORECAST (Bonavita et al. 2022), BINARYS (Leclerc et al. 2023), and kepmodel (Delisle & Ségransan 2022) have been developed to analyse imaging, RV, and astrometric data, they are not primarily designed for analysing Gaia catalogue data as performed by orvara (Brandt et al. 2021a) and htof (Brandt et al. 2021b). Compared with orvara and htof, our method has the following features:

• we simultaneously optimize all model parameters through MCMC posterior sampling;

• we utilize multiple Gaia DRs by fitting multiple five-parameter models to Gaia epoch data simulated by gost;⁴

• instead of conducting calibration a priori, we employ jitters and offsets to model systematics a posteriori;

• we model time-correlated noise in the RV data.

While the use of multiple Gaia DRs does not significantly improve the constraint on long-period orbits by increasing the temporal baseline between Gaia and Hipparcos, it does provide additional information about the raw abscissae. This, in turn, leads to stronger orbital constraints at Gaia epochs. In other words, Gaia DR2 provides a 5D astrometric data point, and when combined with the 5D data point from Gaia DR3, the orbital constraint becomes stronger. The additional information obtained from multiple data releases enhances the accuracy and reliability of our orbital solutions.

Through sensitivity tests, we find that our orbital solutions are not strongly sensitive to whether or not we calibrate frame rotation and zero-point parallax a priori by adopting values from previous studies (Brandt 2018; Kervella et al. 2019; Lindegren 2020). It is important to note that there is uncertainty in the estimation of these calibration parameters, as indicated by studies such as Brandt (2018) and Lindegren (2020). Additionally, uncertainties can be amplified during the transformation from the Gaia frame to the Hipparcos frame. Considering these factors, it is more appropriate to consider astrometric systematics on a case-by-case basis, taking into account individual target characteristics. These issues have been discussed in our previous work (Feng et al. 2021).

To validate the accuracy of gost emulations, we perform a comparison between the astrometric epochs generated by gost and the G-band transit times, provided in Gaia DR 3 (GDR3) for a randomly selected sample of 1000 stars that have epoch photometry. We assess the number of mismatched epochs between gost and the GDR3 epoch-photometry catalogue, as well as the distribution of these mismatches. We define N_gp as the number of epochs predicted by gost but not present in the GDR3 epoch-photometry catalogue, and N_pg as the number of epochs present in the photometry catalogue, but not predicted by gost. The distribution of the sample over G magnitude and the number of mismatched epochs is depicted in Fig. 1. Notably, it is apparent that for bright stars with G magnitudes less than 10, the number of photometric epochs mismatched between gost and GDR3 is at most 1.

Figure 1.

Distribution of 1000 stars over the G-band magnitude and the number of mismatched epochs between gost and the GDR3 epoch-photometry catalogue. The left-hand panel is for N_pg, while the right-hand panel is for N_gp.

To quantify the fraction of missing epochs relative to the total photometric epochs, we calculate η_gp = N_gp/N_p and η_pg = N_pg/N_p, where N_p represents the total number of G-band photometric transits. The median of N_p is 22. The mean values of η_gp and η_pg for bright stars are found to be 2.1 and 0.7 per cent, respectively. For all stars, the mean values of η_gp and η_pg are 7.0 and 0.9 per cent, respectively. This analysis serves as a validation of the accuracy of gost emulations, further supporting the reliability of our approach.

4 RESULTS

We find the optimal orbital solutions for ε Ind A b and ε Eridani b based on the MCMC samplings of posteriors and show the solutions in Fig. 2, and the posterior distributions of orbital parameters in Figs A1 and A2. To optimize the visualization of Fig. 2, we project the synthetic abscissae along the R.A. and decl. direction, encode the orbital phase with colours, and represent the orbital direction using circles with arrows. In the panels for Gaia synthetic abscissae (third column), we use segments and shaded regions to visualize the catalogue astrometry (⁠|$\Delta \vec \iota -\Delta \hat{\vec \iota} ^b$|⁠), and fitted astrometry (⁠|$\Delta \hat{\vec \iota} -\Delta \hat{\vec \iota} ^b$|⁠) after subtracting the astrometry of TSB, respectively. The centre of the segment is determined by the R.A. and decl. relative to the TSB, the slope is the ratio of proper motion offsets in the decl. and R.A directions, and the length is equal to the proper motion offset multiplied by the time span of Gaia DR2 or DR3. The fitted astrometry is determined through a five-parameter linear fit to the synthetic data, which is represented by colourful dots in the panels of the third column of Fig. 2. We also predict the position of ε Ind A b and ε Eridani b on 2024 January 1st. Their angular separations are 2.1 ± 0.1 and 0.64 ± 0.10 arcsec, and position angles are 243 ± 14 and 258 ± 30°, respectively. In addition, we present the fit to the five-parameter astrometry of Gaia DR2 and DR3 in Fig. 3. It is apparent that the deviation caused by ε Ind A b from the barycentric astrometry is more pronounced compared to that induced by ε Eridani b. The model fit consistently remains within 1σ of the five-parameter catalogue astrometry.

Figure 2.

Optimal orbital solutions for ε Ind A b (top) and ε Eridani b (bottom). The panels from the left to the right, respectively show the best fit to RV, Hipparcos IAD, and Gaiagost data, and the predicted planetary position on 2024 January 1st. The first column shows the binned RV data sets encoded by colours and shapes. For optimal visualization, each RV set is binned with a 100 d time window, while the un-binned RVs and residuals are shown in grey. The second column shows the post-fit Hipparcos abscissa residual projected along the R.A. and decl. directions. The multiple measurements for each epoch are binned to present the binned data encoded by colours. The darker colours encode earlier phases while the brighter colours represent later phases. The directions of the error bars indicate the along-scan direction. The third column shows the optimal fit to the Gaiagost data and the comparison between best-fitting and catalogue proper motions and positions at Gaia DR2 (GDR2) and GDR3 reference epochs. The shaded regions represent the uncertainty of position and proper motion. The dot and slope of each line, respectively represent the best-fitting position and proper motion offsets induced by the reflex motion at certain reference epoch. The fourth column shows the predicted planet position on 2024 January 1st. The 1σ contour line is shown to indicate the prediction uncertainty.

Figure 3.

Model fit to the five-parameter astrometry of GDR2 and GDR3 for ε Ind A b (left) and ε Eridani b (right). In both cases, the barycentric astrometry is subtracted from the five-parameter solutions for both the data (represented by a dot with error bar) and the model prediction (represented by a boxplot). The subscripts of the labels on the x-axes indicate the Gaia DR number.

The inferred parameters are presented in Table 1, B1, and B2. In the table, the following orbital parameters are directly inferred through posterior sampling: the orbital period (P), RV semiamplitude (K), eccentricity (e), argument of periastron (ω) of the stellar reflex motion,⁵ inclination (I), longitude of ascending node (Ω), and mean anomaly at the minimum epoch of RV data (M₀). The table also presents the astrometric offsets, which need to be subtracted from the catalogue astrometry of Gaia DR3 to derive the TSB astrometry. The bottom three rows present the derived parameters, including the semimajor axis of the planet-star binary orbit (a), planet mass (m_p), and the epoch when a planet crosses through its periastron (T_p).

Table 1.

Parameters for ε Ind A b and ε Eridani b.

Parametera	Unit	Meaning	ε Ind A b	ε Eridani b	Priorc	Minimum	Maximum
P	d	Orbital period	|$15676.48_{-1492.59}^{+2331.54}$|	|$2688.60_{-16.51}^{+16.17}$|	Log-Uniform	−1	16
K	m s−1	RV semiamplitude	|$31.69_{-4.25}^{+4.55}$|	|$9.98_{-0.38}^{+0.43}$|	Uniform	10−6	106
e	–	Eccentricity	|$0.42_{-0.04}^{+0.04}$|	|$0.26_{-0.04}^{+0.04}$|	Uniform	0	1
ωb	deg	Argument of periapsis	|$98.73_{-11.46}^{+9.46}$|	|$166.48_{-6.66}^{+6.63}$|	Uniform	0	2π
I	deg	Inclination	|$84.41_{-9.94}^{+9.36}$|	|$130.60_{-12.62}^{+9.53}$|	CosI-Uniform	−1	1
Ω	deg	Longitude of ascending node	|$243.38_{-13.41}^{+14.36}$|	|$206.07_{-17.48}^{+15.14}$|	Uniform	0	2π
M0	deg	Mean anomaly at the reference epoch	|$127.64_{-24.09}^{+31.05}$|	|$352.80_{-10.00}^{+10.80}$|	Uniform	0	2π
Δα	mas	α offset	|$-7.51_{-1.72}^{+1.59}$|	|$0.26_{-0.20}^{+0.19}$|	Uniform	−106	106
Δδ	mas	δ offset	|$-3.83_{-2.16}^{+2.40}$|	|$-0.04_{-0.20}^{+0.18}$|	Uniform	−106	106
Δϖ	mas	ϖ offset	|$-0.59_{-0.33}^{+0.28}$|	|$-0.31_{-0.15}^{+0.14}$|	Uniform	−106	106
|$\Delta \mu ^{*}_\alpha$|	mas yr−1	|$\mu ^{*}_\alpha$| offset	|$1.00_{-0.08}^{+0.11}$|	|$0.36_{-0.02}^{+0.02}$|	Uniform	−106	106
Δμδ	mas yr−1	μδ offset	|$0.95_{-0.19}^{+0.21}$|	|$0.85_{-0.01}^{+0.01}$|	Uniform	−106	106
P	yr	Orbital period	|$42.92_{-4.09}^{+6.38}$|	|$7.36_{-0.05}^{+0.04}$|	–	–	–
a	au	Semimajor axis	|$11.08_{-0.74}^{+1.07}$|	|$3.53_{-0.06}^{+0.06}$|	–	–	–
mp	MJup	Planet mass	|$2.96_{-0.38}^{+0.41}$|	|$0.76_{-0.11}^{+0.14}$|	–	–	–
Tp − 2400000	JD	Periapsis epoch	|$43293.68_{-1455.27}^{+1202.96}$|	|$44411.54_{-81.95}^{+76.60}$|	–	–	–

Notes. ^aThe first 12 rows show parameters that are inferred directly through Markov Chain Monte Carlo (MCMC) posterior sampling, while the last five rows show the parameters derived from the directly sampled parameters. The semimajor axis a and planet mass m_p are derived from the orbital parameters by adopting a stellar mass of 0.81 ± 0.04 M_⊙ for ε Eridani, and 0.74 ± 0.04 M_⊙ for ε Ind A, from the Gaia Final Luminosity Age Mass Estimator (Creevey et al. 2023).

^bThis is the argument of periastron of stellar reflex motion and ω + π is the argument of periastron of planetary orbit.

^cThe rightest three columns show the prior distribution, and the corresponding minimum and maximum values for a parameter. ‘Log-Uniform’ is the logarithmic uniform distribution, and ‘CosI-Uniform’ is the uniform distribution over cos I.

As the nearest cold Jupiters, ε Ind A b and ε Eridani b have been intensively studied in the past. Here we compare our solutions with previous ones. Based on combined analyses of RV, Hipparcos and Gaia DR2, Feng et al. (2019c) determine a mass of 3.25|$_{-0.65}^{+0.39}$| M_Jup , a period of |$45.20^{+5.74}_{-4.77}$| yr, and an eccentricity of |$0.26_{-0.03}^{+0.07}$| for ε Ind A b. Recently, Philipot et al. (2023) analyse the RV, Hipparcos and Gaia early data release 3 (EDR3; Gaia Collaboration 2023) data using the htof package (Brandt et al. 2021b), and estimate a mass of 3.0 ± 0.1 M_Jup , a period of |$29.93_{-0.62}^{+0.73}$| yr, and an eccentricity of 0.48 ± 0.01. In addition to the Gaia EDR3 (equivalent to DR3 for five-parameter solutions) data used by Philipot et al. (2023), we use Gaia DR2 to constrain the orbit. We find a mass of |$2.96_{-0.38}^{+0.41}$| M_Jup , a period of |$42.92^{+6.38}_{-4.09}$| yr, and an eccentricity of |$0.42_{-0.04}^{+0.04}$|⁠. The orbital period estimated in this work and in Feng et al. (2019c) is significantly longer than the one given by Philipot et al. (2023), probably due to the following reasons: (1) we optimize all parameters simultaneously instead of marginalizing some of them;⁶ (2) as suggested by Zechmeister et al. (2013), we subtract the perspective acceleration (about 1.8 m s⁻¹ yr⁻¹) from both the CES LC and VLC sets; (3) we model the time-correlated RV noise as well as instrument-dependent jitters (Feng et al. 2019c); (4) we use astrometric data from both Gaia DR2 and DR3. A detailed comparison of our methodology with others is given by Feng et al. (2021).

Although the existence of a cold Jupiter around ε Eridani was disputed due to consideration of stellar activity cycles in RV analyses (Hatzes et al. 2000; Anglada-Escudé & Butler 2012), ε Eridani b has been gradually confirmed by combined analyses of RV, astrometry, and imaging data (Benedict et al. 2017; Mawet et al. 2019; Llop-Sayson et al. 2021; Benedict 2022; Roettenbacher et al. 2022). Based on RV constraint and direct imaging upper limit of ε Eridani b, Mawet et al. (2019) determine a mass of |$0.78_{-0.12}^{+0.38}$| M_Jup , a period of 7.37 ± 0.07 yr, and an eccentricity of |$0.07_{-0.05}^{+0.06}$|⁠. Llop-Sayson et al. (2021) analyse RV, absolute astrometry from Hipparcos and Gaia DR2, and imaging data, and estimate a mass of |$0.66_{-0.09}^{+0.12}$| M_Jup , a period of |$2671_{-23}^{+17}$| d, and an eccentricity of |$0.055_{-0.039}^{+0.067}$|⁠. With additional RV data from EXPRES, Roettenbacher et al. (2022) get a solution consistent with that given by Mawet et al. (2019). Using only the astrometric data obtained by the Fine Guidance Sensor 1r of HST, Benedict (2022) determines a mass of |$0.63_{-0.04}^{+0.12}$| M_Jup , a period of 2775 ± 5 d, and an eccentricity of 0.16 ± 0.01. In this work, the combined analyses of RV, Hipparcos, Gaia DR2 and DR3 astrometry determine a mass of |$0.76_{-0.11}^{+0.14}$| M_Jup , a period of |$2688.60_{-16.51}^{+16.17}$| d, and an eccentricity of 0.26 ± 0.04. Modelling the time-correlated RV noise using the MA model, our solution gives an eccentricity significantly higher than the value given by Mawet et al. (2019), Llop-Sayson et al. (2021), and Roettenbacher et al. (2022), who use Gaussian processes to model time-correlated RV noise. While Gaussian process interprets the extra RV variations at the peaks and troughs of the Keplerian signal as time-correlated noise (see the bottom left panel of Fig. 2 in this paper, and fig. 2 in Roettenbacher et al. 2022), we interpret them as part of the signal because these excessive variations always amplify the peaks and troughs at the corresponding epochs with a similar periodicity. In previous studies where Gaussian process is not employed in RV modelling, high eccentricities have typically been estimated, such as 0.608 ± 0.041 (Hatzes et al. 2000), 0.25 ± 0.23 (Butler et al. 2006), and 0.40 ± 0.11 (Anglada-Escudé & Butler 2012). Furthermore, astrometry-only analyses have also resulted in high-eccentricity solutions, as demonstrated by previous studies such as e = 0.702 ± 0.039 (Benedict et al. 2006) and e = 0.16 ± 0.01 (Benedict 2022). Therefore, it is plausible to conclude that the orbital eccentricity of ε Eridani b is significantly higher than zero, as indicated by the findings of this study. In our study, the inclination is |$130.60_{-12.62}^{+9.53}$| °. This is equivalent to an inclination of |$49.40_{-9.53}^{+12.62}$| °, which is calculated in the absence of knowledge regarding the longitude of the ascending node. This value is largely consistent with the debris disc inclination ranging from 17 to 34 °(MacGregor et al. 2015; Booth et al. 2017), while Llop-Sayson et al. (2021) determines a significantly high inclination of |$78.81_{-22.41}^{+29.34}$| °.

5 CONCLUSION

Using gost to emulate Gaia epoch data, we analyse both Gaia DR2 and DR3 data in combination with Hipparcos and RV data to constrain the orbits of the nearest Jupiters. For ε Ind A b and ε Eridani b, the orbital periods are |$42.92_{-4.09}^{+6.38}$| and |$7.36_{-0.05}^{+0.04}$| yr, the eccentricities are |$0.42_{-0.04}^{+0.04}$| and |$0.26_{-0.04}^{+0.04}$|⁠, and the masses are |$2.96_{-0.38}^{+0.41}$| and |$0.76_{-0.11}^{+0.14}$| M_Jup , respectively. It is the first time that both Gaia DR2 and DR3 catalogue data are modelled and analysed by emulating the Gaia epoch data with gost. Compared with previous studies, our approach avoids approximating instantaneous astrometry by catalogue astrometry at the reference epoch, and thus enable robust constraint of orbits with period comparable or shorter than the DR2 and DR3 observation time span.

The orbital period of ε Ind A b in this work is much longer than the value given by Philipot et al. (2023), but is consistent with the solution based on combined analyses of RV and Gaia DR2 (Feng et al. 2019c). While we estimate an orbital eccentricity of ε Eridani b much higher than value given by Mawet et al. (2019), Llop-Sayson et al. (2021) and Roettenbacher et al. (2022), who use Gaussian process to model time-correlated noise, our solution is largely consistent with previous studies (Hatzes et al. 2000; Benedict et al. 2006; Anglada-Escudé & Butler 2012; Benedict 2022), which don’t use Gaussian process to model stellar activity. It is possible that Gaussian process may interpret part of the Keplerian signal as time-correlated noise (Feng et al. 2016; Ribas et al. 2018).

The combined RV and astrometry model presented in this work shares similarities with the successful model used to detect and confirm numerous sub-stellar and planetary companions in previous studies (Feng et al. 2019c, 2021, 2022). The orbital solution obtained for ε Ind A using our new method aligns closely with the results determined in Feng et al. (2019c), which highlights the reliability of our approach in detecting cold and massive companions. The upcoming imaging observations of ε Eridani and ε Ind A with JWST will provide further validation of our method, particularly in utilizing multiple Gaia DRs to identify cold Jupiters. While future Gaia DR4 will release epoch data, the modelling framework developed in our study offers a general methodology for incorporating both catalogue data and epoch data to constrain orbital parameters. This approach can be extended to combined analyses of astrometric data from various sources, including Gaia, photographic plates (Cerny et al. 2021), Tycho-2 (Høg et al. 2000), and future space missions such as the China Space Station Telescope (Fu et al. 2023) and the Nancy Roman Space Telescope (Yahalomi et al. 2023).

ACKNOWLEDGEMENTS

This work is supported by Shanghai Jiao Tong University 2030 Initiative. We would like to extend our sincere appreciation to the Scientific Editor and the anonymous referee for their valuable comments and feedback on our manuscript. We also thank Xianyu Tan for helpful discussion about atmosphere models of our targets. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has also made use of the services of the portal exoplanet.eu of The Extrasolar Planets Encyclopaedia, the ESO Science Archive Facility, NASA’s Astrophysics Data System Bibliographic Service, and the SIMBAD data base, operated at CDS, Strasbourg, France. We also acknowledge the many years of technical support from the UCO/Lick staff for the commissioning and operation of the APF facility atop Mt. Hamilton. All analyses were performed using r Statistical Software (v4.0.0; R Core Team 2020). This paper is partly based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes: 072.C-0488,072.C-0513,073.C-0784,074.C-0012,076.C-0878,077.C-0530,078.C-0833,079.C-0681, and 192.C-0852,60.A-9036.

DATA AVAILABILITY

The new RV data are available in the appendix, while the Gaia and Hipaprcos data are publicly available.

Footnotes

1

gost website: https://gaia.esac.esa.int/gost/ and gost user manual: https://gaia.esac.esa.int/gost/docs/gost_software_user_manual.pdf.

2

Note that μ_α is defined as |$\dot{\alpha }\cos \delta$| and is equivalent to |$\mu ^{*}_\alpha$|⁠.

3

It is actually the motion of system photocentre rather than the mass centre of the target star. Considering that the companion is far smaller than its host star in this work, the two centres are almost identical.

4

Although htof also simulates Gaia epoch data using gost, it does not utilize multiple Gaia DRs.

5

Note that the argument of periastron for the planetary orbit is ω_p = ω + π.

6

As an example, orvara performs marginalization of the RV offsets before fitting the orbital parameters. It assumes that there is no correlation between the offset and the orbital period.

APPENDIX A: POSTERIOR DISTRIBUTION OF ORBITAL PARAMETERS

Figure A1.

Corner plot showing the 1D and 2D posterior distribution of the orbital parameters for ε Ind A b. In each histogram, the blue lines indicate the 1σ confidence intervals, while the red line represents the median of the posterior distribution. The contour lines depict the 1σ and 2σ confidence intervals. Grey dots represent the posterior samples for each pair of parameters. It is important to note that ω_b denotes the argument of periastron for the stellar reflex motion, while the argument of periastron for the planetary orbit, ω_p, is equal to ω_b + π. Furthermore, M_0b corresponds to the mean anomaly at JD 2447047.96844.

Figure A2.

Similar to Fig. A1, but for ε Eridani b. M_0b is the mean anomaly at JD 2448930.56223.

APPENDIX B: OTHER MODEL PARAMETERS

Table B1.

Other parameters for ε Ind A b.

Parametera	Unit	Meaning	ε Ind A b	Prior	Minimum	Maximum
γHARPSpost	m s−1	RV offset for HARPSpost	|$-28.17_{-3.31}^{+2.97}$|	Uniform	−106	106
JHARPSpost	m s−1	RV jitter for HARPSpost	|$1.11_{-0.08}^{+0.09}$|	Uniform	0	106
|$w_{1}^{\rm HARPSpost}$|	–	Amplitude of component 1 of MA(1) for HARPSpost	|$0.97_{-0.02}^{+0.02}$|	Uniform	−1	1
lnτHARPSpost	–	Logarithmic time-scale of MA(1) for HARPSpost	|$1.01_{-0.21}^{+0.21}$|	Uniform	−12	12
γHARPSpre	m s−1	RV offset for HARPSpre	|$-25.62_{-3.33}^{+3.01}$|	Uniform	−106	106
JHARPSpre	m s−1	RV jitter for HARPSpre	|$1.00_{-0.01}^{+0.01}$|	Uniform	0	106
|$w_{1}^{\rm HARPSpre}$|	–	Amplitude of component 1 of MA(5) for HARPSpre	|$0.40_{-0.01}^{+0.01}$|	Uniform	−1	1
|$w_{2}^{\rm HARPSpre}$|	–	Amplitude of component 2 of MA(5) for HARPSpre	|$0.25_{-0.02}^{+0.01}$|	Uniform	−1	1
|$w_{3}^{\rm HARPSpre}$|	–	Amplitude of component 3 of MA(5) for HARPSpre	|$0.15_{-0.02}^{+0.01}$|	Uniform	−1	1
|$w_{4}^{\rm HARPSpre}$|	–	Amplitude of component 4 of MA(5) for HARPSpre	|$0.06_{-0.01}^{+0.01}$|	Uniform	−1	1
|$w_{5}^{\rm HARPSpre}$|	–	Amplitude of component 5 of MA(5) for HARPSpre	|$0.12_{-0.01}^{+0.01}$|	Uniform	−1	1
lnτHARPSpre	–	Logarithmic time-scale of MA(5) for HARPSpre	|$0.16_{-0.12}^{+0.11}$|	Uniform	−12	12
γLC	m s−1	RV offset for LC	|$-5.00_{-4.31}^{+3.98}$|	Uniform	−106	106
JLC	m s−1	RV jitter for LC	|$4.71_{-2.07}^{+1.74}$|	Uniform	0	106
γUVES	m s−1	RV offset for UVES	|$-17.88_{-3.36}^{+2.96}$|	Uniform	−106	106
JUVES	m s−1	RV jitter for UVES	|$0.61_{-0.09}^{+0.09}$|	Uniform	0	106
γVLC	m s−1	RV offset for VLC	|$-14.95_{-3.58}^{+3.13}$|	Uniform	−106	106
JVLC	m s−1	RV jitter for VLC	|$0.94_{-0.60}^{+0.92}$|	Uniform	0	106
lnJhip	–	Logarithmic jitter for hip	|$0.20_{-0.13}^{+0.13}$|	Uniform	−12	12
lnJgaia	–	Logarithmic jitter for gaia	|$2.82_{-0.69}^{+0.79}$|	Uniform	−12	12

Note. ^aThe MA model of order q or MA (q) is Parametrized by the amplitudes of q components and the logarithmic time-scale lnτ, where τ is in units of days. The MA model is introduced by Tuomi et al. (2013) and is frequently used in red-noise modelling for RV analyses (e.g. Tuomi & Anglada-Escudé 2013 and Feng et al. 2016). The superscripts of MA parameters and offsets represent the names of data sets. The definition of Gaia and Hipparcos jitter (lnJ_gaia and lnJ_hip) can be found in Feng et al. (2019c).

Table B2.

Other parameters for ε Eridani b.

Parameter	Unit	Meaning	ε Eridani b	Prior	Minimum	Maximum
γAPFp	m s−1	RV offset for APFp	|$0.68_{-0.37}^{+0.25}$|	Uniform	−106	106
JAPFp	m s−1	RV jitter for APFp	|$7.24_{-0.34}^{+0.35}$|	Uniform	0	106
γAPFh	m s−1	RV offset for APFh	|$5.74_{-0.31}^{+0.34}$|	Uniform	−106	106
JAPFh	m s−1	RV jitter for APFh	|$4.93_{-0.27}^{+0.28}$|	Uniform	0	106
|$w_{1}^{\rm APFh}$|	–	Amplitude of component 1 of MA(-Inf) for APFh	|$0.87_{-0.06}^{+0.05}$|	Uniform	−1	1
lnτAPFh	–	Logarithmic time-scale of MA(-Inf) for APFh	|$0.06_{-0.33}^{+0.24}$|	Uniform	−12	12
γEXPRES	m s−1	RV offset for EXPRES	|$-6.96_{-0.21}^{+0.24}$|	Uniform	−106	106
JEXPRES	m s−1	RV jitter for EXPRES	|$2.87_{-0.14}^{+0.15}$|	Uniform	0	106
|$w_{1}^{\rm EXPRES}$|	–	Amplitude of component 1 of MA(2) for EXPRES	|$0.87_{-0.07}^{+0.07}$|	Uniform	−1	1
|$w_{2}^{\rm EXPRES}$|	–	Amplitude of component 2 of MA(2) for EXPRES	|$0.14_{-0.08}^{+0.08}$|	Uniform	NA	NA
lnτEXPRES	–	Logarithmic time-scale of MA(2) for EXPRES	|$1.19_{-0.20}^{+0.23}$|	Uniform	−12	12
γHARPN	m s−1	RV offset for HARPN	|$7.83_{-0.27}^{+0.37}$|	Uniform	−106	106
JHARPN	m s−1	RV jitter for HARPN	|$13.57_{-0.55}^{+0.51}$|	Uniform	0	106
γHARPSpost	m s−1	RV offset for HARPSpost	|$-2.56_{-0.30}^{+0.21}$|	Uniform	−106	106
JHARPSpost	m s−1	RV jitter for HARPSpost	|$4.87_{-0.28}^{+0.30}$|	Uniform	0	106
γHARPSpre	m s−1	RV offset for HARPSpre	|$-12.08_{-0.37}^{+0.36}$|	Uniform	−106	106
JHARPSpre	m s−1	RV jitter for HARPSpre	|$1.47_{-0.05}^{+0.05}$|	Uniform	0	106
|$w_{1}^{\rm HARPSpre}$|	–	Amplitude of component 1 of MA(3) for HARPSpre	|$0.78_{-0.04}^{+0.04}$|	Uniform	−1	1
|$w_{2}^{\rm HARPSpre}$|	–	Amplitude of component 2 of MA(3) for HARPSpre	|$0.17_{-0.05}^{+0.05}$|	Uniform	−1	1
|$w_{3}^{\rm HARPSpre}$|	–	Amplitude of component 3 of MA(3) for HARPSpre	|$0.05_{-0.04}^{+0.04}$|	Uniform	−1	1
lnτHARPSpre	–	Logarithmic time-scale of MA(3) for HARPSpre	|$0.10_{-0.20}^{+0.20}$|	Uniform	−12	12
γKECK	m s−1	RV offset for KECK	|$5.05_{-0.31}^{+0.34}$|	Uniform	−106	106
JKECK	m s−1	RV jitter for KECK	|$6.86_{-0.21}^{+0.22}$|	Uniform	0	106
|$w_{1}^{\rm KECK}$|	–	Amplitude of component 1 of MA(1) for KECK	|$0.95_{-0.07}^{+0.03}$|	Uniform	−1	1
lnτKECK	–	Logarithmic time-scale of MA(1) for KECK	|$1.56_{-0.19}^{+0.20}$|	Uniform	−12	12
γKECKj	m s−1	RV offset for KECKj	|$-0.22_{-0.19}^{+0.20}$|	Uniform	−106	106
JKECKj	m s−1	RV jitter for KECKj	|$6.43_{-0.33}^{+0.37}$|	Uniform	0	106
γLC	m s−1	RV offset for LC	|$5.17_{-0.55}^{+0.43}$|	Uniform	−106	106
JLC	m s−1	RV jitter for LC	|$5.27_{-0.33}^{+0.33}$|	Uniform	0	106
|$w_{1}^{\rm LC}$|	–	Amplitude of component 1 of MA(1) for LC	|$0.75_{-0.11}^{+0.11}$|	Uniform	−1	1
lnτLC	–	Logarithmic time-scale of MA(1) for LC	|$2.65_{-0.27}^{+0.29}$|	Uniform	−12	12
γLICK13	m s−1	RV offset for LICK13	|$5.67_{-0.42}^{+0.56}$|	Uniform	−106	106
JLICK13	m s−1	RV jitter for LICK13	|$4.88_{-0.34}^{+0.48}$|	Uniform	0	106
γLICK6	m s−1	RV offset for LICK6	|$0.38_{-0.22}^{+0.19}$|	Uniform	−106	106
JLICK6	m s−1	RV jitter for LICK6	|$8.14_{-0.38}^{+0.60}$|	Uniform	0	106
γLICK8	m s−1	RV offset for LICK8	|$-0.42_{-0.13}^{+0.10}$|	Uniform	−106	106
JLICK8	m s−1	RV jitter for LICK8	|$9.22_{-0.28}^{+0.28}$|	Uniform	0	106
|$w_{1}^{\rm LICK8}$|	–	Amplitude of component 1 of MA(2) for LICK8	|$0.62_{-0.10}^{+0.11}$|	Uniform	−1	1
|$w_{2}^{\rm LICK8}$|	–	Amplitude of component 2 of MA(2) for LICK8	|$0.30_{-0.13}^{+0.13}$|	Uniform	−1	1
lnτLICK8	–	Logarithmic time-scale of MA(2) for LICK8	|$2.48_{-0.33}^{+0.31}$|	Uniform	−12	12
γVLC	m s−1	RV offset for VLC	|$3.15_{-0.37}^{+0.33}$|	Uniform	−106	106
JVLC	m s−1	RV jitter for VLC	|$1.41_{-0.55}^{+0.40}$|	Uniform	0	106
|$w_{1}^{\rm VLC}$|	–	Amplitude of component 1 of MA(1) for VLC	|$0.80_{-0.12}^{+0.10}$|	Uniform	−1	1
lnτVLC	–	Logarithmic time-scale of MA(1) for VLC	|$3.98_{-0.43}^{+0.31}$|	Uniform	−12	12
lnJhip	–	Logarithmic jitter for hip	|$-5.61_{-0.40}^{+0.38}$|	Uniform	−12	12
lnJgaia	–	Logarithmic jitter for gaia	|$2.00_{-0.38}^{+0.36}$|	Uniform	−12	12

APPENDIX C: NEW APF RVS FOR ε Eridani

Table C1.

APF data for ε Eridani.

BJD	RV [m s−1]	RV error [m s−1]	S-index
2456582.93034	26.64	2.73	0.524
2456597.91368	6.40	2.36	0.528
2456606.68427	16.52	0.75	0.531
2456608.10376	4.69	0.78	0.530
2456610.76250	16.04	1.18	0.512
2456618.88476	−2.11	0.78	0.530
2456624.72004	4.20	1.11	0.519
2456626.81421	24.46	0.75	0.521
2456628.72976	24.14	0.70	0.540
2456631.42746	−2.26	0.88	0.502
2456632.80921	14.46	0.62	0.523
2456644.75696	8.20	2.30	0.522
2456647.81171	14.44	0.63	0.535
2456648.59184	12.62	1.10	0.538
2456662.63738	9.77	0.73	0.536
2456663.75415	10.43	1.11	0.531
2456667.52792	18.00	0.78	0.535
2456671.68695	19.96	1.05	0.604
2456675.75647	7.84	1.12	0.519
2456679.83732	17.70	1.05	0.529
2456682.56608	17.80	0.82	0.550
2456689.76638	26.34	0.75	0.500
2456875.02028	7.12	2.18	0.501
2456894.88054	8.28	1.30	0.470
2456901.06193	9.95	1.54	0.479
2456909.10279	−4.71	1.21	0.476
2456922.07953	12.25	2.13	0.461
2456935.94021	−2.43	1.27	0.479
2456937.92403	−0.55	1.35	0.468
2456950.03798	3.82	1.44	0.472
2456985.64755	−1.80	2.28	0.441
2456988.63095	5.93	1.29	0.478
2456999.76434	8.84	1.37	0.459
2457015.72916	−2.17	1.10	0.465
2457026.78021	−1.44	1.34	0.464
2457058.45996	−3.69	1.89	0.435
2457234.08236	7.73	1.39	0.525
2457245.86234	−4.19	1.41	0.519
2457249.93007	−3.94	1.31	0.500
2457253.11257	5.63	1.33	0.511
2457257.15719	−1.02	1.15	0.506
2457258.94437	−12.69	1.23	0.517
2457261.02221	−2.76	1.32	0.501
2457262.94505	−7.81	1.36	0.496
2457265.95783	9.67	1.24	0.516
2457275.01304	−1.91	1.23	0.515
2457283.96368	1.88	1.29	0.507
2457287.02735	−1.11	1.35	0.524
2457290.95635	3.19	1.42	0.534
2457305.83659	−5.63	1.23	0.515
2457308.90844	13.30	1.26	0.534
2457318.83435	8.72	1.26	0.557
2457321.79157	6.64	1.36	0.540
2457325.84352	2.87	1.41	0.543
2457331.10764	9.90	1.36	0.552
2457332.78237	9.64	1.25	0.558
2457334.82998	5.22	1.30	0.548
2457337.78910	5.41	1.59	0.545
2457340.95644	−1.99	1.27	0.553
2457347.86896	4.10	1.29	0.556
2457348.77993	4.65	1.27	0.556
2457350.72611	5.83	1.20	0.558
2457354.70613	−0.88	1.65	0.548
2457361.64656	17.26	1.43	0.549