© The Authors 2025

1 Introduction

The star β Pictoris (hereafter β Pic) is a young (${18.5}_{-2.4}^{+2.0}$ Myr, Miret-Roig et al. 2020) and nearby (19.63 ± 0.06 pc, Gaia Collaboration 2020) southern star that has consistently captured the attention of astronomers since the discovery of its circumstellar dust disk by Smith & Terrile (1984), the first of its kind ever imaged. Earlier that same year, the presence of circumstellar dust was also inferred for the first time through infrared observations of Vega (Aumann et al. 1984). Since these groundbreaking discoveries, it has been established that at least ∼20% of FGK stars host such a circumstellar dust disk (Eiroa et al. 2013; Sibthorpe et al. 2018).

Initially, the observed dust was believed to be of primitive origin. However, Backman & Paresce (1993) demonstrate that this dust is too short-lived for that paradigm to be valid. Indeed, the lifetime of dust grains, when considering removal mechanisms such as radiation pressure, collisions, and Poynting- Robertson drag (see, e.g., the review by Krivov 2010), is shorter by orders of magnitude than the ages of the stars around which they are detected. Consequently, this dust is understood to be second-generation material, widely believed to be continuously replenished by an underlying population of kilometer-sized bodies. These planetesimals serve as a reservoir, capable of sustaining dust production throughout the lifetime of stars, either through slow evaporation (Lecavelier Des Etangs et al. 1996) or collisions (Backman & Paresce 1993; Artymowicz 1997).

By observing the dust and understanding the mechanisms behind particle production and movement, it is possible to deduce the structure of the planetesimal reservoir from which the dust originates. However, this connection is rarely straightforward. Radiation pressure causes small dust particles to deviate from the orbits of their parent bodies, often reaching apoastrons far from their initial production sites (Lecavelier Des Etangs et al. 1996). Some of these particles are even observed in scattered light at significant distances from β Pic (see, e.g., Janson et al. 2021). To address this complexity in the β Pic’s system, models have been developed to constrain the distribution of parent bodies based on the observed dust. Several authors have made observations and/or proposed models suggesting that the inner edge of the disk is at ∼50 au (Augereau et al. 2001; Dent et al. 2014; Apai et al. 2015; Ballering et al. 2016). However, the edge-on orientation of β Pic’s disk causes irregularities at its inner edge in these models due to projection effects, introducing significant uncertainty into the estimates. This problem also occurs at millimetre wavelengths. Millimeter emission is a valuable tool for directly revealing the structure of planetesimal disks. It is dominated by large dust particles, which are only minimally affected by radiation pressure, follow the same orbits as their parent bodies, and remain close to them. However, excellent resolution at these wavelengths is crucial for obtaining accurate information about the inner edge of the disk in this edge-on configuration. For example, Matrà et al. (2019) reports observations with at least twice the resolution of Dent et al. (2014), confirming an inner edge of the disk at ∼50 au.

The debris disk around β Pic is one of the most extensively studied, offering valuable insights into the dynamical and evolutionary processes of disks interacting with planets. To date, at least one giant planet has been hypothesized to explain various phenomena observed in the disk. For example, the misalignment between its inner and outer regions (Mouillet et al. 1997), the asymmetries observed between the two arms of the disk (Kalas & Jewitt 1995; Heap et al. 2000), and the recurring detection of exocomets in the star’s spectrum (Beust & Morbidelli 2000), all suggest the gravitational influence of at least one planet (Augereau et al. 2001). Indeed, a planet was observed. β Pic b, a gas giant with a mass of ∼12 M_Jup, was detected using high-contrast imaging techniques (Lagrange et al. 2009). Its orbit was since refined through regular monitoring, revealing a moderately eccentric orbit with an eccentricity of ∼0.1 at a distance of ∼10 au from the star (Lacour et al. 2021).

Since then, a number of intriguing features have been identified in the disk, potentially linked to dynamical gravitational interactions with planets. Mid-infrared observations reveal a large dust clump on the western side of the disk (Telesco et al. 2005). Submillimeter observations by Atacama Large Millime- ter/submillimeter Array (ALMA) also detect a concentration of CO coinciding spatially with this dust clump (Dent et al. 2014; Matrà et al. 2017). However, its nature remains uncertain. If the clump is in motion, it could indicate the presence of a gas vortex, potentially caused by an unseen planet trapping the dust (Skaf et al. 2023). Conversely, a stationary clump might result from a past massive collision (Han et al. 2023), possibly consistent with the “cat’s tail” recently observed by James Webb Space Telescope (JWST), which may represent debris from a collision between two large planetesimals (Rebollido et al. 2024). It should be noted that these features are located much further out in the disk than the region studied in this work, although there may be a link with the dynamical gravitational interactions of planets closer to the star.

Furthermore, a second planet, another gas giant with a mass of ∼9 M_Jup , β Pic c, was detected. It was discovered through radial velocity monitoring (Lagrange et al. 2019) and subsequently confirmed by high-contrast imaging (Nowak et al. 2020). This planet, while slightly less massive than β Pic b, orbits closer to the star at ∼3 au, following a more eccentric orbit with an eccentricity of ∼0.3 (Lacour et al. 2021).

Our goal in this paper is to establish a connection between the β Pic’s planetary system and the structure of the planetesimal disk, particularly the inner cavity at ∼50 au. We assumed this cavity is shaped by the gravitational influence of the planets, neglecting other dynamical effects that remain poorly constrained today. We adopted an age of 20 Myr for the β Pic system and assumed that the inner edge of the disk lies precisely at 50 au. However, it is important to recognise the uncertainties associated with this model, as the results can be refined with the arrival of new, more precise observations. In Sect. 2, we first conduct a numerical exploration of the gravitational effects of the two known planets, concluding that they cannot sculpt the disk up to 50 au. In Sect. 3, we explore the possibility of an additional, yet-to-be-discovered planet orbiting outside β Pic b that could fulfill this role, and we constrain the characteristics of this hypothetical planet. In Sect. 4, we consider alternative scenarios, such as a more eccentric planet or the presence of two additional planets instead of one. We present our conclusions in Sect. 5.

2 Exploration with the two known planets

Here, we focus on the dynamical effects of the two known planets in the β Pic system on the debris disk’s inner edge and demonstrate that they are insufficient to create this edge at 50 au.

2.1 Semi-analytical theory

According to the theory proposed by Wisdom (1980) and Mustill & Wyatt (2012), a planet can create a gap in a disk when its first- order mean-motion resonances j : j ± 1 overlap. This overlap is inevitable for a certain value of j, because as j increases, the resonances are progressively located closer to the planet and closer to each other. Wisdom (1980) provides an estimate for the critical value of j at which this overlap occurs as $$j_{\text{overlap }}=0.51{\mu }^{-2/7},$$(1)

where µ is the ratio of planetary to stellar mass. The location of this critical resonance provides an estimate of the chaotic zone, and consequently, the size of the gap in terms of semi-major axis, given by $${\left(\frac{\delta a}{a}\right)}_{\text{chaos }}=1.3{\mu }^{2/7},$$(2)

where the full gap size, including both the inner and outer regions relative to the planet’s orbit, should be about twice this value. Mustill & Wyatt (2012) highlight that this result depends on the planet’s eccentricity. Equation (2) remains valid for small eccentricities, up to a critical eccentricity estimated to $$e_{\text{crit }}=0.21{\mu }^{3/7}.$$(3)

For larger planet eccentricities, a new regime applies where $${\left(\frac{\delta a}{a}\right)}_{\text{chaos }}=1.8e^{1/5}{\mu }^{1/5}.$$(4)

Let us assume that the inner edge of the planetesimal disk is shaped by β Pic b, as β Pic c orbits much closer to the star. A numerical application using β Pic b’s parameters gives a critical eccentricity of e_crit ≃ 0.02. In this case, the large eccentricity regime described by Eq. (4) applies, leading to a chaotic zone width of (δa/a)_chaos ≃ 0.4. This places the outer edge of the gap at ∼14 au. This implies that the two known planets β Pic b and β Pic c alone lack the capacity to carve the disk up to 50 au, suggesting the presence of an additional planet.

However, this result should be interpreted with caution, as the formulae defined by Wisdom (1980) and Mustill & Wyatt (2012) are asymptotically valid for large j. Here, due to the high mass of β Pic b, the overlap of the mean-motion resonances near β Pic b occurs at j_overlap ≃ 2, which is not particularly large.

Another possible semi-analytical prediction that accounts for the high mass of β Pic b is based on the work of Morrison & Malhotra (2015). These authors showed that the size of the outer chaotic zone, when the planet-to-star mass ratio µ exceeds 10⁻⁴, is described by $${\left(\frac{\delta a}{a}\right)}_{\text{chaos }}=1.7{\mu }^{0.31}.$$(5)

In the case of β Pic b, with µ = 0.0065, this formula predicts an outer gap edge at ~13 au. This value is only slightly different than the ~14 au obtained by the theory of Mustill & Wyatt (2012). However, it is important to note that the formalism of Morrison & Malhotra (2015) does not account for the planet’s eccentricity.

Thus, while existing semi-analytical approaches can accurately characterize the chaotic zone for planets that are either massive or eccentric, none of them simultaneously takes into account both regimes, nor consider small j_overlap values. Given these limitations, numerical simulations become indispensable for studying such regimes corresponding to the β Pic case.

Fig. 1 Preliminary simulation of the dynamics of the β Pic planetary system with only the two planets known to date, β Pic b (red) and β Pic c (blue), and the disk of planetesimals (black). The initial orbital parameters of β Pic b and c are taken from Table 1. Left: upper view of the system. The planetesimals are depicted by small black dots, and the orbits of the planets by colored dashed lines. Middle: view of the system in terms of semi-major axis and eccentricity. The planets are here represented by colored points. Right: radial profile of the planetesimal disk (solid black line), superimposed to the models of Augereau et al. (2001) (blue dotted line) and Dent et al. (2014) (red dotted line). Note the clear mismatch in the location of the inner edge between this simulation and observational models.

2.2 Simulation

We present a simulation that includes the two known planets and a disk of planetesimals. The physical and orbital parameters of the planets are taken from the data set presented in Lacour et al. (2021) and are summarized in Table 1. The initial disk of planetesimals consisted of 400 000 massless particles that do not interact with each other. Their initial semi-major axes were randomly selected between 20 and 80 au, with initial eccentricities ranging from 0 to 0.05 and initial inclinations between 0 and 2° relative to the invariant plane of the two-planets system. The remaining initial orbital elements of the disk particles, including the longitudes of ascending nodes, arguments of periastron, and mean longitudes, were also randomly assigned values between 0 and 360°.

Calculations were performed using the Regularized Mixed Variable Step Size Symplectic (RMVS3) integrator (Levison & Duncan 1994), a modified version of the original Mixed Variable Symplectic (MVS) scheme by Wisdom & Holman (1991), which includes a first-order but rapid treatment of close encounters. This approach is particularly relevant here, as we aim to examine the location of the disk’s inner edge as sculpted by the planets, specifically focusing on particles near the instability region. The integration was carried out over 20 Myr, corresponding to the adopted age of the system.

The results are presented in Fig. 1. As shown in the figure, the planets carve the disk out to ~28 au. This is inconsistent with observations that suggest an inner disk cavity extending out to 50 au. There are still far too many particles remaining between ~28 and 50 au. Additional tests, varying the orbits of the planets within the uncertainties provided by Lacour et al. (2021), did not alter this outcome.

Given this result, we can therefore conclude that the semi- analytical theories of Mustill & Wyatt (2012) and Morrison & Malhotra (2015) seem to underestimate the range of gravitational interactions of β Pic b, which is either too massive or too eccentric. However, the conclusion remains the same: the inner edge of the planetary disk at 50 au cannot be attributed solely to the perturbative action of β Pic b and β Pic c. In the following sections, we explore how this could be achieved with the introduction of a hypothetical additional planet orbiting outside β Pic b’s orbit.

3 An additional planet

We now hypothesize the presence of an additional planet in the β Pic system, which we shall refer to as β Pic d. We present simulations using the same initial planetesimal disk as before, but this time assuming a three-planets system. The initial orbital parameters of β Pic b and β Pic c were still set to the values listed in Table 1.

3.1 Constraints and initial orbital parameters

Since this planet is proposed to be part of the same planetary system as the two known planets, it is reasonable to assume that it is more or less coplanar with them. Therefore, we assumed that its inclination relative to the midplane of the two-planet system does not exceed 2°. Angular parameters such as the longitude of the ascending node, the longitude of periastron, and the initial mean longitude were selected randomly, as these are expected to secularly process under the gravitational influence of β Pic b and β Pic c and thus play a minor role in the long term.

The most critical parameters to determine for this hypothetical β Pic d are its mass m_d, orbital semi-major axis a_d, and eccentricity e_d . Several constraints must be satisfied for each of these parameters.

The first constraint is observational. Current observations, such as those made with Spectro Polarimetric High contrast Exoplanet REsearch (SPHERE), rule out the presence of superJupiter planets exterior to β Pic b (Lagrange et al. 2020).

Second, based on the semi-analytical theories discussed in Sect. 2.1, reproducing a disk with an inner edge at 50 au requires β Pic d to be located at a semi-major axis no greater than ~40 au. Although these semi-analytical theories are not always applicable, this constraint provides an indicative range for the planet’s semi-major axis.

Finally, the planet must remain dynamically stable against perturbations from other planets. Assuming that β Pic c is too close to the star, we focused solely on the two-planet system comprising β Pic b and β Pic d. According to Petrovich (2015), a stability criterion between these two planets, with mass ratios relative to the central star such that 10⁻² ≥ µ_b > µ_d > 10⁻⁴, can be defined by $$\frac{a_d\left(1-e_d\right)}{a_b\left(1+e_b\right)}>2.4{\mu }_b{}^{\frac{1}{3}}{\left(\frac{a_d}{a_b}\right)}^{\frac{1}{2}}+\mathrm{1.15.}$$(6)

This stability condition between β Pic b and β Pic d directly imposes a lower limit on β Pic d’s semi-major axis as a function of its eccentricity. This relationship is depicted in Fig. 2. As the eccentricity increases, the permissible mass range for β Pic d, constrained by this stability condition with β Pic b and the inner edge of the disk, narrows progressively. For eccentricities exceeding ~0.5, no viable solutions remain.

Therefore, our study focused on a planet β Pic d with a mass m_d ≤ 1 M_Jup, a semi-major axis a_d between ~30 and ~40 au, and an eccentricity e_d ≲ 0.5. For the eccentricity, we first assumed a moderate value in line with β Pic b and initially set it at e_d = 0.05.

3.2 Results

Figure 3 illustrates an example of a simulation with three planets that successfully reproduces the desired disk profile. As opposed to the previous simulation, we now find that the inner cavity aligns with the observational models.

Of course, the solution is not unique. If we assume a more massive β Pic d, it would likely create a larger gap in the disk, necessitating its placement closer to the star to reproduce a disk edge at 50 au. For each assumed value of m_d , we should be able to identify a corresponding a_d that yields the desired inner edge of the disk. To explore suitable combinations of m_d and a_d , that can reproduce the desired disk profile, we conducted several dozen simulations. The overall results of this investigation are summarized in Fig. 4, which displays the acceptable ranges of a_d for various values of m_d . As a_d depends, to within one translation, on the position of the simulated disk inner edge at 50 au, the acceptable range for ad corresponds to the measurement uncertainty of this simulated inner edge. A wide acceptable range of semi-major axes therefore indicates a very unclear inner edge.

In addition, Fig. 4 includes an attempt to adapt to a power law, similar to that of the semi-analytic theories discussed in Section 2.1, according to $$\frac{\delta a_d}{a_d}=c_1{\mu }_d{}^{c_2}\iff m_d={\left(\frac{\delta a_d}{c_1{a}_d}\right)}^{\frac{1}{c_2}}{M}_{*}.$$(7)

To achieve an inner edge of the disk at 50 au, we set $$\delta a_d=50-a_d[\text{au}],$$

where c₁ and c₂ are free parameters. For a low eccentricity of 0.05, the resulting c₁ ≃ 5.1 is significantly higher than that predicted by the theories discussed in Sect. 2.1. Here again it appears that these theories underestimate the amplitude of the chaotic zone of β Pic d because its value of j_overlap ≃ 7 is still relatively small. However, the phenomenon behavior indicated by c₂ ≃ 0.27 is reminiscent of the predictions of 0.29 and 0.31 from the theories of Mustill & Wyatt (2012) and Morrison & Malhotra (2015), respectively.

As expected, a more massive planet must be positioned closer to the star. This imposes an independent upper limit on m_d . For m_d exceeding a certain value, β Pic would need to be placed too close to the stability threshold with β Pic b, making it impossible to achieve a stable and suitable configuration. Such values of m_d are already excluded here by the observational constraints established by Lagrange et al. (2020). Nonetheless, this consideration remains crucial for the subsequent analysis.

Our study also allows us to derive a lower limit for m_d . A low-mass β Pic d generates a small gap in the disk and must therefore be positioned closer to 50 au to effectively carve the disk to this distance. However, this simultaneously creates a wider region between β Pic b and β Pic d where planetesimals could potentially orbit safely. This is illustrated in Fig. 5, which shows the output of a simulation similar to that of Fig. 3, but with md ≲ 0.15 M_Jup. In agreement with the overall results shown in Fig. 4, this combination successfully reproduces the inner edge of the disk at 50 au. However, it leaves enough space between β Pic b and β Pic d for an additional ring of stable particles to exist around ~30 au. We indeed recover the stability limit of particles at ~28 au due to β Pic b, as β Pic d is now too far away and not massive enough to destabilize the planetesimals orbiting in this region. This additional ring of particles is, of course, incompatible with the observational models of Augereau et al. (2001), Dent et al. (2014), and Ballering et al. (2016). Therefore, this configuration must be excluded.

Fig. 2 Accessible semi-major axis of the hypothetical additional planet β Pic d, prior to simulations, as a function of its eccentricity. Shaded areas are not accessible by β Pic d. Its semi-major axis must be smaller than the inner edge of the disk at 50 au (solid black line) to sculpt it up to this distance, while being greater than the stability limit with β Pic b (dashed black line) to avoid close encounters, based on the criterion of Petrovich (2015).

Fig. 3 Example of a simulation of the dynamics of the β Pic planetary system with three planets, β Pic b (red), β Pic c (blue), and an additional planet β Pic d (green), and the disk of planetesimals (black). The initial orbital parameters of β Pic b and c are taken from Table 1. In this example, the initial orbital parameters of β Pic d are: md = 0.2 MJup, ad = 35 au, and ed = 0.05. The plotting conventions are identical to Fig. 1.

Fig. 4 Combinations of β Pic d’s mass and semi-major axis within observational constraints that successfully reproduce the disk profile at 50 au. For several values of md, the corresponding acceptable ranges of ad are displayed in green. A power-law fit, according to Eq. (7), is overlaid in blue, taking into account the ranges of the semi-major axis. Shaded areas are not accessible byβ Pic d. The limits in semi-major axis are similarly presented as in Fig. 2. Additionally, above a mass of 1 MJup, corresponding to the observational limit (solid red line), β Pic d would have already been detected in previous observations (Lagrange et al. 2020). And below 0.15 MJup, corresponding to the ring limit (dashed red line), β Pic d does not completely clear the inner zone.

4 Alternative scenarios

The solutions presented in Sect. 3, which successfully reproduce the desired inner edge of the planetesimal disk, rely on specific initial parameter choices. But, other combinations of initial parameters might also yield acceptable solutions. In Fig. 5, we demonstrate that a low eccentricity and low-mass additional planet is insufficient to efficiently carve the disk up to 50 au, resulting in a stable, unperturbed ring of planetesimals around ~30 au, which is inconsistent with observational models. This observation led us to derive a minimum mass of ~0.15 M_Jup for β Pic d. This result is directly linked to our consideration of only one additional planet with low eccentricity. To address this limitation, in Section 4.1, we explore the possibility that the additional planet could have a higher eccentricity than initially assumed. Additionally, in Section 4.2, we consider the potential presence of two additional planets instead of just one.

4.1 A more eccentric additional planet

The simulations presented in Sect. 3 all assumed an initial low eccentricity of 0.05 for the additional planet. Although this eccentricity undergoes secular evolution due to gravitational perturbations from β Pic b and β Pic c, the simulations indicate that β Pic d’s eccentricity varies only slightly, never exceeding ~0.1. It is therefore of interest to study the potential for larger initial eccentricities of β Pic d.

A preliminary estimate of the outcome for such a configuration can be obtained by requiring that the apoastrons of β Pic d coincide in both scenarios. We anticipated that an eccentric planet would sculpt the disk in the same way as a low eccentricity planet with the same apoastron. Since eccentric planets dig the disk inwards and outwards more efficiently than a noneccentric planet, this hypothesis suggests that they would have to move closer to the star to maintain a similar disk structure.

We numerically explored various configurations with increasing initial values of β Pic d’s eccentricity, similarly to the approach presented in Sect. 3. The results are summarized in Fig. 6, which presents plots similar to the one in Fig. 4 for the increasing eccentricity regimes of β Pic d. For each eccentricity, we again fitted the semi-major axis by a power law according to Eq. (7).

As eccentricity increases, the phenomenon evolves, as shown by the corresponding values of c₂ in the Table 2. It remains consistent with the predictions of Mustill & Wyatt (2012) but diverges from those of Morrison & Malhotra (2015). This is expected, as the theory of Mustill & Wyatt (2012) accounts for high eccentricity of ~0.3, whereas Morrison & Malhotra (2015) does not. However, for extreme eccentricities of ~0.5, our values no longer align with either theory, as they do not apply to such eccentricities.

As expected, increasing eccentricity effectively clears the remaining dust ring when present, rapidly lowering the lower mass limit and extending the range of acceptable parameters to smaller masses. However, with higher eccentricity, β Pic d becomes more likely to enter the region of close encounters with β Pic b, falling below the stability limit determined in Sect. 3.1. The higher the eccentricity, the smaller the range of acceptable masses for β Pic d. Thus, for each eccentricity regime, a maximum mass is defined, as summarized in Fig. 7. As eccentricity increases, this maximum mass limit decreases, becoming more restrictive than the observed mass limit for e_d ≳ 0.25, and tends towards zero as eccentricity continues to increase. Beyond an eccentricity threshold of ~0.4, it becomes impossible for a planet, regardless of its mass and semi-major axis, to match the observational constraints while remaining stable within the system.

Fig. 5 Example of a simulation of the dynamics of the β Pic planetary system with three planets, β Pic b (red), β Pic c (blue) and an additional planet β Pic d (green), and the disk of planetesimals (black). The initial orbital parameters of β Pic b and c are taken from Table 1. In this example, the initial orbital parameters of β Pic d are: md = 0.05 MJup, ad = 40 au, and ed = 0.05. The plotting conventions are identical to Fig. 1. This combination still generates a disk inner edge at 50 au, but there is now enough space between β Pic b and β Pic d to allow an additional ring of particles to exist at ~30 au.

Fig. 6 Combinations of β Pic d’s mass and semi-major axis for various eccentricities, within observational constraints that successfully reproduce the disk profile at 50 au. The plotting conventions are identical to Fig. 4. The ring limit quickly tends to zero as the eccentricity increases. The minimum semi-major axis stability limit with β Pic b imposes a maximum mass limit (dot-dashed black line), determined by the intersection (black dot) with the fit’s right extremum of the acceptable range of β Pic d’s semi-major axis. This maximum mass limit decreases as eccentricity increases, becoming more restrictive than the observational mass limit for ed ≳ 0.25, and tends towards zero as eccentricity continues to increase (see Fig. 7).

Fig. 7 Accessible mass of the potential β Pic d as a function of its eccentricity ed, based on the results of Fig. 6, whose plot conventions are similar.

4.2 Several additional planets

In this section, we briefly explore the possibility of having two additional planets in the β Pic’s planetary system instead of one. The available parameter space for such a configuration is significantly limited. Clearly, these two planets must remain dynamically stable. Any configurations where a single additional planet is sufficient to efficiently clear the disk up to 50 au should not be considered here, as the insertion of an other additional planet would inevitably lead to instability. We focused on configurations like those illustrated in Fig. 5, where a single additional planet is insufficient. Furthermore, Figs. 6 and 7 show that as the eccentricity of the additional planet increases, the space available for any other additional planets is greatly reduced. Therefore, we worked on a model involving two additional planets, both with low mass (≤0.15 M_Jup) and a low eccentricity of ~0.05.

Figure 8 illustrates a simulation involving two additional planets which successfully carves the disk up to 50 au without leaving a stable ring at ~30 au. We explored various similar configurations by varying the masses and locations of the planets. Considering all requirements, the available parameter space around the configuration shown in Fig. 8 is quite limited. Notably, it appeared that among the two additional planets, the inner one must have a mass at least equal to that of the outer one to ensure its dynamical stability.

The existence of two additional low-mass, low-eccentricity planets presents therefore an interesting alternative to the model of a single, more massive planet. However, this hypothesis was not explored further, as these planets are relatively small and well below current detection limits, making it unnecessary to immediately pursue this avenue.

5 Conclusion and discussion

This paper investigates how the β Pic planetary system can dynamically clear out the disk of planetesimals up to the ~50 au threshold, as modeled by Augereau et al. (2001), Dent et al. (2014), and Ballering et al. (2016). This analysis was conducted under the assumption that the planetesimal disk initially extends from the star to well beyond 80 au, without accounting for potential planetary migration or other dynamical effects.

The first result is that the currently known planetary system is unable to clear the disk as needed. With only β Pic b and β Pic c, the disk would not be carved beyond ~28 au. Therefore, the presence of additional, yet undiscovered planets, such as β Pic d, at greater distances can be hypothesized.

Our simulations indicate that an additional planet is sufficient to achieve the desired outcome. Various combinations of mass, semi-major axis, and eccentricity are possible, but the higher the eccentricity, the more restricted the parameter space becomes. In any case, the planet’s eccentricity cannot exceed ~0.4. Alternatively, we show that a model involving two planets of low mass and low eccentricity, instead of a more massive one, could also be a viable option.

However, it is important to note that our dynamical analysis could be further refined. We did not attempt to reproduce the complete radial profile of the disk, and simplified the problem to a single radius representative of the inner edge, which is itself uncertain due to the difficulty of deprojecting the observations. As some of the simulated profiles have a complex structure, this approximation introduces uncertainties about the semi-major axis of the additional planets.

Also, we did not explore specific dynamical configurations, such as a planet trapped in an outer mean-motion resonance with β Pic b. In such a case, a phase-protection mechanism could prevent β Pic d from being ejected while still orbiting within the instability zone of β Pic b, similar to the relationship between Neptune and Pluto in the Solar System (Greenberg 1977). Future modeling efforts should focus on these specific scenarios, but we emphasize that a more accurate understanding of the exact orbital parameters of β Pic c and β Pic b is necessary first. Indeed, Beust et al. (2024) demonstrated that these two planets could be temporarily locked in a 7:1 mean-motion resonance. If this were the case, any additional planets in resonance with β Pic b would form a chain of resonances involving three planets, as discussed for example in the case of HR8799 (Wang et al. 2018). This could lead to a significantly different dynamical result.

As it stands, the observed disk structure could therefore be explained by the presence of at least one additional planet, with a mass still well below the current detection limit for high- contrast imaging, estimated at ~1 M_Jup (Lagrange et al. 2020). Recent observations with JWST (Kammerer et al. 2024) could potentially improve upon this limit. However, as this depends on the projected separation between the planet and the star, rather than the semi-major axis of the orbit, it is difficult to draw definitive conclusions about the presence of a planet. In the most favorable scenario, a planet with low eccentricity, and an apparent separation equal to its semi-major axis of 35 au at the time of observation, should be detectable if it has a mass of at least 0.3 M_Jup (the equivalent of Saturn). No planet of this type has been observed, but it is possible that the planet is in an unfavourable configuration, perhaps occulted by β Pic, or that it is completely absent. So far, the available JWST data do not yet allow us to reduce the detection limit below the current limit of ~1 M_Jup (Lagrange et al. 2020). We hope that future observations will enable us to lower this limit and strengthen the constraints on our model, or even detect the suspected planet directly.

It would now be interesting to examine the impact of a system with three or even four planets on the disk, in order to try to understand the dynamical origin of the various intriguing features observed so far, particularly the asymmetries (Kalas & Jewitt 1995; Heap et al. 2000) and the clump (Telesco et al. 2005; Dent et al. 2014; Matrà et al. 2017), through gravitational interactions with the planets.

Fig. 8 Example of a simulation of the dynamics of the β Pic planetary system with four planets, β Pic b (red), β Pic c (blue) and two additional planets (green and orange), and the disk of planetesimals (black). The initial orbital parameters of β Pic b and c are taken from Table 1. In this example, the additional planets have a mass of 0.05 MJup, an initial eccentricity of 0.05, and initial semi-major axes of 30 and 40 au respectively. The plotting conventions are identical to Fig. 1.

Acknowledgements

We are grateful to an anonymous referee for feedback that helped us improve this manuscript. All computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (COBREX; grant agreement nº885593). V. Faramaz acknowledges funding from the National Aeronautics and Space Administration through the Exoplanet Research Program under Grants No. 80NSSC21K0394 (PI: S. Ertel) and No. 80NSSC23K0288 (PI: V. Faramaz).

References

All Tables

All Figures

Fig. 1 Preliminary simulation of the dynamics of the β Pic planetary system with only the two planets known to date, β Pic b (red) and β Pic c (blue), and the disk of planetesimals (black). The initial orbital parameters of β Pic b and c are taken from Table 1. Left: upper view of the system. The planetesimals are depicted by small black dots, and the orbits of the planets by colored dashed lines. Middle: view of the system in terms of semi-major axis and eccentricity. The planets are here represented by colored points. Right: radial profile of the planetesimal disk (solid black line), superimposed to the models of Augereau et al. (2001) (blue dotted line) and Dent et al. (2014) (red dotted line). Note the clear mismatch in the location of the inner edge between this simulation and observational models.

In the text

Fig. 2 Accessible semi-major axis of the hypothetical additional planet β Pic d, prior to simulations, as a function of its eccentricity. Shaded areas are not accessible by β Pic d. Its semi-major axis must be smaller than the inner edge of the disk at 50 au (solid black line) to sculpt it up to this distance, while being greater than the stability limit with β Pic b (dashed black line) to avoid close encounters, based on the criterion of Petrovich (2015).

In the text

Fig. 3 Example of a simulation of the dynamics of the β Pic planetary system with three planets, β Pic b (red), β Pic c (blue), and an additional planet β Pic d (green), and the disk of planetesimals (black). The initial orbital parameters of β Pic b and c are taken from Table 1. In this example, the initial orbital parameters of β Pic d are: md = 0.2 MJup, ad = 35 au, and ed = 0.05. The plotting conventions are identical to Fig. 1.

In the text

Fig. 4 Combinations of β Pic d’s mass and semi-major axis within observational constraints that successfully reproduce the disk profile at 50 au. For several values of md, the corresponding acceptable ranges of ad are displayed in green. A power-law fit, according to Eq. (7), is overlaid in blue, taking into account the ranges of the semi-major axis. Shaded areas are not accessible byβ Pic d. The limits in semi-major axis are similarly presented as in Fig. 2. Additionally, above a mass of 1 MJup, corresponding to the observational limit (solid red line), β Pic d would have already been detected in previous observations (Lagrange et al. 2020). And below 0.15 MJup, corresponding to the ring limit (dashed red line), β Pic d does not completely clear the inner zone.

In the text

Fig. 5 Example of a simulation of the dynamics of the β Pic planetary system with three planets, β Pic b (red), β Pic c (blue) and an additional planet β Pic d (green), and the disk of planetesimals (black). The initial orbital parameters of β Pic b and c are taken from Table 1. In this example, the initial orbital parameters of β Pic d are: md = 0.05 MJup, ad = 40 au, and ed = 0.05. The plotting conventions are identical to Fig. 1. This combination still generates a disk inner edge at 50 au, but there is now enough space between β Pic b and β Pic d to allow an additional ring of particles to exist at ~30 au.

In the text

Fig. 6 Combinations of β Pic d’s mass and semi-major axis for various eccentricities, within observational constraints that successfully reproduce the disk profile at 50 au. The plotting conventions are identical to Fig. 4. The ring limit quickly tends to zero as the eccentricity increases. The minimum semi-major axis stability limit with β Pic b imposes a maximum mass limit (dot-dashed black line), determined by the intersection (black dot) with the fit’s right extremum of the acceptable range of β Pic d’s semi-major axis. This maximum mass limit decreases as eccentricity increases, becoming more restrictive than the observational mass limit for ed ≳ 0.25, and tends towards zero as eccentricity continues to increase (see Fig. 7).

In the text

	Fig. 7 Accessible mass of the potential β Pic d as a function of its eccentricity ed, based on the results of Fig. 6, whose plot conventions are similar.
In the text

Fig. 8 Example of a simulation of the dynamics of the β Pic planetary system with four planets, β Pic b (red), β Pic c (blue) and two additional planets (green and orange), and the disk of planetesimals (black). The initial orbital parameters of β Pic b and c are taken from Table 1. In this example, the additional planets have a mass of 0.05 MJup, an initial eccentricity of 0.05, and initial semi-major axes of 30 and 40 au respectively. The plotting conventions are identical to Fig. 1.

In the text