© The Authors 2024

1 Introduction

The southern A6V (Gray et al. 2006) star β Pictoris (β Pic) has been the focus of intense research for decades, initially for its wide edge-on debris disk (Smith & Terrile 1984). The exact age of this young main-sequence star has been the subject of controversy. The most robust determinations (22 ± 6 Myr, Binks & Jeffries 2014; Shkolnik et al. 2017; ${18.5,}_{-2.4}^{+2.0}\text{\hspace{0.17em}Myr}$, Miret-Roig et al. 2020) are based on the young moving group it belongs to. In the following we assume an age of 20 Myr as a typical standard value.

The presence of hidden planets in the circumstellar disk was suspected very early, as a natural outcome of the disk evolution. Specific clues suggesting the presence of at least one giant planet were identified. First, the warped profile of the disk as well as its various asymmetries (Kalas & Jewitt 1995; Heap et al. 2000) were tentatively attributed to the perturbing action of a planet (Mouillet et al. 1997; Augereau et al. 2001). Similarly, the photometric variations observed in 1981 (Lecavelier Des Etangs et al. 1995) were also tentatively attributed to a planetary transit (Lecavelier Des Etangs et al. 1997). Finally, the repeated transient spectral variations monitored in the β Pic absorption spectrum, interpreted as resulting from the evaporation of star-grazing exocomets (Beust et al. 1996) (see details below), were also attributed to the perturbing action of a giant planet (Beust & Morbidelli 1996,2000; Thébault & Beust 2001). These independent studies all agreed on a giant planet orbiting β Pic at ~10 au.

This planet, today known as β Pic b, was indeed detected in high-contrast imaging first with NaCo (Lagrange et al. 2010), and later with GPI (Wang et al. 2016) and SPHERE (Lagrange et al. 2019a). It was thus immediately tempting to identify β Pic b with the previously suspected planet. Before being able to give a robust conclusion, an accurate knowledge of β Pic b’s orbit was necessary. This was done thanks to a regular astrometric follow-up of the planet (Chauvin et al. 2012; Bonnefoy et al. 2014; Wang et al. 2016; Lagrange et al. 2019a). Currently, these studies all converge towards a semi-major axis distribution peaking close to 9–10 au and an orbital period of ~20 yr. The eccentricity is harder to constrain, but it is presumably small.

The mass of β Pic b is less easy to constrain. Fitting the relative astrometric data, it is only possible to fit the total dynamical mass of the system, which falls in the range 1.75–1.85 M_⊙ (Wang et al. 2016; Dupuy et al. 2019). This is in fact no more than a fit of the central star’s mass. Photometric constraints deduced from direct imaging and models, combined with radial velocity data, helped to derive an estimate around 10 M_Jup. Based on a joint analysis of Gaia and HIPPARCOS data, Snellen & Brown (2018) report a mass of 11 ± 2 M_Jup, while Dupuy et al. (2019) derive ${13.1}_{-3.2}^{+2.8}{M}_{\text{Jup }}$.

All these characteristics were closely compatible with those deduced from the former theoretical studies listed above for the hypothetical planet orbiting β Pic, so that the whole picture appeared clear. Meanwhile, a careful analysis of the HARPS radial velocity of the star led to the identification of another massive planet termed β Pic c orbiting the star at 2.7 au (Lagrange et al. 2019b) (i.e. well inside β Pic b’s orbit). Thanks to high-precision astrometric data from the GRAVITY interferometer, that planet was directly confirmed (Nowak et al. 2020), so that its reality is now certain. All orbital determinations aim now at fitting the two planetary orbits together (Lagrange et al. 2020; Nowak et al. 2020; Lacour et al. 2021), making use of available astrometric data for both planets and radial velocity data of the star. Currently, β Pic c appears slightly less massive (~8 M_Jup) than β Pic b, but significantly more eccentric, with an eccentricity around 0.3.

Gas in the circumstellar disk of β Pic was detected very early in the stellar spectrum (Hobbs et al. 1985), thanks to the edge-on orientation of the disk. Repeated observations of various metallic species (Ca II, Mg II, Al III, Vidal-Madjar et al. 1994) rapidly revealed frequent and transient changes appearing as additional absorption components in the bottom of the lines, most of the time redshifted by tens to hundreds of km s⁻¹ with respect to the stellar velocity, and with variation timescales of a few days or even less (Ferlet et al. 1987; Lagrange-Henri et al. 1992; Petterson & Tobin 1999; Kennedy 2018; Tobin et al. 2019). These events have been convincingly modelled as the result of the sublimation of star-grazing exocomets in the vicinity of the star, each transient event corresponding to the passage across the line of sight of a single body (Beust et al. 1990, 1996; Crawford et al. 1998). This scenario has been termed falling evaporating bodies (FEBs), and it implies an activity of hundreds of such FEBs per year having star-grazing periastron passages. These bodies do not actually graze the stellar surface at each periastron passage. They only need to be close enough to the star to allow their refractory compounds to sublimate and produce the observed absorption components. As shown by Beust et al. (1998), in the environment of an A-type star like β Pic, this occurs below a threshold distance of ~0.5 au (i.e. far above the stellar surface). For comparison purposes, this threshold distance around the Sun is no more than a few solar radii. Any exocomet orbiting β Pic that gets closer to the central star than this distance is susceptible to becoming a FEB. The various subsets or families of variable components (in terms of velocity, depth, and variation timescales) have been successfully interpreted as FEBs crossing the line of sight at different distances within the quoted threshold, down to a few stellar radii only (Beust et al. 1996; Kiefer et al. 2014). The reality of this scenario was recently reinforced by the detection of similar events in transit photometry with the TESS satellite (Zieba et al. 2019; Lecavelier des Etangs et al. 2022). A detailed analysis of the size distribution of the exocomets even revealed a differential distribution close to that of the collisional equilibrium.

In this framework, the redshift (or blueshift) velocity of each transient component corresponds to the projection of the FEB’s orbital velocity onto the line of sight at the time it crosses it, as what we observe is basically a cloud of metallic ions surrounding it. If we assume that these bodies move on very elongated orbits that can be assimilated to parabolic around periastron passage, this velocity mainly depends on the periastron distance of the FEBs and on the longitude of their periastron with respect to the line of sight. For instance, any redshifted (resp. blueshifted) component corresponds to a FEB crossing the line of sight before (resp. after) periastron. The large number of FEB events recorded (see e.g. Tobin et al. 2019) subsequently allows one to draw statistics on the FEB population (Kiefer et al. 2014). The most obvious outcome is the statistical difference between redshifts and blueshifts. Although a few blueshifted events have been reported, the variable events appear most of the time on the red wing of the spectral lines. This shows up in all observational campaigns quoted above. This indicated a non-uniform distribution of the longitudes of periastra. The FEBs actually seem to come from one specific side of the disk (Beust et al. 1990). Further detailed studies were able to identify several subfamilies with different orbital characteristics (Beust et al. 1998; Kiefer et al. 2014), but this main result still holds.

These constraints challenge any dynamical model willing to explain this high number of star-grazers. The common feature of all models is that planetary perturbations are needed to trigger the FEB phenomenon. For instance the Kozai-Lidov resonance, often presented as a source of star-grazers in the Solar System (Bailey et al. 1992), cannot be invoked here because of its natural rotational invariance. In this context, there should be as numerous redshifted as blueshifted events. Conversely, mean-motion resonances (MMRs) represent a potential anisotropic source of stargrazers. This idea was first investigated by Yoshikawa (1989), who showed that fictitious bodies trapped in some inner MMRs such as 4:1, 3:1, 5:2, 7:3, and 7:4 with a moderately eccentric giant planet can undergo large resonant eccentric increases. This idea was further developed by Beust & Morbidelli (1996, 2000) and Thébault & Beust (2001) in the context of the FEB model for β Pic. Semi-analytical and numerical studies showed that this model was indeed able to efficiently trigger the FEB phenomenon and reproduce nearly all their observed dynamical characteristics, provided the planet assumes an eccentricity e ≳ 0.05, and has a suitable value $\varpi =-{70}^{°}\pm {20}^{°}$ of longitude of periastron with respect to the line of sight. Detailed modelling (Beust & Morbidelli 2000; Thébault & Beust 2001) led to the conclusion that the planet was presumably orbiting the star at ~ 10 au within 50% uncertainty. The mass of the perturbing planet was less easy to constrain, as the topology of the secular dynamics involved is nearly independent of it. Recently, this dynamical mechanism was theoretically re-investigated by Pichierri et al. (2017) who confirmed its general efficiency, showing that it is still active even if the perturbing planet is significantly more eccentric.

This resonant mechanism could apply in other systems than β Pic. The Kirkwood gaps in the solar asteroid belt have been emptied by a similar process in its early evolution (Morbidelli & Moons 1993; Moons & Morbidelli 1995; Farinella et al. 1994). The β Pic system presumably represents a likely analogue to the young Solar System.

As mentioned above, the first detected planet β Pic b closely matches the characteristics of the putative planet of 2001. However, the presence of the second planet β Pic c orbiting inside complicates the picture. Any FEB progenitor initially trapped in an inner MMR with β Pic b¹ orbits initially outside β Pic c’s orbit. During the resonant eccentricity increase process that drives it to the FEB state, it inevitably comes to regularly cross that orbit as the semi-major axis only undergoes small changes (Pichierri et al. 2017). Subsequent close encounters with β Pic c could presumably eject it well before reaching the FEB state. Beust & Morbidelli (2000) already investigated the possible perturbing action on FEB progenitors of a second planet orbiting inside the first one, but that planet was supposed to be Earth-sized. Such a small planet was indeed able to divert some of the FEB progenitors from their resonant dynamical route, but most of them were still able to safely reach the FEB state. Here β Pic c is much more massive and could actually kill the resonant mechanism.

In the present paper we re-investigate the resonant FEB scenario based on the present-day knowledge of the dynamical configuration of the β Pic planetary system. In Sect. 2 we first investigate the stability and dynamical perturbations of the two-planet system, taking as input basis the latest orbital determination by Lacour et al. (2021). Then in Sect. 3 we add a disk of test particles to the simulation to see where FEB progenitors can be stable. Then we focus on the surviving part of the disk inside β Pic c’s orbit, showing that it can be a valuable source of FEBs in the presence of both planets, leading to a revision of the scenario. In Sect. 4 we discuss the consequences of the revised model. Our conclusions are presented in Sect. 5.

Fig. 1 Secular evolution over 108 yr of the β Pic b+c two-planet system corresponding to solution #2 from Table 1. From left to right: semi-major axes, eccentricities, and mutual inclination. The logarithmic time axes highlight the short-term as well as the long-term behaviour. Red always stands for β Pic c and green for β Pic b.

2 Two-planet secular dynamics

The two-planet system orbiting β Pic is assumed to be dynamically stable. Its secular evolution is briefly mentioned in Lacour et al. (2021), revealing stability despite non-negligible eccentricity oscillations. Interestingly, Lacour et al. (2021) points towards a possible 7:1 MMR between the two planets. To investigate this issue, we integrated the system starting from the orbital solutions given by Lacour et al. (2021). More specifically, we integrated two configurations, listed in Table 1. Solution #1 is the best one from Table 3 of Lacour et al. (2021) computed using all available data; solution #2 is the one labelled ‘truth’ from Table B.1 of Lacour et al. (2021), which is basically an average of all the possible values. Many other configurations could be tested, but these can be considered representative. The two solutions are obviously close to each other, but the main difference between them appears if we compute the ratio of the two orbital periods P_b/P_c in each case. We find P_b/P_c = 7.0097 for solution #1 and P_b/P_c = 7.466 for solution #2. Hence, solution #1 is presumably locked in 7:1 MMR, whilst solution #2 is not.

As the solutions of Lacour et al. (2021) are given in Jacobi coordinates, the integrations were carried out over 10⁸ yr (i.e. far above the current age of the star) using the symplectic N-body code SWIFT_HJS (Beust 2003), which naturally works in that coordinate system. We note that prior to all computations, all angular orbital elements have been converted relative to the invariable plane of the three-body system as a more relevant reference system. Figure 1 shows the evolutions of the semi-major axes, the eccentricities of the two planets, and their mutual inclination for solution #2 from Table 1. We note a remarkable stability of the semi-major axes, indicating a stable system. Both eccentricities and mutual inclination exhibit the regular oscillations characteristic of a non-resonant system with regular dynamics. The longitudes of ascending nodes and arguments of periastra (not shown here) undergo very regular precession motions.

Figure 2 shows the same evolution, but starting with solution #1 from Table 1. We do not show the general evolution of the semi-major axes as this appears identical to Fig. 1. We note however a more irregular evolution of the eccentricities and of the mutual inclination. This is a consequence of the 7:1 MMR. To check this point, we also show the evolution of the critical argument σ_7:1 of that resonance, which is defined as $${\sigma }_{7:1}=\frac{7}{6}{\lambda }_b-\frac{1}{6}{\lambda }_c-{\varpi }_c,$$(1)

where λ_c = is β Pic c’s mean longitudes (and similarly for β Pic b) and ${\varpi }_c={\text{Ω}}_c+{\omega }_c$ is its longitude of periastron (see e.g. Beust 2016, for a more general definition). The quantity σ_7:1 is presumably a slowly variable quantity thanks to the MMR. Resonant configurations are expected to be characterized by a libration of σ_7:1, whilst non-resonant ones should exhibit circulation of σ_7:1.

The evolution of σ_7:1 over the first 50 000 yr of the simulation shows a hybrid configuration. We first note a circulation up to ~18 000 yr, and then a libration motion that lasts approximately 20 000 yr before the circulation starts again. What we see here is a temporary capture in 7:1 MMR. This situation appears to repeat itself all along the simulation, with a 20 000 yr resonant capture every ~40 000 yr. On average, the two-planet system appears to spend nearly half of its time within the resonance and half out of it. This hybrid behaviour could still be artificial and due to some numerical aliasing in Fig. 2 as a consequence of an insufficient time resolution of the plot. This is actually not the case. Orbital data were saved every 100 yr in the simulations and plotted with the same resolution (i.e. much higher than the frequency of the σ_7:1 oscillations).

Moreover, the hybrid behaviour is not specific to the configuration tested. We checked many other solutions close to solution #1 from Table 1, in particular changing β Pic c’s semi-major axis to match the exact 7:1 MMR value. We found no configuration where the system permanently remains locked in the MMR. In some cases σ_7:1 always circulates, which denotes a non-resonant configuration, and thus similar to solution #2. In all other cases we have temporary resonant captures as in Fig. 2, so that the two secular evolutions we detail here are representative for all possible situations.

In Fig. 2 we also show the detailed evolution of the semi-major axes of the two planets over the first 50 000 yr. We note rapid very small amplitude changes when the system is out of the 7:1 MMR, and more pronounced oscillations in the resonant phase. This behaviour is characteristic of a temporary resonant capture. Out of any MMR, the semi-major axes undergo no secular variations apart from minor short-term changes. These rapid changes would also show up in Fig. 1 if we focused on the detailed evolution of each semi-major axis. Conversely, resonant configurations are known to generate secular oscillations of semi-major axes coupled with eccentricity and σ oscillations (Beust 2016). This is exactly what we observe in Fig. 1.

Hence, we confirm that solution #1 corresponds to a mode of regular temporary captures of the two-planet system in 7:1 MMR. Given the high masses of the two planets, this resonance (a sixth-order MMR) is actually too weak to sustain itself permanently.

Fig. 2 Secular evolution over 108 yr of the β Pic b+c two-planet system corresponding to solution #1 from Table 1. Upper row, from left to right: eccentricities, mutual inclination, and (over the first 50 000 yr) the characteristic resonant argument σ7:1 (see text). Our definition of σ7:1 implies a 60° degeneracy. In our calculation, we always keep the values between 0 and 60°. Hence the variations from 0 to 60° appearing here actually mean circulation of σ7:1. Lower row: short-term semi-major axis variations of the two planets over the first 50 000 yr.

3 Simulations with a disk of planetesimals

3.1 Global simulation

In this section we deal with presenting simulations involving the two planets β Pic b and β Pic c and a disk of planetesimals. The main goal is to investigate where planetesimals can dynamically survive in this environment, and whether the suspected FEB-generating mechanism is still active. The subsequent simulations assume a starting disk of 400000 planetesimals treated as massless particles with initial semi-major axes randomly drawn between given bounds that depend on the simulation. Eccentricities are drawn between 0 and 0.05, and inclinations between 0 and 2° with respect to the invariable plane of the two-planet planetary system. Other orbital elements such as longitudes of ascending nodes, longitudes of periastra, and initial mean longitudes are randomly drawn between 0 and 360°. Computations are carried out using the RMVS3 symplectic integrator (Levison & Duncan 1994), a version of the initial mixed-variable sym-plectic integrator (Wisdom & Holman 1991) that better treats close encounters. We note that from a technical point of view, RMVS3 works in heliocentric coordinates, whilst the orbits listed in Table 1 are given in Jacobi coordinates. Hence, prior to launching any simulation, we first convert the planetary orbital elements of Table 1 to heliocentric. This actually only affects β Pic b, as for β Pic c (the innermost planet) the two sets of elements coincide; however, even for β Pic b the difference is small as the stellar mass dominates.

Figure 3 shows the result of such a simulation after 10 Myr assuming orbital solution #2 (non-resonant) from Table 1 for the planets, and a disk of planetesimals ranging initially between 0.5 au and 45 au. The remaining disk is shown as a planar upper view and in (semi-major axis a, eccentricity e; herefter a, e) space. The planetesimals appear as black dots. The disk appears strongly depleted interior to ~25 au. Most of the planetesimals initially located closer in have been removed by perturbations and close encounters due to the planets. The inner edge of the disk appears strongly carved by various MMRs with β Pic bas all subsisting rings inside ~27 au actually correspond to these MMRs: the 2:3 and 2:5 at a = 13.03 au and a = 18.33 au (marginally); the 1:3 at a = 20.7 au, the 2:7 at a = 22.94 au, and the 1:4 at a = 25.07 au. Particles trapped in these resonances subsist here thanks to a phase-protection mechanism that prevents close encounters with β Pic b. The resonant dynamics excites the eccentricities of these particles, so that their mean eccentricity is significantly larger than those from the main surviving non-resonant disk beyond 25 au. This particularly shows up in the (a, e) map of Fig. 3. A ring of more eccentric particles, actually corresponding to the 1:5 MMR with β Pic b at a = 29.09 au, even shows up in the middle of the unperturbed disk.

It is important to note that no particle was able to survive between the two planets. In the context of the FEB-generating mechanism by MMRs with β Pic b, this is problematic as the model described by Beust & Morbidelli (1996, 2000) and Thébault & Beust (2001) pointed out the major inner (4:1, 3:1) MMRs with that planet as their probable source. These MMRs all fall in the region between β Pic b and β Pic c. For instance, the 4:1 MMR, which appeared to be the most powerful source of FEBs, is located close to 4 au where no planetesimal appears dynamically stable over 10 Myr (and actually much less). It thus turns out that even if β Pic b appears to almost perfectly match the giant planet suspected more than 20 yr ago to be responsible for the FEB mechanism, the presence of β Pic c orbiting inside invalidates this picture.

The outer MMRs with β Pic b quoted above that carve the outer planetesimal disk could not also constitute valuable sources of FEBs. It actually shows up in Fig. 3 that even if the eccentricity of the resonant particles is enhanced by the resonant mechanism, none of them is able to reach a high enough eccentricity to become a FEB (e ≳ 0.96 would be required here). This may also be due to the presence of β Pic c orbiting inside that acts as a barrier.

Fig. 3 Outcome after a 10 Myr simulation of a disk of planetesimals initially between 0.5 and 45 au in the environment of β Pic b and β Pic c taken from solution #2 from Table 1, viewed from above (left), and displayed in (a, e) space (right). The red bullets show the location of the planets and the star, and the black dots are the remaining planetesimals.

3.2 Focus on the inner ring

3.2.1 Technical issue

In Fig. 3, however, another ring of surviving particles is present inner to ßPic c’s orbit. We focus here on this inner ring. To achieve a better spatial resolution in this area, a new simulation was initiated, now taking 400 000 particles between 0.5 and 1.5 au, as the previous simulation had revealed that no particle was able to survive beyond 1.5 au due to the presence of β Pic c. The inner bound of 0.5 au was fixed according to the evaporation limit of refractory compounds determined by Beust et al. (1998). Hence no planetesimal was expected to sustainably orbit the star closer to this threshold. This is of course an approximate limit that depends actually on the real composition of the bodies under consideration. Beust et al. (1998) shows for instance that more carbonaceous FEBs may be able to survive closer to the star. The simulation was also extended up to 20 Myr (i.e. the present-day estimated age of β Pic).

Technically speaking, this new simulation was initially carried out with the same integrator (RMVS3) as before, but although the results were satisfying at first glance, they also appeared inaccurate. Many particles around ~1 au underwent a secular drift in semi-major axis that seemed spurious. Other attempts with reduced time-steps revealed that this secular drift actually depends on the time-step, proving its numerical origin. Nonetheless, it was not possible to find any convenient time-step value able to eliminate this effect. After analysis, the origin of this phenomenon was identified as inherent to the way close encounters are treated in RMVS3. Whenever a test particle gets sufficiently close to a massive planet, RMVS3 automatically reduces the time-step for that particle during the encounter to better resolve it. This temporarily breaks the symplecticity of the integration, but does not affect the global relevance of the integration, as individual particles undergo close encounters only occasionally. This statement applies to the simulation described in Fig. 3, in particular for the outer part of the disk. Due to the high mass of β Pic c and to its fairly high eccentricity, particles orbiting around ~ 1 au regularly get sufficiently close to β Pic c (near conjunction) to be considered as having a close encounter with it by RMVS3 without being ejected. However, as this occurs quite often, the integration for such particles is never symplectic. This is the origin of the semi-major axis drift observed.

To overcome these difficulties, we changed the integrator. We note that adopting other integrators that also temporarily break the symplecticity of the integration, such as MERCURY (Chambers 1999) and REBOUND (Rein & Liu 2012), would not help much more. We tried SYMBA (Duncan et al. 1998), a more sophisticated integrator that remains strictly symplec-tic while handling close encounters. This also turned out to be inaccurate, but for another reason. Due to the use of the democratic heliocentric method (DH), SYMBA hardly resolves orbits with very small periastron, which is exactly what we have here with FEBs. Therefore, the computation became extremely long, although the semi-major axis drifts were no longer present. Conversely, the classical second-order Wisdom-Holman (WH) mapping naturally handles such orbits, but is not able to compute close encounters. We finally moved to a higher-order symplec-tic integrator than WH with no change in time-step during close encounters, and modified WH accordingly. We adopted the sixth-order S6B scheme described by Chambers & Murison (2000), assuming a time-step of 1.33 × 10⁻² yr (i.e. 1/20 of the smallest orbital period). This way semi-major axis drifts were no longer present, and thus we can trust our integrations.

Fig. 4 Outcome in (a, e) space of a disk of planetesimals initially between 0.5 and 1.5 au, under the same conditions as in Fig. 3, after 2 Myr (left) and 20 Myr (right). Both planets are outside the panel.

3.2.2 Results

Figure 4 presents the result of this new simulation, a view of the inner ring in (a, e) space after 2 Myr and 20 Myr. It appears severely carved down to ~1 au by perturbations arising mainly from β Pic c located at 2.6 au. Up to ~ 1 au, most disk particles appear to still be present after 20 Myr, but many of them have reached high eccentricity values, thus becoming FEB candidates. We note that for a planetesimal starting at ~0.8 au, reaching the FEB state with a periastron ≲0.4 au only requires an eccentricity e ≳ 0.5, which is easily realized for many particles here. As can be seen in Fig. 4, the eccentricity increase is preferably done along vertical lines (i.e. at specific semi-major axis values) that can be easily identified with many inner MMRs with β Pic c, for example the 4:1 at a = 1.030 au, the 9:2 at a = 0.952 au, the 5:1 at a = 0.888 au, the 6:1 at a = 0.786 au, the 7:1 at a = 0.709 au, the 8:1 at a = 0.649 au, the 9:1 at a = 0.600 au. The process is highly active at t = 2 Myr, but it is still present at t = 20 Myr.

Figure 5, to be compared to Fig. 4, presents the result of another simulation, but where solution #1 from Table 1 (partially resonant) for the initial orbital configurations of β Pic b and β Pic c is assumed instead of solution #2 in Fig. 4. There are actually minor differences between the two simulations, but in both cases the disk appears carved with very similar shapes down to 1–1.1 au. The number of particles reaching high eccentricities at 20 Myr is comparable, with slight differences in the distribution between different MMRs. The disk appears slightly more carved close to 1 au in the non-resonant case. This is due to the semi-major axis difference for β Pic c between the two solutions (Table 1), which causes a more important carving in the non-resonant configuration. In the resonant case as well the resonances appear somewhat more distinctly. Both simulations nevertheless show that MMRs with β Pic c constitute potential sources of FEBs. Figure 6 displays the location of all these MMRs. We note that some very high-order MMRs, such as 13:2 and 16:3, are active FEB sources.

This study drives us towards a new picture of the FEB generation mechanism in the β Pic disk. MMRs are still the preferred route as, like Kozai resonance, they can constitute a long-term source of stargrazers and, contrary to Kozai resonance, they trigger an asymmetric infall (Beust & Morbidelli 1996, 2000; Thébault & Beust 2001) that better matches the observational statistics (Beust et al. 1996). We were previously considering inner MMRs with a single giant planet similar to β Pic b, with FEB progenitors starting around 4–5 au. The presence of β Pic c invalidates this picture. We are now considering planetesimals trapped in MMRs with β Pic c, thus starting from a region much further in, between ~0.5 au and ~1 au. Moreover, the MMRs involved in the previous model were mainly 4:1,3:1, and 5:2. We consider here much higher-order MMRs with β Pic c, for example 7:1 and 8:1. This is due to β Pic c’s high mass and eccentricity. Under such conditions, lower-order MMRs located further out, and thus closer to β Pic c fall in the chaotic unstable region. The 4:1 MMR, located slightly outside 1 au at the outer edge of the surviving disk of planetesimals, is marginally active at 2 Myr, but no longer at 20 Myr. Conversely, higher-order MMRs located closer to the star appear now as sources of FEBs thanks to β Pic c’s high eccentricity. All our previous studies (Beust & Morbidelli 1996, 2000; Thébault & Beust 2001; Beust & Valiron 2007) considered a low-eccentricity perturber with e ≲ 0.1. Pichierri et al. (2017) showed that the mechanism remains active in higher-eccentricity regimes. Actually, the higher the planet’s eccentricity, the more efficient the mechanism. The high-order MMRs considered here do not constitute efficient sources of FEBs for e ≲ 0.1, but do with e = 0.33, as for β Pic c here.

Fig. 5 Same as Fig. 4, with the same initial disk of particles between 0.5 and 1.5 au, but starting from solution #1 from Table 1 for β Pic b and β Pic c.

3.3 Taking evaporation into account

Figures 4 and 5 could appear misleading as they actually show potential FEB candidates moving in MMRs with β Pic c at t = 2 Myr and at t = 20 Myr, but as long as bodies get trapped in the resonant eccentricity increase process and get closer to the star, they start to evaporate and may be quickly fully destroyed. Evaporation was not taken into account in our dynamical simulation as this would have resulted in an excessively long computation time and, in first approximation, the evaporation process does not affect the dynamics. It would nevertheless be relevant to determine whether particles moving in MMRs at 20 Myr in the preceding simulations are virtually already evaporated or not as they may account for FEB statistics in the latter case only.

Computing gradual evaporation of FEBs along successive periastron passages is not straightforward as this may depend highly on its chemical composition, porosity, and other parameters. To date, the most complete study of evaporation process of exocomets in the vicinity of β Pic has been made by Karmann et al. (2003). It shows that, on average, most of the sublimation of refractories occurs closer to ~0.2 au. The process actually starts below ~0.4 au (Beust et al. 1998). To simplify the treatment, we assume a surface evaporation rate ∝ 1/r¹⁶ for r < 0.4 au. This actually represents the best power-law fit of the production rate curves of Karmann et al. (2003) as long as the simulated FEBs are not too destroyed (i.e. for r ≳ 0.2 au). More specifically, we consider a reference mass evaporation rate z₀ for a reference body of radius s₀ = 10 km at reference distance r₀ = 0.2 au, defined as $$\alpha =\frac{z_0{r}_0^{16}}{4\pi s_0^2}.$$(2)

Then, the mass loss rate of any body of radius s at distance r reads $$\frac{\text{d}m}{ \text{d}t}=\frac{4\pi \alpha s^2}{r^{16}}$$(3)

and the radius decrease of the body is $$\frac{\text{d}s}{ \text{d}t}=\frac{1}{4\pi s^2\rho }\frac{\text{d}m}{ \text{d}t}=\frac{\alpha }{\rho r^{16}},$$(4)

where ρ is the bulk density of the material. Integrating this expression over one periastron passage, we get $$\text{Δ}s=2{\displaystyle {\int }_0^{v_{\max }}\frac{\text{d}s}{ \text{d}t}}\times \frac{\text{d}v}{ \text{d}v/\text{d}t},$$(5)

where Δs is the full radius loss over one periastron passage, v is the true anomaly along the orbital path, and v_max is the value of v corresponding to r = r_max = 0.4 au, the maximum evaporation distance of refractories, $$1+e\cos v_{\max }=\frac{a\left(1-e^2\right)}{r_{\max }}$$(6)

where a is the semi-major axis and e the eccentricity. Using then Kepler’s second law to obtain dν/dt, the integral becomes $$\text{Δ}s=\frac{2\alpha }{n\rho a^{16}}\frac{1}{{\left(1-e^2\right)}^{29/2}}{\displaystyle {\int }_0^{v_{\max }}{(1+e\cos v)}^{15}} \text{d}v,$$(7)

where n is the mean motion of the orbit related to a via Kepler’s third law. The integral can be calculated in closed form as a polynomial expression involving the eccentricity and ν_max. We assume z₀ = 10⁷ kg s⁻¹ to be in agreement with the results of Karmann et al. (2003), and integrated the radius loss (when relevant) over successive periastron passages for all particles involved in our simulations. This is a simplified model that deserves better treatment, but it is a convenient way to derive an estimate of the importance of evaporation from the computed dynamical evolution.

Figure 6 shows the result of this computation. Both plots are actually identical to the 20 Myr plots in Figs. 4 and 5, but dead particles (i.e. those that have been computed to be already fully evaporated at that time) are highlighted here in blue. In addition, we show the location of various MMRs with β Pic c with red vertical bars. This shows how the inner disk is sculpted by those MMRs. We note three major facts. First, many particles close to the inner edge of the disk at 0.5 au are already dead at 20 Myr. This should not be surprising, as the average eccentricity they reach by secular perturbations by the two planets lets them constantly get below 0.4 au and evaporate. Second, many particles that have been involved in the resonant eccentricity process are also dead. This also is not surprising, as they regularly enter the evaporation zone that eventually drives them to destruction. This is the real FEB phenomenon simulated here. Third, even if many resonant particles have already disappeared at t = 20 Myr, some of them are still present (black dots in the upper part of the plots). This shows that the FEB phenomenon is still active at t = 20 Myr.

The survival of bodies after 20 Myr could appear surprising as we could expect evaporation to act more quickly. This is actually not true close to the inner edge of the disk. A numerical application of Eq. (7) gives ∆s ≃ 10⁻⁵−10⁻⁴ m at 0.4 au. Therefore, assuming an orbital period of 0.27 yr at a = 0.5 au, the whole radius loss after 20 Myr ranges between 800 m and 8 km. This is enough to allow some of the bodies concerned to be fully evaporated, but not all, which is what Fig. 6 confirms. Evaporation is conversely much more efficient for resonant bodies when they reach high eccentricities (∆s ≃ 1 m for a periastron of 0.2 au). None of these bodies should be able to survive 20 Myr. Hence, the presence of living bodies at high eccentricities in the MMRs at that time should be surprising.

The reason why some living bodies appear at high eccentricity even after 20 Myr can be seen in Fig. 7, which displays the temporal evolution of the periastron of two particles out of the simulations. We adopt here the same convention as in Fig. 6 (i.e. the orbital path appears in blue when the bodies are fully evaporated). These parts of the curves should therefore be considered virtual, as the corresponding bodies are already dead. Both particles are involved in the resonant process, but the decrease in periastron, a consequence of the eccentricity increase process, does not occur at the same time. The first particle (left plot) gets quickly involved in that process and disappears. The second particle (right plot) remains at low eccentricity (i.e. in the main disk) for almost all the simulation, and gets trapped in the same process only close to 20 Myr and subsequently disappears. Hence fresh bodies permanently enter the FEB regime at any time and contribute to sustaining the FEB activity even even after 20 Myr.

There are several reasons why some particles enter the FEB process after a long time. First, Beust & Morbidelli (2000) showed (with the 4:1 resonance) that not all resonant particles are subject to the eccentricity increase process, but only those having a small enough resonant libration amplitude. Particles that have initially large libration amplitudes may stay at low eccentricity in the MMRs and never get involved in the FEB process, or possibly after some delay. Second, the location of the β Pic c resonances themselves is subject to secular changes, thanks to mutual perturbations between the two planets. Consequently, particles that are initially not resonant may then be captured and enter the FEB regime, possibly after several million years.

We also note from Fig. 6 that some of the highest-order MMRs involved in the FEB mechanism seem to still be at the beginning of this process, even after 20 Myr. This is for instance the case of the 13:2 and 16:3 MMRs, which seem not to have fully carved the corresponding zone in the disk and to be sending particles at high eccentricity that are still active. The reason is that the eccentricity increase process associated with those very high-order and weak MMRs is very slow compared to stronger ones, such as 7:1 and 8:1. Hence, after 20 Myr, the associated resonant particles have still not evaporated.

Finally, in addition to those various mechanisms, mutual scattering between planetesimals can also contribute to continuously replenishing the MMRs. This additional effect is not simulated here as our planetesimals are treated as massless particles.

Fig. 6 Same as Figs. 4 and 5 at t = 20 Myr, but taking into account gradual evaporation of FEBs, with solution #2 (left plot) and #1 (right plot) from Table 1. All particles appearing in blue are virtually already evaporated and should not account for the FEB statistics at that time. The locations of the main MMRs with β Pic c are indicated with red vertical bars.

Fig. 7 Temporal evolution of the periastron value of two individual bodies out of the simulations from Figs. 4 and 5. The colour convention follows that of Fig. 6: the evolution appears in blue when the corresponding body is already evaporated.

4 FEBs velocity statistics

4.1 Theoretical background and semi-analytical study

In Beust & Morbidelli (2000) and Thébault & Beust (2001) we showed that the resonant mechanism supposed to be the source of FEBs towards β Pic was able to trigger a non-axisymmetric infall of FEBs, thus generating asymmetric statistics between blueshifted and redshifted spectral absorption events. In the single-planet model, this was fully constrained by the orientation of the planet’s periastron. In comparison with the Doppler velocity statistics of observed events, which appear strongly biased towards redshifts (Petterson & Tobin 1999), it was even possible in the framework of that model to deduce a suitable value for the planet’s longitude of periastron with respect to the line of sight that satisfactorily reproduced that statistics. This ability to trigger non-axisymmetric FEB velocity statistics was the main argument in favour of the resonant origin for FEBs in contrast with Kozai-Lidov resonance.

With the new model involving both planets, β Pic b and β Pic c, and FEBs arising from higher-order MMRs with β Pic c, the picture is somewhat more complex, and it is worth wondering about the FEB velocity statistics it should generate (i.e. the radial velocities of all potential FEBs appearing in Figs. 4 and 5 at the time they cross the line of sight). Following previous studies and Pichierri et al. (2017), we first performed a semi-analytical analysis of the phenomenon. In Beust & Morbidelli (1996), portrait phases of the secular resonant Hamiltonian indeed showed how resonant particles reach high eccentricities with constrained longitudes of periastron. We then tried to build similar maps adapted to the present model, and considered a massless particle orbiting β Pic, perturbed by both planets, β Pic b and β Pic c.

The instantaneous interaction Hamiltonian reads $$H=-\frac{GM}{2a}+H_b+H_c$$(8)

with $$\begin{array}{l}{H}_c=-G{m}_c\left(\frac{1}{\left|r-r_c\right|}-\frac{r\cdot r_c}{r_c^3}\right)\text{\hspace{1em}and}\\ H_b=-G{m}_b\left(\frac{1}{\left|r-r_b\right|}-\frac{r\cdot r_b}{r_b^3}\right),\end{array}$$(9)

where M is the mass of the star; m_c and m_b are the masses of the two planets; and r, r_c, and r_b are respectively the heliocentric radius vectors of the particles and the planets. The secular Hamiltonian is obtained taking the time average of this Hamiltonian over the various orbital motions. Let us consider that the particle is trapped in (p + q : p) MMR with β Pic c and not in any MMR with β Pic b (unless a very high-order weak one, which we neglect). Hence the interaction Hamiltonian with β Pic b H_b can be averaged over the orbital motions of the particle and β Pic b independently. This reads $$H_b=-\frac{G{m}_b}{4{\pi }^2}{\displaystyle \oint {\displaystyle \oint \left(\frac{1}{\left|r-r_b\right|}-\frac{r\cdot r_b}{r_b^3}\right)}}\text{d}l \text{d}{l}_b\text{, }$$(10)

where l and l_b are the mean anomalies of the particle and β Pic b, respectively. To the lowest (quadrupolar) order in semi-major axis a/a_b ratio (we implicitly assume a ⩽ a_b), and assuming full coplanarity, this reads in closed form $$H_{\text{sec},b,0}=-\frac{G{m}_b}{4}\frac{a^2}{a_b^3}\frac{1+(3/2)e^2}{{\left(1-e_b^2\right)}^{3/2}}$$(11)

where e and e_b are the eccentricities of the two bodies. As explained in Beust (2016) and Batygin & Brown (2016), this expression is only valid in the case of a fixed orbit for β Pic b. Considering β Pic b’s orbital precession under perturbations by β Pic c, the full expression is $$H_{\text{sec},b}=H_{\text{sec},\text{b},0}+\sqrt{aGM}\left(\sqrt{1-e^2}-1\right){\dot{\varpi }}_b,$$(12)

where ${\dot{\varpi }}_b$ is the periastron precession velocity of β Pic b.

A similar treatment must be done with the interaction H_c with β Pic c, but an independent averaging over the orbital motions of the particle and of β Pic c cannot be done because of the MMR configuration. A canonical transformation must first be done involving the critical argument of the resonance $$\sigma =\frac{p+q}{p}{\lambda }_c-\frac{p}{q}\lambda -\varpi .$$(13)

Here again, the precession of β Pic c’s longitude of periastron introduces additional terms in the canonical transformation (see Beust 2016, for details). The Hamiltonian then becomes $$H_{\text{res}}=H_c+{\check{\varpi }}_c\sqrt{aGM}\left(\frac{p+q}{p}-\sqrt{1-e^2}\right)-\frac{p+q}{p}{n}_c\sqrt{aGM},$$(14)

with n_c being β Pic c’s mean motion. The final averaging of H_c is then done over λ_c only at constant σ. As we are considering orbits with rather small semi-major axes, we must also add, following Pichierri et al. (2017), the first order post-Newtonian correction due to General Relativity: $$H_{\text{GR}}=\frac{1}{c^2}{\left(\frac{GM}{a}\right)}^2\left(\frac{15}{8}-\frac{3}{\sqrt{1-e^2}}\right).$$(15)

It should be noted that that expression is already averaged over the particle’s orbital motion. Following again Pichierri et al. (2017), we neglect the extra contribution due to the rotational bulge of the star. The resulting secular Hamiltonian is therefore $$H_{\text{sec}}=-\frac{GM}{2a}+H_{\text{sec},b}+H_{\text{GR}}+{\displaystyle \oint H_{\text{res}}}\text{d}{\lambda }_c.$$(16)

Now, to first order in perturbations between β Pic c and β Pic b, we derive $${\dot{\varpi }}_b=\frac{3n_b}{8}\frac{a_c^2}{a_b^2}\frac{m_c}{M+m_b}\frac{3e_c^2+2}{{\left(1-e_b^2\right)}^2},$$(17) $${\dot{\varpi }}_c=\frac{3n_c}{4}\frac{a_c^3}{a_b^3}\frac{m_b}{M+m_b}\frac{\sqrt{1-e_c^2}}{{\left(1-e_b^2\right)}^{3/2}}.$$(18)

Assuming coplanarity, the first three terms of H_sec alone would result in a secularly constant semi-major axis and constitute a one degree of freedom Hamiltonian. However the presence of the fourth term does not easily reduce to one degree of freedom, as the averaging process is done for constant σ and not constant a. As detailed by Morbidelli & Moons (1993) and Pichierri et al. (2017), the secular evolution is characterized by another adiabatic invariant J corresponding roughly to the amplitude of the σ-libration. As pointed out by Beust & Morbidelli (1996) and Pichierri et al. (2017), the study is considerably simplified if we consider orbits with negligible σ-libration amplitude. In this case, the semi-major axis remains secularly constant, and the Hamiltonian reduces to one degree of freedom. Portrait phases in $\left(\varpi -{\varpi }_c,e\right)$ space can be drawn, as was done in Beust (2016), for instance.

Figure 8 shows two such portrait phases for particles trapped in 7:1 and 8:1 MMRs with β Pic c, as they are the main contributors of FEBs in our simulations. The maps were made assuming solution #2 of Table 1 (non-resonant), but assuming the other solution does not reveal significant changes. Basically, these maps appear more complex than those corresponding to the 4:1 and 3:1 MMRs presented in Beust & Morbidelli (1996). This is due to the high order or the MMRs under consideration here, and to the significant eccentricity of β Pic c. We note that the maps were computed numerically with no specific assumption about that eccentricity. They may actually be compared to similar maps in Beust (2016). Nonetheless, particles starting at low eccentricity have many possible routes to reach high eccentricities and become FEBs following the level curves drawn here. Hence, the quoted MMRs are plausible theoretical FEB generators. Moreover, the eccentricity increase occurs near specific values of $\varpi -{\varpi }_c,$, suggesting a non-axisymmetric infall of FEBs as depicted above. However, the relative complexity of the phase portraits shown here, combined with the fact that many MMRs may be simultaneous FEB contributors, and taking into account the intrinsic chaotic nature of the real resonant motion suggests that we should return to the numerical study.

Fig. 8 Phase portraits of the full secular Hamiltonian (16) in $\left(\varpi -{\varpi }_c,e\right)$ space for particles trapped in 7:1 (left) and 8:1 (right) MMRs with β Pic c having negligible resonant libration amplitude, taking into account perturbations by β Pic b and General Relativity (see text). The horizontal bars denote the location of orbits with specific periastron values, down to the stellar surface (1 R*).

Fig. 9 Statistical description of the simulated FEB spectral events generated from the simulation described in Fig. 4, and corresponding to solution #2 from Table 1. Left: Doppler velocities of the FEB events occurring in the first 105 yr after the end of the simulation. The velocity scale is on the left side of the plot; positive velocities correspond to redshifts. The simultaneous evolution of β Pic c’s argument of periastron ωc,sky (see text) is superimposed in blue, with the corresponding scale in blue on the right side of the plot. The red bars highlight epochs where ωc,sky matches the fitted value from Table 1. Right: histogram of FEB velocities occurring ±100 yr around epochs corresponding to the red bars.

4.2 Numerical exploration

In this section we return to the simulations of the inner ring of particles depicted in Sect. 3.2. Figures 4 and 5 show that high-order MMRs with β Pic c may constitute active sources of FEBs, and the above semi-analytical study gives it a theoretical background, suggesting in addition that the FEB infall may be non-axisymmetric. We therefore investigate the statistics of FEB Doppler velocities that should be expected from those simulations; however, this is not as straightforward as it could first seem. First, the simulations were run in the invariable plane of the two-planet system, and did not assume any specific direction for the line of sight. Second, as the simulations extended over 20 Myr, the orbital elements of all particles were stored in the output every 10 000 yr to save disk space. This sampling was sufficient to reveal the reality of the mechanism, as can be seen from Figs. 4 and 5, but not to build the relevant statistics of all periastron passages of all particles. Hence, we decided to let the simulations run 200 000 yr after their end at 20 Myr, but we began storing orbital elements every 100 yr and interpolating all periastron passages occurring in between.

Next, we needed to fix a direction for the line of sight to be able to compute FEB velocities at the time they cross it. We first decided to assume for simplicity that the line of sight lies in the orbital plane of β Pic c (i.e. assimilating its orbital inclination to 90° with respect to the sky plane). As can be seen from Table 1 and Lacour et al. (2021), the actual value is close to 89°, so that the error on the projected velocities should not be significant.

Then we decided to fix the line of sight so as to coincide with the OX-axis of our referential frame. We note that this arbitrary choice does not mean assuming anything on the orbital configuration of β Pic c relative to the line of sight, as its orbit is expect to precess under the perturbations of β Pic b. In this context, and remembering that the orbital conventions of Lacour et al. (2021) and Blunt et al. (2020) assume an OZ-axis pointing away from the line of sight, it is straightforward to derive the following relation between β Pic c’s longitude of periastron ${\varpi }_c,$ as computed from our simulations, and its argument or periastron ω_c,sky following the quoted conventions: $${\omega }_{c,\text{ sky }}={\varpi }_c-\frac{\pi }{2}.$$(19)

Lagrange et al. (2020) assume an opposite convention with an OZ-axis pointing towards the observer. This results in a 180° shift in the value of the argument of periastron ω, as can be seen by comparing the orbital fits between Lagrange et al. (2020) and Lacour et al. (2021).

Figures 9 and 10 show the result of this computation for the simulations described in Figs. 4 and 5, corresponding respectively to solution #2 (non-resonant) and #1 (partly resonant) from Table 1. In each case the expected velocities of all FEB events simulated are plotted as a function of time in the first 10⁵ yr of the extended simulation. In this computation, however, we only took the still-not-fully-evaporated planetesimals into account (black dots in Fig. 6). We also kept computing the evaporation of planetesimals during the extended simulation.

The velocity statistics at any epoch can therefore be viewed directly from the plots. A clear temporal modulation over a few 10⁴ yr appears, which exactly corresponds the precession of β Pic c’s periastron. To highlight this, the temporal evolution of ω_c,sky is superimposed in blue to the plots. Then, red bars corresponding to epochs when ω_c,sky matches the fitted values from Table 1 are overplotted. Finally, the right plots of both figures display the extracted histograms of Doppler velocities of FEB events occurring ±100 yr around epochs, corresponding to the red bars in the left plots. These histograms simulate the expected velocity statistics of FEBs events for the present-day configuration of β Pic c with respect to the line of sight.

The first thing we note is a global agreement with the observed statistics of FEB velocities (i.e. a predominance of red-shifts, although blueshifts are present). The histograms in Figs. 9 and 10 actually compare very well with the observational statistics of Kiefer et al. (2014) and Tobin et al. (2019). The only noticeable difference is perhaps the absence of simulated high-velocity events (≳ ±100 km s⁻¹) compared to the observations. We point out here that the number of such events is probably highly underestimated in our simulation. Higher-velocity events correspond to FEBs crossing the line of sight very close to the star (≲0.1 au). Our evaporation rules (Sect. 3.3) are based on the simulations of Karmann et al. (2003) for average silicate material. However, the same work showed that, depending on the chemical composition of the body, the sublimation occurs at different distances. In particular, FEBs made of graphite should be able to resist as close as 0.1 au. Fernández et al. (2006) and Brandeker (2011) indeed showed that carbonaceous FEBs should be more able to retain their evaporated metallic compounds and generate observable signatures in absorption than others.

These trends are clear in both solutions, but the number of FEB events is obviously larger in the non-resonant configuration. The reason is that in the resonant case more FEBs miss the line of sight, due to larger chaos in their inclination evolution. Beust & valiron (2007) actually showed that the resonant FEB evolution is characterized by large inclination oscillations in the high-eccentricity phase. Here the fact that the FEBs are sometimes locked in MMR configuration with both planets causes larger chaos in that evolution, resulting in fewer events.

The statistics of FEB velocities is expected to change with time as β Pic c’s periastron precesses. The present-day configuration favours redshifts, but blueshifts should be dominant at other epochs. There are also epochs where FEBs events are much less numerous. This corresponds to configurations where the perias-tron of β Pic c is oriented in such a way that most FEB periastron passages occur behind the star as seen from the Earth.

It should also be noted here that all the cycles over 10⁵ in Figs. 9 and 10 are not identical. This highlights the discrete role of β Pic b. Although the periodic configuration of β Pic c relative to the line of sight dominates the process, the configuration relative to β Pic b does not evolve with the same periodicity, resulting in small changes in the successive cycles.

Fig. 10 Same as Fig. 9, but for the simulation described in Fig. 5, and corresponding to solution #1 from Table 1.

5 Discussion

The presence of both β Pic b and β Pic c orbiting β Pic drives us to a totally renewed view of the FEB generation mechanism. All previous studies pointed out the role of a distant planet that closely matches β Pic b as being responsible for the whole process via its low-order inner MMRs. The unexpected presence of β Pic c changes this picture. MMRs with β Pic b can no longer be sources of FEBs as the corresponding regions are dynamically unstable. Conversely, much higher-order MMRs between ~0.6 au and ~ 1.5 au with β Pic c now appear to be good candidates.

As in the previous model involving one perturbing planet, the FEB infall towards the star is not axisymmetric. The result is asymmetric statistics of FEB velocities at the time they cross the line of sight, favouring redshifs or blueshifts depending on the configuration of β Pic c’s periastron with respect to the line of sight. Our simulations show that the present-day fitted configuration actually favours redshifts, in agreement with observations. The simulated histograms in Figs. 9 and 10 confirm this general trend. Several MMRs actually contribute to the FEB population, and many dynamical routes can lead to a FEB state inside individual MMRs (see Fig. 8). This could provide a dynamical origin for the different families of FEB events quoted by Kiefer et al. (2014), thanks to a careful analysis of the events statistics.

The fact that planetesimals may spend a significant time in the disk before being captured into a MMR and involved in the FEB generation mechanism explains the duration of the mechanism. Thébault & Beust (2001) showed in the previous model with only β Pic b, that the 4:1 and 3:1 MMRs with β Pic b quite rapidly cleared out, causing the FEB process to nearly stop after 1–2 Myr. Collisions among planetesimals next to those MMRs had been invoked as a way to replenish the MMRs and sustain the process. Here thanks to its distant perturbations added to the global action of β Pic c, β Pic b helps the resonant FEB generation mechanism to still be active at the present age of β Pic. This activity is not expected to last for ever, however. At some point the whole planetesimal disk further out than ~0.6 au will be cleared out, and the FEB activity will cease. But the point here is that this is expected to happen later than the present age of β Pic. Letting the simulation run (far) beyond 20 Myr should help in specifying this issue. The perturbing action of β Pic b also influences the observability of the FEBs in absorption. β Pic b causes β Pic c’s periastron to secularly precess, and subsequently the statistics of FEB velocities to periodically change. The present-day configuration of β Pic c’s periastron with respect to the line of sight is compatible with redshifts, but this is about to evolve with time.

The new model also questions the physical nature of the FEBs themselves. In the preceding picture, FEBs were supposed to originate from regions located between ~3.5 and ~4.5 au (4:1 and 3:1 MMRs with β Pic b). Depending on the exact location of the ice line in the β Pic system, those FEB were likely to be at least partly icy. Ices were invoked in previous studies (Beust et al. 1990; Karmann et al. 2003) as a necessary source of volatile material in the FEB comas that prevents the metallic ions from being immediately blown away by the stellar radiation pressure, letting them first expand radially around the nucleus. This is required because it allows the ion clouds to cover a significant part of the stellar surface, and subsequently renders the FEBs spectral components visible. Here our revised model puts the FEB reservoir between ~0.6 au and ~ 1.5 au (i.e. much closer to the star). At that distance, the planetesimals are not expected to be icy at all. Ices may not be necessary, however. In early studies (Beust et al. 1990, 1996), volatile compounds were introduced as a braking agent for metallic ions as they are not affected by any noticeable radiation pressure from the star, contrary to metallic ions. Collisional activity between various species is sufficient to render the gas self-braking where the density is high enough (i.e. close to the nucleus). Fernández et al. (2006) showed in a detailed analysis that if the gas surrounding the FEB nucleus (arising from dust grain sublimation) contains enough carbonaceous compounds, the gas can be self-braking in the same way as before, but with no volatiles. The presence of solid material orbiting inside β Pic c’s orbit is also supported by recent independent observations. In a careful analysis of β Pic’s infrared spectrum and its associated silicate features, Lu et al. (2022) note the presence of a weak excess at 5µm tentatively attributed to a hot dust population at ~600 K, located within 0.7 au from the star. This hot dust could obviously be related to our FEB reservoir. Similarly, the stable gas component observed at rest with respect to the star could also originate from this reservoir. As shown by Lagrange et al. (1998) in an early study, this component is likely to arise from the FEBs themselves as the corresponding gas amount needs to be continuously replenished.

Another open issue is the survival of this planetesimal population to collisional activity. This is particularly relevant here as the mean eccentricities generated by the planetary perturbations on this population are large (see Fig. 6). Actually, the main question is not whether the population survives collisions, but how long the population should survive. The only requirement is that it should be able to survive until the current age of the star. Investigating this issue is beyond the scope of the present paper, and will require a dedicated study. There are nevertheless a few clues indicating that its survival could be possible. First, as noted above, there is observational evidence for an existing dust population inner to 1 au (Lu et al. 2022); second, the total population required to generate the observed FEB rate might not be that high as (i) given the large number of MMRs that are active source of FEBs, up to ~10% of the total population may be affected by the resonant process and (ii) thanks to the eccentricities reached (up to ~0.3), even non-resonant FEBs can participate to some extent in the FEB activity; and third, at the same time the FEB progenitor population reaches a significant eccentricity dispersion, a similar spread happens to the inclinations. This involves in particular resonant FEBs, which can reach inclinations up to several tens of degrees. This behaviour was already noted in Beust & valiron (2007), and can be explained as a combination of Kozai resonance and MMRs. Nonetheless, if a given planetes-imal population is able to reach a large inclination dispersion, its collisional activity should be reduced compared to the same population confined to a flattened disk. As a final argument, it can be seen that collisions among planetesimals also naturally contribute to replenishing the MMRs from regions with adjacent semi-major axes, and thus sustain the FEB activity.

Finally, the last open issue is the presence of additional planets in the β Pic system that could potentially affect the secular dynamics of the already-known planets, and subsequently the dynamics of FEBs. Our feeling is that the present model should be fairly robust to the presence of other planets. The simulations presented here show that there is not much space for additional planets closer to ~20 au as these would presumably be dynamically unstable. Planets located further out could conversely be present in the β Pic disk. Lagrange et al. (2020) indeed give observational upper limits showing that planets up to ~ 1 M_Jup in the 20–40 au region could still remain undetected. However, an additional planet, presumably less massive than β Pic b and β Pic c, and located much further out, is likely to have very little impact on the dynamics of planetesimals located around 1 au. Even β Pic c s periastron precession velocity should be only marginally affected. Hence, the presence of additional planets in the β Pic system is not expected to significantly affect the present model.

6 Conclusions

We presented a renewed model for the dynamical origin of the FEBs in the β Pic system. In the presence of both β Pic b and β Pic c, MMRs with β Pic b located at 3–5 au can no longer be a valuable source of FEBs as this region is intrinsically dynamically unstable. Higher-order MMRs with β Pic c located much closer to the star constitute conversely a new potential source, their efficiency being furthermore enhanced by the more distant action of β Pic b. In this context, a significant portion of the planetesimal disk between 0.6 and 1.5 au is potentially capable of becoming FEBs, in contrast to the previous model where only bodies initially trapped in one or two MMRs under consideration were able to do so. This causes the FEB activity in the β Pic system to potentially last much longer than initially estimated, and to remain observable today at β Pic’s present age. Moreover, the present-day configuration of the two orbits relative to the line of sight is able to reproduce the observed statistics of FEB velocities that shows a clear trend favouring redshifts, and thus reinforcing the plausibility of the scenario. The main parameter controlling this process is actually the orientation of β Pic c s periastron relative to the line of sight, and the current fitted orientation is compatible with a predominance of redshifts.

Many pending issues are still present in the new model and will require more dedicated studies. The statistics of FEB velocities needs to be more precisely simulated in the framework of the new model to more specifically address the issue of differentiated FEB families. This will be the purpose of future modelling work involving simulations with many more particles. Moreover, FEB progenitors arising from regions around 1 au are presumably not icy. As quoted above this is not a huge difficulty as carbon is likely to play a braking role similar to that played by volatiles (previously invoked) for metallic ions in FEB comas. However, the exact physics involved here needs to be modelled in the same way as was done with volatiles in early studies (Beust et al. 1990, 1996).

Finally, an observational test could be proposed. If the FEBs are actually not icy, it should be possible in principle to distinguish them from icy FEBs via specific spectroscopic studies. Interferometric observations have the capability to resolve the inner regions at sub-au scales where the dust generated by the FEB lies (Defrère et al. 2012) to shed light on its composition and water content, for instance using the low- to medium-resolution spectroscopic capability at L- and M-band (3–5 µm) of the VLTI/MATISSE instrument (e.g. Kokoulina et al. 2021).

Acknowledgements

All (or most of) the 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 work used astrometric measurements collected at the European Southern Observatory under ESO large programme ExoGRAVITY (ID 1104.C-0651). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COBREX; grant agreement n°885593). S.L. acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-21-CE31-0017 (project ExoVLTI).

References

All Tables

All Figures

	Fig. 1 Secular evolution over 108 yr of the β Pic b+c two-planet system corresponding to solution #2 from Table 1. From left to right: semi-major axes, eccentricities, and mutual inclination. The logarithmic time axes highlight the short-term as well as the long-term behaviour. Red always stands for β Pic c and green for β Pic b.
In the text

Fig. 2 Secular evolution over 108 yr of the β Pic b+c two-planet system corresponding to solution #1 from Table 1. Upper row, from left to right: eccentricities, mutual inclination, and (over the first 50 000 yr) the characteristic resonant argument σ7:1 (see text). Our definition of σ7:1 implies a 60° degeneracy. In our calculation, we always keep the values between 0 and 60°. Hence the variations from 0 to 60° appearing here actually mean circulation of σ7:1. Lower row: short-term semi-major axis variations of the two planets over the first 50 000 yr.

In the text

	Fig. 3 Outcome after a 10 Myr simulation of a disk of planetesimals initially between 0.5 and 45 au in the environment of β Pic b and β Pic c taken from solution #2 from Table 1, viewed from above (left), and displayed in (a, e) space (right). The red bullets show the location of the planets and the star, and the black dots are the remaining planetesimals.
In the text

	Fig. 4 Outcome in (a, e) space of a disk of planetesimals initially between 0.5 and 1.5 au, under the same conditions as in Fig. 3, after 2 Myr (left) and 20 Myr (right). Both planets are outside the panel.
In the text

	Fig. 5 Same as Fig. 4, with the same initial disk of particles between 0.5 and 1.5 au, but starting from solution #1 from Table 1 for β Pic b and β Pic c.
In the text

	Fig. 6 Same as Figs. 4 and 5 at t = 20 Myr, but taking into account gradual evaporation of FEBs, with solution #2 (left plot) and #1 (right plot) from Table 1. All particles appearing in blue are virtually already evaporated and should not account for the FEB statistics at that time. The locations of the main MMRs with β Pic c are indicated with red vertical bars.
In the text

	Fig. 7 Temporal evolution of the periastron value of two individual bodies out of the simulations from Figs. 4 and 5. The colour convention follows that of Fig. 6: the evolution appears in blue when the corresponding body is already evaporated.
In the text

Fig. 8 Phase portraits of the full secular Hamiltonian (16) in $\left(\varpi -{\varpi }_c,e\right)$ space for particles trapped in 7:1 (left) and 8:1 (right) MMRs with β Pic c having negligible resonant libration amplitude, taking into account perturbations by β Pic b and General Relativity (see text). The horizontal bars denote the location of orbits with specific periastron values, down to the stellar surface (1 R*).

In the text

Fig. 9 Statistical description of the simulated FEB spectral events generated from the simulation described in Fig. 4, and corresponding to solution #2 from Table 1. Left: Doppler velocities of the FEB events occurring in the first 105 yr after the end of the simulation. The velocity scale is on the left side of the plot; positive velocities correspond to redshifts. The simultaneous evolution of β Pic c’s argument of periastron ωc,sky (see text) is superimposed in blue, with the corresponding scale in blue on the right side of the plot. The red bars highlight epochs where ωc,sky matches the fitted value from Table 1. Right: histogram of FEB velocities occurring ±100 yr around epochs corresponding to the red bars.

In the text

	Fig. 10 Same as Fig. 9, but for the simulation described in Fig. 5, and corresponding to solution #1 from Table 1.
In the text