1. Introduction
The first planets discovered orbiting Sun-like stars (e.g. Mayor & Queloz Reference Mayor and Queloz1995; Butler & Marcy Reference Butler and Marcy1996; Marcy & Butler Reference Marcy and Butler1996; Butler et al. Reference Butler, Marcy, Williams, Hauser and Shirts1997) were found using the radial velocity (RV) technique – observing the line-of-sight wobble of a star to infer the presence of an unseen planet.Footnote a In the three decades since, that technique has proven to be extremely powerful – only outperformed in pure discovery numbers by the transit method, which has the advantage of being able to study thousands or tens of thousands of stars at once.
At the time of writing, the RV technique has yielded 1 186 confirmed exoplanet discoveries.Footnote
b
In addition, the technique has been widely employed in the characterisation of exoplanets discovered through other methods (e.g. Mills et al. Reference Mills2019; Addison et al. Reference Addison2021; Akana Murphy et al. Reference Akana Murphy2023; Weiss et al. Reference Weiss2024; Stefánsson et al. Reference Stefánsson2025), since it offers a direct mechanism by which a planet’s mass can be measured. The more massive the planet, the larger the wobble it will induce in the motion of its host star. The RV method provides strong constraints on orbital periods, eccentricities, and companion minimum masses across a broad range of separations, and as a result remains an excellent tool for the characterisation of planetary systems. However, the RV method only provides a minimum mass (
$M\textrm {sin}i$
), with the true mass and orbital inclination to the line of sight remaining unknown. This can be solved by combining RV measurements with astrometry, which provides information on a star’s reflex motion on the plane of the sky, motion that is perpendicular to that detected by the RV method (Brown Reference Brown2015). This information therefore complements the RV measurements, allowing the direct determination of the inclination of the orbits relative to the plane of the sky, and from this, the true mass of the planets. The Gaia space telescope (Perryman et al. Reference Perryman, Hartman, Bakos and Lindegren2014) has already proven the effectiveness of this technique in determining the true mass of high-mass planets and brown dwarfs (e.g. Stevenson et al. Reference Stevenson, Haswell, Barnes and Barstow2023; Unger et al. Reference Unger2023; Wallace et al. Reference Wallace, Casey, Brown and Castro-Ginard2025) on periods of
$\sim$
1–5 yr.
While RV sensitivity generally decreases with increasing orbital period, the astrometric reflex signal increases approximately linearly with orbital separation. Long-period giant planets are therefore particularly well suited to astrometric characterisation. However, companions with orbital periods substantially exceeding the Gaia observing baseline can remain difficult to detect from epoch astrometry alone, due to the fact that the observing baseline doesn’t cover a full orbit. The orbital motion is therefore difficult to distinguish from the star’s proper motion (Lammers & Winn Reference Lammers and Winn2026). In such systems, the reflex motion instead manifests as a discrepancy between short-term and long-term proper motions. The Hipparcos–Gaia Catalogue of Accelerations (HGCA) (Brandt Reference Brandt2021) provides a homogeneous framework for exploiting this effect over a baseline of more than two decades. In recent years, this has been the widest application of Gaia astrometry to studies of exoplanets (e.g. Damasso et al. Reference Damasso2020; Venner et al. Reference Venner, Vanderburg and Pearce2021; Llop-Sayson et al. Reference Llop-Sayson2021; Barbato et al. Reference Barbato2025; Pinamonti et al. Reference Pinamonti2026).
HD 65216 is a nearby G5V star (Houk Reference Houk1978) hosting at least two giant planets identified through RV monitoring (Wittenmyer et al. Reference Wittenmyer2019). The system also contains two M7–8 stellar companions 253 AU away from the primary (Mugrauer et al. Reference Mugrauer, Seifahrt and Neuhäuser2007), forming a triple system. The study by Wittenmyer et al. (Reference Wittenmyer2019) determined the inner planet, HD 65216 b, has a minimum mass of 1.30
$\pm$
0.06 M
$_{\mathrm{J}}$
and a period of 578
$\pm$
1 days, while the outer companion, HD 65216 c, has a minimum mass of 2.03
$\pm$
0.11 M
$_{\mathrm{J}}$
and a period of 5 370
$\pm$
20 days. Owing to its wide orbit and comparatively large minimum mass, HD 65216 c is an excellent candidate for astrometric detection through the HGCA proper motion anomaly.
In this paper, we combine HARPS RV observations with astrometric constraints from the HGCA to refine the orbital architecture of the HD 65216 system. Section 2 presents an updated RV analysis with new orbital parameters for both planets. In Section 3, we model the HGCA proper motion anomaly to derive constraints on the inclination and true mass of HD 65216 c. Finally, in Section 4, we investigate the detectability of the system with Gaia DR4 epoch astrometry and the expected DR3–DR4 internal proper motion anomaly.
2. Radial velocity analysis
A previous RV analysis of the system (Wittenmyer et al. Reference Wittenmyer2019) identified the outer planet and constrained the minimum masses and periods of the two planets. This used data from CORALIE (Mayor et al. Reference Mayor2004) comprising 70 observations between March 1999 and March 2003, and HARPS comprising 24 observations between December 2003 and December 2013. The minimum mass of HD 65216 b was found to be 1.295
$\pm$
0.062 M
$_{\mathrm{J}}$
and the period was found to be 577.6
$\pm$
1.3 days with an eccentricity of
$0.27\pm 0.02$
, while HD 65216 c had a minimum mass of 2.03
$\pm$
0.11 M
$_{\mathrm{J}}$
, period of 5 370
$\pm$
20 days and eccentricity of
$0.17\pm 0.04$
. Since that study, more RV data has been collected with HARPS, comprising 44 observations between October 2021 and December 2022, covering almost a full period of the inner planet. This new dataset was taken after a fibre link upgrade of the HARPS spectrograph (Lo Curto et al. Reference Lo Curto2015). We analyse this updated dataset using a Bayesian framework implemented with numpyro (Phan et al. Reference Phan, Pradhan and Jankowiak2019) using JAX (Bradbury et al. Reference Bradbury2018) for gradient computation and accelerated sampling with the the No-U-Turn Sampler (Hoffman & Gelman Reference Hoffman and Gelman2011).
For our RV model, we adopt a two-planet model where the total RV is given by
\begin{equation} v(t) = \sum\limits ^{2}_{n=1}K_{n}[{\cos}{(\nu _{n}(t)+\omega _{n})}+e_{n}\cos \omega _{n}]\space +\space \gamma _{\mathrm{inst}} \end{equation}
where
$\nu _{n}(t)$
is the true anomaly of planet
$n$
, and
$K_{n}$
is the RV semi-amplitude given by
\begin{equation} K_{n} = \frac {2\pi a_{n}}{P_{n}}\frac {M_{n}}{M_{\star }+M_{n}}\frac {\sin i}{\sqrt {1-e_{n}^{2}}}, \end{equation}
where
$M_{n}$
and
$M_{\star }$
are the masses of planet
$n$
and the star, respectively, in M
$_{\odot }$
,
$a_{n}$
, and
$P_{n}$
are the planet’s semi-major axis and period, respectively. In our model, we solve for the period (
$P$
), mass ratio (
$M$
/
$M_{\star }$
), eccentricity (
$e$
), argument of periastron (
$\omega$
) and initial phase (
$\phi _{0}$
) for both planets, as well as the instrument offset (
$\gamma$
) and the RV jitter (
$\sigma$
) for each dataset.
The eccentricity is sampled from a
$\beta$
prior (Johnson et al. Reference Johnson, Kotz and Balakrishnan1995) with parameters
$a=0.867$
and
$b=3.03$
given by Kipping (Reference Kipping2013), which favours low-to-moderate eccentricities while still allowing highly eccentric solutions. The reference time for
$\phi _{0}$
is set at a Julian date (JD
$_{0}$
) of 2457388.5 which corresponds to the epoch time of Gaia DR3 (2016.0) in order to simplify combining this study with astrometry in Section 3. The initial phase is then converted into a time of periastron by
\begin{equation} T_{0} = \mathrm{JD}_{0}-\frac {\phi _{0}}{2\pi }P. \end{equation}
For convenience, we report the earliest periastron passage occurring within the RV observing window. The priors and resultant values for orbital parameters of both planets are shown in Table 1, along with the RV offset and jitter. Note the mass ratio is sampled from a log uniform prior, between ratios of
$10^{-4}$
and
$10^{-1}$
. Since we keep the host mass fixed at 0.9 M
$_{\odot }$
, this corresponds to a prior of
$\log \mathcal{U}(0.09,90)$
for Msin
$i$
in Jupiter masses, which is what we report in the table. The semimajor axis and periastron time are not fitted parameters, but quantities derived from period and initial phase respectively. Hence, no priors are listed for these parameters. The offset and jitter is labelled with the subscript COR for the CORALIE dataset and H03 and H15 for the HARPS datasets before and after the fibre link upgrade, respectively. As shown in Table 1, the period of HD 65216 b agrees with the previous study, while HD 65216 c is found to have a significantly longer period with higher uncertainty. The most notable change is that HD 65216 c now has a significantly higher eccentricity, with a peak value of 0.58. This is probably due to the sparse RV measurements between JD
$\sim$
2454000 and 2459000, which weakens constraints on curvature, and allows for long-period, high-eccentricity solutions. This gap in data also allows for a wider range of solutions for RV amplitude, which could account for the larger uncertainty in mass compared to Wittenmyer et al. (Reference Wittenmyer2019). The fitted RV curves and the phase-folded curves for both planets are shown in Figure 1. The best-fit values used for the solid curves are shown in Appendix B in Table B1.
Results from RV analysis.
2.a. Long description
A table with 12 rows and 5 columns. The columns are labeled Parameter, Prior, b, and c. The table presents the results from radial velocity analysis for two planets, labeled b and c. Each row lists a different parameter, including Msin(Mj), P(days), a(AU), e, ω(°), ϕ_0(°), T_0(JD), γ_COR(m/s), γ_H03(m/s), γ_H15(m/s), σ_COR(m/s), σ_H03(m/s), and σ_H15(m/s). The values for each parameter are provided for both planets b and c, along with their respective priors. The table includes various orbital parameters and their derived values, with some parameters not having priors listed as they are derived quantities.
Priors and results from HGCA analysis of HD 65216 c, removing RV contribution from the inner planet. There are two solutions for
$i$
and
$\Omega$
: anti-clockwise (
$i\lt 90$
$^\circ$
) and clockwise (
$i\gt 90$
$^\circ$
). Both are listed, separated by commas.
2.a. Long description
A table with two columns labeled Prior and Value, and rows labeled with different parameters. The table compares the priors and results from HGCA analysis of HD 65216 c, removing RV contribution from the inner planet. Row 1: M(M), log L(0.09, 9.0), 3.32. Row 2: P(days), log L(10^2, 10^4), 7 007. Row 3: a(AU), -, 6.92. Row 4: e, -, 0.42. Row 5: i(degrees), beta(0.867, 3.03), 35.6, 147.1. Row 6: OMEGA(degrees), L(0, 360), 107, 154. Row 7: omega(degrees), L(0, 360), 120. Row 8: phi(degrees), L(0, 360), 169. Row 9: T0(JD), -, 2 447 086. Row 10: gamma(COR)(m/s), L(-100, 100), 23.6. Row 11: gamma(H03)(m/s), L(-100, 100), 15.8. Row 12: gamma(H15)(m/s), L(-100, 100), 17.9. Row 13: sigma(COR)(m/s), L(0, 20), 4.0. Row 14: sigma(H03)(m/s), L(0, 20), 4.5. Row 15: sigma(H15)(m/s), L(0, 20), 2.1. Row 16: mu_alpha*(mas/yr), L(-128.74, -118.74), 123.88. Row 17: mu_delta(mas/yr), L(-151.60, -141.60), 146.71.
Radial velocity fits and phase-folded curves for each planet. Faint lines are from a random draw of 500 samples and dark lines are the best-fit values.
Figure 1. Long description
Panel A: A line graph shows radial velocity fit with the offset subtracted. The x-axis is labeled JD-2450000 and the y-axis is labeled RV (m/s). The graph includes data points from CORALIE, HARPS03, and HARPS15, represented in orange, green, and red respectively. Panel B: Two phase-folded curves are shown. The top graph represents HD 65216 b and the bottom graph represents HD 65216 c. Both graphs have the x-axis labeled with phase (φ) and the y-axis labeled RV (m/s). Faint lines represent a random draw of 500 samples, while dark lines represent the best-fit values.
Figure 1(a) shows that, although our two-planet solution can explain the measured RV time series, there is a long (
$\sim$
2 000 days) gap in the data which makes precise measurement of the outer planet difficult. The posterior for HD 65216 c is slightly asymmetric, with a tail extending towards periods exceeding 8 000 days, allowing for solutions in which the outer planet has not yet completed a full orbital cycle within the RV baseline. Additional RV monitoring over the coming decade will therefore be important for refining the orbital architecture of the system.
However, all our solutions for HD 65216 c’s period place it within the ideal period range for a long-term astrometric study, which is detailed in the next section.
3. Astrometric analysis
3.1. Detection threshold
The gravitational pull of a planet causes an offset of its host star’s photocentre which is given by
\begin{equation} a^{*} = \frac {a}{d}\left (\frac {M_{p}}{M_{\star }+M_{p}}\right )\!, \end{equation}
where
$a$
is the semimajor axis of the planet and
$d$
is the distance from Earth. For HD 65216 b and c, this offset is equal to 0.051 and 0.357 mas, respectively, based on our mass and period results from RV analysis, assuming an edge-on (
$i$
= 90
$^\circ$
) orbit. This implies that HD 65216 b is likely below the detection threshold for astrometry. HD 65216 c would produce a significant offset, but its long period hinders direct detection through Gaia DR3 alone. This can be solved by using the long Hipparcos-Gaia baseline of 25 yr.
The detectability of planets with Gaia DR3 alone can be quantified by the re-normalised unit weight error (RUWE) a planet’s motion could induce. This is a dimensionless quantity that represents the error in Gaia’s fitted position, parallax, and proper motion. Reflex motion due to a companion can disrupt this fit and increase the RUWE value. A value greater than 1.25 could be indicative of a companion (Penoyre et al. Reference Penoyre, Belokurov and Evans2022; Castro-Ginard et al. Reference Castro-Ginard2024). We therefore take this as our Gaia DR3 detection threshold. From this, we can calculate the minimum detectable companion mass as a function of period.
The minimum detectable mass as a function of period is shown in Figure 2 for RV, Gaia’s RUWE, and the Hipparcos-Gaia Catalogue of Accelerations (HGCA). We define as our detection thresholds, RV semi-amplitude
$ \gt $
10 m/s, RUWE
$ \gt $
1.25, and proper motion anomaly magnitude greater than the combined error of 0.064 mas/yr. From Figure 2, it can be seen that, although HD 65216 b has a period in the ideal range for detection by Gaia DR3, its mass is below the detection threshold, if we assume an edge-on orbit with mass equal to the minimum mass. This is also reflected in the RUWE value of 1.09, which is below the accepted companion detection threshold for DR3. However, HD 65216 c’s orbital period sits within the ideal range for detection with the Hipparcos–Gaia proper motion anomaly. The sensitivity of the HGCA arises because the
$\sim\!25$
yr Hipparcos–Gaia baseline samples a substantial fraction of the outer companion’s orbit, allowing the orbital acceleration to manifest as a measurable proper motion discrepancy. For HD 65216, the proper motion anomaly given by the HGCA (Brandt Reference Brandt2021) is
\begin{align} \Delta \mu _{\alpha *} = 0.14\pm 0.05\,\mathrm{mas/yr}, \\[-8pt] \nonumber \end{align}
\begin{align} \Delta \mu _{\delta } = -0.11\pm 0.04\,\mathrm{mas/yr}, \\[12pt] \nonumber \end{align}
which is given by the long-baseline proper motion between Hipparcos and Gaia (
$\boldsymbol{\mu }_{HG}$
) subtracted from the short-term proper motion recorded by Gaia (
$\boldsymbol{\mu }_{G}$
). This indicates a significant acceleration which hints at the presence of a detectable long-period companion.
Minimum detectable mass as a function of period for three methods. The detected planets from the radial velocity analysis are shown for comparison.
3.2. Proper motion modelling
In order to calculate the planets’ true masses, we rerun our numpyro analysis, this time including the inclination (
$i$
), longitude of ascending node (
$\Omega$
), and true proper motion as free parameters. As input data, in addition to the RV measurements, we include the proper motion calculated by Gaia, Hipparcos and by the long-term change in position (
$\boldsymbol{\mu }_{G}$
,
$\boldsymbol{\mu }_{H}$
and
$\boldsymbol{\mu }_{HG}$
, respectively) and associated uncertainties taken from the HGCA. These uncertainties are used to construct a 6
$\times$
6 covariance matrix.
The long baseline proper motion is modelled by calculating the offset of the star’s photocentre from the centre of mass at Hipparcos and Gaia’s epochs times, due to the orbital motion of both planets. We model the star’s position due to planet
$n$
with:
\begin{equation} r_{n}(t) = a^{*}_{n}\left [\frac {1-e_{n}^{2}}{1+e_{n}\cos \nu _{n}(t)}\right ]\!, \end{equation}
where
$a^{*}_{n}$
is the average angular offset of the photocentre due to planet
$n$
, defined by Equation (4). Other parameters are defined as in Equations (1) and (2). The offset in the R.A. (
$\alpha$
) and Dec. (
$\delta$
) directions are given by:
\begin{align} x_{\alpha } = \sum \limits _{n=1}^{2}r_{n}[\! \sin \Omega \cos (\nu _{n}+\omega _{n})+\cos \Omega \sin (\nu _{n}+\omega _{n})\cos i], \\[-10pt] \nonumber \end{align}
\begin{align} x_{\delta } = \sum \limits _{n=1}^{2}r_{n}[\! \cos \Omega \cos (\nu _{n}+\omega _{n})-\sin \Omega \sin (\nu _{n}+\omega _{n})\cos i]. \\[10pt] \nonumber \end{align}
The long-baseline proper motion is then given by:
\begin{equation} \boldsymbol{\mu }_{HG} = \frac {\boldsymbol{x}_{G}-\boldsymbol{x}_{H}}{t_{G}-t_{H}}+\boldsymbol{\mu }_{0}, \end{equation}
where
$\boldsymbol{\mu }_{0}$
is the proper motion of the barycentre (the ‘true’ proper motion) and
$\boldsymbol{x}_{G}$
and
$\boldsymbol{x}_{H}$
are the astrometric offsets of the photocentre at times
$t_{G}$
and
$t_{H}$
(2 016.0 and 1 991.25, respectively.)
The short-term proper motion from Gaia is modelled using a simplified version of the Astrometric Global Iterative Solution (Lindegren et al. Reference Lindegren2012) in which the photocentre’s offsets in Gaia’s scanning direction are modelled in order to then model the observed proper motion and parallax. We define the photocentric offset in the scanning direction due to planetary motion as
\begin{equation} x_{AL,\mathrm{planets}} = x_{\alpha }\sin \psi +x_{\delta }\cos \psi , \end{equation}
where
$\psi$
is the scanning angle. The total offset due to parallax, proper motion, and planetary orbits is then given by
\begin{equation} x_{AL} = \boldsymbol{{A}}\begin{bmatrix} \Delta \alpha * \\ \Delta \delta \\ \varpi \\ \mu _{0\alpha *} \\ \mu _{0\delta } \end{bmatrix}+x_{AL,\mathrm{planets}}, \end{equation}
where
$\Delta \alpha *$
and
$\Delta \delta$
are constant offsets in R.A. and Dec., respectively (assumed to be 0),
$\varpi$
is the parallax, and
$\mu _{0\alpha *}$
and
$\mu _{0\delta }$
are the proper motions of the barycentre.The matrix
$\boldsymbol{{A}}$
is given by
\begin{equation} \boldsymbol{{A}} = [\! \sin \psi \quad \cos \psi \quad f_{\varpi } \quad (t-t_G)\sin \psi \quad (t-t_G)\cos \psi ], \end{equation}
where
$\psi$
is the scanning angle and
$f_{\varpi }$
is the ‘parallax factor’, which is the projection of the barycentric parallax displacement vector onto the Gaia along-scan direction (Lindegren et al. Reference Lindegren2012). The expected observations are then calculated using the inverse of the first part of Equation (12):
\begin{equation} \begin{bmatrix} \Delta \alpha * \\ \Delta \delta \\ \varpi \\ \mu _{\alpha *} \\ \mu _{\delta } \end{bmatrix} = (\boldsymbol{{A}}^{T}\boldsymbol{{A}})^{-1}\boldsymbol{{A}}^{T}x_{AL} \end{equation}
For simplicity, we do not solve for the true parallax, and keep this fixed at the value reported by Gaia. This is because a planetary mass companion is unlikely to have a significant effect on the observed parallax. The Hipparcos proper motion (
$\boldsymbol{\mu }_{H}$
) is modelled the same way, using the Hipparcos scanning angles and observation times.
Although HD 65216 is a triple system, we do not consider the effect of the two M-dwarf companions, as their separation of 253 AU from the primary (Mugrauer et al. Reference Mugrauer, Seifahrt and Neuhäuser2007), with a period of
$\sim$
4 000 yr, makes the reflex motion negligible over the Hipparcos-Gaia baseline.
3.3. Results from joint RV-astrometry analysis
Given the negligible astrometric contribution expected from HD 65216 b, the current data do not meaningfully constrain mutual inclination. We therefore adopt a coplanar configuration to reduce the dimensionality of the astrometric fit and assume the same
$i$
and
$\Omega$
values for both planets. We remove the RV contribution from the inner planet and conduct a combined RV-astrometry analysis on the residual signal from the outer planet. We adopt a uniform prior in
$\Omega$
spanning 0 to 2
$\pi$
and a uniform prior in cos
$i$
spanning −1 to 1. We use uniform priors for the true proper motion of the system, spanning a space within 5 mas/yr of the Gaia-reported value.
The resultant posterior distributions of the outer planet’s orbital parameters are shown in Appendix A in Figure A2.
We find a bimodal distribution in inclination, which peaks at
$\sim$
30
$^\circ$
and 157
$^\circ$
, corresponding to anticlockwise and clockwise motion on the sky plane, respectively. The resultant true parameters of HD 65216 c, and the barycentric proper motion are summarised in Table B2. We recover a true mass for HD 65216 c of 3.32
$^{+1.40}_{-1.16}$
M
$_{\mathrm{J}}$
, placing this companion comfortably within the planetary regime. Assuming the planets are co-planar, we apply our best-fit minimum mass of HD 65216 b (1.25 M
$_{\mathrm{J}}$
) to obtain a true mass of 2.23
$^{+0.96}_{-0.77}$
M
$_{\mathrm{J}}$
. We note however, that the co-planarity assumption is not necessarily realistic and is only meant to provide a first-order approximation of the inner planet’s mass. Visualisations of the planet orbits for the anti-clockwise and clockwise cases are shown in Figure 3.
Planet orbit solutions for anti-clockwise (a) and clockwise (b) cases. Dark lines are best-fit values and faint lines are a random draw of 500 samples.
4. Prospects for DR4
Gaia’s fourth data release (DR4), expected in December 2026, will contain along-scan epoch astrometry for all sources, greatly improving our ability to detect and characterise planetary systems with astrometry. In addition to improving the precision of measurements of parallax and proper motion, the individual observations are expected to have improved astrometric precision. Based on its apparent G magnitude, HD 65216 has an astrometric uncertainty of 0.22 mas in DR3 (Lindegren et al. Reference Lindegren2021). This is expected to decrease to
$\sim$
0.16 mas in DR4, due to the longer observing baseline of 66 months, up from DR3’s 34 months. The value of this uncertainty can help assess the detectability of a companion.
4.1. Epoch astrometry
Expected epoch astrometry from parameter distributions shown in Figure A2.
Using the posterior distributions of
$i$
and
$\Omega$
, we predicted the expected along-scan offset as a function of time in DR4. The epoch times and scanning angles are obtained from Gaia’s Observation Forecast Tool. For each posterior sample from our joint fit, we simulate the along-scan offset from Equation (12) (
$x_{AL}$
), now including both planets and assuming co-planarity. We then calculate the fitted parameters using Equation (14). These are then fed back into the first half of Equation (12) to calculate the expected offsets without planets (
$x_{AL0}$
). The difference between these two offsets (
$\Delta x_{AL}$
) can be used as an indicator of a planet’s astrometric signal. In Figure 4, we show the distribution of
$\Delta x_{AL}$
for each DR4 epoch. Figure 4 shows that, while many solutions produce small offsets within the expected along-scan error of 0.16 mas, there are some solutions which produce significantly higher offsets. This suggests that epoch data may be capable of further constraining inclination and, from this, true mass. The ability of the epoch data to constrain orbital parameters can be quantified by the predicted RUWE for DR4, given by
\begin{equation} RUWE_{4} = \sqrt {\frac {\sum \nolimits _{i}^{N}(\Delta x_{i})^{2}}{\sigma ^{2}(N-5)}}, \end{equation}
where
$\Delta x_{i}$
are the along-scan offsets shown in Figure 4,
$N$
is the number of epochs and
$\sigma$
is the along-scan error. In DR3, the threshold for companion detection was a RUWE of 1.25. In DR4, this is expected to be 1.15 (Guerriero et al. Reference Guerriero, Penoyre and Brown2026). In Figure 5, we plot the distribution of inclination and mass of HD 65216 b, coloured by the predicted RUWE
$_{4}$
. As shown in Figure 5, only solutions with inclination
$\lt\!35^\circ$
or
$\gt\!\!145^\circ$
produce RUWE
$_{4}$
greater than the expected detectability threshold for DR4, which includes 43.5% of our posterior samples. These inclinations correspond to planet masses of 2.00 M
$_{\mathrm{J}}$
and 2.98 M
$_{\mathrm{J}}$
for HD 65216 b and c, respectively. With the release of DR4, if the astrometric offset is insignificant, these will be the upper limits of the planet masses. If the offset is significant, we will be able to place tighter constraints on inclination and planet masses. Additionally, it is worth noting that the inner planet HD 65216 b provides the greatest contribution to RUWE
$_{4}$
as DR4 only covers a fraction of the outer planet’s
$\sim$
19 yr orbit. This indicates that we can use epoch astrometry to constrain the properties of the inner planet and Hipparcos-Gaia proper motion anomalies for the outer planet. In this work, we assume the planets are co-planar, but a combination of these two methods will allow us to constrain the mutual inclination of the system. It may also be possible to further constrain the properties of HD 65216 b by examining the proper motion anomaly over Gaia’s changing baseline.
Distribution of inclinations and mass of HD 65216 b, coloured by predicted RUWE
$_{4}$
. Grey points are values of RUWE
$_{4}$
$ \lt $
1.15 and thus considered undetectable. Posterior distributions for inclination and planet mass are shown with medians and 1-
$\sigma$
values. In the case of inclination, the medians and 1-
$\sigma$
values are shown for both clockwise and anticlockwise orbits.
Figure 5. Long description
Panel A: A scatter plot showing the relationship between inclination (i in degrees) and mass (M b in M J) of HD 65216 b. The data points are colored by predicted RUWE, with a color bar ranging from 1.2 to 2.8. Grey points represent values of RUWE 1.15 and are considered undetectable. Panel B: Histograms showing the posterior distributions for inclination and planet mass. The medians and 1-sigma values are shown for both clockwise and anticlockwise orbits.
4.2. DR3-DR4 proper motion anomaly
With an extra
$\sim$
2.7 yr of observing baseline, Gaia DR4 will cover a larger portion of the outer planet’s orbit. As a result, the measured proper motion would be different than that measured by DR3. This ‘internal’ proper motion anomaly can also be used to constrain orbital parameters.
We simulate the expected proper motion anomaly by first simulating the epoch astrometry in DR3 and DR4 for all posterior samples, using Equation (12), and then calculating the proper motion fits using Equation (14). The proper motion anomaly is then simply the difference between these fits. In DR3, the uncertainty in
$\mu _{\alpha *}$
is 0.029 mas/yr for HD 65216 and the uncertainty in
$\mu _{\delta }$
is 0.023 mas/yr. In DR4, these uncertainties are expected to decrease by a factor of
$2\sqrt {2}$
(Brown Reference Brown2025). We adopt these new values as the DR4 uncertainties, and calculate the uncertainties in the proper motion anomaly
$\Delta \boldsymbol{\mu }$
by combining the DR3 and DR4 uncertainties. This is produces uncertainties of 0.030 and 0.025 mas/yr for
$\Delta \mu _{\alpha *}$
and
$\Delta \mu _{\delta }$
, respectively.
In Figure 6, we show the posterior distribution of inclination (
$i$
) and longitude of the ascending node (
$\Omega$
), coloured by the proper motion anomaly in both the R.A. and Dec. dimensions. From this study, we find that 94.8% of our posterior samples produce a detectable proper motion anomaly in at least one of the R.A. and Dec. dimensions. Furthermore, Figure 6 demonstrates how the sign of the DR3–DR4 proper motion is affected by the direction of the orbit. Which suggests that, with a significant proper motion anomaly, we will go beyond the capabilities of the Hipparcos-Gaia anomalies and, in addition to calculating the true mass, we will distinguish between clockwise and anti-clockwise orbits.
Distribution of
$i$
and
$\Omega$
, coloured by DR3–DR4 proper motion anomaly. Grey points are solutions with anomalies below the expected uncertainties of 0.030 and 0.025 mas/yr for
$\Delta \mu _{\alpha *}$
and
$\Delta \mu _{\delta }$
, respectively. Dashed lines show the 1-
$\sigma$
level of the posterior.
5. Conclusions
We have presented a joint RV and astrometric analysis of the HD 65216 planetary system using updated HARPS RV observations and the HGCA. Our revised RV analysis confirms the presence of the long-period outer companion HD 65216 c and favours a longer-period, more eccentric orbit than previously reported. In particular, we find evidence for a highly eccentric orbit with
$e = 0.61^{+0.11}_{-0.10}$
and a period of
$6\,455^{+273}_{-84}$
days.
By combining the RV observations with the HGCA proper motion anomaly, we constrained the orbital inclination and derived true dynamical masses for both companions. The posterior distribution is bimodal, corresponding to anticlockwise and clockwise motion, with preferred inclinations near 30
$^\circ$
and 157
$^\circ$
, respectively. The inferred true masses of HD 65216 b and c remain firmly within the planetary regime, supporting their classification as giant planets rather than brown dwarfs.
We further investigated the detectability of the system with upcoming Gaia DR4 astrometry. Simulations of both epoch astrometry and the DR3–DR4 internal proper motion anomaly indicate that future Gaia measurements should provide meaningful additional constraints on the orbital geometry of the inner planet. In particular, low-inclination (high mass) solutions are expected to produce along-scan astrometric signals comparable to or exceeding the anticipated DR4 single-epoch precision. As the Hipparcos-Gaia proper motion anomaly is more sensitive to the outer planet and Gaia epoch astrometry is more sensitive to the inner planet, these methods can be used in combination to allow constraints on the mutual inclination of the system.
HD 65216 therefore represents an excellent benchmark system for studying the synergy between long-baseline RV monitoring and Gaia astrometry in the characterisation of long-period giant planets. Continued RV observations, combined with future Gaia releases, will enable substantially improved constraints on the orbits of both planets and may ultimately permit a fully dynamical three-dimensional reconstruction of the system architecture.
Acknowledgements
This work was funded by the Australian Research Council Discovery Grant DP250101273.
This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange, and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch.
Data availability statement
The data underlying this analysis are publicly available in the DACE and HGCA databases and also stored with the analysis programs on the repository https://github.com/awallace142857/HD65216
Appendix A. Posterior corner plots
The posteriors of the RV fits of both planets are shown in Figure A1. The posteriors of the joint RV-astrometric fit of HD 65216 c are shown in Figure A2.
Posterior distribution of planet parameters from RV analysis. Median and 68% confidence intervals are shown for each parameter.
Figure A1. Long description
Panel A: Parameters of HD 65216 b. This panel contains multiple scatter plots and histograms. Each scatter plot shows the relationship between two parameters, with contour lines indicating density regions. The histograms display the distribution of individual parameters. The parameters include mass (M sub p) in Earth masses, period (P sub b) in days, eccentricity (e sub b), argument of periastron (omega sub b) in degrees, and time of periastron passage (T sub p) in Julian days. Panel B: Parameters of HD 65216 c. Similar to Panel A, this panel also contains multiple scatter plots and histograms. The parameters include mass (M sub c) in Earth masses, period (P sub c) in days, eccentricity (e sub c), argument of periastron (omega sub c) in degrees, and time of periastron passage (T sub p) in Julian days. Each scatter plot and histogram provides insights into the posterior distribution of these parameters from the RV analysis.
Posterior distribution of the properties of HD 65216 c, determined from joint RV-astrometric fit. Median and 68% confidence intervals are shown for each parameter. Due to the bimodal inclination distribution, two medians and confidence intervals are shown for
$i$
and
$\Omega$
distributions.
Appendix B. Additional tables
Tables B1 and B2 show the best fit results of the RV and joint fits, respectively.
Maximum likelihood values from RV analysis. These values were used for the solid curves in Figure 1.
Table B1. Long description
A table with two columns labeled ‘Parameter’ and ‘Prior’ and two rows labeled ‘b’ and ‘c’. The table compares various orbital parameters and their values for two planets. Row 1: Parameter, Prior, b, c. Row 2: Msin(M), 1.27, 2.00. Row 3: P(days), 577.9, 6689. Row 4: a(AU), 1.31, 6.71. Row 5: e, 0.25, 0.58. Row 6: ω(deg), 225, 116. Row 7: T0(JD), 2451747, 2454211. Row 8: γCOR(m/s), 16.9. Row 9: γH03(m/s), 15.4. Row 10: γH15(m/s), 14.6. Row 11: σCOR(m/s), 6.3. Row 12: σH03(m/s), 6.3. Row 13: σH15(m/s), 7.1.
Maximum likelihood values from RV-astrometry joint analysis of HD 65216 c.
Table B2. Long description
A table summarizing the resultant true parameters of HD 65216 c and the barycentric proper motion. The table has 16 rows and 2 columns. The first column lists the parameters, and the second column lists their corresponding values. The parameters include M2 (MJ), Pc (days), ec, ic, Ωc (degrees), ωc (degrees), T0c (JD), γCOR (m/s), γH03 (m/s), γH15 (m/s), σCOR (m/s), σH03 (m/s), σH15 (m/s), μH15 (mas/yr), and μc (mas/yr). The values are specific numerical results from the analysis.
i
Ω
i<90
∘
i>90
∘
4
4
<
σ
σ
i
Ω
Δμα∗
Δμδ
σ
i
Ω
Open access
