Hint of curvature in the orbital motion of the exoplanet 51 Eridani b using 3 years of VLT/SPHERE monitoring

TL;DR

This study uses three years of VLT/SPHERE data to suggest a potential curvature in the orbit of exoplanet 51 Eridani b, providing improved orbital constraints and hints of orbital dynamics.

Contribution

First direct astrometric monitoring over three years indicating possible orbital curvature and refined orbital parameters of 51 Eridani b.

Findings

Detected orbital motion with a decreasing separation of ~10 mas.

Suggested an orbital period of approximately 32 years.

Indicated possible alignment or misalignment within 18 degrees of the star-planet system.

Abstract

Context. The 51 Eridani system harbors a complex architecture with its primary star forming a hierarchical system with the binary GJ 3305AB at a projected separation of 2000 au, a giant planet orbiting the primary star at 13 au, and a low-mass debris disk around the primary star with possibly a cold component and a warm component inferred from the spectral energy distribution. Aims. We aim to better constrain the orbital parameters of the known giant planet. Methods. We monitored the system over three years from 2015 to 2018 with the VLT/SPHERE exoplanet imaging instrument. Results. We measure an orbital motion for the planet of ~130 mas with a slightly decreasing separation (~10 mas) and find a hint of curvature. This potential curvature is further supported at 3$\sigma$ significance when including literature GPI astrometry corrected for calibration systematics. Fits of the SPHERE and GPI data using three complementary approaches provide broadly similar results. The data suggest an orbital period of 32$^{+17}_{-9}$ yr (i.e. 12$^{+4}_{-2}$ au in semi-major axis), an inclination of 133$^{+14}_{-7}$ deg, an eccentricity of 0.45$^{+0.10}_{-0.15}$, and an argument of periastron passage of 87$^{+34}_{-30}$ deg [mod 180 deg]. The time at periastron passage and the longitude of node exhibit bimodal distributions because we do not detect yet if the planet is accelerating or decelerating along its orbit. Given the inclinations of the planet's orbit and of the stellar rotation axis (134-144 deg), we infer alignment or misalignment within 18 deg for the star-planet spin-orbit. Further astrometric monitoring in the next 3-4 years is required to confirm at a higher significance the curvature in the planet's motion, determine if the planet is accelerating or decelerating on its orbit, and further constrain its orbital parameters and the star-planet spin-orbit.

Tables

Table 1. Table 1: Observing log of SPHERE observations of 51 Eridani.

UT date	$ϵ$ (′′)	${\tau }_0$ (ms)	AM start/end	Mode	Bands	DIT (s) $\times$ Nfr	FoV rot. (∘)	SR
2015/09/25	0.5–1.0	4–9	1.10–1.09	IRDIFS_EXT	$YJH$+$K12$	16$\times$256	41.5	0.68–0.88
2015/09/26	0.7–1.3	6–12	1.10–1.09	IRDIFS	$YJ$+$BBH$	4(64)$\times$918(64)	42.6	0.66-0.90
2016/01/16	1.6–2.3	1	1.09–1.10	IRDIFS	$YJ$+$H23$	16(64)$\times$256(64)	41.8	0.63–0.86
2016/12/12	1.6–2.8	2	1.09–1.16	IRDIFS	$YJ$+$H23$	64$\times$54	25.2	0.55–0.62
2016/12/13	0.6–1.0	4–8	1.10–1.09	IRDIFS	$YJ$+$H23$	64$\times$72	44.4	0.78–0.92
2017/09/28	0.4–0.7	5–12	1.10–1.09	IRDIFS_EXT	$YJH$+$K12$	24(32)$\times$192(144)	44.1	0.85–0.91
2018/09/18	0.7–1.2	2–5	1.21–1.08	IRDIFS_EXT	$YJH$+$K12$	24(32)$\times$200(160)	38.5	0.64–0.87

Table 2. Table 2: SPHERE astrometry relative to the star of 51 Eridani b.

Epoch	Filter	$\rho$	PA	$\Delta$RA	$\Delta$Dec	Pixel scale	North correction angle
		(mas)	(∘)	(mas)	(mas)	(mas/pix)	(∘)
2015.74	$K1$	453.4$\pm$4.6	167.15$\pm$0.56	100.8$\pm$2.9	$-$442.0$\pm$3.6	12.267$\pm$0.009	$-$1.813$\pm$0.046
2015.74	$H$	453.9$\pm$16.3	166.1$\pm$2.0	108.7$\pm$8.8	$-$440.7$\pm$13.7	12.251$\pm$0.009	$-$1.813$\pm$0.046
2016.04	$H2$	456.7$\pm$6.9	165.50$\pm$0.84	114.3$\pm$4.5	$-$442.2$\pm$5.2	12.255$\pm$0.009	$-$1.82$\pm$0.06
2016.95	$H2$	453.6$\pm$5.7	160.30$\pm$0.72	152.9$\pm$3.4	$-$427.1$\pm$4.6	12.255$\pm$0.009	$-$1.808$\pm$0.043
2017.74	$K1$	449.0$\pm$2.9	155.67$\pm$0.38	185.0$\pm$2.0	$-$409.2$\pm$2.1	12.267$\pm$0.009	$-$1.735$\pm$0.043
2018.72	$K1$	443.3$\pm$4.2	150.23$\pm$0.55	220.2$\pm$2.8	$-$384.8$\pm$3.1	12.267$\pm$0.009	$-$1.796$\pm$0.068

Table 3. Table 3: Preliminary orbital parameters of 51 Eridani b

Parameter	Unit	Median	Lower	Upper	${\chi }_{\min }^2$
$P$	yr	32	23	49	30
$a$	au	12	10	16	12
$e$		0.45	0.30	0.55	0.43
$i$	∘	133	126	147	135
$\Omega$	∘	103	13	156	17
$\omega$	∘	87	57	121	69
$T_0$		2011	2006	2027	2006

Table 4. Table 4: SPHERE astrometry relative to the star of 51 Eridani b obtained with SDI+TLOCI and SDI+ANDROMEDA.

Epoch	Spectral band	$\rho$ (mas)	PA (∘)
SDI+TLOCI
2015.74	$K1$	453.4$\pm$4.4	167.15$\pm$0.55
2015.74	$H$	453.9$\pm$15.9	166.1$\pm$2.0
2016.04	$H2$	456.7$\pm$6.6	165.50$\pm$0.83
2016.95	$H2$	453.6$\pm$5.7	160.30$\pm$0.72
2017.74	$K1$	449.0$\pm$2.9	155.67$\pm$0.37
2018.72	$K1$	443.3$\pm$4.2	150.23$\pm$0.54
SDI+ANDROMEDA
2015.74	$K1$	448.6$\pm$1.4	167.45$\pm$0.06
2015.74	$H$	467.4$\pm$2.9	167.09$\pm$0.07
2016.04a𝑎aa𝑎afootnotemark: $a$	$H2$	–	–
2016.95	$H2$	456.1$\pm$1.6	160.06$\pm$0.06
2017.74	$K1$	447.9$\pm$1.3	155.80$\pm$0.04
2018.72	$K1$	439.0$\pm$1.2	150.09$\pm$0.03

Table 5. Table 5: Orbital parameters of 51 Eridani b derived from the MCMC analysis.

Parameter	Unit	Median	Lower	Upper	${\chi }_{\min }^2$
$P$	yr	30	23	46	53
$a$	au	12	10	15	17
$e$		0.49	0.36	0.58	0.39
$i$	∘	135	126	146	124
$\Omega$	∘	76	14	139	5
$\omega$	∘	85	54	111	121
$T_0$		2007	2005	2027	2009

Table 6. Table 6: Orbital parameters of 51 Eridani b derived from the OFTI analysis.

Parameter	Unit	Median	Lower	Upper	${\chi }_{\min }^2$
$P$	yr	31	23	47	34
$a$	au	12	10	16	13
$e$		0.47	0.33	0.57	0.39
$i$	∘	135	126	147	132
$\Omega$	∘	82	15	148	134
$\omega$	∘	86	59	118	101
$T_0$		2012	2006	2027	2026

Full text

*Aims. *We aim to better constrain the orbital parameters of the known giant planet.

*Methods. *We monitored the system over three years from 2015 to 2018 with the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) instrument at the Very Large Telescope (VLT).

Results. We measure an orbital motion for the planet of $\sim$130 mas with a slightly decreasing separation ($\sim$10 mas) and find a hint of curvature. This potential curvature is further supported at 3$\sigma$ significance when including literature Gemini Planet Imager (GPI) astrometry corrected for calibration systematics. Fits of the SPHERE and GPI data using three complementary approaches provide broadly similar results. The data suggest an orbital period of 32${}^{+17}_{-9}$ yr (i.e., 12${}^{+4}_{-2}$ au in semi-major axis), an inclination of 133${}^{+14}_{-7}$ deg, an eccentricity of 0.45${}^{+0.10}_{-0.15}$, and an argument of periastron passage of 87${}^{+34}_{-30}$ deg [mod 180∘]. The time at periastron passage and the longitude of node exhibit bimodal distributions because we do not yet detect whether the planet is accelerating or decelerating along its orbit. Given the inclinations of the orbit and of the stellar rotation axis (134–144*∘), we infer alignment or misalignment within 18∘* for the star–planet spin-orbit. Further astrometric monitoring in the next 3–4 years is required to confirm at a higher significance the curvature in the motion of the planet, determine if the planet is accelerating or decelerating on its orbit, and further constrain its orbital parameters and the star–planet spin-orbit.

1 Introduction

The giant planet 51 Eridani b is the first of its kind discovered in the Gemini Planet Imager (GPI) exoplanet imaging survey (Macintosh et al., 2015). The methane-rich planet is a bound companion to the young star 51Eridani, which is a member of the 24-Myr $\beta$ Pictoris moving group (Zuckerman et al., 2001; Torres et al., 2008; Bell et al., 2015). The star is located at 29.78$\pm$0.15 pc, inferred from the inverse of the parallax measured by Gaia (Gaia Collaboration et al., 2018). It is in good agreement with the value derived with an optimized approach (29.76$\pm$0.12 pc, Bailer-Jones et al., 2018). Our uncertainty of 0.15 pc includes in addition to the statistical error of 0.12 pc an uncertainty term of 0.1 mas to account for potential parallax systematics (https://www.cosmos.esa.int/web/gaia/dr2). The star forms a hierarchical system with the M-dwarf binary GJ3305AB with separation $\sim$10 au located at a projected separation of 2000 au (Feigelson et al., 2006; Montet et al., 2015). Simon & Schaefer (2011) measured a stellar radius of 1.63$\pm$0.03 $R_{\odot}$ with the Center for High Angular Resolution Astronomy (CHARA) interferometer and inferred a stellar mass of 1.75$\pm$0.05 $M_{\odot}$. The primary star also harbors a debris disk inferred from the spectral energy distribution (SED; Riviere-Marichalar et al., 2014; Patel et al., 2014). Riviere-Marichalar et al. (2014) estimated a low infrared (IR) fractional luminosity $L_{\rm{IR}}$/$L_{\odot}$ = 2.3$\times$10*-6* from Herschel photometry. Since their analysis is based on fitting a three-parameter model of a modified blackbody to three data points with excess IR emission at wavelengths $\geq$70 $\muup$m, the resulting value for the inner edge of the cold dust belt is largely uncertain with 82${}^{+677}_{-75}$ au. They also estimated an upper limit for the dust mass of 1.6$\times$10*-3* $M_{\oplus}$ and did not report gas detection ([OI], [CII]). Patel et al. (2014) observed the target with WISE as part of a survey for warm debris disks and inferred a warm disk with a temperature of 180 K (upper limit 344 K) and a radius of 5.5 au (lower limit 1.5 au) assuming the disk radiates as a blackbody. The 51 Eridani system could therefore harbor a two-belt debris disk architecture, a feature observed in other young systems with giant planets, such as HR 8799 (Marois et al., 2008, 2010; Su et al., 2009) and HD 95086 (Rameau et al., 2013b; Moór et al., 2013).

The planet 51 Eridani b has a projected separation of $\sim$13 au from the primary star. Macintosh et al. (2015) could not confirm the companionship with a proper motion test because of the very short time baseline of their GPI measurements ($\sim$1.5 months, between December 2014 and January 2015). Instead, the planetary nature hypothesis is based on the spectrum showing methane absorption. De Rosa et al. (2015) presented a new GPI astrometric epoch obtained in September 2015, confirming that the planet is gravitationally bound, and detected orbital motion. They also carried out a preliminary assessment of its orbital elements using Bayesian rejection sampling and Markov-chain Monte Carlo methods. Their analysis suggests most probable values with 1$\sigma$ error bars for the semi-major axis of 14${}^{+7}_{-3}$ au, for the period of 41${}^{+35}_{-12}$ yr, and for the inclination of 138${}^{+15}_{-13}$ deg. The other parameters are marginally constrained. The authors also noted that the orbital inclination of the planet is different from the inclination of the orbital plane of the binary GJ 3305AB ($i$ = 92.1$\pm$0.2*∘*, Montet et al., 2015), implying that they cannot be coplanar.

We present in this paper astrometric follow-up observations of 51 Eridani b obtained with the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) instrument (Beuzit et al., 2019) at the Very Large Telescope (VLT) as part of the SpHere INfrared survey for Exoplanets (SHINE; Chauvin et al., 2017). We describe the observations and the data reduction (Sect. 2). We then use the new astrometric data of the planet to analyze its orbital motion (Sect. 3). We subsequently fit the SPHERE astrometry in combination with GPI data to derive its orbital parameters (Sect. 4).

2 Observations and data analysis

We observed 51 Eridani eight times from September 2015 to September 2018 with the IRDIFS mode of SPHERE. In this mode, the near-infrared (NIR) camera IRDIS (Dohlen et al., 2008; Vigan et al., 2010) and integral field spectrograph IFS (Claudi et al., 2008) are operated in parallel, either in the $YJ$ bands for IFS and the $H23$ filter pair for IRDIS (standard IRDIFS mode) or in the $YJH$ bands for IFS and the $K12$ filter pair for IRDIS (IRDIFS_EXT mode). Four datasets were published in an analysis of the SED of the planet in Samland et al. (2017). Table 1 lists the published observations used for astrometry and the new observations. We only considered the IRDIS data in this work because the planet astrometry could be extracted from a higher number of datasets due to limitations in signal-to-noise ratio. The challenging contrast of the planet meant that it could be detected and its astrometry measured in six datasets only (Table 2).

For all sequences, an apodized pupil Lyot coronagraph (Carbillet et al., 2011; Martinez et al., 2009) was used. For calibrating the flux and the centering of the images, we acquired unsaturated non-coronagraphic images of the star (hereafter reference point-spread function or reference PSF) and coronagraphic images with four artificial crosswise replicas of the star (Langlois et al., 2013) at the beginning and end of the sequences. For all datasets obtained starting from December 2016, the science images were recorded with the stellar replicas simultaneously to minimize the frame centering uncertainties in the astrometric error budget. Night-time sky background frames were taken and additional daytime calibration performed following the standard procedure at ESO.

The data were reduced with the SPHERE Data Center pipeline (Delorme et al., 2017), which uses the Data Reduction and Handling software (v0.15.0, Pavlov et al., 2008) and custom routines. It corrects for the cosmetics and instrument distortion, registers the frames, and normalizes their flux. Subsequently, we sorted the frames using visual inspection to reject poor-quality frames (adaptive optics open loops, low-wind effect) and an automatic criterion to reject frames with low flux in the coronagraphic spot (semi-transparent mask). After this step, we were left with 77–97% of the frames, depending on the sequence. Finally, the data were analyzed with a consortium image processing pipeline (Galicher et al., 2018). Figure 1 shows the IRDIS images obtained for the best epochs with a two-step process333The 2015 September 26 dataset was obtained with the broad $H$-band filter, so it was processed with angular differential imaging only.: simultaneous spectral differential imaging (SDI; Racine et al., 1999) and angular differential imaging (Marois et al., 2006a) with the Template Locally Optimized Combination of Images algorithm (TLOCI; Marois et al., 2014).

For all epochs, the planet astrometry and photometry was measured in the SDI+TLOCI images using the fit of a model of planet image built from the reference PSF and processed with TLOCI (Galicher et al., 2018). The position and flux of the model of planet image was optimized to minimize the image residuals within a circular region of radius 1.5 full width at half maximum centered on the measured planet location. The values reported in Table 2 were calibrated following the methods in Maire et al. (2016). We also compared them with the astrometry extracted using SDI in combination with the ANgular DiffeRential Optimal Method Exoplanet Detection Algorithm (ANDROMEDA, Mugnier et al., 2009; Cantalloube et al., 2015) and found most values to agree within the TLOCI measurement uncertainties (Appendix A). We use the SDI+TLOCI astrometry for the astrometric and orbital analyses in the following sections, because TLOCI was tested and validated on a larger number of SPHERE datasets to retrieve the astrometry and photometry of detected companions (Galicher et al., 2018).

3 Orbital motion

The astrometry of the planet is given in Table 2. The data are represented in Fig. 2 with the GPI measurements reported by De Rosa et al. (2015), who also revised the astrometry published in Macintosh et al. (2015). The SPHERE data over three years confirm the orbital motion of the planet at a high significance: $\sim$119 mas at $\sim$30$\sigma$ in right ascension and $\sim$57 mas at $\sim$12$\sigma$ in declination. While there is a hint of a decrease in separation by $\sim$10 mas, the position angle clearly decreases at a rate of 5.7$\pm$0.2*∘/yr. The trend in position angle is similar to the trend seen in the GPI data (De Rosa et al., 2015). The position angle variation is not compatible with the expectations for a face-on circular orbit ($\sim$10∘*/yr, assuming a semi-major axis for the planet of 13 au and a stellar mass of 1.75 $M_{\odot}$), suggesting an inclined and/or eccentric orbit. The data also show signs of curvature, hinting at orbital inflexion (see below).

The GPI data obtained in December 2014 and January 2015 show a discrepant increasing trend in separation (in particular, the separations measured on two consecutive nights in January 2015 are not included within the measurement uncertainties of each other hence disagree at the 1-$\sigma$ level: 454.0 $\pm$ 6.4 mas on January 30 and 461.8 $\pm$ 7.1 mas on January 31). Macintosh et al. (2015) noted that the conditions for this last observation were average. We also note a small systematic offset in position angle between the SPHERE and GPI data using two measurements obtained close in time in September 2015 (the GPI point has PA = 166.5 $\pm$ 0.6*∘, which is smaller by 1.1 $\sigma$ than the position angle of the SPHERE point which was obtained more than three weeks later; we should also expect a further displacement of about -0.38 $\pm$ 0.02∘* in position angle due to the orbital motion of the planet over this elapsed time). This offset is likely due to systematic uncertainties related to differences in the astrometric calibration of the instruments (see Appendix B). It is well accounted for by an offset in the measured position angles of 1.0$\pm$0.2*∘*. The inclusion of the recalibration uncertainty in the GPI measurement uncertainties has negligible effects.

We show the recalibrated GPI data and the SPHERE measurements in Fig. 3, which represents the temporal evolution of the relative right ascension, relative declination, separation, and position angle. The curvature of the SPHERE data suggested from Fig. 2 is better seen in the top panels of Fig. 3. We also show in these panels second-order polynomial fits to all the SPHERE and GPI data except for the GPI data point taken on 2015 January 31 due to its discrepant separation with respect to the other GPI data points. The separation measured for this epoch is 461.8 mas whereas the other data points have separations smaller than 455 mas (upper limit of 460.4 mas at 1$\sigma$). Even when increasing the error bars on this data point to include the other GPI measurements, the LSMC orbital fit is affected with respect to a fit where this data point is excluded and shows a stronger paucity in low-eccentricity orbits. Second-order polynomial fits provide significantly better unreduced chi-square goodness-of-fit parameters (1.4 for $\Delta$RA vs. time and 1.0 for $\Delta$Dec vs. time) with respect to linear fits (7.4 and 12.2, respectively). We followed the approach of Konopacky et al. (2016a) to test if acceleration is detected (the measured acceleration plus its uncertainty at 3$\sigma$ shall stay negative). From the above-mentioned second-order polynomial fits, we estimated the cartesian components of the acceleration and converted them into radial and tangential components. The radial acceleration component is $-$4.03$\pm$1.34 mas.yr*-2*, which implies that acceleration is detected at the 3.0$\sigma$ level. New measurements should help to confirm the acceleration estimate, and measure it with better accuracy.

4 Orbital analysis

4.1 Determination of the orbital parameters

We assumed for the system the distance estimated from the Gaia parallax and a total mass of 1.75 $M_{\sun}$ (Simon & Schaefer, 2011). We also make the assumption that the GJ 3305AB binary does not dynamically disturb the orbit of the planet. Montet et al. (2015) showed that given the wide separation of the binary and the young age of the system, it appears unlikely that Kozai-Lidov oscillations (Kozai, 1962; Lidov, 1962) could have had the time to disturb the semi-major axis of the planet (typical timescale of 200 Myr for a perturber on a circular orbit). Nevertheless, moderate changes in the inclination and eccentricity cannot be excluded (Fabrycky & Tremaine, 2007). Quicker dynamical effects (secular precession due to, e.g., an unseen inner companion) could suppress Kozai-Lidov oscillations. On the other hand, assuming the criterion in Holman & Wiegert (1999), we could also expect that the current orbit remains unchanged despite the wide binary.

We first used a least-squares Monte Carlo (LSMC) procedure to fit the SPHERE and GPI astrometry (Esposito et al., 2013; Maire et al., 2015, Appendix C). We also performed complementary analyses using a Markov-chain Monte Carlo (MCMC) procedure (Chauvin et al., 2012, Appendix D) and the Bayesian rejection sampling approach OFTI (Blunt et al., 2017, Appendix E).

We used the LSMC and OFTI methods to analyze the effects of the inclusion of the GPI astrometric points on the parameter distributions. Adding the GPI data strengthens short-period and eccentric orbits. We also used these two approaches to test the effect of a different initial eccentricity distribution on the resulting eccentricity distribution, given that the analyses favor eccentric orbits over circular orbits. We used a distribution that gives more weight to low-eccentricity orbits, similar to the fit to the eccentricity distribution of radial-velocity planets in Nielsen et al. (2008). The resulting eccentricity distributions have similar shapes using both types of initial distributions. Finally, we checked that the uncertainty on the stellar mass ($<$3%) has negligible effects on the parameter distributions.

Further constraints could also be obtained with detection limits from archival high-contrast imaging data by rejecting orbits predicting overly large separations for the planet in the past. Unfortunately, the current detection limits are not deep enough to provide useful constraints (Heinze et al., 2010; Biller et al., 2013; Rameau et al., 2013a; Hagan et al., 2018; Stone et al., 2018).

4.2 Parameter intervals and correlations

The LSMC distributions and intervals of the orbital parameters are shown in Fig. 4 and Table 3. The shapes of the distributions and the parameter intervals are broadly similar to those obtained with the MCMC and OFTI approaches. The three $T_{0}$ distributions show two peaks around $\sim$2005 and $\sim$2025. The three $\Omega$ distributions display two broad peaks around $\sim$10*∘* and $\sim$130*∘. Nevertheless, we note some differences in the detailed shape of the distributions. The LSMC eccentricity distribution does not show a high-eccentricity tail (correlated with a long-period tail) as seen in the MCMC and OFTI distributions. We verified that this is due to the correction for the LSMC fitted orbits of the bias on the time at periastron passage for eccentric orbits (Konopacky et al., 2016a). The LSMC distribution for $i$ has a more pronounced peak toward values smaller than $\sim$135∘* with respect to the MCMC and OFTI distributions. This is caused by a larger number of low-eccentricity orbits with $\Omega$ around $\sim$160∘ in the LSMC fitted orbits. The LSMC distribution for $\Omega$ exhibits a deeper dip between the two broad peaks mentioned above with respect to the MCMC and OFTI distributions. This feature is related to a paucity of orbits with $\Omega$ around $\sim$65*∘* and low to moderate eccentricities (up to $\sim$0.4) in the LSMC fitted orbits. The MCMC distribution for $\omega$ appears sharper than the LSMC and OFTI distributions.

When comparing our results obtained with the OFTI approach with those of De Rosa et al. (2015), who employed a similar method, we note that most of the parameters are better defined in our analysis. The 68% interval for the longitude of node in De Rosa et al. (2015) was derived by wrapping the values to the range 30–120*∘, whereas we consider the full [0;180∘] range. These authors also wrapped their distribution of the time of periastron passage to the range 1995–1995+$P$, whereas our distributions extend to previous epochs. Nevertheless, we verified that applying a similar wrapping has a negligible effect on the derived $T_{0}$ ranges in our analysis. We derive a 68% interval for the inclination of 126–147∘, which is slightly better constrained with respect to the range of 125–153∘* in De Rosa et al. (2015). Thus, we confirm after De Rosa et al. (2015) that the orbital plane of the planet cannot be coplanar with the orbital plane of the wide-separated binary GJ 3305AB ($i$ = 92.1$\pm$0.2*∘, Montet et al., 2015). Our eccentricity interval is 0.33–0.57 at 68%, which suggests eccentric orbits. This feature is related to our more extended dataset with respect to De Rosa et al. (2015) (Appendix E). Nevertheless, orbits with low to moderate eccentricities are not formally excluded. Low-eccentricity orbits are characterized by parameters around $\Omega$ $\sim$ 160∘, $i$ $\sim$ 130∘*, and $P$ $\sim$ 35–40 yr. New data covering a larger fraction of the orbit of the planet are needed to verify any suggestion of large eccentricities.

Figure 5 shows a random sample of fitted orbits from the LSMC analysis. In particular, low-eccentricity orbits predict steeper decreases of the separation and position angle in the coming years. If new data taken in the next couple of years show a weak decrease of the separation and still follow the linear trend in position angle seen with the current data, this will rule out low-eccentricity orbits.

4.3 An unseen inner companion that could bias the eccentricity of the planet?

Finally, we used the methods in Pearce et al. (2014) to test the scenario of an unseen inner low-mass companion which could bias the eccentricity of 51 Eridani b toward large values due to the orbital motion that the unseen companion induces on the host star around the center of mass of the system. We considered the case where this putative inner companion lies on a circular orbit. For this, we used the period and eccentricity distributions in Fig. 4. Figure 6 shows the minimum mass of a putative inner companion as a function of the eccentricity of the planet. Such a companion would lie at an angular separation of $\sim$0.21*′′* ($\sim$6 au). By comparing these masses to the SPHERE/IRDIS mass limit of $\sim$3 $M_{\rm{J}}$ measured at this separation (Samland et al., 2017) according to the atmospheric and evolutionary models of Baraffe et al. (2015, 2003), we can conclude that if 51 Eridani b has a nonzero eccentricity, this eccentricity is genuine and does not result from an unseen low-mass inner companion. The companion could in theory lie at a larger separation but its mass would be larger, hence its detection would be even easier.

5 Conclusions

We presented VLT/SPHERE observations over three years of the young giant exoplanet 51 Eridani b to further characterize its orbital motion and parameters. The planet moved by $\sim$130 mas over this elapsed time with hints of orbital curvature and a decreasing trend in its separation to the star. We compared the results of three orbital fitting approaches based on LSMC, MCMC, and Bayesian rejection sampling and found similar distribution shapes for all parameters. With respect to the study of De Rosa et al. (2015), our orbital analysis based on a similar Bayesian rejection sampling approach provides narrower ranges for the orbital parameters. The time at periastron passage and the longitude of node exhibit bimodal distributions, the ambiguity being related to the non-detection of changes in the orbital speed of the planet. We derived an inclination range of 126–147*∘, which is slightly narrower than the 125–153∘* range derived by De Rosa et al. (2015). We note that the orbital inclination of the planet is compatible with an orbit lying in the stellar equatorial plane (the stellar rotation axis has an inclination of 134–144*∘; see Appendix F) or offset by less than 18∘*. Given that the star is expected to host a debris disk, this might suggest a coplanar planet–disk configuration and dynamical interactions. Further astrometric monitoring will help to solve for the ambiguities in the time at periastron passage and the longitude of node as well as refine the orbital inclination of the planet and the analysis of the system spin-orbit.

Our orbital analysis suggests an eccentric orbit for the planet with a 68% interval of 0.30–0.55. If the eccentricity of 51 Eridani b is indeed genuine, this may hint at dynamical interactions between the planet and another body in the system to produce such a large eccentricity. This putative additional body could be an unseen inner or outer planet, although we note that the current imaging detection limits are relatively deep ($>$4 $M_{\rm{J}}$ beyond 5 au and $>$2 $M_{\rm{J}}$ beyond 9 au assuming hot-start models; Samland et al., 2017). Another possibility would be gravitational perturbations from GJ 3305AB, such as Kozai-Lidov oscillations, but this scenario may face timescale issues because of the large separation of the binary and the system youth. Fabrycky & Tremaine (2007) predict that a close-in giant planet experiencing Kozai-Lidov oscillations from a distant binary companion to its host star will typically have an orbit that is misaligned with the stellar equatorial plane. Our analysis favors alignment or misalignment within 18*∘* for the spin-orbit of the 51 Eridani star–planet system. Interactions of the planet with the circumstellar disk could also be a possible mechanism. Although this kind of interaction is usually thought to dampen the eccentricity of a planet, simulations have shown that for massive giant planets (¿4–5 $M_{\rm{J}}$) lying near the disk, midplane (¡10*∘*) interactions with a protoplanetary disk increase their eccentricity (Papaloizou et al., 2001; Kley & Dirksen, 2006; Bitsch et al., 2013). The current mass estimate of 51 Eridani b is $\sim$2–4 $M_{\rm{J}}$ assuming hot-start models, but it could be as large as 12 $M_{\rm{J}}$ assuming warm-start models (Samland et al., 2017).

Further astrometric monitoring in the next 3–4 years will be critical to confirm at a higher significance the curvature in the motion of the planet, determine if it is accelerating or decelerating on its orbit, and further constrain its orbital parameters. It will also be critical for preparing future observations. If the angular separation of 51 Eridani b strongly decreases, this might prevent follow-up investigations near its periastron passage with SPHERE, GPI, and the James Webb Space Telescope to better constrain its orbital and atmospheric properties and leave such observations feasible with Extremely Large Telescope instruments only.

In addition to further orbital follow-up studies, resolved images of the host star debris disk will be valuable to determine if the planet orbits in the disk plane by providing the disk orientation and, if the data confirm the two disk belts inferred from the stellar SED, if it dynamically shapes these belts by providing their radial extent. Such information would also help to better characterize the orbital period and eccentricity of 51 Eridani b.

Acknowledgements.

We thank an anonymous referee for a constructive report that helped to improve the manuscript. The authors thank the ESO Paranal Staff for support in conducting the observations and Philippe Delorme, Eric Lagadec, and Nadège Meunier (SPHERE Data Center) for their help with the data reduction. We also thank Remko Stuik for discussions about the analysis of the MASCARA data and Sarah Blunt, Henry Ngo, Jason Wang, Isabel Angelo, Devin Cody, Robert De Rosa, and Logan Pearce for making the orbitize code publicly available. We acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU. This work has also been supported by a grant from the French Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). D.M. was supported by the South Carolina Space Grant Consortium. J.C. was supported by the Fulbright Colombia program as well as the South Carolina Space Grant Consortium. This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice) and Observatoire de Paris/LESIA (Paris). This research made use of the SIMBAD database and the VizieR Catalogue access tool, both operated at the CDS, Strasbourg, France. The original description of the VizieR service was published in Ochsenbein et al. (2000, A&AS 143, 23). This research has made use of NASA’s Astrophysics Data System Bibliographic Services. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF – Osservatorio di Padova (Italy), Observatoire astronomique de l’Université de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France), and ASTRON (Netherlands), in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from the CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland), and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programs as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012), and grant number 312430 for FP7 (2013–2016).

Appendix A Comparison of astrometric measurements from different algorithms

Appendix B Comparison of the SPHERE and GPI astrometry without GPI/SPHERE recalibration

Figure 7 shows further comparisons of all the SPHERE astrometry and the GPI astrometry reported in De Rosa et al. (2015) without applying any recalibration of the latter on the SPHERE data. While we do not see any clear GPI/SPHERE offset in the separations (Fig. 3), we note that the position angle of the GPI point taken on 2015 September 1 is 0.65*∘* (1.1 times the measurement uncertainty) smaller than the position angle of a SPHERE point taken more than three weeks later, on 2015 September 25. However, we expect an additional decrease in position angle of $\sim$0.38$\pm$0.02*∘* between the GPI and SPHERE epochs due to the orbital motion of the planet, meaning that the actual offset between these two measurements is possibly $\sim$1*∘* (1.7 times the measurement uncertainty).

We did not consider for the comparison the SPHERE point taken on 2015 September 26 because the position angle measured at this epoch has significantly larger uncertainties and deviates from a linear fit matching all the other SPHERE data points well (purple line in the bottom panel of Fig. 7). Contrary to the other SPHERE datasets, this dataset was not obtained with the dual-band imaging mode of IRDIS that allows for simultaneous imaging in two spectral bands in and out of a methane absorption band. Therefore, SDI could not be used in the image post-processing to attenuate fast quasi-static stellar speckles which are not attenuated with angular differential imaging, resulting in a poorer detection of the planet. Furthermore, we did not consider other GPI and SPHERE data points because they were taken at large time intervals. We note that the GPI position angle measurements exhibit a steeper slope with respect to the SPHERE data although the measured uncertainties are large ($-$6.7$\pm$1.3*∘) and include the slope value derived from the SPHERE data ($-$5.7$\pm$0.2∘*).

In order to further analyze potential position angle systematics between the SPHERE and GPI data, we reduced all the GPI H-band data of 51 Eridani available in the Gemini archive (eight datasets taken from December 2014 to November 2017) using the GPI data reduction pipeline v1.4.0 (Perrin et al. 2014, 2016), which applies an automatic correction for the north offset of $-$1.00$\pm$0.03*∘* measured by Konopacky et al. (2014). We then post-processed them using SDI+ANDROMEDA. The SDI step was necessary to enhance the S/N of the planet. For this, we selected spectral channels where the planet is not expected to show large fluxes due to strong methane absorption (Macintosh et al. 2015; Samland et al. 2017; Rajan et al. 2017). We were able to recover the planet in five of the datasets. The GPI position angle measurements shown in Fig. 8 display a negative slope ($-$5.9$\pm$0.3*∘) in agreement with the slope measured with the SPHERE data but are offset by 0.65$\pm$0.17∘* toward smaller values. We also note that our GPI measurements in common with De Rosa et al. (2015) are offset by $\sim$0.35$\pm$0.05*∘* toward larger values.

From these analyses, we applied a recalibration in position angle of 1.0$\pm$0.2*∘* to the GPI measurements in De Rosa et al. (2015) before fitting the SPHERE and GPI data. Due to the location of 51 Eridani b, the offset in position angle produces an offset mainly in relative right ascension as seen in the top-left panel of Figure 3.

The large GPI/SPHERE position angle offset that we found in our analysis is likely related to differences in the astrometric calibration of the instruments. It is currently unclear if this position angle offset should be considered systematically when combining SPHERE and GPI astrometry in orbital fits because it is not seen for other targets observed with both instruments and with published observations close in time (HD 95086, HR 2562, $\beta$ Pictoris, HR 8799; Rameau et al. 2016; Konopacky et al. 2016b; Wang et al. 2016, 2018; Chauvin et al. 2018; Maire et al. 2018; Lagrange et al. 2019; Zurlo et al. 2016). Further analysis is needed to conclude on this point but is considered to be beyond the scope of this paper.

Appendix C LSMC orbital fitting

We drew 1 000 000 random realizations of the astrometric measurements assuming Gaussian distributions around the nominal values. We then fit the six Campbell elements simultaneously using the downhill simplex AMOEBA procedure provided in the EXOFAST library (Eastman et al. 2013): orbital period $P$, eccentricity $e$, inclination $i$, longitude of node $\Omega$, argument of periastron passage $\omega$, and time at periastron passage $T_{0}$. Initial guesses of the parameters were computed assuming uniform distributions in log $P$, $e$, cos $i$, $\Omega$, $\omega$, and $T_{0}$. We considered no restricted ranges except for the period ($P$=10–1000 yr). We also included the correction for the bias on the eccentricity and time at periastron passage due to the small orbital arc covered by the data following the method in Konopacky et al. (2016a) ($\sim$34% of fitted orbits are rejected when applying this method, because of long periods, large eccentricities, and/or times at periastron passage close to the epochs of the data).

For the corner plot and the 68% intervals of the parameters shown in Sect. 4, we retained for the analysis all the derived solutions with $\chi_{\rm{red}}^{2}$ ¡ 2. The longitude of node and the argument of periastron passage are restrained in the [0;180*∘*] range to account for the ambiguity on the longitude of node inherent to the fitting of imaging data alone.

We note two broad peak features around $\sim$10*∘* and 130*∘* in the distribution of the longitude of node, which are associated with eccentricities of $\sim$0.2–0.5 and show correlations with the time at periastron passage. The peak around $\sim$10*∘* appears to produce more orbits with $T_{0}$ around $\sim$ 2005, whereas the peak around $\sim$130*∘* appears to produce more orbits with $T_{0}$ around $\sim$ 2025.

Appendix D MCMC orbital fitting

Here we provide the parameter distribution obtained using an MCMC approach (see details in Chauvin et al. 2012). We assumed uniform priors in log $P$, $e$, cos $i$, $\Omega+\omega$, $\omega-\Omega$, and $T_{0}$. Ten chains of orbital solutions were conducted in parallel, and we used the Gelman–Rubin statistics as convergence criterion (see the details in Ford 2006). We picked a random sample of 500 000 orbits from those chains following the convergence.

Appendix E OFTI orbital fitting

Here we provide the parameter distribution obtained using a custom Interactive Data Language (IDL) implementation of the Orbits For The Impatient (OFTI) approach described in Blunt et al. (2017). Briefly, we drew random orbits from uniform distributions in $e$, cos $i$, $\omega$, and $T_{0}$ and adjusted their semi-major axis and longitude of node by scaling and rotating the orbits to match one of the measured astrometric points. As explained in Blunt et al. (2017), the scale-and-rotate method to adjust the semi-major axis and longitude of node imposes uniform priors in log $P$ and $\Omega$. Subsequently, the $\chi^{2}$ probability of each orbit was computed assuming uncorrelated Gaussian errors before performing the rejection sampling test. To speed up the procedure in order to obtain a meaningful number of orbits (29 870), we applied the procedure at each iteration over 8 000 trial orbits simultaneously and we also restrained the prior ranges using the statistics of the first 100 accepted orbits.

We cross-checked our code with the OFTI procedure available as part of the Python orbitize package (Blunt et al. 2019). For verification, we also performed a fit using only GPI data points (2014 December 18, 2015 January 30, and 2015 September 1) from De Rosa et al. (2015) and assuming the same priors for the parameters; we found parameter distributions and intervals very similar to those obtained by these latter authors.

Appendix F Stellar rotation axis

Koen & Eyer (2002) estimated a rotation period of 0.65 d for 51 Eridani from Hipparcos photometric data without giving an uncertainty. This rotation period was used by Feigelson et al. (2006) with a stellar projected rotational velocity of $v$ sin $i_{\star}$= 71.8 $\pm$ 3.6 km s*-1* (Reiners & Schmitt 2003) and a stellar radius of 1.5 $R_{\odot}$ to estimate an inclination of 45*∘* for the rotation axis of the star.

In order to better constrain the stellar rotation axis and estimate in particular an uncertainty on this parameter, we reanalyzed the Hipparcos photometric data (Perryman et al. 1997; van Leeuwen et al. 1997). Figure 11 shows the results. Both the Lomb-Scargle periodogram (Scargle 1982) and the CLEAN periodogram (Roberts et al. 1987) show a peak at $P_{\star}$ = 0.65 $\pm$ 0.03 d. We also analyzed as a cross-check analysis archival data from the Multi-site All-Sky CAmeRA (MASCARA; Talens et al. 2017) and found a rotation period in good agreement with the value derived from the Hipparcos data ($P_{\star}$ = 0.66 d, Fig. 12). We considered only the Hipparcos results in the remainder of the analysis.

Using a $V$ magnitude $V$ = 5.20 mag, a distance $d$ = 29.78 $\pm$ 0.15 pc, a bolometric correction $BC_{V}$ = 0 mag, and an average effective temperature from the literature $T_{\star}$ = 7250 K, we infer a stellar radius $R_{\star}$ = 1.53 $\pm$ 0.04 $R_{\odot}$. Combining the rotation period, the stellar radius, and an average projected rotational velocity $v$ sin $i_{\star}$ = 83 $\pm$ 3 km s*-1* (estimated from an average of the measurements in Royer et al. 2007; Luck 2017, 84 km s*-1* and 81.2 km s*-1*, respectively), we infer an inclination of the stellar rotation axis with respect to the line of sight 39*∘* $<$ $i_{\star}$ $<$ 51*∘* or 129*∘* $<$ $i_{\star}$ $<$ 141*∘. Given the derived orbital inclination of the planet (126–147∘), this suggests alignment for the spin-orbit of the star-planet system. Using the stellar radius measurement of 1.63 $\pm$ 0.03 $R_{\odot}$ in Simon & Schaefer (2011) gives 36∘* $<$ $i_{\star}$ $<$ 46*∘* or 134*∘* $<$ $i_{\star}$ $<$ 144*∘*, which also suggests spin-orbit alignment.