© The Authors 2025

1 Introduction

The young, bright, and nearby β Pictoris system presents an iconic planet formation laboratory. Two massive gas giants (Lagrange et al. 2009, 2019) orbiting at intermediately close distances of 2.7 and 9.9 AU (Lacour et al. 2021) create a dynamically multifaceted system, with a structured outer debris disk. The disk was immediately hypothesized to be associated with planet formation upon its discovery (Aumann 1985; Smith & Terrile 1984), long before the eventual detections of β Pictoris b and c, and has been a target of intense study ever since. Observations of the infrared excess emission of the disk have shown that multiple dust components exist, both far away from the star as well as within the inner few AU (Okamoto et al. 2004; Chen et al. 2007; Lu et al. 2022). Recent observations by JWST/MIRI MRS have shown that the dust emission is variable over decadetimescales, providing evidence for large planetesimal collisions, including those occurring close to the star (Chen et al. 2024). The star is also surrounded by gas-phase atoms and ions at various distances, visible in numerous absorption lines of calcium, sodium, iron, and other metals (e.g., Vidal-Madjar et al. 1986; Hobbs et al. 1988).

An emblematic property of β Pictoris is the existence of planetesimals on highly eccentric orbits that violently sublimate on close approach to the star. These objects, classically referred to as falling evaporating bodies (FEBs) but now commonly referred to as exocomets, were first observed in the mid-1980s owing to their strong line-absorption near the wavelengths of the Ca II H&K lines (e.g., Ferlet et al. 1987) as well as other metal ion lines (e.g., Lagrange et al. 1987), which vary on short timescales and occur at large Doppler shifts. Throughout the 1990s, many targeted spectroscopic observations and theoretical studies were carried out to further develop and cement the exocomet model (see Beust et al. 2024, for a review). Final confirmation of the exocomet hypothesis has recently come in the form of detection of exocomet transits in broadband photometry from the Transiting Exoplanet Survey Satellite (TESS) (Zieba et al. 2019; Pavlenko et al. 2022; Lecavelier des Etangs et al. 2022). These transit light curves have distinct shapes consistent with material arranged in a cometary tail, as originally predicted by Lecavelier Des Etangs et al. (1999).

The initial detections of massive close-in exoplanets using high-resolution spectroscopic monitoring of bright stars (Mayor & Queloz 1995) spurred the development of dedicated, highly stabilized, high-resolution spectrographs. Since its inception, the HARPS high-resolution spectrograph on ESO’s 3.6-m telescope in La Silla observatory (Pepe et al. 2002) has been used to systematically observe the β Pictoris system, both for the purpose of exocomet monitoring as well as later establishing the masses of β Pictoris b and c (e.g., Lagrange et al. 2012). By 2011, over 1000 HARPS spectra of β Pictoris were obtained, which were analyzed statistically for the first time by Kiefer et al. (2014). Assuming the obscuring cloud of a comet to be homogeneous, fitting for the line ratio of the Ca II H&K lines allows the optical depth as well as the fraction of the stellar disk that is covered by the cloud to be constrained (e.g., Lagrange-Henri et al. 1989). Kiefer et al. (2014) found that two populations of exocomet absorption features exist: the first, events with a broad distribution in radial velocities, relatively wide absorption lines, and small fill fractions; and second, lower redshift events with narrower absorption lines and larger fill fractions. This was interpreted as evidence that these populations tend to transit at different distances from the star, as larger clouds can be sustained at larger distances from the star where radial velocities are typically lower. Kiefer et al. (2014) also investigated the region of the neutral sodium (Na I) doublet near 589 nm and confirmed the existence of a static absorption line attributed to the circumstellar disk (Vidal-Madjar et al. 1986).

This paper is the first in a series in which we investigate exocomets in the β Pictoris system. Here, we introduce a new analysis methodology of the HARPS spectra and highlight a number of general aspects of the observed exocomets that we consider novel. We also focus on two accelerating exocomet features in 2019 and apparent exocomet activity in the 589 nm sodium doublet, while leaving a comprehensive statistical analysis of the large sample of observed Ca II features (similar to e.g., Kiefer et al. 2014) for a later study.

2 Methods

2.1 Data

The β Pictoris system has been observed extensively since 2003 with the HARPS high-resolution spectrograph on ESO’s 3.6m telescope in La Silla observatory. The public data archive lists 9133 spectra obtained until April 2020. We download these spectra¹, rejecting a small number of spectra that are incorrectly associated with β Pictoris or have excessive noise, leaving a total of 9071 spectra taken as part of 26 programs on 198 unique nights, with a total exposure time of 159.4 hours. At present, we do not aim to carry out a comprehensive statistical analysis, but we provide annual statistics of these observations in Table 1.

Previous studies (e.g., Kiefer et al. 2014) have used the blaze-corrected, and stitched 1D spectra produced by the HARPS Data Reduction Software (DRS). In this study, we follow our previous practice (e.g., Hoeijmakers et al. 2020; Prinoth et al. 2022; Hoeijmakers et al. 2024) by using the individual extracted spectral orders (e2ds files) instead, to avoid unnecessary reinterpolation that the pipeline carries out to produce data on a uniform wavelength grid and to average wavelength regions where orders overlap. Using the individual spectral orders is operationally nearly identical to using the 1D spectra, with the caveat that some wavelength ranges are covered by two orders (this is the case for the Ca II H&K lines), and the fact that the stellar continuum is modified by the blaze function – which does not affect our methodology (see Methods). As a result, we are able to directly use the count values in each of our spectral channels as representative of the flux that was recorded by the detector. Assuming that our spectra are in the photon-noise limit because β Pictoris is a very bright (V=3.9) star, the uncertainties are estimated as the square root of the recorded count levels – mitigating the fact that the reduced 1D spectra are not provided with estimates of the uncertainty at each wavelength.

Using the individually extracted orders (e2ds files) requires the application of the barycentric velocity correction for each spectrum in the dataset. For this, we use the value of the barycentric velocity provided by the pipeline. This leads to a unique wavelength axis for each spectrum in the dataset, even if the wavelength solution was constant in each night of data. Calculation of time-average spectra (e.g., the reference spectrum; see Section 2.3) still requires an interpolation onto a common wavelength grid, but for all other purposes in this paper, we retain a separate wavelength axis for each order of each of the spectra in the dataset.

For comparing spectra obtained at different times (e.g., to compute a reference spectrum; see Section 2.3), we need to account for possible variations in the broadband continuum, such as those due to atmospheric dispersion or flat-field changes. Although Kiefer et al. (2014) makes use of pipeline-calibrated spectra, such variations may have been largely removed. However, the e2ds spectra can attain end-to-end variations in the flux level of the continuum by 5% to 10%, which we correct following a methodology similar to our previous work (Hoeijmakers et al. 2020, e.g.,): We fit fourth-order polynomials to the residuals of the time-average spectrum using a least absolute deviation (LAD) optimization, while iteratively rejecting outliers and bad regions (telluric residuals, varying exocomet absorption lines).

2.2 Telluric correction

Besides Ca II, we intend to search for exocomet activity in the Na I doublet near 589 nm. This region is contaminated by telluric water lines, hindering the identification of exocomet sodium lines. We removed telluric contamination systematically from all spectra in the dataset using molecfit (Smette et al. 2015; Kausch et al. 2015). Because the line-spread function of HARPS varies from the short to the long wavelength side of each order (e.g., Zhao et al. 2021), and because we are interested only in correcting the water lines at the wavelengths of the sodium doublet, we restricted molecfit to ten water lines immediately adjacent to Na I lines. We avoided all water lines inside the broad photospheric Na I lines because strong exocomet sodium lines could bias telluric fitting. A typical correction is shown in Figure 1.

Fig. 1 Example of telluric correction near the neutral sodium doublet. The doublet from the circumstellar disk is readily visible.

2.3 Determining a reference spectrum

Similar to Kiefer et al. (2014), we proceeded to use the fact that exocomet absorption is highly variable in strength and wavelength to identify a reference spectrum of the system. The principle purpose of this was to correct for the absorption lines of the circumstellar gas, which were assumed to be static but partially blend with exocomet features, and to reliably identify exocomet absorption in the steep wings of the stellar absorption lines. To do this, we followed an approach similar to that of Kiefer et al. (2014), by estimating the mean flux at each wavelength channel using only their highest flux values, rejecting flux values that are potentially affected by stochastic exocomet absorption. The method adopted by Kiefer et al. (2014) assumes that all flux values in a spectral channel are drawn from the same random distribution. However, the large number of spectra available today span a wide range of signal-to-noise ratios (from a few dozen to over 200 in the center of the order covering Ca II K). Therefore, we adopted the following method that allows for the variable uncertainties between spectra:

Consider that the data consists of 9071 exposures 0...i...9071. Assume that in the absence of exocomet absorption, the flux X_ij (of exposure i in each wavelength channel j) was drawn from a normal distribution centered on a mean flux µ_j with a standard deviation σ_ij, equal to the measurement uncertainty (photon noise) of that particular exposure. Taking the mean of this collection of samples 〈X〉 would normally yield an accurate estimate for µ_j. However, due to the stochastic presence of exocomet absorption, an unknown number of these flux values will be depressed, biasing the estimate of µ_j toward lower values. Let us define a threshold, T_j, above which the flux samples in the wavelength channel are unlikely to be affected by exocomet absorption and are thus drawn from the unbiased distribution. This collection of samples forms a truncated normal distribution, and the true mean, µ_j, can be obtained from the mean of the truncated sample (Tallis 1961) as follows: $$⟨X_{ij}|X>T⟩={\mu }_j+{\sigma }_{ij}\frac{\varphi (\alpha )}{1-\mathrm{\Phi }(\alpha )},$$(1)

where φ and Φ, respectively, are the probability density function (PDF) and cumulative density function (CDF) of the normal distribution, and α is the standardized truncation threshold given by $$\alpha ={\alpha }_{ij}=\frac{T-{\mu }_j}{{\sigma }_{ij}},$$(2)

that is, the number of standard deviations by which the truncation threshold differs from the mean. But as µ_j is unknown in the presence of exocomets, equation 1 is to be solved iteratively: For each sample i, j, assume an estimate for µ_j to compute α. Subtract the term ${\sigma }_{ij}\frac{\varphi }{1-\mathrm{\Phi }}$ from each element in the truncated sample X_ij > T, and take the mean to obtain an estimate ofµ_j, and repeat until convergence.

At this stage, the computation of the reference spectrum still depends on the choice of threshold value T, which we take to be the 80^th percentile of the entire distribution of flux values X_ij. To arrive at this choice, we varied T between the 10^th to 98^th percentile and determined that the reference spectrum obtained by truncating the samples below the 80^th percentile lies within 2% of the reference spectrum that would be obtained when using the 98^th percentile, at all wavelengths beyond 15 km/s from the center of the disk line. This provides an accurate estimate of the reference spectrum based on a statistically large (∼2000) number of samples (see Fig. A.1).

However, similar to the reference spectrum obtained by Kiefer et al. (2014), the estimated reference spectrum is affected by noise. Because of this, we restricted the above formalism only to the 50% of spectra that have a higher than average signal-to-noise ratio (in the range of 100 to 200 measured in the center of order 6). In addition, we used the following Tikhonov regular-ization² to obtain a weakly smoothed approximation to the noisy reference spectrum by solving the system: $$\mu =(I-\lambda D^{T}D){\mu }_{\mathrm{s}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}},$$(3)

where I is the identity matrix and λD^T D is a regularization term, in which D is a sparse matrix populated by the finite difference approximation of the fourth-order derivative (1, −4, 6, −4, 1) on the diagonal and λ set to unity. The regularization term effectively penalizes roughness in the fourth derivative of the reference spectrum, acting to smooth variations in adjacent wavelength channels while preserving narrow spectral lines. We determined that the smoothed reference spectrum is within 0.5% of the un-smoothed reference spectrum at all wavelengths greater than 7 km/s away from the center of the disk lines. The reference spectrum thus obtained is shown in Fig. 2.

2.4 Exocomet line fitting

We followed Kiefer et al. (2014) in constructing an analytical model of the exocomet absorption lines. The absorbing cloud was approximated as a homogeneous medium with an optical depth at the line center of τ₀. The line shape was assumed to be Gaussian, making the relative absorption equal to $$A=1-e^{-{\tau }_0\mathrm{\Phi }},$$(4)

where Φ is the Gaussian line profile given as $$\mathrm{\Phi }=e^{-\frac{(\lambda -{\lambda }_0{)}^2}{2{\sigma }^2}},$$(5)

with λ₀ equal to the center wavelength of the line, and σ the line width.

For both calcium and sodium, we studied the line doublets with ratios in their oscillator strengths equal to 2. For small optical depth τ₀ < 1, the line ratio approaches 2. For large optical depths τ₀ ≫ 1, both lines saturate, resulting in equal line depths. However, if the cloud does not fully cover the stellar disk, the absorption lines do not extinguish the entirety of the stellar flux at the line center despite the large optical depth. This means that observations of line doublets are sensitive to the fill factor f of the stellar disk – a proxy for the radial extent of the cloud if assuming point-symmetry around the line of sight. We thus model the fraction of transmitted flux in one line i as follows, equivalent to Kiefer et al. (2014) and more recently Vrignaud et al. (2024) using a large collection of FeII lines simultaneously: $$F_i(\lambda )=1-f(1-e^{-{\tau }_0\mathrm{\Phi }}).$$(6)

The fractional line depth measured at the center of the line is $D=f(1-e^{-{\tau }_0})$.

Because we used non-normalized spectra, the continuum flux F_c (λ) is arbitrarily scaled and slowly varying with wavelength because our target lines are situated in the cores of the broad stellar lines, and the spectra extracted by the pipeline are not flux-calibrated. We therefore assume that the continuum is described by a low-order polynomial function F_c = c₀ + c₁ (λ -λ₀) + ... + c_n(λ - λ₀)ⁿ, where we typically set n = 2 or n = 3 to get a satisfactory fit without a dramatic increase in the number of free parameters. As there may be multiple (overlapping) exocomet features in any spectrum, we allow multiple components to be fit simultaneously using $$F(\lambda )=F_0(\lambda )F_1(\lambda )...F_N(\lambda )F_c(\lambda ).$$(7)

The doublet lines we targeted are widely separated in wavelength compared to the typical line widths and velocities involved. We wished to restrict the range of wavelengths over which we fit Eq. (7) for the continuum polynomial approximation to hold. Moreover, the Ca II H&K lines are covered by multiple spectral orders of HARPS simultaneously, with the K-line appearing in orders 6 and 7, and the H-line in orders 7 and 8. This means that for fitting the Ca II doublet, we included four narrow wavelength slices, each with its own continuum. The H line has a line strength that is half of that of the K line, so we fixed its optical depth accordingly. The same situation applies to the Na doublet, but these lines are covered by a single spectral order, meaning that only two wavelength slices needed to be included in a fit.

For any one spectrum, the model therefore has the following free parameters when fitting N line components:

• N center wavelengths λ₀ of each line component.

• N line widths σ.

• N center optical depths τ₀ of one of the doublet lines (the other being fixed via the known ratio of oscillator strengths).

• N surface fill fractions f.

• n + 1 polynomial components for each of the slices. For fitting the Ca II lines, we have four slices, giving 12 to 16 free parameters.

• A factor β that uniformly scales the uncertainties to capture unmodeled variance (e.g., small exocomet features that were not fit explicitly with a separate line component or telluric residuals). Using the same notation as in Section 2.3, this effectively modifies the log-likelihood of spectrum j as log $\log L_j\propto -\frac{1}{2}{\mathrm{\Sigma }}_i{\left(\frac{X_{ij}-F_{ij}}{\beta {\sigma }_{ij}}\right)}^2$. This use of β to scale uncertainties is common in exoplanet model fitting applications, including the fitting of transit light curves that are prone to correlated noise (e.g., Winn et al. 2008). A commonly used alternative method to capture unmodeled variance is by adding a constant term to the reported uncertainties, often referred to as a red noise or jitter term (see Pont et al. 2006, for an early application).

• Optionally a wavelength shift dλ between one slice and the others (i.e., the wavelength shift between the two lines of the doublet or a shift between the wavelength solutions of overlapping orders) to allow for any errors in the wavelength solution of the spectrograph. For the calcium lines, this gives two extra free parameters, for the sodium lines one.

Because this model is analytical, we implemented the model as an auto-differentiable function in Jax, and use NumPyro to sample from the prior distributions (see below) with a no U-turn sampler (Betancourt 2017; Bradbury et al. 2018; Bingham et al. 2018; Phan et al. 2019). An example of a best-fit model to a single exocomet feature with two components to describe the circumstellar disk line is shown in Fig. 3. This is a model with 27 free parameters.

A key choice for any spectrum is the wavelength range to fit over, the number of line components to include, and the definition of prior parameter distributions. There are nearly 10 000 spectra that we may wish to fit in this way, which would be prohibitively labor intensive. Using a relatively small subset of the present data, Kiefer et al. (2014) used an automated algorithm that chooses the number of exocomet components by evaluating and comparing the Bayesian information criterion. Our model would allow for a similar workflow, but in the present study, we target particular exocomet signatures in two nights of observation in 2019, and the sodium lines, which, as we show, typically have only one coherent absorption signature shifted away from the disk line (see Results). Secondly, it is worth noting that exocomet features typically do not vary strongly between the spectra taken in any one night of observation. In theory, this means that the definition of components and priors can, in most cases, be done on a night-by-night basis. The computational efficiency offered by Jax makes it feasible to execute model fitting on large numbers of spectra, and we aim to use this framework to carry out more rigorous statistical analysis of exocomet absorption signatures in these data in a future study.

A final choice is the assumption regarding the systemic velocity, which sets the zero-point for any measurements of radial velocities of exocomet spectral lines or disk lines. An independent measurement of the radial velocity of the β Pictoris host star is complicated by its high rotation velocity and δ Scuti pulsations, which cause time-variability in the photospheric line shapes. The systemic velocity is often measured using the narrow disk line components. For the systemic velocity, we adopt a value of 20.18 km/s as the average of FeI disk lines measured recently by Kiefer et al. (2019).

Fig. 2 Reference spectrum obtained following the methodology described in Section 2.3 (red), compared to the simple time-averaged spectrum (green) and the HARPS data (black). Smoothing using Tikhonov regularization reduces the noise in the reference spectrum at the 0.5% level (inset). As expected for an unbiased estimator, the reference spectrum tends to the time-averaged spectrum at wavelengths where no exocomet activity is present.

Fig. 3 Example of a fit to the Ca II doublet (K-line in the top panel and H-line in the bottom panel) of a blueshifted exocomet observed on December 10, 2019. The exocomet is modeled as a single component, and the disk line as a blend of two. The narrow fitting region is restricted to the blue side of the disk line to avoid additional unrelated exocomet components. Because the CaII lines are covered by three spectral orders of HARPS, the fit is carried out over four wavelength regions simultaneously. The posterior distribution of models is shown by the orange lines. Posterior distributions for model parameters are shown in Fig. 4.

3 Results

3.1 Rare exocometary Na I absorption

After telluric correction, we find a general absence of strong, variable Na I absorption similar to the strongly absorbing stochastic activity in the Ca II lines. Out of 198 nights of observation, we identify only two nights in which clear additional (>~2%) sodium absorption is present in both lines of the Na I doublet, in the nights of March 16, 2008, and March 22, 2009, strongly redshifted from the central disk line (see Fig. 5). In addition, we discern redshifted Na I absorption during an epoch spanning approximately 13 nights in February 2004, where it is present persistently with a depth of approximately 1% (see below).

In addition, we find strong, narrow Na I absorption at blueshifts between 12 and 28 km/s away from the NaI disk lines in approximately one third of the observed nights. These velocities coincide with the barycentric velocity at the time of observation (see Fig. A.2 for data obtained in 2017–2018), and we identify this as uncorrected telluric sodium absorption at the 1% level.

We used our fitting method (see Methods) to fit the two exocomet features of March 2008 and March 2009, as well as one epoch from the sequence of February 2004 in which the weak exocomet absorption is most clearly visible. The resulting fits are shown in Fig. 6, and compared with the Ca II K-line at the same epochs. In each fit, we include the Na I disk line as two blended components, a strong central line with a secondary red or blueshifted line (see Fig. 5). Because both the disk lines and the variable exocomet features are observed to have line depth ratios approximately equal to two, we infer the absorption to be unsaturated. At the same time, the lines are relatively weak, with line-depths on the order of 1% of the stellar flux. The optical depth of the absorbing cloud must thus be small (τ₀ ≪ 1). In this regime, f and τ₀ are strongly degenerate, which prevents us from determining confidently whether the exocometary Na I clouds fill the entire stellar disk. For this reason, we fix the fill fraction f to 100%, making τ₀ essentially a lower limit on the true optical depth, and a proxy for the line depth D (see above), which is tightly constrained in all three cases. In the 2004 sequence, the D1 line might actually be deeper than half of the depth of the D2 line (see Fig. 6), which would indicate mild saturation. Leaving the fill fraction f free leads the algorithm to prefer small fill fractions, but these are only inconsistent with 100% at a level between 3 and 4 σ, so we considered this constraint to be only tentative and continued to work with effective line depths in the remainder. We also systematically fit the blueshifted telluric features, where it is clearly visible by eye in 71 of the 198 nights of observation, and constrain the distribution of telluric sodium line depths (see Fig. A.3), demonstrating that the Earth’s atmosphere routinely contributes absorption at the 1–5% level.

All fits were carried out with the same priors (see Table A.1), apart from the prior on the center wavelength of the target component, which we shifted according to its approximate blue or redshift. An example of the posterior distributions of the model parameters of the exocometary Na I lines in the night of February 17, 2004, which is the weakest signal that we decided to fit. The line depth is constrained to 0.89 ± 0.13%, and thus the line is detected with high confidence. A comparison between the wavelengths of the Na I and the Ca II absorption lines reveals that exceptionally deep exocometary Ca II are present at the rare times during which significant Na I is detectable, at the same radial velocity (see Fig. 6). Although more detections of Na I would be needed to prove a more statistically robust correlation, we infer that the Na I we detect is sourced from the same exo-cometary bodies that are producing large clouds of Ca II. It is not unexpected that rocky exocometary material would source both calcium and sodium, but because of the short photo-ionization timescale of sodium, its column density may be expected to be insignificant. Indeed, most calcium absorption events are not associated with significant Na I absorption. The present observations suggest that Na I can still be detectable during episodes of strong Ca II production.

The observation of lines from multiple species can, in theory, be used to determine the elemental ratios of the evaporating exocometary material, and this has recently been convincingly demonstrated by Vrignaud et al. (2025). However, we expect that in the case of neutral sodium, the simple models used until now are unlikely to be able to provide accurate abundance measurements. In this work and in previous literature (e.g., Vrignaud et al. 2025), exocomet clouds have been modeled as obscuring the star homogeneously, while in reality they may exhibit radially varying density profiles, depending on their rates and balances of ionization and their dynamics. As Na I is very easy to ionize, we expect it to be present only close to the nucleus (potentially consistent with low values for f tentatively found earlier), while Ca II may be present much further away. In cases of very strong outgassing, there may be regions that are shielded from ionizing photons, and this may explain why Na I occurs simultaneously with very deep Ca II events. We conclude that determining the Na I/Ca II ratio may be possible, but that this would require a more detailed physical modeling of the particle dynamics, resulting morphology, and radiative properties of the exocometary cloud.

Fig. 4 Posterior distributions for five relevant model parameters of the fit shown in Fig. 3, out of a total of 27 free parameters (see Methods). We note that there is a degeneracy between the peak optical depth τ0 and the fill factor f, which is typical for fits of exocomet features that are not optically thick. However, in this example, f is sufficiently tightly constrained to conclude that the star is not fully covered by the cloud. Also note that the radial velocity of the cloud is constrained to within an uncertainty of about 100 m/s, and that the 1σ line width is approximately 5 km/s.

Fig. 5 Telluric-corrected HARPS observations at the sodium doublet, obtained on two nights nearly one year apart, vertically offset for clarity. On both nights, significant disk line absorption is visible, as well as redshifted lines attributed to a form of exocomet activity. Labels in black indicate the nights of observation, the number of spectra taken, and the total exposure time. The night-average spectrum prior to telluric correction is overplotted in orange, demonstrating the relative positions of the telluric water lines.

3.2 Two anomalous exocomet events

On two sequential nights of observation in 2019 (program ID 0104.C-0418, P.I. Lagrange, on December 11 and 12, 2019, hereafter nights 1 and 2), strongly blueshifted exocomet signatures were observed in the Ca II H&K lines, which experience significant acceleration during their transit (see Fig. 8), and occur at similar radial velocities. Assuming that their orbital periods are much greater than 24 hours (see below), we conclude that these are two objects on the same “comet train” on a similar orbit, transiting nearly 24 hours apart.

Observations of accelerating exocomets are relatively rare (Ferlet et al. 1987; Kennedy 2018), as these require observations spanning multiple hours and a comet that is very close to the star for its orbital curvature to cause significant acceleration. Accelerating blueshifted features are especially rare because exocomets are mostly observed while moving toward the star, away from the observer (Kiefer et al. 2014). blueshifted exocomets are significant because they are retreating from the star, implying that they have survived at least one close-in periastron passage, and because such a trajectory is only a rare outcome of the meanmotion resonance formation model (Beust et al. 2024), meaning that these comets may instead have a different dynamical origin (see Jaworska et al. in preparation for more discussion about the dynamics of blueshifted exocomets) from the population that is observed predominantly at redshifted velocities. Previous efforts to determine the statistics of orbital elements have been hampered by the fact that most exocomet absorption lines are stable on ∼1 hour time-scales, and a singular radial velocity measurement leaves strong degeneracies in possible orbital parameters. Observable acceleration of a transiting exocomet allows for some of this degeneracy to be lifted (Kennedy 2018).

We proceeded to use our line-fitting algorithm (see Methods) to determine the radial velocities and line widths of the accelerating features on both nights, with priors shown in Table A.2. Singular examples of these fits were already shown in Figures 4 and 3. On the second night of observation, a clear signature of a static or weakly accelerating blueshifted exocomet appears in the vicinity of −39 km/s, which partially blends with the stronger accelerating feature, overlapping with it toward the end of the observations (see Fig. 8). To fit the accelerating exocomet in the second night reliably, we therefore added a second exocomet component in the fit with priors included in Table A.2, with a typical result shown in Fig. 9. Best-fit parameters for both nights are plotted in Fig. 10, and we note that the uncertainty on fitting parameters increases at the end of the second night because of increased blending between the two components. Peak optical depths τ were typically found to be close to unity, with fill factors near 0.5 and consistent with full coverage of the stellar disk in some of the spectra, though typically with ∼10% uncertainty. These properties may be anomalous compared to the exocomet families distinguished by Kiefer et al. (2014), who showed that small line-widths and large fill factors are generally associated with exocomets with small positive radial velocities. The features detected here have relatively large fill factors, consistent with D-family comets, but also strong negative blueshift consistent with S-family comets.

On both nights of observation, the absorption lines accelerated linearly near –32 km/s for the first three hours and two hours, respectively. This is expected from a Keplerian orbit, because the radial velocity is, to first order, a narrow linear segment of a pseudo-sinusoid, covering a few hours of an orbit with a much longer period (Kennedy 2018). On the first night, we determined that this acceleration is nearly 1.1 km/s/h. To determine which class(es) of orbit are consistent with these velocities (initially chosen in preparation for defining Bayesian fitting priors), we computed the mean in-transit velocity of a 3D space varying semimajor axis a, eccentricity e, and argument of periastron ω, leaving the inclination i fixed to 0 (because this is a transiting geometry) and fixing the longitude of ascending node Ω to 180 degrees, as it is fully degenerate with ω in this geometry. We then determined which orbits have a mean radial velocity during transit between –27.5 and –35 km/s, while experiencing an acceleration between 0.9 and 1.2 km/s/h. The parameter combinations satisfying these criteria are shown in Fig. 11, revealing that a wide range of orbits with e > 0.3 could be consistent with the observed motion. Not shown in Fig. 11 is a matching class of orbits of arbitrary eccentricity, where the periastron is along the line-of-sight (ω = 90°), if only the first half of the transit (blueshifted) is covered by the observations. Also not shown are hyperbolic orbits. This serves to illustrate that robustly determining the orbital parameters of these objects is not possible without additional data or assumptions. A common assumption is to adopt a parabolic orbit, which removes two degrees of freedom from the problem. However, recent dynamical simulations by Beust et al. (2024) place the origin of exocomet progenitors within 1 AU from the star. Given that Ca II sublimation sets on within approximately 0.4 AU (Beust et al. 1998), assuming very high eccentricity trajectories may not be justified. In addition, we remark that such fitting procedures assume that the observed cloud of calcium follows the Keplerian velocity of the sourcing body, while in reality the ion tail should be expected to be blown in the anti-stellar direction, introducing an additional and unknown (range of) blueshift(s).

Instead of attempting to further constrain the orbital parameters of the transiting exocomet features in the two nights of 2019 data, we turn our attention to the clear departures from Keplerian velocities occurring later during the observations. In the fourth hour of observation on the first night, the acceleration of the comet increases, before decreasing one hour later. This behavior cannot be reconciled with Keplerian motion. More dramatic is the acceleration of the body in the second night—presumably on the same orbit as the first body if both are on the same “comet train” nearly 24 hours apart. After two hours into the observation, the acceleration of this object appears to reverse, before fully blending with the more static second component at the end of the sequence.

An apparent blue-ward acceleration of an exocomet may be expected if the observations cover the transit egress. In this scenario, the leading part of the cometary tail is removed from view, while the rest remains in-transit. Since material further along the tail is presumably more blueshifted due to radiation pressure than closer to the nucleus, this could cause an apparent blueshift of the disk-integrated absorption line. However, this scenario would predict that in the early part of the observations, the line would be wide, before narrowing as the least blueshifted components egress. In these observations (particularly in night one), the opposite happens. Secondly, the reversal of the acceleration in the second night is sudden, taking place in approximately an hour, while extended clouds have ingress and egress times of several hours. Finally, we note that the behavior of both objects is different. Assuming that these objects are indeed moving on the same orbit, their behavior during egress could be expected to be similar.

The presence of the blended component in night two complicates the fit, and this is evident from the increased uncertainties on the best-fit parameters. In theory, the presence of weaker, blended exocomet features on similar or different orbits could act to skew the centroid velocities, but we do not see clear evidence of additional absorbing features and note that blueshifted exocomets are altogether relatively rare (Kiefer et al. 2014).

Instead, we propose an alternative scenario to interpret the sudden departure of Keplerian motion seen in night two as providing evidence for the final fragmentation of the exocomets’ nucleus and subsequent blow-out of the cometary tail, briefly after periastron passage. Destruction of stargrazing comets may be sudden if the core of the nucleus contains unexposed, previously unevaporated volatiles (ice) that evaporate in response to the intense heating during periastron passage. This evaporation may be delayed until after periastron passage, as the heat pulse first needs to propagate through the rocky exterior before evaporating the icy interior. This scenario is comparable to Kreutz-family Sungrazers in the Solar system, for example, the post-periastron destruction of the nucleus of comet Lovejoy in 2011 (Sekanina & Chodas 2012). As this scenario may favor an icy interior composition, this further strengthens the possibility that the origin of these two exocomets is not among the more commonly observed redshifted exocomets presumed to be generated from the devolatilized inner planetesimal belt via mean-motion resonance with β Pictoris c but, instead, may have originated from the outer system exterior to the snow-line relatively recently. Rodet & Lai (2024) as well as our companion paper (Jarworska et al. in preparation) suggest that secular evolution of icy exocomet progenitors from the outer system is a viable pathway to producing observable exocomets. Our dynamical simulations (Jaworska et al., in prep.) indicate that exocomet progenitors likely spend less than 10² to 10³ years within the water sublimation limit before reaching within 0.4 AU of the star where they are spectroscopically observed as exocomets. Given a sufficiently large nucleus, an icy core may survive throughout the comets’ passage through the inner system before becoming an exocomet, and this has been explored previously by Karmann et al. (2001).

Fig. 6 Best-fit models to the NaI D lines (left and middle columns) on three nights of data (top, middle, and bottom rows), compared to the CaII K line (right column). The night of March 17, 2008, (middle row) shows the strongest absorption of exocometary sodium of all 198 nights of data. The night of February 17, 2004, is part of a sequence of 13 nights in which weak exocometary sodium absorption persists, and this is the night with the best detectability of this sequence but the weakest exocometary signal detected in the entire dataset. The posterior parameter distributions of this fit are shown in Fig. 7 as a representative example. In all three cases, the exocometary sodium absorption occurs simultaneously with strong Ca II absorption at the same relative radial velocity – indicating that the sodium is sublimated together with calcium by the same exocomets.

Fig. 7 Posterior distributions of the parameters that describe the exo-cometary sodium lines on the night of February 17, 2004, shown in Fig. 6. Because the observed Na I lines are highly optically thin, τ0 and f were not independently constrained. Instead, we fixed f to 1 and converted τ0 into line depth D, which is constrained to 0.89 ± 0.13%. Note that β may be elevated due to unmodeled variance in the spectrum caused by residuals from the telluric correction process.

Fig. 8 Spectroscopic time series of two nights of observations blueward of the Ca II K-line, calibrated using the reference spectrum. The observing sequences both lasted for about 6 hours and took place 24 hours apart. The disk line is near 0 km/s, and a strong blueshifted exocomet feature absorbs near –30 km/s on both nights, decelerating toward smaller blueshifts. In the second night, a third component appears statically near –40 km/s that blends with the decelerating feature. This component is indicated with the- dashed line. The line parameters of these exocomets are shown in Fig. 10.

Fig. 9 Example of a fit of the blended components during night two, indicated by the vertical blue markers. The fit confidently distinguishes the two components.

Fig. 10 Centroid velocities (top panel) and Gaussian line widths (bottom panel) of the exocomet features shown in Fig. 8. Both bodies follow a near-constant acceleration at the start of the observations, as expected from a Keplerian orbit (see Fig. 11). However, they then deviate from Keplerian acceleration toward the second half of the observing sequence. Simultaneously, the observed line width increases dramatically in the first night and may increase on the second night as well. In some cases, in particular toward the end of Night two when the two components blend, the fit was unstable, resulting in prior-bound parameter estimates. These points are omitted.

3.3 Coherence of exocomet lines and disk line variability

As previously shown for a subset of this data (Kiefer et al. 2014), exocometary Ca II absorption is extremely variable and occurs over a wide range of (mostly redshifted) radial velocities up to 100 km/s. Some exocomet absorption lines vary on timescales of hours, while some absorption signatures are coherently present for days. Variability on hour timescales is naturally expected from transiting objects on eccentric orbits at close distances to the star, while longer-term coherence supports the notion of “comet families”: trains of objects on similar orbits but separated from one another in time. This is naturally expected from a formation scenario in the mean-motion resonance model (Beust et al. 2024). However, we hypothesize that planetesimals from the inner system that turn into exocomets via mean-motion resonance with β Pictoris c may make close passes with the planet and be susceptible to tidal disruption. Tidal disruption may also explain the prevalence of long timescale coherence in the observed Ca II exocomet signatures of planetesimals originating from the outer system, as these typically undergo chaotic dynamical evolution after being scattered into the inner system.

However, the data also provide examples of much longer time-stability: From early 2017 to mid-2018, there is a ∼5 km/s blueshifted feature that appears to last continuously (see Fig. A.6). Dynamical simulations should be carried out to investigate whether such long coherence timescales are readily explained by the exocomet model, or whether this feature is due to some other component of the circumstellar material.

Finally, from investigating the large volume of observations at the Ca II lines by eye, and by analyzing the best-fit line parameters of the Na I D-lines, we conclude that the disk lines show significant variability that may hinder the usage of percentilebased reference spectra (as per Kiefer et al. 2014). Figure A.4 shows that the relative line-depth of the Na I disk lines varies by more than a factor of two from ∼5% up to 13%. The line center varies by a few hundred meters per second, with a mean measurement uncertainty of 50 m/s. Although the Ca II doublet appears to be more stable, it also shows variability. Fig. A.5 shows an example of a night of observation in 2014 where the observed disk line is significantly narrower and shifted compared to the reference spectrum (see Methods). Kiefer et al. (2019) also observed variability in the Fe I disk line, indicating that the tenuous circumstellar gas is generally observably variable.

Fig. 11 Family of orbital solutions with a mean radial velocity during the transit between –35 and –27.5 km/s and a mean radial acceleration between 0.9 and 1.3 km/s/h, consistent with the linear part of the exocomet radial velocity observed during the first three hours of observation on the first 2019 night, and assuming an inclination i equal to 0.0 and a fixed longitude of ascending node Ω (degenerate with argument of periastron ω when i = 0). A linearly accelerating absorption line still permits a wide combination of orbital eccentricities, semimajor axes, and orbital orientations.

4 Conclusions

In this paper, we investigated the growing archive of HARPS observations of the β Pictoris system in search for new exocomet phenomenology. This paper is the first in a series that seeks to investigate the nature of the exocomets of β Pictoris, leveraging the large number of high-quality spectra obtained for nearly two decades. We started by focusing on the region of the sodium D-line doublet, which overlaps with significant telluric water absorption. We systematically corrected these spectra using Molecfit and concluded that exocometary sodium absorption is relatively rare, with only three distinct epochs showing detectable signals. These events are contemporaneous with exceptionally strong Ca II absorption occurring at the same relative radial velocity, and we conclude that the observed sodium and calcium components are sourced by the same bodies. We also investigated two anomalous exocomet events at the Ca II lines observed in 2018, which show a strong and sudden onset of non-Keplerian acceleration that we hypothesize to be evidence for the destructive fragmentation of the exocomet nucleus and the dispersal of the Ca II tail. Lastly, after creating a qualitative overview of the phenomenology in the β Pictoris spectra but without doing a thorough statistical analysis, we conclude that there is significant disk line variability in Ca II and Na I, as well as long time-scale coherence in some presumed exocomet absorption signals that as yet remain unexplained.

Data availability

The telluric-corrected spectra are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/700/A239

Acknowledgements

This work was supported by grants from eSSENCE (grant number eSSENCE@LU 9:3), the Swedish National Research Council (project number 2023-05307), The Crafoord foundation and the Royal Physiographic Society of Lund, through The Fund of the Walter Gyllenberg Foundation. We are grateful to Brett Morris for useful discussions regarding the inference method used extensively in this work, and the attendees of the 2024 ISSI workshop on Exocomets in Bern, Switzerland, for important discussions regarding the telluric Na I lines. We are also grateful by the insightful comments and advice provided by an anonymous referee, who is to be credited for remarking the existence of the blended component in the second night of the 2019 observations.

Appendix A Additional tables and figures

Fig. A.1 Comparison of the reference spectrum (red) computed for varying values of the threshold T as the ninety-eighth, ninetieth, eightieth, fiftieth, twenty-fifth, or tenth percentile (red lines from top to bottom). The time-average spectrum is shown in green, corresponding to the zeroth percentile reference spectrum.

Fig. A.2 Nightly averages of telluric-corrected spectra at the sodium doublet, vertically offset for clarity. These spectra are obtained between January 2017 and June 2018 and show how a telluric NaI feature is continuously present with a radial velocity varying within the range of radial velocities of the Earth around the barycenter (shaded orange).

Fig. A.3 Best-fit depths (D) of the telluric D2 line, for all 71 nights in which telluric Na I was discerned.

Fig. A.4 Best-fit radial velocities (top panel), line widths (middle panel), and depths (bottom panel) of the main component of the D2 disk line.

Fig. A.5 One night of observation in which the Ca II disk line is significantly narrower and offset compared to the reference spectrum, indicating intrinsic variability in the disk line.

Fig. A.6 Spectra taken from early 2017 until mid-2018 compared to the reference spectrum. In May 2017 a blueshifted feature appeared that lasted for at least a year, indicated with the orange line at −5 km/s.

References

All Tables

All Figures

	Fig. 1 Example of telluric correction near the neutral sodium doublet. The doublet from the circumstellar disk is readily visible.
In the text

Fig. 2 Reference spectrum obtained following the methodology described in Section 2.3 (red), compared to the simple time-averaged spectrum (green) and the HARPS data (black). Smoothing using Tikhonov regularization reduces the noise in the reference spectrum at the 0.5% level (inset). As expected for an unbiased estimator, the reference spectrum tends to the time-averaged spectrum at wavelengths where no exocomet activity is present.

In the text

Fig. 3 Example of a fit to the Ca II doublet (K-line in the top panel and H-line in the bottom panel) of a blueshifted exocomet observed on December 10, 2019. The exocomet is modeled as a single component, and the disk line as a blend of two. The narrow fitting region is restricted to the blue side of the disk line to avoid additional unrelated exocomet components. Because the CaII lines are covered by three spectral orders of HARPS, the fit is carried out over four wavelength regions simultaneously. The posterior distribution of models is shown by the orange lines. Posterior distributions for model parameters are shown in Fig. 4.

In the text

Fig. 4 Posterior distributions for five relevant model parameters of the fit shown in Fig. 3, out of a total of 27 free parameters (see Methods). We note that there is a degeneracy between the peak optical depth τ0 and the fill factor f, which is typical for fits of exocomet features that are not optically thick. However, in this example, f is sufficiently tightly constrained to conclude that the star is not fully covered by the cloud. Also note that the radial velocity of the cloud is constrained to within an uncertainty of about 100 m/s, and that the 1σ line width is approximately 5 km/s.

In the text

Fig. 5 Telluric-corrected HARPS observations at the sodium doublet, obtained on two nights nearly one year apart, vertically offset for clarity. On both nights, significant disk line absorption is visible, as well as redshifted lines attributed to a form of exocomet activity. Labels in black indicate the nights of observation, the number of spectra taken, and the total exposure time. The night-average spectrum prior to telluric correction is overplotted in orange, demonstrating the relative positions of the telluric water lines.

In the text

Fig. 6 Best-fit models to the NaI D lines (left and middle columns) on three nights of data (top, middle, and bottom rows), compared to the CaII K line (right column). The night of March 17, 2008, (middle row) shows the strongest absorption of exocometary sodium of all 198 nights of data. The night of February 17, 2004, is part of a sequence of 13 nights in which weak exocometary sodium absorption persists, and this is the night with the best detectability of this sequence but the weakest exocometary signal detected in the entire dataset. The posterior parameter distributions of this fit are shown in Fig. 7 as a representative example. In all three cases, the exocometary sodium absorption occurs simultaneously with strong Ca II absorption at the same relative radial velocity – indicating that the sodium is sublimated together with calcium by the same exocomets.

In the text

Fig. 7 Posterior distributions of the parameters that describe the exo-cometary sodium lines on the night of February 17, 2004, shown in Fig. 6. Because the observed Na I lines are highly optically thin, τ0 and f were not independently constrained. Instead, we fixed f to 1 and converted τ0 into line depth D, which is constrained to 0.89 ± 0.13%. Note that β may be elevated due to unmodeled variance in the spectrum caused by residuals from the telluric correction process.

In the text

Fig. 8 Spectroscopic time series of two nights of observations blueward of the Ca II K-line, calibrated using the reference spectrum. The observing sequences both lasted for about 6 hours and took place 24 hours apart. The disk line is near 0 km/s, and a strong blueshifted exocomet feature absorbs near –30 km/s on both nights, decelerating toward smaller blueshifts. In the second night, a third component appears statically near –40 km/s that blends with the decelerating feature. This component is indicated with the- dashed line. The line parameters of these exocomets are shown in Fig. 10.

In the text

	Fig. 9 Example of a fit of the blended components during night two, indicated by the vertical blue markers. The fit confidently distinguishes the two components.
In the text

Fig. 10 Centroid velocities (top panel) and Gaussian line widths (bottom panel) of the exocomet features shown in Fig. 8. Both bodies follow a near-constant acceleration at the start of the observations, as expected from a Keplerian orbit (see Fig. 11). However, they then deviate from Keplerian acceleration toward the second half of the observing sequence. Simultaneously, the observed line width increases dramatically in the first night and may increase on the second night as well. In some cases, in particular toward the end of Night two when the two components blend, the fit was unstable, resulting in prior-bound parameter estimates. These points are omitted.

In the text

Fig. 11 Family of orbital solutions with a mean radial velocity during the transit between –35 and –27.5 km/s and a mean radial acceleration between 0.9 and 1.3 km/s/h, consistent with the linear part of the exocomet radial velocity observed during the first three hours of observation on the first 2019 night, and assuming an inclination i equal to 0.0 and a fixed longitude of ascending node Ω (degenerate with argument of periastron ω when i = 0). A linearly accelerating absorption line still permits a wide combination of orbital eccentricities, semimajor axes, and orbital orientations.

In the text

	Fig. A.1 Comparison of the reference spectrum (red) computed for varying values of the threshold T as the ninety-eighth, ninetieth, eightieth, fiftieth, twenty-fifth, or tenth percentile (red lines from top to bottom). The time-average spectrum is shown in green, corresponding to the zeroth percentile reference spectrum.
In the text

	Fig. A.2 Nightly averages of telluric-corrected spectra at the sodium doublet, vertically offset for clarity. These spectra are obtained between January 2017 and June 2018 and show how a telluric NaI feature is continuously present with a radial velocity varying within the range of radial velocities of the Earth around the barycenter (shaded orange).
In the text

	Fig. A.3 Best-fit depths (D) of the telluric D2 line, for all 71 nights in which telluric Na I was discerned.
In the text

	Fig. A.4 Best-fit radial velocities (top panel), line widths (middle panel), and depths (bottom panel) of the main component of the D2 disk line.
In the text

	Fig. A.5 One night of observation in which the Ca II disk line is significantly narrower and offset compared to the reference spectrum, indicating intrinsic variability in the disk line.
In the text

	Fig. A.6 Spectra taken from early 2017 until mid-2018 compared to the reference spectrum. In May 2017 a blueshifted feature appeared that lasted for at least a year, indicated with the orange line at −5 km/s.
In the text