# Dynamical Mass Estimates of the *β* Pictoris Planetary System through Gaussian Process Stellar Activity Modeling

Thomas Vandal, Julien Rameau, René Doyon

► To cite this version:

Thomas Vandal, Julien Rameau, René Doyon. Dynamical Mass Estimates of the *β* Pictoris Planetary System through Gaussian Process Stellar Activity Modeling. The Astronomical Journal, 2020, 160, 10.3847/1538-3881/abba30. insu-03705170

HAL Id: insu-03705170

https://insu.hal.science/insu-03705170v1

Submitted on 3 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Dynamical Mass Estimates of the *β Pictoris* Planetary System through Gaussian Process Stellar Activity Modeling

Thomas Vandal<fnref n="1" />,<fnref n="2" /> ID, Julien Rameau<fnref n="1" />,<fnref n="3" /> ID, and René Doyon<fnref n="1" /> ID

## Abstract

Nearly 15 yr of radial velocity (RV) monitoring and direct imaging enabled the detection of two giant planets orbiting the young, nearby star *β Pictoris*. The *δ Scuti* pulsations of the star, which overwhelm planetary signals, need to be carefully suppressed. In this work, we independently revisit the analysis of the RV data following a different approach than available in the literature to model the activity of the star. We show that a Gaussian process (GP) with a stochastically driven damped harmonic oscillator kernel can model the *δ Scuti* pulsations. It provides similar results to parametric models but with a simpler framework, using only three hyperparameters. It also enables us to model poorly sampled RV data that were excluded from previous analyses, hence extending the RV baseline by nearly five years. Altogether, the orbit and mass of both planets can be constrained from RV only, which was not possible with the parametric modeling. To characterize the system more accurately, we also perform a joint fit of all available relative astrometry and RV data. Our orbital solutions for *β Pic b* favor a low eccentricity of $0.029_{-0.024}^{+0.061}$ and a relatively short period of $21.1_{-0.8}^{+2.0}$ yr. The orbit of *β Pic c* is eccentric with $0.206_{-0.063}^{+0.074}$ with a period of $3.36 \pm 0.03$ yr. We find model-independent masses of $11.7 \pm 1.4$ and $8.5 \pm 0.5 M_{Jup}$ for *β Pic b* and *c*, respectively, assuming coplanarity. The mass of *β Pic b* is consistent with the hottest start evolutionary models, at an age of $25 \pm 3$ Myr. A direct detection of *β Pic c* would provide a second calibration measurement in a coeval system.

*Unified Astronomy Thesaurus concepts: Exoplanets (498); Radial velocity (1332); Gaussian Processes regression (1930); Astrometry (80); Orbit determination (1175); Delta Scuti variable stars (370)*

## 1. Introduction

Directly detected exoplanets must rely on theoretical models to determine their masses from their luminosity and age. Not only is age difficult to evaluate (Soderblom et al. 2014), but evolutionary models are also uncalibrated at young ages (≤1 Gyr) and low masses (≤70 $M_{Jup}$). Model-independent mass measurements must be obtained to test these models, which ultimately will provide confidence in the inferred masses for all directly imaged planets. Absolute astrometry with exquisite precision enables such mass measurements. Doppler spectroscopy of planet host stars can also determine the planet masses, biased with the orbit inclination which can be constrained from relative astrometry to overcome this bias. However, directly imaged planets are usually found at long periods around young stars because they are bright and spatially resolved with current instruments (see Bowler et al. 2014 and references therein). Radial velocity (hereafter RV) data must therefore be obtained over a long timescale. These data are however severely plagued by strong stellar activity signal due to spots, plages, and pulsations (Galland et al. 2005; Lagrange et al. 2013; Hillenbrand et al. 2015). Planet detection is thus very challenging on both observational and data analysis sides. However, these limitations are less severe for low-mass star and high-mass brown-dwarf companions to stars, enabling independent mass and luminosity measurements using a combination of relative, absolute astrometry, and/or RV data. Evolutionary model calibration is therefore being undertaken for these objects. (Crepp et al. 2012; Ryu et al. 2016; Dupuy & Liu 2017). In the planetary mass regime, reduced signal amplitudes (RV, astrometry, contrast) make this task much harder.

*β Pictoris* (hereafter *β Pic*) is a nearby ($19.44 \pm 0.05$ pc; van Leeuwen 2007), A6V (Gray et al. 2006), intermediate-mass ($1.80_{-0.04}^{+0.03} M_{⊙}$) Wang et al. (2016)) star. *β Pic* has been classified as a *δ Scuti* variable (Koen 2003) characterized by pulsation frequencies in the range $22.81-75.68$ day$^{-1}$ (Mékarnia et al. 2017; Zwintz et al. 2019; Zieba et al. 2019). *β Pic* is a member of the *β Pictoris* moving group (Zuckerman et al. 2001), with recent age estimates of $24 \pm 3$ Myr (Bell et al. 2014) and $26 \pm 3$ Myr (Nielsen et al. 2016) using isochrone fitting, or $25 \pm 3$ Myr by measuring the lithium depletion boundary on group members (Messina et al. 2016). *β Pic* is surrounded by a complex environment. The large circumstellar disk of gas and debris (Smith & Terrile 1984; Dent et al. 2013) has a primary component seen almost edge on and an inner warp component titled by $+4.0 \pm 0.6$ and further inclined by $6^\circ \pm 1^\circ$ (Burrows et al. 1995; Mouillet et al. 1997; Golimowski et al. 2006; Ahmic et al. 2009; Lagrange et al. 2012a). Falling evaporating bodies or exocomets were evidenced (Lagrange et al. 1996; Kiefer et al. 2014; Zieba et al. 2019). Both warp and comets were modeled with gravitational perturbations by planet(s) within the system (Beust et al. 1990; Beust & Morbidelli 1996; Beust et al. 1998; Mouillet et al. 1997). Two giant planets were discovered within the system: *β Pic b* by direct imaging (Lagrange et al. 2010) and *β Pic c* with RV (Lagrange et al. 2019b). Monitoring the stellar RVs from 2008 (2003) up to 2018 combined with a multiparametric sinusoidal fit for the *δ Scuti* pulsations led to a residual noise of $10$ m s$^{-1}$, allowing for the detection of *β Pic c*. The system is therefore a rare case for

<footnote n="1">¹ Institut de Recherche sur les Exoplanètes, Département de physique, Université de Montréal, CP 6128 Succ. Centre-ville, H3C 3J7, Montréal, QC, Canada
vandal@astro.umontreal.ca</footnote>

<footnote n="2">² Department of Physics, McGill University, 3600 rue University, H3A 2T8, Montréal, QC, Canada</footnote>

<footnote n="3">³ Université Grenoble Alpes/CNRS, Institut de Planétologie et d'Astrophysique de Grenoble, F-38000 Grenoble, France
*Received 2020 June 30; revised 2020 September 11; accepted 2020 September 17; published 2020 November 3*</footnote>
which independent mass measurement can be foreseen, as well as a deep understanding of a complex planetary system.

Simultaneous imaging of β Pic b and the disk proved that β Pic b is not coplanar with the main disk but consistent with being responsible for the inclined warp (Lagrange et al. 2012a). Monitoring the motion of β Pic b with direct imaging has led to various constraints on its orbit. The inclination is now well constrained, excluding a perfectly edge-on orbit, at $i = 89.04^\circ \pm 0.04^\circ$ (Wang et al. 2016; Lagrange et al. 2019a; Nielsen et al. 2020; Nowak et al. 2020). Period and eccentricity were also constrained with direct imaging. Lagrange et al. (2019a) reported $P_b = 20.29_{-1.35}^{+0.86}$ yr and $e_b = 0.01_{-0.01}^{+0.03}$ when using relative astrometric data only from VLT instruments (NaCo, SPHERE) while Nielsen et al. (2020) obtained $P_b = 21.3_{-1.0}^{+2.2}$ yr and $e_b = 0.038_{-0.029}^{+0.063}$ with all relative astrometric measurements except those from VLT/SPHERE. However, the inclusion of absolute astrometry from Hipparcos and Gaia in the analysis by Nielsen et al. (2020) lead to significantly higher period and eccentricity values, namely $P_b = 24.4_{-1.5}^{+1.0}$ yr and $e_b = 0.12_{-0.03}^{+0.04}$. An analysis of the absolute astrometric motion of the star allowed Snellen & Brown (2018) to directly measure a planetary mass of $M_b = 11 \pm 2$ $M_{Jup}$; Dupuy et al. (2019) found a mass of $M_b = 13 \pm 3$ $M_{Jup}$ with additional corrections to the same data and further including relative astrometry and RV data and fitting as well for the mass of the star. Analysis of the stellar RVs using VLT relative astrometric priors yielded $M_b \sim 10$ $M_{Jup}$ (Lagrange et al. 2019b). Lastly, a joint analysis of direct imaging, RV, and astrometric data accounting for both planets resulted in a mass of $8.0 \pm 2.6$ $M_{Jup}$ for β Pic b (Nielsen et al. 2020). This is lower than previous estimates, but still consistent within $1\sigma$. These last two mass measurements account for the presence of a second giant planet in the system. The addition of a VLTI/GRAVITY relative astrometric point with an exquisite precision of 80 µas, but without the RVs and without accounting for a second planet, led Nowak et al. (2020) to measure $M_b = 12.7 \pm 2.2$ $M_{Jup}$ while excluding circular orbit ($e_b = 0.15_{-0.04}^{+0.05}$).

From stellar RVs only, the mass and period of β Pic c were estimated to be $M_c \sim 8.9$ $M_{Jup}$ (assuming the same inclination as β Pic b) and $P_c \sim 1220$ days (Lagrange et al. 2019b), while Nielsen et al. (2020) found $M_c = 9.18_{-0.87}^{+0.96}$ $M_{Jup}$, $P_c = 1238_{-11}^{+7}$ days, and $e_c = 0.24_{-0.10}^{+0.11}$.

RV data are therefore of primary importance when fitting for the masses of both planets, the absolute astrometry being still plagued by imperfect knowledge, bias, and required corrections of Gaia DR2 data onto bright stars (Brown et al. 2018). For β Pic, the $\delta$ Scuti pulsations result in RV noise with a semiamplitude of $\sim 500$ m s$^{-1}$, while β Pic b and c have expected semiamplitudes of $K_b \sim 70$ m s$^{-1}$ and $K_c \sim 120$ m s$^{-1}$, respectively. The main challenge is therefore to properly model and subtract the stellar activity, thus revealing planetary signals in the RV data.

A useful and increasingly used method to mitigate the stellar noise in the RV data of an active star is to use Gaussian processes (GPs). There are several successful attempts to use this method in the literature (e.g., Haywood et al. 2014; Grunblatt et al. 2015; Cloutier et al. 2017a, 2017b). However, these examples are generally limited to less massive stars from the field for which the stellar activity has fewer components and less temporal variations than for young stars like β Pic.

Here, we investigate possible improvements of the modeling of stellar activity in the RV data of β Pic in a Bayesian framework. We explore GPs for correlated Gaussian noise modeling over traditional parametric approaches from the

<figure><img src="image_3.png" /><figcaption>**Figure 1.** ASTEP photometric light curve of β Pic. The black dots are the observations between 2017 June 7 (top left) and 13 (bottom right). The data are split into segments of length $\Delta$JD = 0.7 day, and a constant offset is added to each segment to show the pulsations clearly. The solid blue line is the fit from an SHO GP, and the shaded area is the $1\sigma$ confidence interval.</figcaption></figure>

literature (Lagrange et al. 2019b). This investigation leads to a measurement of the orbits and masses of both β Pic b and c from RV data only, an improvement over the literature that was unable to constrain the orbital parameters and mass of β Pic b from solely the same data. We then present a new joint analysis of all RV and astrometric data available, including RV points between 2003 and 2008 that were previously ignored in the literature. This analysis leads to new mass estimates for both planets.

The paper is structured as follows. Section 2 gives an overview of all archival data used for the analysis. In Section 3, we outline the main elements of the GP method used to define a new stellar activity model for β Pic, followed in Section 4 by a detailed orbit-fitting analysis using RV alone, then adding astrometry data. A discussion of these results is presented in Section 5 followed by a summary and conclusions in Section 6.

## 2. Literature Data

In this work, we made use of published data that monitored the system over nearly the last 15 yr. We were agnostic to the type and origin of the data but deliberately excluded absolute astrometric measurements from Gaia. Their analysis requires multiple known correction factors to compensate for the star being very bright for the instrument (Brown et al. 2018) but also still unknown analysis steps to provide accurate measurements. We defer their inclusion into this type of work once the Gaia data release will be secured for bright stars.

### 2.1. ASTEP Photometry Time Series

We had access to 6358 photometric observations of β Pic from the Antarctica Search for Transiting Extrasolar Planets (ASTEP) 400 mm telescope (Crouzet et al. 2018), obtained between 2013 June 7 and 13 (JD = 2457910 to JD = 2457917), with a sampling rate of $\sim 1000$ day$^{-1}$ (Mékarnia et al. 2017, private communication). This light curve is shown in Figure 1, along with a GP fit explained in Section 3.1. These observations were part of a longer monitoring campaign of β Pic (2017 March to September, Mékarnia et al. 2017). For this longer light curve, a periodogram analysis allowed Mékarnia et al. (2017) to detect 31 $\delta$ Scuti pulsation modes and to model the photometric stellar activity as a sum of sine waves. Because all of the 31 detected frequencies were in the 34.76–75.68 day$^{-1}$ interval, the fact that we only had access to a small portion of the photometric time
<figure><img src="image_4.png" /><figcaption>**Figure 2.** HARPS RV time series of β Pic. The inset on the bottom shows an example of the pulsations causing the high-amplitude noise in the RV data. The top inset shows how the sparse pre-2008 sampling hides these pulsations.</figcaption></figure>

series (6.8 days) should not prevent us from characterizing the activity of the star.

## 2.2. HARPS RV Time Series

Lagrange et al. (2012b, 2019b) presented 5362 RV measurements of β Pic obtained with the High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003). These RV measurements, shown in Figure 2 and not corrected for stellar activity, were computed with the SAFIR package (Galland et al. 2005). These measurements, taken between 2003 October and 2018 May, were sampled with two different methods. The sampling is much denser after 2008 March and reveals high-frequency variations in the RV data, with a similar timescale to the δ Scuti pulsations detected in the photometry (Section 2.1), that were not visible before that date (Lagrange et al. 2012b). The two insets in Figure 2 show an example of these two different sampling methods (pre-2008 and post-2008). The poor sampling from before 2008 needs to be considered carefully when analyzing this data, especially if we attempt to model the high-frequency stellar activity in the radial velocities. One option is to ignore pre-2008 data and use only the well-sampled post-2008 data (Lagrange et al. 2019b; Nielsen et al. 2020). However, the observation baseline of the RV data is shorter than the period of β Pic b, so five additional years of data, leading to more than 14 yr of RV coverage, would provide valuable information about the system. This is addressed in Section 3.2, where we discuss the inclusion of pre-2008 data.

## 2.3. Relative Astrometry

We use all relative astrometry data available for β Pic. This includes 9 epochs from VLT/NaCo (Lagrange et al. 2010; Currie et al. 2011; Chauvin et al. 2012), 7 epochs from Gemini-South/NICI (Nielsen et al. 2014), 2 epochs from Magellan/MagAO (Morzinski et al. 2015), 15 epochs from Gemini-South/GPI (Wang et al. 2016; Nielsen et al. 2020), and 12 epochs from VLT/SPHERE (Lagrange et al. 2019a). These measurements are all shown in Figure 10.

We also include the single epoch of relative astrometry obtained with VLTI/GRAVITY by Nowak et al. (2020) on 2018 September 22. This data point, obtained by averaging 17 exposure files, is given by the bivariate normal distribution $(\Delta R.A., \Delta Decl.) \sim N(\mu, \Sigma)$. $\mu$ is the mean relative planet-to-star position and $\Sigma = \text{Cov}(\Delta R.A., \Delta Decl.)$, the covariance matrix of all exposure files, gives the $1\sigma$ confidence interval.

The mean coordinates and their covariance, as reported by Nowak et al. (2020), are

$$ \left\{ \begin{aligned} \mu &= \begin{bmatrix} \Delta R.A. \\ \Delta Decl. \end{bmatrix} = \begin{bmatrix} 68.48 \\ 128.31 \end{bmatrix} \text{mas} \\ \Sigma &= \begin{bmatrix} 0.0027 & -0.0035 \\ -0.0035 & 0.0045 \end{bmatrix} \text{mas}^2. \end{aligned} \right. \qquad (1) $$

This yields a precision more than an order of magnitude better than the other relative astrometry measurements mentioned above. The GRAVITY data point is shown in the inset of Figure 10.

# 3. Stellar Activity Modeling

## 3.1. ASTEP Photometry

Mécarnia et al. (2017) detected 31 δ Scuti pulsation frequencies between 34.76 and 75.68 days⁻¹ in the ASTEP photometric light curve of β Pic, with a peak at 47.43 day⁻¹ (period of 30.4 minutes). Similar high-frequency pulsations appear in the RV data (Lagrange et al. 2012b, 2019b), and their amplitude is far greater than the expected signals from β Pic b and c. Efficiently modeling and removing this stellar activity in the RV data is crucial if we intend to detect Keplerian planetary signals. Because the RV and photometric data sets exhibit similar pulsation patterns, we first train a null-mean GP on the ASTEP photometry. We use the celerite Python package (Foreman-Mackey et al. 2017) to compute the GP models used in the following analysis. To obtain posterior distributions of the GP hyperparameters, we use the affine-invariant Markov Chain Monte Carlo (MCMC) sampler from emcee (Foreman-Mackey et al. 2013). We use 128 walkers and ensure that all chains are longer than 40 times the integrated autocorrelation time. Once we have a distribution of the GP hyperparameters from photometry, we can apply it as a prior when fitting stellar activity in the RV data.

We first explore the use of a quasi-periodic (QP) kernel typically employed to model stellar rotation for stars quieter than β Pic (Haywood et al. 2014; Grunblatt et al. 2015; Rajpaul et al. 2015; Cloutier et al. 2017a; Angus et al. 2018). We refer to Foreman-Mackey et al. (2017) for an explicit definition of the QP kernel, both in the celerite framework (their Equation (56)) and in its more classical form (their Equation (55)). We find that this kernel is able to properly model the pulsations of β Pic, recovering the main period of $P = 30.4$ minutes.

Foreman-Mackey et al. (2017) introduced a covariance function that represents a stochastically driven damped simple harmonic oscillator (SHO). This kernel is defined as

$$ k_{\text{SHO}}(\tau; S_0, Q, \omega_0) = S_0 \omega_0 Q \exp\left(-\frac{\omega_0 \tau}{2Q}\right) \times \begin{cases} \cosh(\eta \omega_0 \tau) + \frac{1}{2\eta Q} \sinh(\eta \omega_0 \tau), & 0 < Q < 1/2 \\ 2(1 + \omega_0 \tau), & Q = 1/2 \\ \cos(\eta \omega_0 \tau) + \frac{1}{2\eta Q} \sin(\eta \omega_0 \tau), & Q > 1/2 \end{cases}, \quad (2) $$

where $\eta = |1 - (4Q^2)^{-1}|^{1/2}$, $\omega_0$ is the frequency of the undamped oscillator, $Q$ is the quality factor, and $S_0$ is proportional to the power at $\omega = \omega_0$, $S(\omega_0) = \sqrt{2/\pi} S_0 Q^2$. $\tau = |t_n - t_m|$ is the absolute time difference between any two observations taken at times $t_n$ and $t_m$. When varying the quality
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Posterior distribution of the SHO GP hyperparameters from the fit on the ASTEP photometry. Broad, uniform priors were used on all three hyperparameters. The dashed lines and gray shaded areas show the median and the $1\sigma$ interval for each hyperparameter. The solid blue line shows the main frequency detected by Mékarnia et al. (2017).</figcaption></figure>

factor, the SHO GP can provide an effective model for a large variety of phenomena. In the large $Q$ limit, it can mimic QP pulsations. Moreover, the spectrum of an SHO covers a certain range of frequencies that depend on its hyperparameters. This ability to account for a wide range of frequencies provides some physical justification to adopt such a kernel in the context of a pulsating star like $\beta$ Pic.

We find that the SHO GP models the $\delta$ Scuti pulsations of $\beta$ Pic adequately, yielding an $\omega_0$ value corresponding to $P \sim 30$ minutes. The posterior distribution of hyperparameters is shown in Figure 3. We note that the BIC favors the SHO kernel over the QP kernel very strongly, with $\Delta\text{BIC} = 227$.

We consider both kernels in the following analysis of the RV data because both sampling and baseline are different from those of the light curve, possibly affecting their responses to the fit.

## 3.2. HARPS RV

Typically, when a GP is used to model stellar activity in RV data, the GP and the Keplerian orbits are jointly fit, i.e., the mean function of the GP is the Keplerian model, and they are modeled all at once on the entire data set. This method has been successfully used in the past (e.g., Haywood et al. 2014; Rajpaul et al. 2015; Grunblatt et al. 2015; Cloutier et al. 2017a). However, the stars studied were typically much quieter than $\beta$ Pic. For this reason, we investigate an alternative approach similar to the one used by Lagrange et al. (2019b): we model the activity locally by splitting the data into subsets corresponding to continuous observing sequences (hereafter referred to as groups). Within these groups, the maximum time separation between observations is 8 days. For each group, we use a GP model trained on the photometry, i.e., using posterior distributions from photometry as priors, with a constant mean. By letting all parameters vary in our MCMC, i.e., the GP

<figure><img src="image_6.png" /><figcaption>**Figure 4.** RV offsets for each group. The multisine offsets from Lagrange et al. (2019b) are shown in red, and the SHO GP offsets from this work are shown in blue. All values for post-2008 offsets are consistent within $1\sigma$, and the uncertainties on the GP offsets tend to be smaller than on the multisine ones by a factor of 1.5.</figcaption></figure>

hyperparameters and the constant mean, we can obtain an offset with an uncertainty (from the MCMC posterior distribution) for each group. These offsets thus represent the RV of the star at each epoch, and we can then fit them with a Keplerian model.

As mentioned in Section 2.2, there is RV data available from 2003 to 2018. However, before 2008, the sampling is too sparse to exhibit the pulsations in the RV data. For this reason, Lagrange et al. (2019b) used only the data from after 2008. We note that using a multisine parametric model fitted on the well-sampled post-2008 groups to fit groups from before 2008 is not possible because of temporal variations in the phase and amplitude of the pulsation modes. A GP, on the other hand, does not retain any phase information and is not strongly affected by small changes in the amplitude of individual modes. Therefore, only the covariance properties of the pulsations matter for a given GP to effectively model various epochs. We use a GP to model all available RV data between 2003 and 2018. After splitting the data as described above, we obtain 36 groups: the same 28 post-2008 groups as Lagrange et al. (2019b) and 8 additional pre-2008 groups. We then fit each group with a GP model previously trained on the photometry. To define the prior on the GP hyperparameters in the RV analysis, we use a kernel density estimate (KDE) of the posterior distributions from the photometry model (Figure 3 for the SHO GP). We do not use such a prior for the offset and the amplitude, which are both expected to vary between the two data sets.

We first test the QP kernel, $k_{\text{QP}}$, which properly accounts for the pulsations but yields a large uncertainty on the RV offset. The uncertainties on the RV offsets obtained with this kernel are approximately four times greater than those from the multisine parametric model used by Lagrange et al. (2019b).

The SHO kernel, on the other hand, results in offsets consistent with those from the multisine model, as shown in Figure 4. It yields uncertainties about 1.5 times smaller, while requiring three hyperparameters instead of up to 90 parameters (a sum of up to 30 sine waves with 3 parameters each; Lagrange et al. 2019b). We also note that the rms of the residuals from the SHO GP fit for all groups combined (with the offset subtracted from each group) is smaller by a factor of 1.4 than the rms of the multisine residuals published by Lagrange et al. (2019b). The SHO GP thus provides an effective alternative to the multisine model. The posterior distribution of the three GP hyperparameters and the constant
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Example fit and posterior distribution for the SHO GP on RV group #15. The dashed lines and gray shaded areas show the median and the $1\sigma$ interval for each parameter (the three GP hyperparameters and a constant offset). The solid blue lines show the median value of the training distribution for hyperparameters that were trained on the photometry. For $S_0$ and the RV offset, broad uniform priors were used. In the top-right corner, the GP fit of the RV data in group #7 is shown as a solid blue line, with a shaded area showing the $1\sigma$ confidence interval. The bottom panel shows the residuals, the dashed blue line being the resulting offset of the group.</figcaption></figure>

offset, as well as the stellar activity fit for group # 15 (2009 December 4), are shown in Figure 5 as an example. Similar fits are obtained for the other groups, and the hyperparameters are approximately constant throughout the groups, always staying within $1\sigma$ of other groups and of the training values. No long-term drift is thus further introduced with the modeling. This is also true for the poorly sampled pre-2008 groups.

Both the QP and SHO GPs were a priori reasonable models for the pulsations of β Pic. However, because the SHO GP is favored by the BIC and also yields a better RV precision, we use RV offsets from this model in the remaining analysis.

## 4. Orbit Fitting

The stellar activity signal has been effectively corrected from the RV, enabling the two planets from the RV offsets to be fitted in order to constrain their orbital parameters and masses. In this section, we first explore the information content of the RV alone. Our goal is to assess whether our extended baseline and increased precision can lead to any constraint on the two planets. Then, we make use of all observables of the two planets in addition to the RVs, including relative astrometry from the ground using adaptive optics and interferometric data, in a joint fit. This complete analysis will provide robust constraints on the planet parameters.

### 4.1. RV Offsets Only

We perform a two-planet fit of the RV offsets using the `radvel` Python package (Fulton et al. 2018), along with `emcee` (Foreman-Mackey et al. 2013) for MCMC sampling. To evaluate the impact of the extended RV baseline and the sensitivity to priors from relative astrometry, we explore five different cases. In

all cases, we use uniform or log-uniform priors for the parameters of $\beta$ Pic c. For each planet, unless specified otherwise, we adjust five free parameters: the period $P$, the time of conjunction $T_c$, $h = \sqrt{e} \cos \omega$, $k = \sqrt{e} \sin \omega$, and the RV semiamplitude $K$. We also fit for a global RV offset $\gamma$ because the RVs were computed with respect to a reference spectrum (Galland et al. 2005; Lagrange et al. 2019b). We use the alternative parameters $h$ and $k$ because they reduce the correlation between the eccentricity $e$ and argument of periastron $\omega^4$ while also reducing the bias of the MCMC toward high eccentricities (Ford 2006). Additionally, we use the time of conjunction $T_c$ instead of the time of periastron $T_p$ because it is well defined for circular and low-eccentricity orbits.

The five cases are the following:

1. Case A. We consider only the post-2008 RV offsets, obtained with the SHO GP. We impose Gaussian priors on the orbital parameters of $\beta$ Pic b based on the results from Lagrange et al. (2019a). This is similar to the analysis from Lagrange et al. (2019b), but we use Bayesian priors instead of repeating the fit with fixed orbital parameters randomly drawn from the prior for $\beta$ Pic b.

2. Case B. This is similar to Case A, but we fix the orbital parameters of $\beta$ Pic b to the best-fit values from Lagrange et al. (2019a) and its mass to $10 M_{\text{Jup}}$. As mentioned above, for each fit performed in Lagrange et al. (2019b), the orbital parameters of $\beta$ Pic are fixed to randomly drawn values from direct imaging results. $M_b$ was also held fixed, when using RV offsets only and not all RV residuals. Case B is thus very close to the analysis from Lagrange et al. (2019b) and enables a direct comparison of the parameters of $\beta$ Pic c between their analysis and ours.

3. Case C. Here, we include all available RV offsets between 2003 and 2018, obtained with the SHO GP. The eight additional RV offsets described in Section 3.2 provide an RV baseline of approximately 14 yr. We use Gaussian priors from Lagrange et al. (2019a) on all parameters of $\beta$ Pic b except its period and its mass. For $M_b$, we use a wide Jeffreys priors, and for $P_b$, we apply a Gaussian prior centered at 22 yr with $\sigma = 4$ yr. This prior includes all period estimates published by Lagrange et al. (2019a; ~20 yr) and Nielsen et al. (2020; ~24 yr) within $1\sigma$.

4. Case D. This is the same as Case C, but the priors from Lagrange et al. (2019a) are replaced with the slightly different priors from Nielsen et al. (2020) to test the sensitivity of the posteriors.

5. Case E. To test whether the extended RV data can reveal a change in eccentricity as it did for the period, we use a different parameterization of the Keplerian. Adjusting the parameters $P$, $T_c$, $e$, $\omega$, and $K$, which is already implemented in `radvel`, simplifies the use a uniform prior on $e_b$ while imposing a Gaussian prior on $\omega$, whose value is consistent in both Lagrange et al. (2019a) and Nielsen et al. (2020).

For each case, the adopted priors and the median of the parameter posterior distributions, along with the $1\sigma$ uncertainties corresponding to the 16th and 84th percentiles, are summarized in Table 1. We also report the Bayesian Information Criterion

<footnote n="4">⁴ In `radvel`, $\omega$ is the argument of periastron of the star's orbit.</footnote>
<table><caption>**Table 1**
Two-planet Fit with RV Data Only</caption><thead><tr><th rowspan="3"></th><th colspan="2">Case A</th><th colspan="2">Case B</th><th colspan="2">Case C</th><th colspan="2">Case D</th><th colspan="2">Case E</th></tr><tr><th colspan="2">2008–2018</th><th colspan="2">2008–2018</th><th colspan="2">2003–2018</th><th colspan="2">2003–2018</th><th colspan="2">(Adopted—RV Only) 2003–2018</th></tr><tr><th>Prior</th><th>Posterior</th><th>Prior</th><th>Posterior</th><th>Prior</th><th>Posterior</th><th>Prior</th><th>Posterior</th><th>Prior</th><th>Posterior</th></tr></thead><tbody><tr><td>P<sub>b</sub> [yr]<sup>a</sup></td><td>N(20.3, 1.1)</td><td>19.4 ± 1.0</td><td>Fixed to 20.3</td><td>N(22, 4)</td><td>18.4<sup>+1.1</sup><sub>-1.0</sub></td><td>See C</td><td>21.1 ± 1.8</td><td>See C</td><td>19.6<sup>+2.4</sup><sub>-1.5</sub></td></tr><tr><td>T<sub>c,b</sub> [JD-2 450 000]</td><td>N(8010, 15)</td><td>8012 ± 15</td><td>Fixed to 8010</td><td>See A</td><td>8013 ± 15</td><td>See A</td><td>8013 ± 15</td><td>See A</td><td>8012 ± 15</td></tr><tr><td>h<sub>b</sub> = √e<sub>b</sub> cos ω<sub>b</sub><sup>b</sup></td><td>N(0.09, 0.10)</td><td>0.14 ± 0.10</td><td>Fixed to 0.09</td><td>See A</td><td>0.10 ± 0.09</td><td>N(0.33, 0.05)</td><td>0.32 ± 0.05</td><td>...</td><td>...</td></tr><tr><td>k<sub>b</sub> = √e<sub>b</sub> sin ω<sub>b</sub><sup>b</sup></td><td>N(0.03, 0.04)</td><td>0.03 ± 0.04</td><td>Fixed to 0.03</td><td>See A</td><td>0.03 ± 0.04</td><td>N(0.11, 0.03)</td><td>0.10 ± 0.3</td><td>...</td><td>...</td></tr><tr><td>K<sub>b</sub> [m s<sup>-1</sup>]</td><td>J*(1, 200)</td><td>58<sup>+15</sup><sub>-17</sub></td><td>Fixed to 69.4</td><td>See A</td><td>81 ± 10</td><td>See A</td><td>81<sup>+10</sup><sub>-11</sub></td><td>See A</td><td>82 ± 10</td></tr><tr><td>e<sub>b</sub></td><td>...</td><td>0.021<sup>+0.035</sup><sub>-0.017</sub></td><td>Fixed to 0.010</td><td>...</td><td>0.013<sup>+0.023</sup><sub>-0.010</sub></td><td>...</td><td>0.115<sup>+0.033</sup><sub>-0.028</sub></td><td>U(0.0, 0.5)</td><td>0.062<sup>+0.075</sup><sub>-0.044</sub></td></tr><tr><td>ω<sub>b</sub> [deg]</td><td>...</td><td>13<sup>+31</sup><sub>-15</sub></td><td>Fixed to 20</td><td>...</td><td>16<sup>+49</sup><sub>-21</sub></td><td>...</td><td>17<sup>+6</sup><sub>-5</sub></td><td>N(20, 10)</td><td>16 ± 11</td></tr><tr><td>a<sub>b</sub> (AU)</td><td>...</td><td>8.7 ± 0.3</td><td>Fixed to 9.1</td><td>...</td><td>8.4<sup>+0.4</sup><sub>-0.3</sub></td><td>...</td><td>9.2 ± 0.5</td><td>...</td><td>8.8<sup>+0.7</sup><sub>-0.5</sub></td></tr><tr><td>M<sub>b</sub> [M<sub>Jup</sub>]</td><td>...</td><td>8.0<sup>+2.1</sup><sub>-2.4</sub></td><td>Fixed to 10.0</td><td>...</td><td>11.0 ± 1.4</td><td>...</td><td>11.5<sup>+1.6</sup><sub>-1.8</sub></td><td>...</td><td>11.4 ± 1.5</td></tr><tr><td>P<sub>c</sub> [yr]<sup>a</sup></td><td>U(2.74, 3.83)</td><td>3.29 ± 0.03</td><td>See A</td><td>3.32 ± 0.03</td><td>See A</td><td>3.32<sup>+0.03</sup><sub>-0.04</sub></td><td>See A</td><td>3.32<sup>+0.04</sup><sub>-0.05</sub></td><td>See A</td><td>3.32 ± 0.04</td></tr><tr><td>T<sub>c,c</sub> [JD-2 450 000]</td><td>U(7500, 9000)</td><td>8118<sup>+32</sup><sub>-33</sub></td><td>See A</td><td>8152<sup>+23</sup><sub>-20</sub></td><td>See A</td><td>8144<sup>+29</sup><sub>-28</sub></td><td>See A</td><td>8141<sup>+31</sup><sub>-35</sub></td><td>See A</td><td>8143<sup>+31</sup><sub>-28</sub></td></tr><tr><td>h<sub>c</sub> = √e<sub>c</sub> cos ω<sub>c</sub><sup>b</sup></td><td>U(-1, 1)</td><td>0.11<sup>+0.18</sup><sub>-0.17</sub></td><td>See A</td><td>0.04 ± 0.12</td><td>See A</td><td>0.03 ± 0.14</td><td>See A</td><td>0.04<sup>+0.16</sup><sub>-0.15</sub></td><td>...</td><td>...</td></tr><tr><td>k<sub>c</sub> = √e<sub>c</sub> sin ω<sub>c</sub><sup>b</sup></td><td>U(-1, 1)</td><td>-0.23<sup>+0.27</sup><sub>-0.16</sub></td><td>See A</td><td>-0.38<sup>+0.10</sup><sub>-0.08</sub></td><td>See A</td><td>-0.34<sup>+0.19</sup><sub>-0.12</sub></td><td>See A</td><td>-0.32<sup>+0.28</sup><sub>-0.13</sub></td><td>...</td><td>...</td></tr><tr><td>K<sub>c</sub> [m s<sup>-1</sup>]</td><td>J*(1, 200)</td><td>101 ± 8</td><td>See A</td><td>110 ± 7</td><td>See A</td><td>96 ± 10</td><td>See A</td><td>96<sup>+10</sup><sub>-11</sub></td><td>See A</td><td>96 ± 10</td></tr><tr><td>e<sub>c</sub></td><td>...</td><td>0.111<sup>+0.082</sup><sub>-0.061</sub></td><td>...</td><td>0.161<sup>+0.064</sup><sub>-0.057</sub></td><td>...</td><td>0.141<sup>+0.084</sup><sub>-0.076</sub></td><td>...</td><td>0.138<sup>+0.102</sup><sub>-0.081</sub></td><td>U(0, 0.99)</td><td>0.130<sup>+0.083</sup><sub>-0.080</sub></td></tr><tr><td>ω<sub>c</sub> [deg]</td><td>...</td><td>-64<sup>+75</sup><sub>-34</sub></td><td>...</td><td>-84<sup>+20</sup><sub>-17</sub></td><td>...</td><td>-84<sup>+34</sup><sub>-23</sub></td><td>...</td><td>-81<sup>+67</sup><sub>-25</sub></td><td>U(-200, 160)</td><td>-84<sup>+38</sup><sub>-25</sub></td></tr><tr><td>a<sub>c</sub> [AU]</td><td>...</td><td>2.68 ± 0.03</td><td>...</td><td>2.69 ± 0.02</td><td>...</td><td>2.69 ± 0.03</td><td>...</td><td>2.69 ± 0.03</td><td>...</td><td>2.69 ± 0.03</td></tr><tr><td>M<sub>c</sub> sin i<sub>c</sub> [M<sub>Jup</sub>]</td><td>...</td><td>7.7 ± 0.6</td><td>...</td><td>8.3 ± 0.5</td><td>...</td><td>7.3 ± 0.8</td><td>...</td><td>7.3<sup>+0.8</sup><sub>-0.9</sub></td><td>...</td><td>7.3 ± 0.8</td></tr><tr><td>γ [m s<sup>-1</sup>]</td><td>U(-100, 100)</td><td>-21 ± 9</td><td>See A</td><td>-31 ± 5</td><td>See A</td><td>-26 ± 9</td><td>See A</td><td>-33<sup>+12</sup><sub>-11</sub></td><td>See A</td><td>-30 ± 10</td></tr><tr><td>BIC<sup>c</sup></td><td></td><td>383</td><td></td><td>367</td><td></td><td>496</td><td></td><td>500</td><td></td><td>498</td></tr></tbody></table>

Notes. $J(a, b)$ denotes a Jeffreys prior between $a$ and $b$, $N(\mu, \sigma)$ denotes a Gaussian prior, $J^*(k, m)$ denotes a modified Jeffreys prior that behaves as a uniform below $k$ and as a Jeffreys prior between $k$ and $m$.

a Periods are given in years for convenience, but the units used for the analysis were days.

b When fitting $h_j$ and $k_j$, we also always ensure that $e_j = h_j^2 + k_j^2 < 1$.

c BIC values can only be compared when the two models fit the same values of the independent variable, i.e., Cases A and B should be compared separately from Cases C, D, and E.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Posterior distribution for the minimum mass of β Pic c, the mass, the period, and the eccentricity of β Pic b from the fit on RV data only. In Case B, only $M_c \sin i_c$ was fitted among these parameters. Fixing the orbit and mass of β Pic b yields a greater, more precise mass for β Pic c (Case B). Using post-2008 data only yields less precise constraints on the mass of β Pic b with a small bimodality near $M_b = 0 M_{\text{Jup}}$. Including pre-2008 data refines the constraints on the mass of β Pic b, but its eccentricity and period are sensitive to the choice of priors. When a wide Gaussian prior is used on $P_b$ and a uniform prior is used on $e_b$ (Case E), low eccentricity and period are favored.</figcaption></figure>

(BIC) for each case to facilitate model comparison. Here, we provide an overview of the main results.

As mentioned above, Lagrange et al. (2019b) used all of the RV residuals for each group between 2008 and 2018 to constrain the mass of β Pic b. When only the multisine offsets are used, the mass of β Pic b is not constrained and has to be held fixed. In Case A, thanks to the increased precision provided by the SHO GP, we constrain the orbital parameters and masses of both planets using the 28 RV offsets only. We find $M_b = 8.0_{-2.4}^{+2.1} M_{\text{Jup}}$ and $M_c = 7.7 \pm 0.6 M_{\text{Jup}}$. These masses are both lower than those found by Lagrange et al. (2019b). As shown in Figure 6, the posterior distribution of $M_b$ includes a small bimodality near $0 M_{\text{Jup}}$. The results from this fit should therefore be interpreted carefully.

In Case B, we obtain results consistent with those of Lagrange et al. (2019b) and Nielsen et al. (2020) for the orbital parameters and the mass of β Pic c. The main difference is with the eccentricity, with $e = 0.161_{-0.057}^{+0.064}$ being lower than values from the literature but still consistent within $2\sigma$. This result suggests that the GP RV offsets might favor slightly lower eccentricities than multisine offsets from Lagrange et al. (2019b). The BIC for Cases A and B yield $\Delta\text{BIC}_{A-B} = 16$, providing very strong evidence in favor of Case B. This is mostly attributable to the difference in the number of parameters. The fact that allowing more free parameters does not yield better agreement with the data is another sign that parameters of β Pic b are not totally well constrained by the 2008–2018 GP offsets only and that results from Case A should be interpreted with caution. However, the fact that they enable constraints on the signal of β Pic b with RV offsets only does show that the smaller error bars of the GP offsets provide more information on the planetary orbits compared to the multisine offsets.

In Cases C and D, we first tried using uniform priors for the orbit of β Pic b, but several parameters were not constrained properly. Hence, as described above, we use Gaussian priors for all parameters of β Pic b except its mass. Performing a two-planet fit, we find $P_b = 18.4_{-1.0}^{+1.1}$ yr for Case C and $P_b = 21.1 \pm 1.8$ yr for Case D. Both cases seem to favor lower periods than the references from which the priors were taken, but are consistent with their respective reference within $2\sigma$. As shown in Figure 6, the eccentricity is also sensitive to the choice of priors with Case C allowing circular orbits and Case D favoring slightly eccentric orbits. The masses, on the other hand, show low sensitivity to the choice of priors. Case C

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Example fit using RV data only. All RV offsets between 2003 and 2018 are shown in black. The fit of β Pic b and c combined is shown as a solid red line, and the corresponding solution for β Pic b only is shown in blue. Both are shown along with 100 random MCMC samples. This fit corresponds to Case E in Table 1.</figcaption></figure>

yields [$M_b$, $M_c$] = [11.0 ± 1.4, 7.3 ± 0.8] $M_{\text{Jup}}$, while Case D gives [$M_b$, $M_c$] = [11.5$_{-1.8}^{+1.6}$, 7.3$_{-0.9}^{+0.8}$] $M_{\text{Jup}}$. We note that the BIC values shown in Table 1 provides evidence in favor of Case C with priors from Lagrange et al. (2019a).

In Case E, we find $P_b = 19.6_{-1.5}^{+2.4}$ yr and $e_b = 0.062_{-0.045}^{+0.072}$. Therefore, the RV data favor low periods and circular orbits when these two parameters are not constrained ($e_b$) or weakly constrained ($P_b$) with priors from direct imaging. Constraints on both masses remain unaffected by this change of parameterization and priors (see Figure 6). Figure 7 shows the RV fit corresponding to Case E, and Figure 8 shows the posterior distribution of the fitted parameters. The distributions are generally well constrained. The apparent correlation between $e_b$ and $P_b$ is consistent with results from the literature where longer periods are generally associated with more eccentric orbits (Dupuy et al. 2019; Lagrange et al. 2019a; Nielsen et al. 2020). We note that the BIC favors Case E to Case D, but not Case C. In both cases, the evidence is only marginal with $\Delta\text{BIC} = 2$. This suggests that priors from astrometry can still provide valuable information about the orbit of β Pic b despite the better constraints provided by the GP offsets and the extended RV baseline.

Throughout the five cases, β Pic c is always detected with a false-alarm probability $\text{FAP} < 2.2 \times 10^{-6}$, and its orbital parameters and its mass show very low sensitivity to the choice of orbital priors for β Pic b. The constraints on $M_b$ are
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Posterior distribution of the Case E RV fit. The dashed lines and gray shaded areas show the median and the 1σ interval for each parameter. Low-eccentricity orbits are preferred for β Pic b. This corresponds to shorter orbital periods due to the correlation between these two parameters.</figcaption></figure>

also unaffected by changes in these priors. As mentioned above, the improved precision of the GP offset plays an important role in our ability to constrain $M_b$ in Case A. However, one can clearly see in Figure 6 that the additional 2003–2008 RV coverage improves the mass estimate, yielding $M_b = 11.4 \pm 1.5 \ M_{\text{Jup}}$ (Case E) and completely ruling out solutions where b is undetected. This mass is consistent with the indicative mass of $M_b \sim 10 \ M_{\text{Jup}}$ reported by Lagrange et al. (2019b). Thanks to the extended RV baseline, we are now able to constrain the mass of both planets from RV data only

with very good precision and low sensitivity to the choice of priors.

## 4.2. RV Offsets and Relative Astrometry

To further assess the impact of the extended RV baseline and of the GP offsets, and to better characterize the system as whole, we perform a joint fit of all available RV and relative astrometry data. To do so, we still use *radvel*, but we include custom models and likelihoods to account for relative astrometry. We still rely on the Keplerian solver from *radvel*,
<table><caption>**Table 2**
Two-planet Joint Fit of RV and Relative Astrometry</caption><thead><tr><th></th><th colspan="2">Joint Fit</th><th colspan="2">Joint Fit With Gravity<br>(Adopted – Final)</th></tr><tr><th></th><th>Prior</th><th>Posterior</th><th>Prior</th><th>Posterior</th></tr></thead><tbody><tr><td>P<sub>b</sub> [yr]</td><td>...</td><td>21.2<sup>+2.4</sup><sub>-0.9</sub></td><td>...</td><td>21.1<sup>+2.0</sup><sub>-0.8</sub></td></tr><tr><td>T<sub>c,b</sub> [JD-2,450,000]</td><td>U(7000, 9000)</td><td>8012.2 ± 2.3</td><td>U(7000, 9000)</td><td>8011.2 ± 0.5</td></tr><tr><td>h<sub>b</sub> = √e<sub>c</sub> cos ω<sub>b</sub><sup>a</sup></td><td>U(-1, 1)</td><td>0.156<sup>+0.137</sup><sub>-0.175</sub></td><td>U(-1, 1)</td><td>0.146<sup>+0.123</sup><sub>-0.160</sub></td></tr><tr><td>k<sub>b</sub> = √e<sub>c</sub> sin ω<sub>b</sub><sup>a</sup></td><td>U(-1, 1)</td><td>0.074<sup>+0.061</sup><sub>-0.101</sub></td><td>U(-1, 1)</td><td>0.077<sup>+0.058</sup><sub>-0.098</sub></td></tr><tr><td>e<sub>b</sub></td><td>...</td><td>0.033<sup>+0.069</sup><sub>-0.027</sub></td><td>...</td><td>0.029<sup>+0.061</sup><sub>-0.024</sub></td></tr><tr><td>ω<sub>b</sub> [deg]</td><td>...</td><td>22<sup>+18</sup><sub>-67</sub></td><td>...</td><td>25<sup>+15</sup><sub>-62</sub></td></tr><tr><td>a<sub>b</sub> [AU]</td><td>J(1, 100)</td><td>9.3<sup>+0.7</sup><sub>-0.3</sub></td><td>J(1, 100)</td><td>9.2<sup>+0.6</sup><sub>-0.2</sub></td></tr><tr><td>M<sub>b</sub> [M<sub>Jup</sub>]</td><td>U(1, 20)</td><td>11.8 ± 1.4</td><td>U(1, 20)</td><td>11.7 ± 1.4</td></tr><tr><td>Ω<sub>b</sub> [deg]</td><td>U(18, 90)</td><td>32.07 ± 0.08</td><td>U(18, 90)</td><td>32.00<sup>+0.06</sup><sub>-0.05</sub></td></tr><tr><td>i<sub>b</sub> [deg]</td><td>sin i</td><td>88.97 ± 0.09</td><td>sin i</td><td>88.88<sup>+0.04</sup><sub>-0.03</sub></td></tr><tr><td>P<sub>c</sub> [yr]<sup>a</sup></td><td>...</td><td>3.36 ± 0.03</td><td>...</td><td>3.36 ± 0.03</td></tr><tr><td>T<sub>c,c</sub> [JD-2,450,000]</td><td>U(7500, 9000)</td><td>8177<sup>+31</sup><sub>-27</sub></td><td>U(7500, 9000)</td><td>8177<sup>+31</sup><sub>-27</sub></td></tr><tr><td>h<sub>c</sub> = √e<sub>c</sub> cos ω<sub>c</sub><sup>a</sup></td><td>U(-1, 1)</td><td>0.024<sup>+0.102</sup><sub>-0.103</sub></td><td>U(-1, 1)</td><td>0.025 ± 103</td></tr><tr><td>k<sub>c</sub> = √e<sub>c</sub> sin ω<sub>c</sub><sup>a</sup></td><td>U(-1, 1)</td><td>-0.442<sup>+0.085</sup><sub>-0.080</sub></td><td>U(-1, 1)</td><td>-0.442<sup>+0.086</sup><sub>-0.081</sub></td></tr><tr><td>e<sub>c</sub></td><td>...</td><td>0.206<sup>+0.074</sup><sub>-0.063</sub></td><td>...</td><td>0.206<sup>+0.074</sup><sub>-0.063</sub></td></tr><tr><td>ω<sub>c</sub> [deg]</td><td>...</td><td>-87<sup>+14</sup><sub>-13</sub></td><td>...</td><td>-87<sup>+14</sup><sub>-13</sub></td></tr><tr><td>a<sub>c</sub> [AU]</td><td>J(2, 3.5)</td><td>2.71 ± 0.02</td><td>J(2, 3.5)</td><td>2.71 ± 0.02</td></tr><tr><td>M<sub>c</sub> sin i<sub>c</sub> [M<sub>Jup</sub>]</td><td>U(1, 20)</td><td>8.5 ± 0.5</td><td>U(1, 20)</td><td>8.5 ± 0.5</td></tr><tr><td>M<sub>c</sub> sin i<sub>c</sub> [M<sub>c</sub>]</td><td>U(1.4, 2.0)</td><td>1.76 ± 0.03</td><td>U(1.4, 2.0)</td><td>1.76 ± 0.03</td></tr><tr><td>π [mas]</td><td>N(51.44, 0.12)</td><td>51.44 ± 0.12</td><td>N(51.44, 0.12)</td><td>51.44 ± 0.12</td></tr><tr><td>ρ<sub>S</sub>/ρ<sub>G</sub></td><td>U(0.5, 1.5)</td><td>0.99 ± 0.6</td><td>U(0.5, 1.5)</td><td>0.99 ± 0.6</td></tr><tr><td>θ<sub>S</sub> – θ<sub>G</sub> [deg]</td><td>N(0.0, 0.5)</td><td>-0.14 ± 0.48</td><td>N(0.0, 0.5)</td><td>-0.15 ± 0.48</td></tr><tr><td>γ [m s<sup>-1</sup>]</td><td>U(-100, 100)</td><td>-42<sup>+7</sup><sub>-8</sub></td><td>U(-100, 100)</td><td>-42 ± 7</td></tr></tbody></table>

**Note.**

<sup>a</sup> When fitting $h_j$ and $k_j$, we also always ensure that $e_j = h_j^2 + k_j^2 < 1$.

but we then derive values of separation (sep.) and position angle (PA) using a model similar to the one used in **orbitize!** (Blunt et al. 2020).

In this joint model, we fit the following orbital parameters for both planets: semimajor axis *a*, time of conjunction *T*<sub>c</sub>, *h* = √e<sub>c</sub> cos ω, and *k* = √e<sub>c</sub> sin ω. To account for the relative astrometry of β Pic b, we also adjust its inclination *i*<sub>b</sub>, the position angle of ascending node Ω<sub>b</sub>, the parallax π, and the total mass M<sub>Tot</sub>. Because the inclination of β Pic b is constrained by relative astrometry, we fit for its mass, M<sub>b</sub>, directly. This is not the case for β Pic c, so we use M<sub>c</sub> sin i<sub>c</sub>. As in Section 4.1, we also include a global RV offset γ. Finally, as pointed out by Nielsen et al. (2020), there is evidence for a systematic position angle offset between GPI and SPHERE. Following their analysis, we include two additional free parameters: a multiplicative correction, ρ<sub>S</sub>/ρ<sub>G</sub>, in sep. and an additive offset in PA, θ<sub>S</sub> – θ<sub>G</sub>, both applied to GPI data points. This adds up to a total of 17 free parameters. We use uniform or log-uniform priors on all parameters except the parallax for which we use a Gaussian prior corresponding to π = 51.44 ± 12 mas, as measured by Hipparcos (van Leeuwen 2007). For Ω<sub>b</sub>, we use a uniform prior between π/10 and π/2 to account for the planetary RV measurement by Snellen et al. (2014). For the inclination, we use a prior uniform in sin(i<sub>b</sub>). A detailed list of adopted priors is provided in Table 2. We perform two joint fits: one with and one without the VLTI/GRAVITY data point from September 2018. The posterior distribution of all parameters for both fits is shown in Appendix A.

We perform a first fit using all available data except the GRAVITY measurement from 2018. In general, values are consistent with the RV-only fit, but have better precision. While constraints on the mass of β Pic b barely change with the inclusion of relative astrometry, yielding M<sub>b</sub> = 11.8 ± 1.4 M<sub>Jup</sub>, the mass of β Pic c is now measured at M<sub>c</sub> = 8.5 ± 0.5 M<sub>Jup</sub>, an increase of about 1 M<sub>Jup</sub> compared to the RV fits of Section 4.1. This constraint on M<sub>c</sub> is consistent with the literature (Lagrange et al. 2019b; Nielsen et al. 2020). However, accounting for β Pic c does not cause a decrease in M<sub>b</sub> as reported by Nielsen et al. (2020, ~8 M<sub>Jup</sub>) and, to a lesser extent, Lagrange et al. (2019b, ~10 M<sub>Jup</sub>). This is likely due to the extended RV baseline, to the increased precision of the GP offsets compared to the multisine offsets, to the fact that we do not include absolute astrometry, or to a combination of these three factors. The tendency toward short-period (20–21 yr), low-eccentricity orbits observed when fitting RV data only is confirmed in the joint fit, and both parameters are constrained more precisely, as shown in Figure 9. We note that the measured eccentricity of β Pic c is now slightly higher at e<sub>c</sub> = 0.206<sup>+0.074</sup><sub>-0.063</sub> and is in better agreement with the literature (Lagrange et al. 2019b; Nielsen et al. 2020). Similarly to Nielsen et al. (2020), we find no evidence of a systematic offset ρ<sub>G</sub>/ρ<sub>S</sub> in separation. In PA, however, the correction stays unconstrained after the MCMC has reached convergence. In an attempt to better constrain this parameter, we use a wide Gaussian prior centered at 0° with σ = 0.5. This prior includes the correction of -0.47 ± 0.14 found by Nielsen et al. (2020).
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Posterior distribution for the minimum mass of β Pic c, the mass, the period and the eccentricity of β Pic b. The parameter distributions are shown for a fit of RV data only (Case E), a joint fit of RV data and relative astrometry, and the same joint fit, but including the high-precision measurement from VLTI/GRAVITY. The joint fit generally provides better precision than the RV-only fit. Low-eccentricity and short-period orbits are also more strongly favored. The inclusion of the GRAVITY measurement slightly improves the precision of most parameters but does not cause a dramatic change in any of them.</figcaption></figure>

within $1\sigma$. We then find a correction $\theta_G - \theta_S = -0.14 \pm 0.48$. This corresponds to a true north offset as reported by Nielsen et al. (2020), but the offset is smaller and is consistent with 0° well within $1\sigma$. In addition, we performed a fit without including any correction and found that the posterior distribution of orbital parameters was unaffected.

When including the GRAVITY data point in our fit, we must add a second custom term to the log-likelihood because the two coordinates (ΔR.A. and ΔDecl.) are correlated. To do so, we first calculate the Pearson correlation coefficient corresponding to the covariance matrix of Equation (1). This allows us to obtain the probability density function (PDF) of a bivariate normal distribution centered on the model values, but with an analytical solution that does not require any matrix inversion, taking advantage of the precomputed correlation coefficient. We then evaluate this PDF at the coordinates of the data point to obtain the likelihood before taking its logarithm and adding it to the other two likelihoods (RV and other relative astrometry). However, the covariance matrix $\Sigma$ from Equation (1), reported by Nowak et al. (2020), is not positive semidefinite. This yields an invalid Pearson correlation coefficient $r < -1$. This issue is likely due to an approximation in the reported covariance matrix. We therefore simply set the off-axis elements of the covariance matrix to $\sigma_{\Delta R.A.}, \Delta Decl. = \sigma_{\Delta Decl.}, \Delta R.A. = -0.00348 \text{ mas}^2$. This small change of 0.00002 mas$^2$ reestablishes the positive semidefiniteness of $\Sigma$ and puts the Pearson coefficient back in the expected interval $-1 < r < 1$.<fnref n="5" /> The principal components of this adjusted covariance matrix give the error bars on the GRAVITY measurement, as shown in the top-left panel Figure 10. Once this fix is applied, we can fit for all available relative astrometry, including the GRAVITY measurement.

As reported in Table 2, the inclusion of the GRAVITY measurement yields results consistent with the joint fit performed without including it. However, in general, the precision is slightly improved when including the GRAVITY measurement, which is more precise than other relative astrometric measurements by more than an order of magnitude. This precision improvement is marginal and can thus hardly be seen by looking at the posterior distribution of parameters, but a careful inspection of Figure 9 does reveal a sharper peak in $P_b$ and a slightly stronger preference for low-eccentricity orbits. Even when using relative astrometry only, Nowak et al. (2020) reported a preference for slightly eccentric orbits ($e_b = 0.15^{+0.05}_{-0.04}$), which we do not observe here.

This could be because they do not account for $\beta$ Pic c in their analysis while we account for both planets using the RV offsets (see Hara et al. 2019 for a detailed discussion of such possible eccentricity biases). Figure 10 shows the joint fit of RV data with all available relative astrometry, including the GRAVITY measurement. The fit is in good agreement with the high-precision GRAVITY measurement, as shown in the top-left panel.

## 5. Discussion

The orbital properties and masses of $\beta$ Pic b and c provide additional knowledge about the full planetary system architecture, including the debris disk, and its evolution with time (see Figure 11). We discuss the impact of our findings on these two aspects and propose future directions to further improve the system characterization.

Our orbital parameters for $\beta$ Pic b confirm that the planet is misaligned with the main disk and consistent with the inner warp (Lagrange et al. 2012b). A direct detection of $\beta$ Pic c would also provide a measurement of the inclination of its orbit, further revealing the inner architecture of the system with respect to the disk. We show the orbit of $\beta$ Pic b to have a very small eccentricity ($<0.06$). This value is still consistent at $1\sigma$ with the planet being responsible for the exocomets ($e > 0.05$; Beust & Morbidelli 1996), but the inclusion of the more eccentric orbit of $\beta$ Pic c should be investigated.

The mass of directly imaged planets cannot be obtained from direct spectrophotometric measurements. Without complementary data like RV or absolute astrometry, the estimate of this parameter solely relies on theoretical models, which predict the temporal evolution of the luminosity of a planet for a variety of masses. The ages of directly imaged planet hosts are usually well constrained (e.g., Mamajek & Bell 2014); the bolometric luminosity of the planets is computed from their spectral energy distributions. However, the evolutionary models are not calibrated at young ages ($<1$ Gyr) and low masses ($<72 M_{Jup}$). The model-based masses of the directly imaged planets might then be incorrect. The independent mass measurement of $\beta$ Pic b, derived from the RV data, offers a rare opportunity to be compared against the predictions of the evolutionary models. In the following, we will compute the predicted mass according to the known age and luminosity of $\beta$ Pic b and see whether the model-derived age is consistent with our dynamical mass.

We considered the classical hot-start evolutionary model grid COND (Baraffe et al. 2003),<fnref n="6" /> the recent hot-start grid

<footnote n="5">⁵ Note that the original covariance matrix does not perfectly match the error bars shown in Figure 2 of Nowak et al. (2020) because both were made separately at different moments of the analysis and not updated for publication (J. Wang, private communication).</footnote>

<footnote n="6">⁶ The most recent grid from Baraffe et al. (2015) does not reach a luminosity that is dim enough at the age of the system.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Joint fit of RV offsets with all available relative astrometry of β Pic b. In all panels, the orbit of β Pic b is shown in blue. In the top-left panel, it is shown along with all relative astrometry measurements. The VLTI/GRAVITY measurement appears in an inset to show how its high precision constrains the fit. All other relative astrometry measurements were fit in separation and position angle rather than RA and DEC, as shown in the top-right panel. The bottom panel shows the 36 RV offsets along with the joint signal of both planets in addition to the signal or β Pic b alone. All fits are shown with 100 random samples from the MCMC.</figcaption></figure>

SONORA (Marley et al. 2018), and the warm-start grid from Spiegel & Burrows (2012). We use grids computed at solar metallicity because the metallicity of β Pic b is not well constrained from existing spectrophotometric data, with estimated values from −0.5 up to 0.63 (Nowak et al. 2020), depending on the data and analysis method. The distributions of the predicted mass for all three evolutionary models were derived with a Monte Carlo process from the normal distributions of the most recent luminosity measurement for β Pic b of $log L/L_\odot = -3.76 \pm 0.02$ (Chilcote et al. 2017) and from the age of 25 ± 3 Myr (Messina et al. 2016; this age

assumes β Pic b formed very quickly after the birth of the host star). The grid luminosity was first interpolated in $log t$, at each randomly selected age, for each mass bin. The log *M* was then interpolated in $log L_{bol}$ at each randomly selected luminosity to derive the model-dependent mass. The process was repeated 10⁵ times to derive the distribution of the model-derived mass. Figure 12 shows the distributions for the three models, compared to our dynamical mass measurement. We found $12.9_{-0.3}^{+0.2} M_{\text{Jup}}$, $12.9 \pm 0.1 M_{\text{Jup}}$, and $13.6 \pm 0.1 M_{\text{Jup}}$ for the COND, SONORA, and warm-start model at the highest available entropy, 13 k<sub>B</sub>/baryon, following the analysis of Chilcote et al. (2017).
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Orbit solutions showing the system architecture in various coordinates. The left panel shows the orbit of both planets in RA and DEC, assuming coplanar orbits for both planets. The planetary orbits are slightly inclined from the main disk in position angle (PA) but are almost aligned with the inner disk (Lagrange et al. 2012a). The middle panel shows both orbits as seen from the pole, with the orbit of β Pic b being nearly circular and the orbit of β Pic c slightly eccentric. The separation between both planets and the star is shown in the left panel. In all three panels, the orbit of β Pic b is shown in blue and that of β Pic c is in red, along with 100 random MCMC realizations. The main uncertainty in the trajectory of β Pic b appears to come from a lack of data in the northeastern part of its orbit, leading to a large range of possible orbital periods.</figcaption></figure>

The results are still consistent with Morzinski et al. (2015), Chilcote et al. (2017), and Nielsen et al. (2020), who used different age and/or luminosity estimates ($23 \pm 3$ Myr, $24 \pm 3$ Myr, and $26 \pm 3$ Myr, respectively, and $\log L_{bol}/L_{\odot} = -3.78 \pm 0.03$). Mostly driven by the tight constraint on the bolometric luminosity, the model-derived masses exhibit a precision 10 times better than our dynamical mass measurement. Yet, the two cold-start model grids are consistent to better than $1\sigma$ with our dynamical mass measurement, and the warm-start model from Spiegel & Burrows (2012) is consistent within $2\sigma$ at the measured age of the system, favoring formation with the highest initial entropy, following the same conclusion reached by Nielsen et al. (2020). Our precision is not high enough to either precisely calibrate the evolutionary tracks or to discriminate which of the hot-start or warm-start scenarios best fit the parameters of β Pic b, or to derive the age of the system assuming the models are accurate.

Alternatively, one can assume the evolutionary models predict correct masses for a given luminosity–age pair. Under such hypothesis, one can conversely use our dynamical mass measurement to determine any significant formation delay of β Pic b, with respect to A, according to its measured luminosity (Currie et al. 2009). However, the large uncertainty on the dynamical mass only translates into an upper limit on this delay of 12 Myr at $1\sigma$ for all three evolutionary models, which is more than any timescale for giant planet formation. A precision of $0.1 M_{Jup}$ is required to test this scenario.

Continuous and long-term monitoring of β Pic b via relative and absolute astrometry, especially with the exquisite precision of VLTI/GRAVITY, and spectroscopy will enable the improvement of the orbital parameters and the mass precision to revisit the calibration of the evolutionary models and/or the formation delay hypothesis, following the same conclusions as Nielsen et al. (2020). As seen in Figure 11, a lot of orbital solutions for β Pic b would be excluded by a better constraint on its period. The fact that both relative astrometry and RV data cover less than a full period for β Pic leads to a greater uncertainty in $P_b$. Continuous monitoring of β Pic b in the next 5 yr will strongly constrain the maximum separation of β Pic b in the northeastern half of the orbit, hence enabling a better precision in $P_b$. Furthermore, continuous RV monitoring of β Pic as well as the use of unbiased absolute astrometric measurements from Gaia in combination with the long baseline

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Distributions of the mass of β Pic b from our dynamical measurement (black) compared to predictions from the evolutionary models COND (Baraffe et al. 2003), SONORA (Marley et al. 2018), and the warm-start model with $S = 13 k_B$/baryon from Spiegel & Burrows (2012) assuming an age of $25 \pm 3$ Myr (Messina et al. 2016) and a bolometric luminosity of $\log L/L_{\odot} = -3.76 \pm 0.02$ (Chilcote et al. 2017).</figcaption></figure>

with Hipparcos would not only help to constrain the period of β Pic b, but it would also improve the accuracy and precision of model-independent mass measurements for both planets.

Similarly, a direct detection of β Pic c together with an adequate sampling of its spectral energy distribution would provide another independent luminosity measurement within the same system, at a lower dynamical mass. This would further strengthen the tests of the evolutionary models. Direct imaging of β Pic c may be possible with the James Webb Space Telescope (JWST) or with next-generation giant ground-based telescopes such as the European Extremely Large Telescope (E-ELT), which should provide an angular resolutions of ~100 mas for the former and ~30 mas for the latter. As shown in Figure 11, the separation between β Pic c and its host is greater than 100 mas for a significant fraction of its orbit
according to the solutions derived in this work. Such a
measurement would independently confirm the existence of the
planet while providing another measurement to calibrate mass-
luminosity models.

**6. Conclusion**

We present a new and independent analysis of the RVs of
β Pic. For the first time, we use an SHO GP framework to
model the stellar activity of the star, training the model from the
photometric light curve, which results in the following results:

1. We show that the SHO GP enables the δ Scuti pulsations with amplitudes as high as several hundreds of meters per second to be modeled;

2. It yields better precision than a QP kernel typically used to model stellar rotation and is favored according to the BIC;

3. It yields RV offsets consistent with parametric models from the literature (Lagrange et al. 2019b);

4. It provides a simpler framework than the parametric models, with three hyperparameters instead of up to 90 parameters;

5. It allows poorly sampled RVs, excluded from previous analyses, to be modeled, hence expanding the baseline by nearly 5 yr, as early as 2003.

This extended RV coverage and the precision of the SHO GP
offsets make it possible to constrain the mass of both planets
using RV data only with low sensitivity to the choice of priors
imposed on the orbit of β Pic b. However, the constraints on
other orbital parameters are sensitive to changes in the priors
from orbital fits using relative and absolute astrometry in the
literature (Lagrange et al. 2019a; Nielsen et al. 2020). To better
characterize the entire system, we also perform a fit of all
available relative astrometry and RV data. The joint fits and the
fits with RV data only yield results that are generally consistent.
We find a model-independent mass of 12.7 ± 1.3 M<sub>Jup</sub> for
β Pic b. A low-eccentricity orbit (e = 0.016<sup>+0.044</sup><sub>-0.012</sub>) with a
relatively short period of P<sub>b</sub> = 20.6<sup>+1.5</sup><sub>-0.4</sub> yr is preferred. For
β Pic c, the mass is constrained to 8.5 ± 0.5 M<sub>Jup</sub>, and the orbit
is consistent with values from the literature (Lagrange et al.
2019b; Nielsen et al. 2020).

The dynamical mass measurement of β Pic b presented
in this work is consistent with hot-start and warm-start
evolutionary models at the highest computed initial entropy
(S = 13 k<sub>B</sub>/baryon) but is not precise enough to discriminate
against any of the models tested, reaching conclusions similar
to Nielsen et al. (2020). Further monitoring with relative and
absolute astrometry, as well as RV, will be required to obtain
stronger constraints on the mass, enabling a better calibration of
mass-luminosity models.

A direct detection of β Pic c could be possible with JWST.
With a minimum mass of 8.5 ± 0.5 M<sub>Jup</sub>, the COND
evolutionary model predicts a contrast with the star around
10.7 mag and 8.9 mag, at the *K* band and *L* band respectively.
Such a detection would provide a second point for the
calibration of the evolutionary model while further improving
the characterization of orbits in the β Pic system.

This work demonstrates that precise modeling of complex
stellar activity with SHO GP is effective. It opens a new path
toward the detection of giant planets with RV only around
young stars as well as the characterization of systems with
directly imaged planets beyond that of β Pic.

We sincerely thank the referee for the thoughtful and careful
review and suggestions, leading to an improvement of the
present work. We thank Anne-Marie Lagrange for providing all
published individual RV measurements and Djamel Mékarnia
for the ASTEP light curve. We thank Lauren Weiss for her
contribution to early phases of this project and Aurélien
Wyttenbach and Raphaëlle Haywood for fruitful discussions
about GP. T.V. acknowledges support from an Institut de
Recherche sur les Exoplanètes (iREx) Trottier Excellence
Grant. J.R. and R.D. acknowledge support from Fonds de
Recherche du Québec as J.R.'s work was performed in part
under contract with the University of Montréal. J.R. is
supported by the French National Research Agency in the
framework of the Investissements d'Avenir program (ANR-15-
IDEX-02), through the funding of the "Origin of Life" project
of the Univ. Grenoble-Alpes."

*Software:* emcee (Foreman-Mackey et al. 2013), celerite (Foreman-Mackey et al. 2017), RadVel (Fulton et al. 2018), Astropy (The Astropy Collaboration et al. 2013), corner (Foreman-Mackey 2016), orbitize! (Blunt et al. 2020).
## Appendix A
### Posterior Distributions for all Fits

Figures 13–16 show the parameter posterior distributions for RV Cases A to D. Figures 17 and 18 show the parameter posterior distribution for the two joint fits of RV and relative astrometry (with and without the VLTI/GRAVITY measurement, respectively).

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Posterior distribution of the Case A RV fit. The dashed lines and gray shaded areas show the median and the 1σ interval for each parameter. A small bimodality in *M*<sub>b</sub> suggests that the orbit of β Pic b is not totally well constrained by the 2008–2018 RV offsets only.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Posterior distribution of the Case B RV fit. The dashed lines and gray shaded areas show the median and the 1σ interval for each parameter. The orbit of β Pic b is fixed to values from Lagrange et al. (2019a).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Posterior distribution of the Case C RV fit. The dashed lines and gray shaded areas show the median and the 1σ interval for each parameter. Gaussian priors from Lagrange et al. (2019a) were applied to the orbital parameter of β Pic b, but with a wide prior of 22 ± 4 yr on $P_b$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Posterior distribution of the Case D RV fit. The dashed lines and gray shaded areas show the median and the 1σ interval for each parameter. Gaussian priors from Nielsen et al. (2020) were applied to the orbital parameter of β Pic b, but with a wide prior of 22 ± 4 yr on $P_b$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Posterior distribution for the joint fit of RV and relative astrometry, without the VLTI/GRAVITY measurement of September 2008. The dashed lines and gray shaded areas show the median and the 1σ interval for each parameter. We note very strong correlations between $a_b$, $h_b = \sqrt{e_b} \cos \omega_b$, and $k_b = \sqrt{e_b} \sin \omega_b$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Posterior distribution for the joint fit of RV and relative astrometry, with the VLTI/GRAVITY measurement of 2008 September. The dashed lines and gray shaded areas show the median and 1σ interval for each parameter. As in Figure 17, $a_b$, $h_b = \sqrt{e_b} \cos \omega_b$, and $k_b = \sqrt{e_b} \sin \omega_b$ are strongly correlated.</figcaption></figure>
<table><caption>**Table 3**
RV Offsets from 2003 to 2018</caption><thead><tr><th>MJD</th><th>RV<br>[m s<sup>-1</sup>]</th><th>σ<sub>RV</sub><br>[m s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>52944.376672</td><td>-77.43</td><td>59.87</td></tr><tr><td>52993.188574</td><td>39.68</td><td>28.94</td></tr><tr><td>53048.602259</td><td>-154.64</td><td>34.06</td></tr><tr><td>53775.641184</td><td>-9.07</td><td>149.54</td></tr><tr><td>54056.325773</td><td>20.93</td><td>106.71</td></tr><tr><td>54112.579919</td><td>-99.04</td><td>247.39</td></tr><tr><td>54441.746852</td><td>-191.08</td><td>148.94</td></tr><tr><td>54483.110698</td><td>-92.49</td><td>135.89</td></tr><tr><td>54544.046717</td><td>-76.74</td><td>16.72</td></tr><tr><td>54780.215069</td><td>-70.24</td><td>24.16</td></tr><tr><td>54799.213508</td><td>-71.54</td><td>25.30</td></tr><tr><td>54829.165935</td><td>-202.14</td><td>27.49</td></tr><tr><td>54850.142609</td><td>-130.78</td><td>22.24</td></tr><tr><td>54913.040806</td><td>-84.29</td><td>18.44</td></tr><tr><td>55131.298795</td><td>2.44</td><td>26.38</td></tr><tr><td>55170.202032</td><td>18.54</td><td>22.49</td></tr><tr><td>55519.296681</td><td>37.59</td><td>30.68</td></tr><tr><td>55566.194698</td><td>148.72</td><td>14.54</td></tr><tr><td>55597.126002</td><td>85.80</td><td>19.29</td></tr><tr><td>56569.311461</td><td>129.29</td><td>41.54</td></tr><tr><td>56582.333626</td><td>95.49</td><td>28.91</td></tr><tr><td>56685.203668</td><td>24.62</td><td>32.80</td></tr><tr><td>56694.067479</td><td>91.60</td><td>27.57</td></tr><tr><td>56706.582371</td><td>198.96</td><td>23.64</td></tr><tr><td>56984.282332</td><td>-15.44</td><td>20.36</td></tr><tr><td>57344.247022</td><td>-148.18</td><td>55.92</td></tr><tr><td>57667.303333</td><td>106.83</td><td>22.29</td></tr><tr><td>57712.283914</td><td>128.33</td><td>21.42</td></tr><tr><td>57771.166012</td><td>81.25</td><td>8.54</td></tr><tr><td>57849.014450</td><td>151.55</td><td>25.00</td></tr><tr><td>58007.347048</td><td>3.42</td><td>21.89</td></tr><tr><td>58037.271945</td><td>37.74</td><td>7.34</td></tr><tr><td>58064.265248</td><td>-62.63</td><td>24.20</td></tr><tr><td>58094.192294</td><td>-35.32</td><td>6.94</td></tr><tr><td>58207.056451</td><td>-32.82</td><td>11.04</td></tr><tr><td>58235.001291</td><td>-70.57</td><td>7.73</td></tr></tbody></table>

**Appendix B**
Table With Newly Derived RV Offset + Error

Table 3 shows the RV offsets from 2003 to 2018 obtained by filtering stellar activity with the SHO GP model, along with the uncertainty on each offset.

**ORCID iDs**

Thomas Vandal https://orcid.org/0000-0002-5922-8267

Julien Rameau https://orcid.org/0000-0003-0029-0258

René Doyon https://orcid.org/0000-0001-5485-4675

**References**

Ahmic, M., Croll, B., & Artymowicz, P. 2009, ApJ, 705, 529

Angus, R., Morton, T., Aigrain, S., Foreman-Mackey, D., & Rajpaul, V. 2018, MNRAS, 474, 2094

Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

Bell, C. P. M., Rees, J. M., Naylor, T., et al. 2014, MNRAS, 445, 3496

Beust, H., Lagrange, A.-M., Crawford, I. A., et al. 1998, A&A, 338, 1015

Beust, H., Lagrange-Henri, A. M., Madjar, A. V., & Ferlet, R. 1990, A&A, 236, 202

Beust, H., & Morbidelli, A. 1996, Icar, 120, 358

Blunt, S., Wang, J. J., Angelo, I., et al. 2020, AJ, 159, 89

Bowler, B. P., Liu, M. C., Shkolnik, E. L., & Tamura, M. 2014, ApJS, 216, 7

Brown, A. G. A., Vallenari, A., Prusti, T., et al. 2018, A&A, 616, A1

Burrows, C. J., Krist, J. E., Stapelfeldt, K. R. & WFPC2 Investigation Definition Team 1995, AAS Meeting, 187, 32.05

Chauvin, G., Lagrange, A.-M., Beust, H., et al. 2012, A&A, 542, A41

Chilcote, J., Pueyo, L., Rosa, R. J. D., et al. 2017, AJ, 153, 182

Cloutier, R., Astudillo-Defru, N., Doyon, R., et al. 2017a, A&A, 608, A35

Cloutier, R., Doyon, R., Menou, K., et al. 2017b, AJ, 153, 9

Crepp, J. R., Johnson, J. A., Howard, A. W., et al. 2012, ApJ, 761, 39

Crouzet, N., Chapellier, E., Guillot, T., et al. 2018, A&A, 619, A116

Currie, T., Lada, C. J., Plavchan, P., et al. 2009, ApJ, 698, 1

Currie, T., Thalmann, C., Matsumura, S., et al. 2011, ApJL, 736, L33

Dent, W. R. F., Thi, W. F., Kamp, I., et al. 2013, PASP, 125, 477

Dupuy, T. J., Brandt, T. D., Kratter, K. M., & Bowler, B. P. 2019, ApJL, 871, L4

Dupuy, T. J., & Liu, M. C. 2017, ApJS, 231, 15

Ford, E. B. 2006, ApJ, 642, 505

Foreman-Mackey, D. 2016, JOSS, 1, 24

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504

Galland, F., Lagrange, A.-M., Udry, S., et al. 2005, A&A, 443, 337

Golimowski, D. A., Ardila, D. R., Krist, J. E., et al. 2006, AJ, 131, 3109

Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006, AJ, 132, 161

Grunblatt, S. K., Howard, A. W., & Haywood, R. D. 2015, ApJ, 808, 127

Hara, N. C., Boué, G., Laskar, J., Delisle, J. B., & Unger, N. 2019, MNRAS, 489, 738

Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517

Hillenbrand, L., Isaacson, H., Marcy, G., et al. 2015, in Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, Vol. 18, ed. G. van Belle & H. C. Harris (Flagstaff, AZ: Lowell Observatory), 759

Kiefer, F., des Etangs, A. L., Boissier, J., et al. 2014, Natur, 514, 462

Koen, C. 2003, MNRAS, 341, 1385

Lagrange, A.-M., Boccaletti, A., Langlois, M., et al. 2019a, A&A, 621, L8

Lagrange, A.-M., Boccaletti, A., Milli, J., et al. 2012a, A&A, 542, A40

Lagrange, A.-M., Bonnefoy, M., Chauvin, G., et al. 2010, Sci, 329, 57

Lagrange, A.-M., De Bondt, K., Meunier, N., et al. 2012b, A&A, 542, A18

Lagrange, A. M., Meunier, N., Chauvin, G., et al. 2013, A&A, 559, A83

Lagrange, A.-M., Meunier, N., Rubini, P., et al. 2019b, NatAs, 3, 1135

Lagrange, A.-M., Plazy, F., Beust, H., et al. 1996, A&A, 310, 547

Mamajek, E. E., & Bell, C. P. M. 2014, MNRAS, 445, 2169

Marley, M., Didier, S., Morley, C., & Fortney, J. 2018, Sonora 2018: Cloud-free, Solar Composition, Solar C/O Substellar Atmosphere Models and Spectra, Zenodo, doi:10.5281/zenodo.1309034

Mayor, M., Pepe, F., Queloz, D., et al. 2003, Msnr, 114, 20

Mékarnia, D., Chapellier, E., Guillot, T., et al. 2017, A&A, 608, L6

Messina, S., Lanzafame, A. C., Feiden, G. A., et al. 2016, A&A, 596, A29

Morzinski, K. M., Males, J. R., Skemer, A. J., et al. 2015, ApJ, 815, 108

Mouillet, D., Larwood, J. D., Papaloizou, J. C. B., & Lagrange, A. M. 1997, MNRAS, 292, 896

Nielsen, E. L., Liu, M. C., Wahhaj, Z., et al. 2014, ApJ, 794, 158

Nielsen, E. L., Rosa, R. J. D., Wang, J., et al. 2016, AJ, 152, 175

Nielsen, E. L., Rosa, R. J. D., Wang, J. J., et al. 2020, AJ, 159, 71

Nowak, M., Lacour, S., Mollière, P., et al. 2020, A&A, 633, A110

Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, MNRAS, 452, 2269

Ryu, T., Sato, B., Kuzuhara, M., et al. 2016, ApJ, 825, 127

Smith, B. A., & Terrile, R. J. 1984, Sci, 226, 1421

Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Natur, 509, 63

Snellen, I. A. G., & Brown, A. G. A. 2018, NatAs, 2, 883

Soderblom, D. R., Hillenbrand, L. A., Jeffries, R. D., Mamajek, E. E., & Naylor, T. 2014, in Protostars and Planets VI (219) ed. H. Beuther et al. (Tucson, AX: Univ. Arizona Press), 219

Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174

The Astropy Collaboration, Roboitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

van Leeuwen, F. 2007, A&A, 474, 653

Wang, J. R., Graham, J. R., Pueyo, L., et al. 2016, AJ, 152, 97

Zieba, S., Zwintz, K., Kenworthy, M. A., & Kennedy, G. M. 2019, A&A, 625, L13

Zuckerman, B., Song, I., Bessell, M. S., & Webb, R. A. 2001, ApJL, 562, L87

Zwintz, K., Reese, D. R., Neiner, C., et al. 2019, A&A, 627, A28