Constraining the Orbit and Mass of epsilon Eridani b with Radial Velocities, Hipparcos IAD-Gaia DR2 Astrometry, and Multiepoch Vortex Coronagraphy Upper Limits

Eridani has been targeted by five multiyear RV planet searches over the past 30 yr. Most of the resulting RV measurements are presented and discussed in Mawet et al. (2019). Since the publication of that paper, we have obtained 76 additional spectra of Eridani with Keck/HIRES (Vogt et al. 1994), and 172 spectra with the Levy Spectrograph on the Automated Planet Finder (APF) telescope (Radovan et al. 2014; Vogt et al. 2014). These new measurements are given in Table 5. As in Mawet et al. (2019), new HIRES observations were taken using the standard iodine cell configuration used for the California Planet Survey (Howard et al. 2010). Observations were taken through either the B5 (087×35) or C2 (087×140) decker, yielding R ≃ 55,000 spectral resolution. Spectra had a median of 293,000 counts per exposure, corresponding to a median extracted flux of 67,000 counts S/N = 260) at 550 nm. The same template spectrum described in Mawet et al. (2019), taken in 2010 using the B3 decker, was used to derive RVs from all new spectra obtained for this work.

New APF measurements were also taken using the standard configuration described in Mawet et al. (2019). Most APF measurements were acquired through the W decker (1'' x 3 "), with a spectral resolution of R ≃ 110,000. The median exposure time across observations was 26 s, yielding a per-pixel SNR of 120 at 550 nm. The same template spectrum described in Mawet et al. (2019), taken in 2014 using the N decker, was used to derive RVs from all new spectra obtained for this work. For our analysis, all exposures taken with a single instrument within 10 hr were binned.

We obtain the astrometry data for Eridani from the Hipparcos catalog Hip2 (van Leeuwen 2007), and the Gaia DR2 catalog (Gaia Collaboration et al. 2018), as well as the intermediate astrometry data (IAD) available from the Hipparcos mission.

2.2.1. Hipparcos Intermediate Astrometry Data

2.2.2. Gaia DR2

2.2.3. Gaia eDR3

Building upon the work presented by Mawet et al. 2019, we present additional Ms-band observations on the Eridani system acquired in 2019 with Keck/NIRC2 and its vector vortex coronagraph (Serabyn et al. 2017). The observations are summarized in Table 2. Except for the night of the 2019 December 08, the new data was acquired using the newly installed infrared pyramid wave front sensor (Bond et al. 2020) instead of the facility Shack–Hartmann wave front sensor. The alignment of the star behind the coronagraph is ensured at the 2 mas level by the quadrant analysis of coronagraphic images for tip-tilt sensing (Huby et al. 2015, 2017). The telescope was used in pupil tracking mode to allow speckle subtraction with angular differential imaging (Marois et al. 2008).

Table 2. Summary of Observations of Eridani with Keck/NIRC2 in Ms Band

Epoch	Tot. time	Exp. time	Coadds	# of exp.	Angle rot.	5σ	5σ	WFS
	(s)	(s)			(deg)	04	10
2017-01-09	6300	0.5	60	210	88.6	4.8e-05	1.1e-05	SH
2017-01-10	7800	0.5	60	260	100.2	4.9e-05	1.1e-05	SH
2017-01-11	4800	0.5	60	160	69.5	6.3e-05	1.0e-05	SH
Combined 2017	5.25h					3.3e-05	7.0e-06
2019-10-20	5190	0.25	120	173	84.9	5.1e-05	2.0e-05	Pyramid
2019-10-21	5850	0.25	120	195	82.0	3.9e-05	1.8e-05	Pyramid
2019-10-22	7170	0.25	120	239	89.3	5.4e-05	2.0e-05	Pyramid
2019-11-04	6660	0.25	120	222	68.2	6.3e-05	1.8e-05	Pyramid
2019-11-05	6510	0.25	120	217	68.0	4.6e-05	1.3e-05	Pyramid
2019-12-08	6510	0.3	100	217	103.7	1.0e-04	2.4e-05	SH
Combined 2019	10.5h					2.3e-05	7.9e-06

Note. From left to right: UT date of the observation, total exposure time, exposure time per integration, number of coadds, number of exposures, total paralactic angle rotation for angular differential imaging, 5σ sensitivity expressed as the planet-to-star flux ratio at 2 projected separation (04 and 10), and finally the wave front sensor (WFS) used (Shack–Hartmann or Pyramid).

Download table as:  ASCIITypeset image

The data was corrected for bad pixels, flat-fielded, sky-subtracted (accounting for detector illumination), and registered following the method detailed in Xuan et al. (2018). Similar to Mawet et al. (2019), the stellar point-spread function (PSF) was subtracted using principal component analysis (PCA) combined with a matched filter (Ruffio et al. 2017) from the open-source Python package pyKLIP (Wang et al. 2015). After high-pass filtering the images with a Gaussian filter with a length scale of twice the PSF full width at half maximum (FWHM), the PSF FWHM is ∼10 pixels (∼01 ) in Ms band. Stellar speckles are subtracted using 30 principal components. For each science frame, the principal components are only calculated from images featuring an azimuthal displacement of the planet on the detector of more than 0.5 FWHM. The final sensitivity of each epoch is plotted in Figure 1 as a function of projected separation. The nightly performance in 2019 was lower than in 2017 resulting in a similar combined sensitivity for each year despite the longer combined exposure time. A weighted mean is used to combine different epochs together. As a way to correct for any correlation between frames, the noise of the combined images is then normalized following the standard procedure of dividing by the standard deviation calculated in concentric annuli. With the exception of the 2019 December 08 observation, the epochs in 2017 and 2019 are close enough in time for the planet to not move significantly compared to the PSF size with a full width at half maximum of ∼10 pixels.

Zoom In Zoom Out Reset image size

Mawet et al. (2019) performed a thorough series of tests to evaluate the possibility that Eridani b is an artifact of stellar activity, finding that the ∼7 yr orbital period is distinct from periods and harmonics of the periodicities in the S_HK activity indicator timeseries. Our aim is not to recapitulate their analysis, but to update their orbital solution using the additional data obtained since the paper's publication, which spans approximately half of one orbital period of Eridani b.

Mawet et al. (2019) identified three peaks in a Lomb–Scargle periodogram of the RVs that rose above the 1% eFAP: one at the putative planet period of 7.3 yr, one at 2.9 yr, and one at 11 day. Applying RVSearch (Rosenthal et al. 2021) to the full data set, we recover this structure of peaks. Like Mawet et al. (2019), we also identify two major periods in the S_HK time series for both HIRES and APF, which coincide with the 11 day and 2.9 yr periods in the RVs. We interpret these as the signatures of rotationally-modulated stellar activity and a long-term activity cycle, respectively. To investigate the effects of these signals on the physical parameters of planet b, we performed two separate RV-orbit fits using RadVel (Fulton et al. 2018) to try different priors on the Gaussian-process (GP) timescale for modeling the stellar activity. In each of these fits, we assume a one-planet orbital solution, parameterized as $\sqrt{{e}_{b}}\cos {\omega }_{b}$, $\sqrt{{e}_{b}}\sin {\omega }_{b}$, T_conj,b, P_b , K_b . We also included RV offset (γ) and white noise (σ) parameters for each instrument in the fit, treating the four Lick velocity data sets independently to account for instrumental upgrades as in Mawet et al. (2019). Finally, we allowed an RV trend $\dot{\gamma }$.

In the first fit, we included a GP noise model to account for the impact of rotationally-modulated magnetic activity on the RVs (Rajpaul et al. 2015). We used a quasiperiodic kernel, following (Mawet et al. 2019). This kernel has hyperparameters η₂, the exponential decay timescale (analagous to the lifetime of active regions on the stellar disk), η₃, the stellar rotation period, η₄, which controls the number of local maxima in the RVs per rotation period, and η₁, the amplitude of the GP mean function, which we treated as independent for each instrument data set. Following López-Morales et al. (2016), we fixed η₄ to 0.5, which allows approximately two local maxima per rotation period. In this first fit, we allowed η₂ and η₃ to vary in the range (0 to 100 day). We calculated a Markov-chain representation of the posterior using emcee (Foreman-Mackey et al. 2013), We visually inspected the chains to ensure appropriate burn in and production periods. In total, the chain contained 450,080 samples. The resulting orbital and nuisance parameters are given in Table 3. Both the orbital parameters and the GP hyperparameters are well constrained; in particular, the marginalized posterior over the rotation period η₃ is Gaussian about the expected period of 11 day. The data allow a trend, although the value is consistent with no trend at the 1σ level, which allows us to conclude that there is no evidence in the current data for a RV trend. One nonintuitive feature of this fit is a preference for extremely small values (10⁻⁶ m s⁻¹) of white noise jitter for the second Lick data set. Even when we ran a fit requiring that all jitter values be at least 0.5 m s⁻¹, the posterior peaked at this lower bound. This may be evidence that much of the noise in this particular data set is correlated, and therefore well modeled by the GP noise model. It could also indicate that the reported observational uncertainties are overestimated for this data set. Whatever the reason, there are fewer than 10 measurements that are affected by the value of this jitter parameter, and neither the orbital parameters nor the GP timescale parameters are affected by its particular value.

Table 3. RV-fit MCMC Posteriors

Parameter	Credible Interval	Maximum Likelihood	Units
Modified MCMC Step Parameters
Pb	${2671}_{-23}^{+17}$	2661	days
Tconjb	${{\rm{2,460,017}}}_{-32}^{+76}$	2460023	JD
Tperib	${{\rm{2,460,054}}}_{-690}^{+680}$	2460235	JD
eb	${0.055}_{-0.039}^{+0.067}$	0.046
ωb	${57.3}_{-154.7}^{+80.2}$	2.1	°
Kb	${10.34}_{-0.93}^{+0.95}$	10.33	m s−1
Other Parameters
γlick4	$-{1.0}_{-2.6}^{+2.7}$	−0.9	m s−1
γlick3	10.0 ± 4.9	9.9	m s−1
γlick2	${7.5}_{-5.7}^{+5.6}$	7.6	m s−1
γlick1	${4.4}_{-6.1}^{+6.2}$	4.5	m s−1
${\gamma }_{{\mathrm{hires}}_{j}}$	2 ± 1	2	m s−1
γharps	${{\rm{16,442.5}}}_{-3.1}^{+3.2}$	16,442.4	m s−1
γces	${{\rm{16,446.6}}}_{-5.5}^{+5.7}$	16,446.6	m s−1
γapf	−0.9 ± 1.3	−0.9	m s−1
$\dot{\gamma }$	$-{0.00026}_{-0.0006}^{+0.00063}$	−0.00026	m s−1 d−1
$\ddot{\gamma }$	≡ 0.0	≡ 0.0	m s−1 d−2
σlick4	${5.17}_{-0.95}^{+1.1}$	5	m s−1
σlick3	${5.6}_{-2.3}^{+2.1}$	5.3	m s−1
σlick2	${2.3}_{-1.4}^{+3.6}$	0.5	m s−1
σlick1	${7.2}_{-3.7}^{+3.5}$	7.4	m s−1
${\sigma }_{\mathrm{hires}{\text{}}j}$	${2.36}_{-0.41}^{+0.46}$	2.26	m s−1
σharps	${4.8}_{-2.7}^{+2.2}$	4.3	m s−1
σces	${8.1}_{-4.5}^{+3.8}$	7.3	m s−1
σapf	${3.64}_{-0.55}^{+0.62}$	3.51	m s−1
η1,apf	${7.82}_{-0.92}^{+0.97}$	7.71	m s−1
${\eta }_{1,\mathrm{hires}{\text{}}j}$	${6.92}_{-0.59}^{+0.64}$	6.78	m s−1
η1,ces	${7.5}_{-4.7}^{+4.1}$	7.0	m s−1
η1,harps	${5.7}_{-3.0}^{+2.3}$	5.3	m s−1
η1,lick,fischer	${8.7}_{-1.4}^{+1.3}$	8.0
η2	${37.6}_{-5.4}^{+6.4}$	36.4	days
η3	${11.68}_{-0.13}^{+0.14}$	11.66	days
η4	≡ 0.5	≡ 0.5

Note. 436,160 links saved

Reference epoch for γ, $\dot{\gamma }$, $\ddot{\gamma }$: 2,457,454.642028.

Download table as:  ASCIITypeset image

To investigate the impact of the long-term activity cycle on the marginalized orbital parameter posteriors, we performed another fit identical to the one described above, except we required that η₂ and η₃ vary between 1 yr and the ∼30 yr observation baseline. The marginalized posteriors for η₂ and η₃ were broad, with power across the entire allowed space, although the period parameter η₃ showed local maxima at both ∼1100 and ∼2000 day. The marginalized 1 day posteriors for both of these parameters did not vary with those of any of the orbital parameters, allowing us to conclude that the long-term activity cycle, while somewhat present in the RVs, did not significantly affect the derived values of the orbital parameters. We therefore adopt the rotation-only GP fit described above.

Our derived orbital parameters, accounting for the effect of rotationally-modulated stellar noise, are very similar to those of Mawet et al. (2019; see their Table 3). We derive an orbital period of ${2671}_{-23}^{+17}$ days, a slightly reduced median semiamplitude of 10.3 m s⁻¹, and a low median eccentricity of 0.067. We show the series of RV measurements, the residuals to the fit and the phase-folded RV curve in Figure 2.

Zoom In Zoom Out Reset image size

We use the RV posterior distributions presented in Section 3.1 that we obtain from the RadVel orbit fit as the prior probabilities for the model fit to the astrometry and direct imaging data. The set of priors from the RV fit consists of: a, e, ω, τ, and M sin i. We use the same orbital parameter convention as in Blunt et al. (2020). However, since the fit to the astrometry and direct imaging data will obtain a separate distribution for both the mass and inclination, we draw a set of correlated values for i and M from the M sin i distribution. We assume a sine distribution as a prior probability for i, which corresponds to an unconstrained prior for i.

Such a set of parameters presents a complex covariance that is practically impossible to reproduce if trying to represent these distributions in a parametric way. Instead, we choose to use a kernel density estimator (KDE) to translate the posteriors to priors. A KDE smooths a probability density by convolving a Gaussian kernel with the discrete points of the Markov-chain Monte Carlo (MCMC). It allows one to preserve the covariance information contained in the original posteriors while allowing for a reliable reproduction of the shapes in the distributions. We use the scipy.stats package implementation of a KDE, which allows for a customized value of the KDE bandwidth parameter. The choice of the bandwidth is critical to maintain the fidelity of the RV fit for the fit of the astrometry and direct imaging data. We explain the process to compute the optimal bandwidth for our data in Appendix C.

The combined astrometric and imaging model consists of thirteen parameters: the six orbital parameters (a, e, i, ω, Ω, τ), the mass of the companion M_b , the total mass of the system M_tot, the parallax, the position of the barycenter of the system at an arbitrary epoch (we choose the Hipparcos epoch at 1992.25, α_1992.25, and δ_1992.25), and the PM vector of the barycenter μ = [μ_α , μ_δ ].

For the astrometric model, we use a similar approach to Rosa et al. (2019a) and Nielsen et al. (2020) while also fitting for the parallax and M_tot. This framework consists on deriving from the model a displacement of the photocenter from the barycenter induced by the presence of a companion. The orbital parameters and the relative mass are used to compute this relative motion of the photocenter about the barycenter. Then, to compare to the astrometry data, we obtain absolute astrometry quantities propagating from a reference R.A./decl. position with the PM of the barycenter, and adding the relative displacement. For this reason we add to our model the reference position α_1992.25 and δ_1992.25, and the PM of the barycenter μ. We assume that the brightness of the planet is negligible compared to the host star (below 10 ppb reflected light in the visible); consequently making the photocenter and the star positions coincide.

We compute the goodness of fit to the astrometric data by deriving two distinct terms for both the Hipparcos and Gaia data: ${\chi }_{\mathrm{astrom}}^{2}={\chi }_{H}^{2}+{\chi }_{G}^{2}$. For the Hipparcos IAD, we calculate the expected position of the barycenter at all IAD epochs (α_H and δ_H ) by propagating α_1992.25 and δ_1992.25 with μ. The displacement of the photocenter with respect to the barycenter (Δα_H and Δδ_H ) is computed using the expected orbit from the model orbital parameters, and the relative mass. We compare this values to the corresponding α_IAD and δ_IAD that we calculate from the 1D scans from IAD. We compute ${\chi }_{H}^{2}$ in the same way as in Nielsen et al. (2020).

Similarly for the Gaia data, we propagate the reference position of the barycenter, α_1992.25 and δ_1992.25, to the Gaia epoch α_G and δ_G , and we compute the displacement of the photocenter with respect to the barycenter (Δα_G and Δδ_G ) with the orbital parameters and the relative mass. We compare this derived quantity to the Gaia DR2 positional data, i.e., the R.A./decl. from the Gaia DR2 catalog α_G,GDR2 and δ_G,GDR2. We can fit to the Gaia DR2 PM data by adding an instantaneous PM disturbance to the reference PM calculated by deriving an average linear velocity of the star due to its orbital motion over a few epochs around 2015.5. The ${\chi }_{G}^{2}$ associated with the Gaia data is then calculated in the same way as described in Rosa et al. (2019a).

The analysis of the Hipparcos–Gaia eDR3 acceleration data is done in the way described in Kervella et al. (2019).

The part of the likelihood function associated with the direct imaging data described in Section 2.3 is computed based on the method described in Mawet et al. (2019). The logarithm of the direct imaging likelihood (${ \mathcal P }$) (Ruffio et al. 2018) for a single epoch can be written as:

$$\begin{eqnarray}&&\mathrm{log}{ \mathcal P }(d| I,x)=-\displaystyle \frac{1}{2{\sigma }_{x}^{2}}({I}^{2}-2I{\tilde{I}}_{x})\end{eqnarray} \tag{ 1 }$$

where I is the planet flux corresponding to the planet mass in the orbital model, ${\tilde{I}}_{x}$ is the estimated flux from the data at the location x, and σ_x is the uncertainty of this estimate. Individual epochs are not combined in the main analysis of this work. The direct imaging epochs are assumed to be independent such that their log likelihoods are simply added together. While this assumption is not perfect, it should not significantly affect the upper limit as the framework marginalizes over the spatial direction which will factor in any brighter speckles. To compute I from the planet mass, we use the COND evolutionary model (Baraffe et al. 2003) and we adopt an upper bound age of 800 Myr for the system (Janson et al. 2015). To compare with this age, we also run model fits with an age of 400 Myr, corresponding to the lower bound. ${\tilde{I}}_{x}$ is the flux measured in the image where the planet is predicted to be based on a given set of orbital parameters.

The fit to the astrometry and direct imaging data using the RadVel posteriors as priors is performed using the orbitize! ²⁰ package (Blunt et al. 2020). We implement a MCMC (Foreman-Mackey et al. 2013) analysis to fit for the six orbital parameters and the mass of planet b. MCMC has been extensively used to do orbit fitting, e.g., first by Hou et al. (2012) and more recently by Blunt et al. (2019), Nowak et al. (2020), Hinkley et al. (2021), and Wang et al. (2021). As discussed in Section 3.2 we also fit for the mass of Eridani, as well as for the astrometric parameters α and δ, and μ_α and μ_δ . Therefore we fit for thirteen parameters: six orbital parameters, the parallax, the photocenter position vector and the proper motion vector, and the masses of the star and planet.

The knowledge of Eridani's mass is accounted with a Gaussian prior 0.82 ± 0.02 M_⊙ (Baines & Armstrong 2011). The priors for the astrometric parameters are set to uniform distributions centered at the Hipparcos values with broad ranges of ±20σ to allow for deviations of the data from the nominal PM and photocenter position. The RV fits do not provide any constraints on the longitude of the ascending node, Ω, for which we set a uniform prior distribution. As discussed in Section 3.2, although the RV fits do not constrain the inclination, i, we assume a sine distribution to draw correlated values for the mass, M_b , and the inclination. The fit is performed with 550 walkers, and 18,000 steps per walker.

The MCMC fits converge to yield the orbital parameters and mass posteriors shown in Table 4. Two fits are presented for both ends of the current bounds on the system's age, i.e., 400 and 800 Myr. A corner plot showing some of these and their correlations, for the 800 Myr fit, is shown in Figure 3. A complete corner plot is shown in Appendix E.

Zoom In Zoom Out Reset image size

Table 4. MCMC Posteriors for Different Sets of Data and Assumptions for the Age of the System

	RV and Direct Imaging	RV, Astrometry,		RV and Astrometry	RV, Astrometry (eDR3 acc.),
		and Direct Imaging			and Direct Imaging
	800 Myr	400 Myr	800 Myr	No Age Assumption	800 Myr
Mb (MJup)	${0.74}_{-0.15}^{+0.37}$	${0.66}_{-0.10}^{+0.11}$	${0.66}_{-0.09}^{+0.12}$	${0.65}_{-0.09}^{+0.10}$	${0.64}_{-0.09}^{+0.09}$
ab (au)	${3.53}_{-0.04}^{+0.04}$	${3.53}_{-0.04}^{+0.03}$	${3.53}_{-0.04}^{+0.03}$	${3.53}_{-0.04}^{+0.04}$	${3.52}_{-0.04}^{+0.04}$
e	${0.07}_{-0.05}^{+0.07}$	${0.07}_{-0.05}^{+0.07}$	${0.07}_{-0.05}^{+0.07}$	${0.07}_{-0.05}^{+0.07}$	${0.07}_{-0.05}^{+0.07}$
ω (° )	$-29.{01}_{-111.95}^{+107.96}$	$-28.{37}_{-99.58}^{+75.90}$	$-19.{15}_{-94.80}^{+88.27}$	$-19.{54}_{-87.06}^{+97.09}$	$-29.{84}_{-115.85}^{+104.59}$
Ω (° )	$181.{27}_{-125.02}^{+124.60}$	$195.{06}_{-72.21}^{+127.20}$	$198.{18}_{-63.14}^{+127.29}$	$195.{08}_{-73.72}^{+122.64}$	$190.{06}_{-151.74}^{+108.72}$
i (° )	$89.{15}_{-43.16}^{+45.03}$	$75.{77}_{-19.94}^{+29.47}$	$78.{81}_{-22.41}^{+29.34}$	$80.{95}_{-21.39}^{+27.55}$	$89.{70}_{-25.08}^{+25.49}$
τperi	$0.{52}_{-0.33}^{+0.33}$	${0.42}_{-0.29}^{+0.34}$	$0.{35}_{-0.24}^{+0.33}$	$0.{37}_{-0.26}^{+0.32}$	$0.{52}_{-0.36}^{+0.34}$
ΔαH (mas)		$-{0.33}_{-0.41}^{+0.54}$	$-{0.25}_{-0.47}^{+0.53}$	$-{0.21}_{-0.48}^{+0.52}$
ΔδH (mas)		$0.{27}_{-0.46}^{+0.54}$	$0.{30}_{-0.47}^{+0.32}$	$0.{27}_{-0.48}^{+0.32}$
${{PM}}_{H}^{\alpha }$		$-{975.20}_{-0.02}^{+0.02}$	$-{975.20}_{-0.02}^{+0.02}$	$-{975.20}_{-0.02}^{+0.02}$
${{PM}}_{H}^{\delta }$		${19.95}_{-0.02}^{+0.02}$	${19.94}_{-0.02}^{+0.02}$	${19.95}_{-0.02}^{+0.02}$

Download table as:  ASCIITypeset image

In order to assess the constraining power of the direct imaging data, given that our data consists of nondetections, we performed an MCMC run to only fit for the astrometry data. Figure 4 show the posterior distribution compared to the prior distribution, i.e., the posteriors from the RV fit.

Zoom In Zoom Out Reset image size

The use of the astrometry data results in a significantly improved constraining of the inclination of the planet with respect to the results presented in Mawet et al. (2019). Starting with a sine distribution as the prior, which is centered at 90°, the walkers converge onto an inclination of $i=78.{81}_{-22\mathop{.}\limits^{^\circ }41}^{+29.34}$. The addition of the astrometry data also yields a different model for the mass of the companion; by using the RV and direct imaging data we obtain a distribution of ${M}_{b}={0.73}_{-0.13}^{+0.34}$ M_Jup (800 Myr), by adding the astrometry, ${M}_{b}={0.66}_{-0.09}^{+0.12}$ M_Jup (800 Myr). The astrometry data from Hipparcos and Gaia seem to favor a lower mass planet and more edge-on inclinations. We compute the most probable perturbation SMA: α_A = 0.89 mas, which is a factor of ∼2 away from the result reported by Benedict et al. (2006), i.e., 1.88 mas. More details on the perturbed orbit can be found in Appendix D.

We add the fit to the Hipparcos–Gaia eDR3 acceleration in the last column of Table 4. The results of this fit are largely consistent with the Hipparcos IAD–Gaia DR2 fits. Two main differences can be appreciated. (1) The inclination median is ∼90°; however, the posterior probabilities converge to a similar inclination with respect to our other fits, i.e., ∼75°, and to its supplementary angle, i.e., ∼105°. This is probably because the rotation information available from the IAD is not accessible when using the PM anomalies. (2) The upper mass bound (see Figure 4) is slightly lower, which brings the mass to a lower-mass solution.

The argument of periapsis, ω, reported in Mawet et al. (2019), is the stellar ω, which is the reason it is ∼180° off with respect to the result presented here.

The interaction between the debris disk and planet in the Eridani system was thoroughly discussed in Mawet et al. (2019). Here we build on this analysis with our updated orbital constraints. Eridani's debris disk is currently known to be composed of three rings (Backman et al. 2008): a main ring from 35–90 au, an inner belt at ∼3 au, and an intermediate belt at ∼20 au. Planet b sits between the two inner belts. Extensive characterization of the disk has been carried out over the years; see Backman et al. (2008), MacGregor et al. (2015), and Booth et al. (2017). In particular, Booth et al. (2017) analysis concluded the presence of a gap in emission in the circumstellar disk between ∼20 and ∼60 au, which would indicate the presence of one or more companions in this range of separations. The inclination of the main ring was constrained to 34° ± 2 (Booth et al. 2017).

The orbital parameters and mass of planet b can set constraints on the edges of the inner belts (Wisdom 1980). We follow the same analysis as Mawet et al. (2019) given that the eccentricity of the planet is expected to be well below 0.3, in which case the chaotic zone structure carved by the planet is independent of the eccentricity (Quillen & Faber 2006). With the results of our MCMC fit (see Table 4), in particular the SMA and relative mass, we expect there to be no particles from 2.97–4.29 au. The edges are slightly moved outwards with respect to Mawet et al. (2019) since the SMA posterior of the planet is now larger, it was then reported no particles from 2.7– 4.3 au. The mass posterior is lower, which reduces the width of the chaotic zone by ∼4%.

It was concluded in Mawet et al. (2019) that the inner belt and outer ring were most likely to be self-stirred, i.e., collisions between disk particles are driven by particle-to-particle gravitational interactions, as opposed to being stirred by the presence of the planet (see Section 5.3.2. and Figure 15 in Mawet et al. 2019). The timescale for the planet to stir and shape the main ring is ∼1 Gyr, which, combined with the distance between the two, makes their dynamical coupling probably not significant. As for the intermediate belt, the timescales of self stirring and planet stirring were computed to be similar in the range of separations of the belt. The new orbital parameters from our MCMC fit (see Table 4), namely a lower mass for the planet and a slightly higher SMA, contribute to larger planet-stirring timescales: an 18% increase in the planet-stirring timescale estimation with respect to Mawet et al. (2019). However, this does not rule out the contribution of planet induced stirring process for the intermediate belt.

The inclination of Eridani b is constrained to $78.{81}_{-22\mathop{.}\limits^{^\circ }41}^{+29.34}$ thanks to the inclusion of the astrometry data (see Section 3.5). This result indicates that the planet orbit is likely inclined with respect to the main ring, for which i = 34° ± 2°, which is ∼2σ away from the most probable inclination. The origin of such a mutual inclination is unknown but could possibly be due to the dynamical effects of a third body.

A mutual inclination could be causing a warping on the main ring. Indeed, the vertical warp in the β Pictoris inner disk is believed to have been produced by its mutual inclination with β Pictoris b (e.g., Dawson et al. 2011). Using a similar analysis as presented in Dawson et al. (2011) based on secular interactions, we find that in the case of Eridani the minimum mass of planet at 3 au to excite the inclination of dust particles in the ring at 70 au after 800 Myr is of ∼0.5 M_Jup. This indicates that planet b could be in the regime of starting to drive a warp in the main ring if it is indeed misaligned with the disk plane. A coplanar solution is still allowed by the data since an inclination of 32° is ∼1σ away from the most probable inclination of the planet. It is worth noting that Benedict et al. (2006) yielded a solution for the inclination of the companion of 301 ± 38.

Rosa et al. (2019b) recently published the prospects for constraining mass of 51 Eridani b with Gaia's final data release. They simulated sets of Gaia data with different astrometric error estimates, and found that the detection was possible only with optimistic mass and astrometric uncertainties. We performed a similar analysis computing the astrometric signal of Eridani b and comparing it to the sensitivity results of Figure 11 on Rosa et al. (2019b). We expect Eridani to have a similar astrometric error due to its brightness since the brightness of the star sets the uncertainty in the scans. We get an estimate for this at Lindegren et al. (2018): ∼50 μas.

We find that the amplitude of the astrometric signal for Eridani b is an order of magnitude stronger than that of 51 Eridani b. Indeed, the shorter distance to the system and the period the planet both favorable factors for a stronger astrometric signal. For the nominal Gaia mission span of 5 yr and for an astrometric uncertainty in the scans of 50 μ a, Eridani b is detectable at ∼1 M_Jup, which falls on the higher end of our mass posterior probability. However, for the extended mission span of 8 yr, for which 51 Eridani is only detectable for a high mass estimate, Eridani b is readily detectable even at ∼0.5 M_Jup, which is on the lower than the median mass of the posterior probability presented in this paper.

The final data release of Gaia's mission is a particularly exciting prospect for the exoplanet science field. Gaia final release will probably have access to constraining the dynamic mass of Eridani b. However, as the work presented in this paper aims at showing, it is by using this data combined with other observations that the best science is attainable.

The results presented in this paper are another example of the power of combining different methods to constrain a system's characteristics. By adding the astrometry data, we have identified new constraints for the inclination and longitude of the ascending node, both of which are inaccessible to an RV-orbit fit. The direct imaging data, although it being a nondetection, sets upper limits on the mass. The distribution of most likely planet positions shown in Figure 5 and how it prefers certain areas and fluxes is no coincidence; the MCMC walkers converge easier where the estimated flux from the coronagraph images is higher.

Zoom In Zoom Out Reset image size

Although the constraints on the position shown in Figure 5 are far from ideal, direct imaging planet hunters will take advantage of any position knowledge however small it may be. Indeed, data reduction techniques in direct imaging greatly benefit from a prior knowledge of the position of the object. For instance, in PCA(Soummer et al. 2012) based methods, a great deal of speckle subtraction power is gained by treating the data by patches; knowing where the planet is more likely to be reduces the computing time and allows the algorithm to focus on a constrained area of interest.

The synergies between RV, astrometry, and direct imaging data are currently being explored, and more work is being done in this direction (GRAVITY Collaboration et al. 2020).

The James Webb Space Telescope (JWST) will provide the community with unprecedented capabilities to do infrared exoplanet science. In this section we discuss the prospect for a detection of Eridani b with JWST's NIRSPEC and MIRI instruments. In Figure 6 we show the probability contour for the position of the planet at an epoch in JWST's Cycle 1.

Zoom In Zoom Out Reset image size

4.4.1. MIRI

The midinfrared instrument (MIRI) aboard JWST is equipped with a coronagraph, and is expected to reach 10⁻⁴ levels of raw contrast at λ = 10 μm (mode F1065C) with conventional star subtraction (Boccaletti et al. 2015). We performed a more sophisticated data reduction as was done for MIRI in Beichman et al. (2020) to get a more accurate representation of the expected performance of MIRI on Eridani. We use the IDL library MIRImSIM, ²¹ and the wave front error drift predictions presented in Perrin et al. (2018). We make use of as much diversity as possible from the jitter of the telescope to perform a reference star subtraction based on PCA (Soummer et al. 2012). A more detailed description of the data processing can be found in Beichman et al. (2020).

In Figure 7 we show some simulation results for MIRI observations. We find that, even for a perfect pointing accuracy of the telescope, MIRI would require ∼75 hr to reach SNR > 5.

Zoom In Zoom Out Reset image size

4.4.2. NIRSpec

Carrión-González et al. (2021) assessed the potential of detecting reflected light from a set of exoplanets in nearby systems with the Roman coronagraph instrument (CGI). In particular, Eridani stands out as a particularly notable, yet risky target; like most targets for the Roman CGI, it would require thousands of hours to get a detection. They conclude that Eridani b would be Roman-accessible at a probability of 57.99%, in the optimistic case, and 51.29% in the pessimistic case. However, they argue that the inclination "is the key factor affecting the detectability of this planet." Indeed, as it can be seen in their Figure C.5, a more edge-on orbit is more favorable to avoiding the outer working angle. Our result for the inclination, $i=78.{81}_{-23\mathop{.}\limits^{^\circ }41}^{+29.34}$, indicates that the orbit should be more favorable for detection, since Carrión-González et al. (2021) assume face on as the most probable orbit.

A recent publication by Pathak et al. (2021) presented new observations of Eridani at 10 μ m with the VLT/VISIR. They obtained comparable sensitivities to the Keck/NIRC2 results presented here. They claim that a sensitivity to 1 M_Jup can be attained with 70 hr of exposure time, assuming an age of 700 Myr and the current setup for the instrument. Unfortunately, the results presented here indicate that the mass is likely lower than 1 M_Jup. However, there are envisioned ways to upgrading VISIR which would improve its sensitivity at smaller separations (Kasper et al. 2019).

Although a formidable task for current ground-based observatories, imaging Eridani b should be much easier with future 30-meter class telescopes. High-contrast imaging at L, M, and N bands with METIS at the ELT will essentially be background limited at the angular separation of Eridani b (Carlomagno et al. 2020). METIS is expected to have access to Earth-like planets around α Centauri A with 5 hr of exposure time N band (Brandl et al. 2021), which would make Eridani b detectable in the order of minutes. Similarly, the TMT with its second generation instrument PSI is expected to reach 10⁻⁸ final contrast at 2 λ/D (Fitzgerald et al. 2019), well above the sensitivity needed to image Eridani b.