# Evidence for an additional planet in the β Pictoris system

A.-M. Lagrange, Nadège Meunier, Pascal Rubini, Miriam Keppler, Franck Galland, Eric Chapellier, Eric Michel, Luis Balona, Hervé Beust, Tristan Guillot, et al.

► To cite this version:

A.-M. Lagrange, Nadège Meunier, Pascal Rubini, Miriam Keppler, Franck Galland, et al. Evidence for an additional planet in the β Pictoris system. Nature Astronomy, 2019, 10.1038/s41550-019-0857-1. hal-02327499

HAL Id: hal-02327499

https://hal.science/hal-02327499v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# Evidence for an additional planet in the β Pictoris system

A.-M. Lagrange<fnref n="1" />*, Nadège Meunier<fnref n="1" />, Pascal Rubini D², Miriam Keppler<fnref n="1" />,³, Franck Galland<fnref n="1" />, Eric Chapellier⁴, Eric Michel⁵, Luis Balona⁶, Hervé Beust<fnref n="1" />, Tristan Guillot D⁴, Antoine Grandjean D<fnref n="1" />, Simon Borgniet D<fnref n="1" />, Djamel Mékarnia⁴, Paul Anthony Wilson⁷,⁸,⁹, Flavien Kiefer⁷, Mickael Bonnefoy<fnref n="1" />, Jorge Lillo-Box D¹⁰, Blake Pantoja¹⁰,¹¹, Matias Jones¹⁰, Daniela Paz Iglesias D¹²,¹³, Laetitia Rodet<fnref n="1" />, Matias Diaz D¹¹, Abner Zapata¹⁴, Lyu Abe⁴ and François-Xavier Schmider⁴

With its imaged debris disk of dust, its evaporating exocomets, and an imaged giant planet, the young (~23 Myr) β Pictoris system is a unique proxy for detailed studies of planet formation processes as well as planet-disk interactions. Here, we study ten years of European Southern Observatory/High Accuracy Radial Velocity Planet Searcher (HARPS) high-resolution spectroscopic data of β Pictoris. After removing the radial velocity (RV) signals arising from the δ Scuti pulsations of the star, a ~1,200-d periodic signal remains, which, within our current knowledge, we can only attribute to a second planet in the system. The β Pic c mass is about nine times the mass of Jupiter; it orbits at ~2.7 AU on an eccentric (e ~ 0.24) orbit. More RV data are needed to obtain more precise estimates of the properties of β Pic c. The current modelling of the planet's properties and the dynamic of the whole system has to be reinvestigated in light of this detection.

After 25 years, our view of planetary system architectures has profoundly evolved thanks to the detection of a large variety of exoplanets. In particular, giant planets that orbit between a few stellar radii (detected by means of indirect techniques such as RV variations or transit) and several hundred AU (detected by direct imaging) from their stars have been found. This diversity suggests different ways in which giant planets or low-mass brown dwarfs may form; it also suggests that, after formation, dynamical evolution (disk migration, third-body interaction) plays an important role in determining the final architecture of the planetary systems<fnref n="1" />.

Young systems aged between ~5 and 50 Myr offer a particularly interesting opportunity to study the early stages of planet formation and evolution, when giant planets are formed and most of the protoplanetary gas has been removed from the disk. The ~23 Myr β Pictoris system is an emblematic prototype of such systems. Its resolved debris disk of dust² is attributed to collisions among planetesimals orbiting between ~50 and 100 AU (refs. ³,⁴). Evaporating exocomets⁵ were detected, and a giant planet⁶,⁷, β Pic b, has been imaged orbiting at about 9 AU from the star. Gravitational interaction between β Pic b and the planetesimals, together with stellar radiation pressure, can explain most of the disk morphological properties, especially the prominent inner warp³,⁸,⁹, and the dust outward extension. β Pic b could be responsible for triggering the infall and evaporation of the exocomets onto the star (the so-called falling evaporating body—FEB—scenario¹⁰).

Thermal data¹¹⁻¹³ revealed gaps and knots that suggest additional planets at separations larger than 20 AU. The apparent void of

material within about 50 AU of the star suggests that planets might be present in this region¹⁴. In addition, Atacama Large Millimeter Array (ALMA) observations revealed side asymmetries that could be due to the presence of other planets in the system⁴. Finally, the slightly inclined orbit of β Pic b could be due to gravitational perturbation by another, as yet unseen, planet.

Recently, HARPS high-precision spectroscopic data were combined with the Very Large Telescope Nasmyth adaptive optics and Conica system (NaCo) high-spatial-resolution imaging data to search for additional giant planets on a wide range of separations from 0.01 to 100+ AU from the star¹⁵, bridging for the first time the gap between RV and direct imaging detection limits in the planetary/brown dwarf domain. It was found that companions with masses 4 M<sub>J</sub> or more orbiting closer than 1 AU or further away than 10 AU are not present in the system. The 2–6 AU region is somewhat less constrained, with limits in the 10 M<sub>J</sub> to brown dwarf mass range. Adding more recent extreme-AO Spectro-Polarimetric High-Contrast Exoplanet Research (SPHERE) data improves the detection limit and allows exclusion of companion planets more massive than about 15 M<sub>J</sub>, and with separations of 2–6 AU.

## The data

We analysed 6,645 HARPS high-resolution ($R=120,000$), high-signal-to-noise (SN, mean SN of 273 at 550 nm) spectra of β Pictoris obtained between 2003 and 2018. Observations before 2008 consist of exposures of a few minutes designed to study the variability of Ca II absorption lines attributed to the exocomets mentioned

<footnote n="1">¹Institut de Planétologie et d'Astrophysique de Grenoble, Université Grenoble Alpes, CNRS, IPAG, Grenoble, France. ²Pixyl, La Tronche, France. ³Max Planck Institute for Astronomy, Heidelberg, Germany. ⁴Université Côte d'Azur, OCA, Lagrange CNRS, Nice, France. ⁵LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, Meudon, France. ⁶South African Astronomical Observatory, Cape Town, South Africa. ⁷CNRS, UMR7095, Institut d'Astrophysique de Paris, Paris, France. ⁸Department of Physics, University of Warwick, Coventry, UK. ⁹Leiden Observatory, Leiden University, Leiden, the Netherlands. ¹⁰European Southern Observatory, Vitacura, Santiago, Chile. ¹¹Departamento de Astronomía, Universidad de Chile, Las Condes, Santiago, Chile. ¹²Instituto de Física y Astronomía, Facultad de Ciencias, Universidad de Valparaíso, Valparaíso, Chile. ¹³Núcleo Milenio de Formación Planetaria—NPF, Universidad de Valparaíso, Valparaíso, Chile. ¹⁴Center of Astro-Engineering, Department of Electrical Engineering, Pontificia Universidad Católica de Chile, Santiago, Chile. *e-mail: anne-marie.lagrange@univ-grenoble-alpes.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1** | RV variations of β Pictoris and generalized Lomb-Scargle periodograms. **a**, RV between 2008 and 2018, computed as described in Methods. The measured RV and the associated uncertainties are listed in Supplementary Table 1. The uncertainties vary depending on the SN with an average of 59 m s⁻¹ and a root mean square (r.m.s.) of 12 m s⁻¹. **b**, A sequence of 6 h of observations of β Pictoris obtained in January 2017. The complex, high-frequency RV variations are due to stellar pulsations. **c**, Generalized Lomb-Scargle periodogram of the RV between 2008 and 2018. Note that, because of the pulsations and of the complex window function, peaks above the FAP do not necessarily indicate the presence of a real periodic signal at the peak position (see text). The orbital period of β Pic b, 7,802 d, is indicated as an orange vertical line, and that of β Pic c by a blue one. **d**, Generalized Lomb-Scargle periodogram of the temporal window, showing the absence of a peak due to the observation window around 1,200 d.</figcaption></figure>

above. After the discovery of β Pic b, the monitoring of the star was intensified to constrain the mass of the planet through RV measurements¹⁶,¹⁷. However, high-frequency stellar pulsations with amplitudes up to 1,200 m s⁻¹ dominate the RV variations and represent a major source of noise. Pulsation modes at 47.44, 38.06 and 47.28 cycles d⁻¹ have been reported¹⁸ using photometric data. These pulsations were also seen in the line profile analysis of about 700 spectra obtained over two weeks¹⁹. They indicate that the star is a member of the δ Scuti group. Recently, high-quality data acquired from the Concordia research station in Antarctica during a long-term monitoring of β Pictoris between April and October 2017 allowed the identification of 31 sinusoidal pulsations²⁰. To mitigate the impact of the stellar pulsations on the RV analysis, we adopted after 2008 a specific observing strategy that consisted in monitoring the stellar variations on uninterrupted sequences of 1–2 h and correcting for them. Furthermore, a few longer (~6 h) sequences were obtained in 2017–2018. After appropriate data selection (Methods), the RVs measured after 2008 are plotted in Fig. 1a with an example of the pulsation pattern in Fig. 1b. We note that the uncertainties associated with the RV are higher than in the case of main-sequence solar-type stars, because β Pictoris is an early-type star (A5V), with fewer spectral lines, and a much higher rotational velocity (120 km s⁻¹) than main-sequence solar-type stars, which have rotational velocities in the range 3–5 km s⁻¹.

The periodogram of the RV data (Fig. 1c) reveals a complex set of peaks, most of them well above the 10% false alarm probability (FAP) level. Most of these peaks are due not to noise but to a combination of pulsation modes and of the complex temporal distribution

of our data. The highest peak of the periodogram is at ~1,200 d. It is not associated with β Pic b's orbital period (~7,800 d), We show in Fig. 1d the periodogram for the temporal window replacing our RV data by a linear function of the Julian date (JD). Peaks are present over the entire range of periods, but not at 1,200 d. This shows that, while a large fraction of the peaks in Fig. 1c are due to the uneven temporal distribution of the observations, this is not the case of the 1,200-d peak.

To further study the signal associated with this ~1,200-d peak, the data must be corrected first for the complex, high-frequency pulsation patterns. To do so, different methods described in Methods were used. The first, straightforward approach was to average the RV obtained from continuous sequences of data obtained over 1–2 h or more. To perform such analysis, only the well sampled and sufficiently long data sequences obtained after 2008 were retained. More sophisticated approaches consisted in fitting each of the longest sets of data, once grouped in bins of less than 8 d, with sine series describing the RV variations due to the pulsations plus an offset, either in an iterative way or simultaneously. The time series of the residuals obtained after subtracting the sine functions and the time series of the offsets can then be analysed.

Altogether, the sine fits confirm the presence of multiple pulsations with frequencies in the range 35–65 cycles d⁻¹ as detected in photometry, but also, possibly, weaker pulsations in the 3–4 cycles d⁻¹ range, characteristic of γ Doradus pulsations (see Methods for more details).

The residual RVs obtained after subtracting the pulsation-induced RV variations are shown in Fig. 2. Their periodogram again
<figure><img src="image_3.png" /><figcaption>**Fig. 2** | Pulsation-corrected RV data with corresponding generalized Lomb–Scargle periodogram. **a**, Pulsation-corrected RV data. **b**, Corresponding generalized Lomb–Scargle periodogram. An FAP probability of 10% is indicated. The orange dashed vertical line indicates the period of β Pic b, and the blue one the period of β Pic c.</figcaption></figure>

shows a prominent peak at ~1,200 d. The offsets also exhibit long-term variations with a total amplitude of about 200 m s⁻¹ (Fig. 2), and their periodogram again shows the same peak.

## Possible origins of the residual RV signal

We first ruled out systematic observing errors due to temporal variations in air mass, proximity of the Moon or variations in the SN level because they show no correlation with the corrected RV residuals.

Instrumental systematic effects can also be ruled out because the observed amplitude of over 300 m s⁻¹ would have been detected on other stars observed with the same set-up. Amplitudes as low as 1–10 m s⁻¹ have been obtained for the RV variation of low-mass planets around solar-type stars. Instrumental effects that would only affect rapid rotators would have been detected in our previous surveys of A–F stars.

An artefact due to the temporal sampling was searched for. The probability of a fortuitous origin of the 1,200-d period was estimated to be less than 10⁻³ (Methods).

Stellar pulsations with periods 1–100 d, which are not well sampled with the present data, cannot explain the ~1,200-d peak. We checked that such pulsations would not create a high peak at 1,200 d. Periods of the order of 1,200 d of stellar origin have not been reported for stars similar to β Pictoris (and would be extremely difficult to measure). We note that radial pulsations that would produce RV variations with amplitudes of more than 100 m s⁻¹ and a period of 1,200 d can definitely be excluded, as they would indicate a radius variation of several hundred stellar radii, which is not possible.

The possibility of an artefact due to an incomplete removal of the stellar pulsations was investigated in depth (Methods). Within the current knowledge of hybrid δ Scuti and γ Doradus stars, a stellar origin is not possible.

No correlation between the offsets and the presence of exo-comet absorptions is found. Examples of such events are given in Supplementary Fig. 3, together with the computed offsets at the same epochs.

We therefore conclude that, under our current knowledge, these long-term RV variations must be caused by an additional companion.

## An additional planet

Hereafter, we assume that the RV residuals are due to the combination of β Pic b and of the additional companion. The pulsation-subtracted RV data series, as well as the offset series, were fitted with

two Keplerian functions, accounting respectively for β Pic b and for the additional companion, plus an offset (hereafter referred to as the global offset) to account for the fact that the RVs are computed relative to a reference spectrum (Methods). An evolutionary algorithm similar to the one used to fit the pulsations was used. A mass of 1.75 M☉ was assumed for the star²¹.

β Pic b's known orbital properties and associated uncertainties have to be taken into account. The astrometric orbit of β Pic b was fitted independently using our MCMC software applied to high-contrast images²². This resulted in a collection of 500,000 orbits that represents the posterior distribution of β Pic b's orbital elements. We selected a subset of 219 orbits from this sample with the following criteria: semi-major axis *a* = 9.0 ± 0.4 AU (1 σ in the main peak), *e* < 0.055 and χ² < 19.2. These criteria ensure that the best possible solutions have been selected. For each of these solutions, we fitted the RV data. The free parameters for each fit were the β Pic b mass, the additional companion mass and orbital parameters, and a global offset. The mass of β Pic b was allowed to vary between 5 and 15 M☉, a range that largely encompasses all published estimated masses, derived primarily from model-dependent brightness-mass relations. For the additional companion, an inclination of 90° was assumed, which is reasonable as the system is seen nearly edge on and β Pic b is only very slightly inclined (<1°) with respect to the star equatorial plane. We then obtained 219 orbital solutions for the companion, as well as 219 masses for both β Pic b and the companion. An example of the two-planet fit of the pulsation-corrected RV data is shown in Fig. 3. This example corresponds to one of the very best solutions found among the 219 nominal ones; this solution is referred to as no. 108.

The distributions of the orbital parameters of the additional companion, of its mass and of the mass of β Pic b are presented in Fig. 4. The additional companion's semi-major axis is about 2.7 AU, corresponding to a ~1,220-d period. Its eccentricity is slightly greater than 0.2, and its mass, 8.9 M☉, falls well in the planet range. It will be referred to as β Pic c in the following. The mass of β Pic b is found to be about 9.9 M☉. Given that our total baseline is still shorter than the β Pic b period (see Fig. 3), this mass is only indicative. We note that it is compatible with recent estimates based on measured Hipparcos and GAIA astrometric measurements²³,²⁴. We also note that the presence of β Pic c will itself require a re-examination of these astrometric constraints on the mass of β Pic b. Finally, we checked that a slightly different mass of β Pic b would not substantially change the mass of β Pic c. A similar conclusion is reached in the case of a slightly different (a few %) mass for the star.
<figure><img src="image_4.png" /></figure>

**Fig. 3** | Keplerian fits of the pulsation-corrected RV data. **a**, Two-Keplerian fit of the pulsation-corrected RV using as an input the best orbital solution of the fit of β Pic b astrometry (no. 108). **b**, Generalized Lomb-Scargle periodogram of the residuals after removing both planets' signals. **c**, Pulsation-corrected and β Pic b-corrected RV, phased to the β Pic c period; the β Pic c orbital solution is shown as a solid blue curve. **d**, Pulsation-corrected and β Pic c-corrected RV, phased to the β Pic b period; the β Pic b orbital solution is shown as a solid blue curve. Note that, because of the limited time baseline, β Pic b phases are not fully covered. **e**, β Pic b-corrected offsets, phased at the β Pic c period. The orbital solution found for β Pic c is shown as the blue curve. Note that the mass of β Pic b was fixed to that found in the Keplerian fit shown in **a,c,d**, because there are not enough data points to allow this mass to be a free parameter.

We conclude that the current data and knowledge on the star's properties provide very strong evidence of a giant planet in the inner system, in addition to β Pic b. Additional HARPS data will be valuable to strengthen this result and improve the planet's characterization. In the following, we discuss the impact of this result.

To our knowledge, β Pictoris is today the only system around which a planet has been detected by direct imaging and a second planet found by indirect techniques. β Pictoris is the second known multiple planetary system with at least one imaged planet. The first

one is HR 8799, with four imaged planets of roughly equal masses and distances from 14 to more than 60 AU<sup>25,26</sup>. The planetary formation mechanism in HR 8799 is still unknown. Given their large orbits, HR 8799 b and c are good candidates for formation by gravitational instability in a disk. Given their small distances from their parent star, β Pic b and c are good candidates for in situ formation by core accretion. Furthermore, we note that β Pic c orbits near the lower limit of the snowline position for a $1.75 M_\odot$ star<sup>27</sup>.

The presence of a planet such as β Pic c raises several dynamical issues for the system as a whole.
<figure><img src="image_5.png" /><figcaption>**Fig. 4 | Properties of the additional planet.** Histograms of the semi-major axis, eccentricity, periastron argument and periastron time of β Pic c and the masses of the two planets.</figcaption></figure>

* The first is whether the two planets are dynamically stable. A quick analysis shows that the system is indeed stable over a time of the order of its age (Methods).

* A second issue concerns the relationship between β Pic c and the FEBs claimed to be responsible for strong, transient absorptions in the lines of ionized elements, monitored in the β Pic spectrum since the 1980s (ref. ² and references therein). Even though no direct correlation is found between the exocomet absorptions (see Supplementary Fig. 3 for examples of such absorptions) and the RV variations, a possible link cannot be excluded (Methods).

* Finally, it was proposed that β Pic b could be responsible for the debris disk inner warp if it is in a slightly inclined orbit¹. Later, it was shown that the planet orbit was indeed inclined²⁸. The presence of β Pic c does not substantially change the conclusions, owing to its smaller distance from the star than β Pic b.

## Future observations

The expected astrometric wobble produced by β Pic c is about 700 µas, substantially smaller than the astrometric excess noise, 2.12 mas, reported in the current GAIA data release, but it will be easily detected when all data become available (~2022). An aggressive monitoring of the β Pictoris RV variations will allow further constraint of the planet orbital properties within the next two or three years.

A direct detection in imaging of β Pic c would allow a unique calibration of the brightness-mass relations for young stars. Such a calibration is required to accurately estimate the masses of directly

imaged planets. It would also allow for a test of the planet formation mechanisms and, in particular, the way the material is accreted onto the planet embryo. Indeed, hot-start models of 9–10 $M_J$ planets predict an H magnitude (at 21 Myr) more than 100 times higher than cold-start models²⁹. The expected contrast between star and planet in H is about 10 mag, assuming a hot-start model. The maximum projected separation depends on the detailed planet orbital properties but is in the range 0.1"–0.15". Achieving such a contrast may then be possible with SPHERE when the planet is at its maximum elongation. Note that, given the present uncertainties on the orbital parameters, it is not possible to provide a precise date for this event. A contrast 100 times larger, on the other hand, would not be achievable with current instruments. Hence, a positive detection would favour a hot-start model and would provide a confirmation of the brightness prediction. A negative detection would be less conclusive, as it would mean that either the hot-start model predictions are overestimated or that the planet formed under a cold-start model.

Finally, because its semi-major axis is substantially smaller than that of β Pic b, β Pic c may transit its star. A transit of β Pic c would allow its atmosphere to be explored, while a transit of its Hill sphere would allow for the search for moons or circumsplanetary rings. Co-orbital bodies trapped in the Lagrangian regions of β Pic c might also be searched for, as they are common outcomes of planet formation hydrodynamical simulations.

## Methods

**Data and RV measurements.** A total of 6,035 high-resolution spectra of β Pictoris were obtained between October 2003 and April 2018. We removed the data taken with air mass larger than 2.0 to avoid any potential shift due to a poor centring
of the star on the fibre. We also removed the few data with the lowest SN and highest SN (to avoid possible saturation effects), to keep the SN at 550 nm in the range 120–500. As explained in the main text, we started in March 2008 to record continuous, typically 1.5–2 h, sequences of spectra; only these spectra were selected for the present analysis.

To compute the RV and associated uncertainties, we used the SAFIR code³⁰, which was developed to measure accurate RVs of rapidly rotating stars such as β Pictoris, and associated uncertainties, which include the photon noise and the instrumental uncertainties. Basically, the code consists in correlating in the Fourier space, each spectrum of the star and a reference spectrum that is built by averaging all the available spectra of the star. The RVs are computed by SAFIR relative to this reference spectrum.

To correct for possible global chromatic effects (for example due to atmospheric absorption, seeing, atmospheric refraction residuals), we implemented in SAFIR the possibility to normalize each spectrum (before computing the RV), to obtain the same evolution of the average flux at the centre of order as a function of the order, based on the reference spectrum. This normalization does not change the measured RV within the uncertainties. Second, as variations of the local continuum within the orders directly impact the RV, particularly in the case of highly rotating stars such as β Pictoris (although moderated by filtering low spatial frequencies), and as there may be variations to the level of a few per cent (for example due to blazes plus calibration lamp residuals³¹), we also normalized the continuum of each spectrum locally (within each order) to that of the reference spectrum before computing the RV. This local normalization allows correction of the continuum from potential variations occurring within the orders. We thus obtained two RV time series (with and without this local normalization). The two corresponding time series appear to be comparable within uncertainties except on JD = 2,456,773. As this indicates a substantial and uncontrolled change in the blaze, we did not take these specific data into account in the rest of the study. The RVs obtained after averaging the two time series are shown in Fig. 1 and listed in the Supplementary Table 1 dataset.

Note that a major instrumental change was made in HARPS in 2015, which included the installation of a new set of fibres and scrambling capabilities. To our knowledge, the impact of this change was estimated only in the case of low-rotation stars³². To check the impact of this change, we first computed the RVs of the data obtained before and after this change separately (using thus two references), and then the RV using all the data (using here one reference only). No change was detected within the uncertainties. This means that the impact is smaller than 10–20 m s⁻¹.

**Correction from stellar pulsations.** The measured RVs were used to compute the generalized Lomb–Scargle periodogram and identify a peak at 1180 ± 130 d, as described in the main text.

As a second step, we corrected the RV from the high-frequency pulsation pattern. To do so, we used long enough (1.5–2 h or more) continuous sequences, presenting at least three or four well sampled extrema, to estimate as accurately as possible the pulsation pattern. Sequences not complying with this criterion were then removed. In rare cases, the sequences were (slightly) shorter but repeated during the same night or over consecutive nights or nights separated by less than eight nights; these sequences were retained. Finally, data separated by less than 8 d were gathered in groups. A total of 28 groups (with 5,108 data points) were selected according to the aforementioned criteria. The corresponding RVs are shown in Fig. 2 in a single figure and detailed in Supplementary Fig. 1.

To correct for the stellar pulsations, we used three different approaches and several methods within these. The first approach (Method 1) consists in averaging the data considered for each group. A second approach (Methods 2 to 4) consists in fitting the pulsations by sine functions and an offset in an iterative way. Finally, our third approach (Method 5) consists in fitting the data simultaneously by sine functions plus an offset. The various methods are described below.

* Method 1. An average of the RVs was computed for each group of data. This method is accurate only when several pulsations are monitored and are properly (regularly enough) sampled. This is not the case at JD = 2,454,780, 2,457,344 and 2,457,849, and to a lesser extent 2,455,170, 2,456,579 and 2,456,583 (Supplementary Fig. 1).

* Method 2. Each data group was independently fitted with a sine function (three free parameters), plus a zero-point offset. The zero-point offset is justified by the fact that long-term (>1,000 d) variations are present, and that the RVs are relative to the RV of a reference spectrum. 1,000 fits were made with frequency estimates between 0.014 and 0.034 d⁻¹ and the best fit was retained. The corresponding sine function was then subtracted. We repeated this fitting process nine times. Hence, each group was fitted by 10 sine functions plus an offset, which corresponds to a total of 31 free parameters. In most cases, the offset and the r.m.s. of the residuals converged rapidly, after few (typically three) iterations. In a few cases, the offset had not fully stabilized at the 10th iteration, but the residuals were much lower than the uncertainties (computed below). The five data groups corresponding to the longest sequences and a large number of data points were fitted with 30 sine functions instead of 10, assuming that the length of the sequence allows us to characterize more pulsations, to be possibly compared with the 31 pulsations derived from the

photometric analysis of Mekarnia and colleagues²⁰ and Dupuy and colleagues²⁴. The convergence was reached within typically 15 iterations. The fits of the pulsations are shown in Supplementary Fig. 1 for each group, and the offset values and associated uncertainties (see below) are listed in Supplementary Table 2. The number of nights in each data group, the number of points and the r.m.s. of the RV residuals after subtraction of the fitted pulsations and offset are also given in this table. To compute the uncertainties associated with the offsets, we proceeded as follows: for each group, we generated a synthetic RV time series using the 30 pulsation frequencies found in the fit of the January 2017 data (755 points) with the time sampling associated with this group. We assumed that this provides a realistic representation of the pulsations observed at both high frequency (in the δ Scuti domain) and low frequency (1–10 d⁻¹). For each point of the data group, we added Gaussian random noise to the synthetic RV time series, with a r.m.s. equal to the RV uncertainties for this point. We then fitted each synthetic RV time series using the same step-by-step procedure as described above. The expected offsets for each of these fits should be zero. This process was repeated 1,000 times, after shifting the temporal sequence to probe different positions in the time series (nodes, noise realizations and so on). The r.m.s. of the 1,000 values of the offsets provides an uncertainty associated with the time sampling and therefore with the RV of the corresponding data group. The median (respectively r.m.s.) of the resulting uncertainties is 31 (respectively 11) m s⁻¹. We then added, for each group, an uncertainty equal to half the difference between the offsets measured with the SAFIR usual method and the additional local normalization described above.

* Method 3a. We fitted a long sequence obtained in January 2017 (seven nights, 755 points, including a long continuous 6.5 h sequence) as in Method 2. We then applied the same iterative method as in Method 2 to each of the 27 other data groups, while keeping the periods equal to those found in the fitting of the January 2017 data. We assumed then here that the frequencies of the pulsations do not vary over our time basis (between 2008 and 2018). The phases were free because our precision on the periods is not enough to fix them, and the amplitudes were also free, to account for possible temporal variations.

* Method 3b. The same as 3a but using the periods obtained in the fit of the October 2017 data (508 points).

* Method 3c. The same as 3a but using the periods obtained in the fit of the Spring 2017 data (782 points).

* Method 4. The same as 3a but using the first 13 periods given by the photometric analysis of Mekarnia and colleagues²⁴.

* Method 5. Each data group was fitted simultaneously with 10 sine functions plus an offset. To do so, an evolutionary algorithm was used. To estimate the uncertainties associated with the offsets, 100 time series were generated where each RV data point was replaced by a random value drawn according to a Gaussian law centred on the original RV with a σ equal to the uncertainty of the measurement, N(RV, uncertainty). These time series were fitted, and a new series was created with, for each date, the average of the resulting offsets, and with the r.m.s. of the offsets as uncertainties.

Supplementary Fig. 2 shows the offsets obtained with these different methods, and Supplementary Table 2 lists the offsets obtained with Method 2 and their associated uncertainties. As expected, the differences are larger when comparing Method 1 and Methods 2–5 when the sequences are short and/or not well sampled (see above). The r.m.s. of the differences between the offsets obtained with Methods 1 and 2 is 20 m s⁻¹ when considering all epochs and is reduced to 10.6 m s⁻¹ if we do not consider the poorly sampled datasets obtained at JD = 2,450,000 = 4,780, 7,344 and 7,849. It is further reduced to 8.2 m s⁻¹ when we do not consider in addition JD = 2,450,000 = 5,170, 6,569 and 6,583. The offsets obtained with Methods 2–5 are in excellent agreement, with averaged differences between Methods 2 and 3a (respectively 2 and 3b, 2 and 3c, 2 and 4, 2 and 5) of less than 2 m s⁻¹ and r.m.s. of 6.5 (respectively 8.1, 8.5, 8.8 and 8.2) m s⁻¹.

Noticeably, most of the high-amplitude pulsations found in Method 2 are in agreement with the analysis of the photometric data²⁴. In the long sequences, at least seven or eight frequencies found in the RV fit are also found from analysis of photometric data by Mekarnia and colleagues²⁴. The fact that Methods 2 and 4 give very similar results indicates that the fit of the lowest-amplitude pulsations does not impact the offsets or the long-term variations of the residuals.

The fits of the longest groups moreover reveal additional, low-frequency (2–6 cycles d⁻¹) pulsation(s) with an amplitude of about 40 m s⁻¹, that is 6 to 11 times lower than the main pulsation peak. Main sequence A–F stars showing pulsations with periods in the range 0.25–3 d are classified as γ Doradus pulsators and the pulsations are attributed to high-order g modes. We conclude that β Pictoris might belong to the small group of ‘hybrid’ δ Scuti/γ Doradus stars, mostly detected using long and continuous monitoring by Corot and Kepler³³. Finally, some low-amplitude variations on timescales of 10–30 d could be present as well, as seen in the longest sequences, but the available data do not allow us to characterize them. However, we checked that such pulsations do not create high peaks at longer periods, even considering beats between several pulsations.

**Could the RV residuals be spurious?** The 1,000 simulated time series described in Correction from stellar pulsations to estimate the uncertainties were used to check whether the pulsations can produce long-term RV variations similar to those
observed. We computed the periodogram for each time series, and the maximum power in the 1,100–1,300-d range. The maximum power (at the 1% level) is slightly below the FAP and far below the peak amplitude of the 1,200-d peak observed in the periodogram of the raw RV time series or that of the pulsation-corrected RV ones. Also, the periodogram of the RV time series with offsets subtracted does not show a peak above the FAP. Hence, the pulsations do not lead to a substantial peak at 1,200 d.

**Dynamical impact of an additional planet.** *Stability of the two-planet system.* The ratio of the β Pic b and β Pic c semi-major axes is about 3.5, so the system is presumably stable. Supplementary Fig. 4 shows an example of secular evolution of the semi-major axis and eccentricity of the two planets over 1 Myr, as computed with the symplectic code HJS³⁴. The real computation was extended to 50 Myr and does not reveal any change with respect to the 1 Myr case. We note that the semi-major axes are extremely stable, while the eccentricities exhibit moderate regular fluctuations. This is characteristic for regular dynamics and allows us to conclude that the system is stable. Interestingly, we note that β Pic b's eccentricity oscillates secularly between ~0 and 0.15. This shows that β Pic c is able to excite β Pic b's eccentricity. Conversely, the eccentricity of β Pic c never reaches zero. This means that at least part of the planet's eccentricity cannot be explained by the sole perturbing action of β Pic b and must be primordial.

*Impact on the FEB scenario.* Dynamical studies, based on the large number of observed FEB events, have shown that the presumed star-grazing asteroids or comets probably originate from internal mean-motion resonances (4:1 or 3:1) with a Jovian planet orbiting β Pictoris at ~10 AU with moderate eccentric orbit ($e \sim 0.1$)³⁵,³⁶. β Pic b was then identified as the planet responsible for the FEB phenomenon, even though its mass ($\sim 9 M_{\oplus}$) is larger than that estimated originally ($\sim 2 M_{\oplus}$). To explain the velocity distribution of the observed blueshifted and redshifted absorptions, an additional, less massive, possibly terrestrial-mass planet orbiting more or less at the current location of β Pic c³⁵ was proposed. The secular evolution of the semi-major axis and eccentricity of the two planets over a period of 1 Myr, computed with the symplectic code HJS³⁴, yields variations of a and e of β Pic b of about 0.015 AU and 0.1 respectively, and less than $5 \times 10^{-3}$ AU and less than 0.15 for β Pic c. Note that the two separate families of FEBs were confirmed recently³⁷ on the basis of a larger set of data. In this scenario, the masses of the planets are much smaller than the masses of the observed planets. Thus, new simulations are needed to see whether the FEBs can be generated under these new conditions. In fact, the presence of β Pic c might actually enhance the FEB generation process. As seen above, β Pic c is able to secularly excite β Pic b's eccentricity to larger values than previously considered, leading to more efficient resonances and enhancement of the FEB generation process. Indeed, with a larger eccentricity, other mean-motion resonances such as 5:2 may become additional sources of FEBs that were not considered before. β Pic c itself is able, to a moderate extent, to excite the eccentricities of FEB progenitors locked in mean-motion resonances with β Pic b and contributing to the FEBs. Of course, the number of progenitors that will actually become FEBs despite encounters with β Pic c is a key issue. The high mass of β Pic b might create wider resonance ranges than previously. Also, with a higher mass, the resonant eccentricity increase that drives the progenitor to the FEB state is faster. Various aspects of this scenario will be the purpose of future work.

**Data availability**

The HARPS spectra are available in the ESO archive, and the measured RV data are given in Supplementary Table 1.

**Code availability**

The codes used for this paper are available from the corresponding author upon reasonable request.

Received: 22 October 2018; Accepted: 4 July 2019;
Published online: 19 August 2019

**References**

1. Baruteau, C., Bai, X., Mordasini, C. & Mollière, P. Formation, orbital and internal evolutions of young planetary systems. *Space Sci. Rev.* **205**, 77–124 (2016).

2. Smith, B. A. & Terrile, R. J. A circumstellar disk around β Pictoris. *Science* **226**, 1421–1424 (1984).

3. Augereau, J. C., Nelson, R. P., Lagrange, A.-M., Papaloizou, J. C. B. & Mouillet, D. Dynamical modeling of large-scale asymmetries in the β Pictoris dust disk. *Astron. Astrophys.* **370**, 447–455 (2001).

4. Dent, W. et al. Molecular gas clumps from the destruction of icy bodies in the β Pictoris debris disk. *Science* **383**, 1490–1492 (2014).

5. Kiefer, F. et al. Two families of exocomets in the β Pictoris system. *Nature* **514**, 462–464 (2014).

6. Lagrange, A.-M. A probable giant planet imaged in the β Pictoris disk. VLT/NaCo deep L'-band imaging. *Astron. Astrophys.* **493**, L21–L25 (2009).

7. Lagrange, A.-M. et al. A giant planet imaged in the disk of the young star β Pictoris. *Science* **328**, 57–59 (2010).

8. Mouillet, D., Larwood, J. D., Papaloizou, J. C. B. & Lagrange, A.-M. A planet on an inclined orbit as an explanation of the warp in the β Pictoris disc. *Mon. Not. R. Astron. Soc.* **292**, 896–904 (1997).

9. Nesvold, E., Kuchner, M. J. & SMACK, A. Model of colliding planetesimals in the β Pictoris debris disk. *Astrophys. J.* **798**, 83–100 (2015).

10. Beust, H. & Morbidelli, A. Falling evaporating bodies as a clue to outline the structure of the β Pictoris young planetary system. *Icarus* **143**, 170–188 (2000).

11. Telesco, C. M. et al. Mid-infrared images of β Pictoris and the possible role of planetesimal collisions in the central disk. *Nature* **433**, 133–136 (2005).

12. Okamoto, Y. K. et al. An early extrasolar planetary system revealed by planetesimal belts in β Pictoris. *Nature* **431**, 660–662 (2004).

13. Wahhaj, Z. et al. The inner rings of β Pictoris. *Astrophys. J.* **584**, L27–L32 (2003).

14. Lagage, P. O. & Pantin, E. Dust depletion in the inner disk of β Pictoris as a possible indicator of planets. *Nature* **369**, 628–630 (1994).

15. Lagrange, A.-M. et al. Full exploration of the giant planet population around β Pictoris. *Astron. Astrophys.* **612**, 108–112 (2018).

16. Lagrange, A.-M. et al. Constraints on planets around β Pic with Harps radial velocity data. *Astron. Astrophys.* **542**, A18–A23 (2012).

17. Bonnefoy, M. et al. Physical and orbital properties of β Pictoris b. *Astron. Astrophys.* **567**, L9–L14 (2014).

18. Koen, C. δ Scuti pulsations in β Pictoris. *Mon. Not. R. Astron. Soc.* **341**, 1385–1387 (2003).

19. Koen, C. et al. Pulsations in β Pictoris. *Mon. Not. R. Astron. Soc.* **344**, 1250–1256 (2003).

20. Mekarnia, D. et al. The δ Scuti pulsations of β Pictoris as observed by ASTEP from Antarctica. *Astron. Astrophys.* **608**, L6–L10 (2017).

21. Crifo, F., Vidal-Madjar, A., Lallemand, R., Ferlet, R. & Gerbaldi, M. β Pictoris revisited by Hipparcos. Star properties. *Astron. Astrophys.* **320**, L29–L32 (1997).

22. Lagrange, A.-M. et al. β Pictoris b post conjunction detection with VLT/SPHERE. *Astron. Astrophys.* **621**, L8–L14 (2019).

23. Snellen, I. A. G. & Brown, A. G. A. The mass of the young planet β Pictoris b through the astrometric motion of its host star. *Nat. Astron.* **2**, 883–886 (2018).

24. Dupuy, T., Brandt, T. D., Kratter, K. M. & Bowler, B. P. A model-independent mass and moderate eccentricity for β Pic b. *Astrophys. J.* **871**, L4–L9 (2019).

25. Marois, C. et al. Direct imaging of multiple planets orbiting the star HR 8799. *Science* **322**, 1348–1352 (2008).

26. Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B. & Barman, T. Images of a fourth planet orbiting HR 8799. *Nature* **468**, 1080–1083 (2010).

27. Kennedy, G. M. & Kenyon, S. J. Planet formation around stars of various masses: the snow line and the frequency of giant planets. *Astrophys. J.* **673**, 502–512 (2008).

28. Lagrange, A.-M. et al. The position of β Pictoris b position relative to the debris disk. *Astron. Astrophys.* **546**, 38–51 (2012).

29. Fortney, J. J., Marley, M. S., Saumon, D. & Lodders, K. Synthetic spectra and colors of young giant planet atmospheres: effects of initial conditions and atmospheric metallicity. *Astrophys. J.* **683**, 1104–1116 (2008).

30. Galland, F. et al. Extrasolar planets and brown dwarfs around A–F type stars. I. Performances of radial velocity measurements, first analyses of variations. *Astron. Astrophys.* **443**, 337–345 (2005).

31. Anglada-Escudé, G. & Butler, P. R. The HARPS-TERRA Project. I. Description of the algorithms, performance, and new measurements on a few remarkable stars observed by HARPS. *Astrophys. J. Suppl.* **200**, 15–34 (2012).

32. Locurto, G. et al. HARPS gets new fibres after 12 years of operation. *Messenger* **162**, 9–15 (2015).

33. Bradley, P. A. et al. Analysis of γ Doradus and δ Scuti stars observed by Kepler. *Astron. J.* **149**, 68–81 (2015).

34. Beust, H. Symplectic integration of hierarchical stellar systems. *Astron. Astrophys.* **400**, 1129–1144 (2003).

35. Beust, H. & Morbidelli, A. Mean-motion resonances as a source for infalling comets toward β Pictoris. *Icarus* **120**, 358–370 (1996).

36. Thebault, P. & Beust, H. Falling evaporating bodies in the β Pictoris system. Resonance refilling and long term duration of the phenomenon. *Astron. Astrophys.* **376**, 621–640 (2001).

**Acknowledgements**

This work has been supported by grants from the Agence Nationale de la Recherche (ANR-14-CE33-0018) and the French Labex OSUG@2020 (Investissements d'avenir—ANR10 LABX56). A.Z. was supported by CONICYT grant no. 2117053. A.-M.L. thanks F. Forbes, K. Zwincks, A. Lecavelier, J. Pepper, P. Kervella and J. C. B. Papaloizou for discussions. T.G., D.M., I.A. and F.-X.S. acknowledge support from Idex UCAJEDI (ANR-15-IDEX-01) and IPEV.

**Author contributions**

A.-M.L. led the monitoring of the variations, the data reduction, the analysis and interpretation of the data, and the paper writing. N.M., P.R., M.K. and F.G.
participated to the data fitting and analysis. E.C., E.M., L.B. and F.-X.S. brought their expertise in stellar variability. H.B. provided analysis of the dynamical stability of the system. T.G., D.M. and L.A. brought expertise on β Pictoris photometric variability. P.A.W. and F.K. brought expertise on β Pictoris spectroscopic variability. M.B., S.B., A.G., J.L.-B., B.P., D.P.I., L.R. and A.Z.S. participated in the observations.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41550-019-0857-1.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Correspondence and requests for materials** should be addressed to A.-M.L.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2019