© The Authors 2023

1 Introduction

Rocky and icy bodies in the size range of ~1 km to a few hundreds of kilometers are very abundant in our Solar System, and they are often concentrated in various locations such as the asteroid belt, the Kuiper belt, the Oort cloud, and the Lagrangian L4 and L5 points of the giant planets. They can be collectively referred to as “planetesimals”, and are thought to be a natural outcome of the planet formation process (e.g., Schlichting et al. 2013; Johansen et al. 2015). Due to their small sizes and low temperatures, individual planetesimals are typically incredibly difficult to detect, even in the Solar System (except in favorable circumstances), let alone in other stellar systems. However, planetesimals are occasionally disrupted, either through collisions with other planetesimals (e.g., Jackson et al. 2014), or through tidal disruption during close encounters with planets (e.g., Cataldi et al. 2018; Janson et al. 2020). The dust produced as part of the disruption has a much larger collective optical surface than the planetesimal itself, and it can therefore produce a visible signal (e.g., Lawler et al. 2015). The dust from an individual disruption disperses over a short timescale due to irradiation from the central star, but a rich planetesimal belt, where disruptions occur sufficiently frequently, can form a continuously visible structure known as a debris disk. Debris disks around other stars typically need to be substantially brighter than the Kuiper belt to be observable, but still, they are relatively common (e.g., Eiroa et al. 2013). For example, ~ 16–17% of Sunlike stars host confirmed debris disks, with a mildly increasing trend toward higher-mass stars (Trilling et al. 2008; Sibthorpe et al. 2018). Additionally, stacking thousands of phase-folded light curves of Kepler planet candidate systems has revealed putative flux deficits around the L4 and L5 points of the planets, which would imply that populations of trojan planetesimals could be common in planetary systems (Hippke & Angerhausen 2015). These observational results imply that a large number of planetesimals must be an essentially ubiquitous feature of stellar systems, as is expected in core accretion-based planet formation theories (Drazkowska et al. 2022). Gaining observational insight into the compositional properties of planetesimals will constrain the chemical makeup of rocky exoplanets in inner planetary systems (Krijt et al. 2022), as the timescale and physical mechanisms of planetesimal differentiation influence the atmospheric composition of mature exoplanets (Lichtenberg et al. 2022).

Through dynamical interactions, for example with planets in the same system (Duncan & Levison 1997), a fraction of the planetesimals can be scattered into highly eccentric orbits. If the periastron is close enough to the parent star (on the order of 1 au or closer), then for a fraction of the orbit, the planetesimal partially evaporates and forms a cloud of dust and gas around it – that is to say, a coma, making it identifiable as a comet. During this phase, the effective optical surface of the object increases drastically due to the coma, and from the heating of the star, the temperature also rises, making the comet many orders of magnitude brighter in both visible and infrared light than in its “normal” quiescent phase.

Even though active comets are much brighter than planetesimals of corresponding sizes, they are still faint relative to the parent stars. Furthermore, in order to receive enough stellar flux to be reasonably bright at visible and infrared wavelengths, potentially observable comets need to reside at a small separation from the parent star, placing them in the same extreme contrast regime as directly imaged planets (with the exception of young, self-luminous planets, e.g. Marois et al. 2010; Janson et al. 2021a; Bonati et al. 2019). Indeed, as discussed in Sect. 3, moderate-sized planets (approximately Earth- to Neptune-sized) overlap with the brightness range of large comets. Exocomets are therefore subject to the same type of high-contrast problem as exoplanets, where the vast majority of them are beyond reach for direct imaging with present-day observational facilities. Hence, for exocomets and exoplanets alike, the first pieces of evidence have come from indirect methods. In the case of exocomets, this evidence has arrived through transit measurements of comets (or material of cometary origin) passing in front of the stellar surface, as summarized in Sect. 3. However, as technology approaches the capability to directly image terrestrial planets, it simultaneously approaches the capacity to directly image exocomets.

The Large Interferometer for Exoplanets (LIFE, see Quanz et al. 2022) is a mission concept for a space-based infrared interferometer capable of nulling interferometry. The primary aim of the conceived mission is to image planets through infrared nulling interferometry, being optimized for sufficient contrasts to detect and characterize several tens of Earth-like planets in the classical habitable zones of their parent stars (e.g., Konrad et al. 2022; Alei et al. 2022). The resolution, sensitivity and contrast provided by such a mission would also open up a wide range of other scientific avenues that could not be addressed with any other existing or planned facility. In this paper, we study observational aspects of exocomets in the context of high-contrast imaging surveys, with a strong focus on LIFE since it is foreseen as being particularly potentially powerful for the detection of comets. There are two recent developments that make a concrete study to this end possible for the first time: firstly, a LIFE simulator has been produced and made available, such that realistic estimations can be performed for the mission performance and required integration times (LIFEsim¹; see Dannert et al. 2022). Secondly, a large number of exocomet transits have now been identified in the β Pictoris system (Zieba et al. 2019; Lecavelier des Etangs et al. 2022), making it possible for the first time to predict the frequency and observational properties of comets in this young debris disk system.

The paper is outlined as follows: in Sect. 2, we discuss the basic parameters assumed for LIFE when used in an exocometary context. In Sect. 3, we outline previous observational evidence for exocomets and what has been learned from them, and how this knowledge is used to predict observational properties for direct imaging purposes. A set of particularly promising potential targets for exocomet searches are described in Sect. 4. We then present the results from a case study for the β Pic system in Sect. 5, which sets the stage for a wider discussion of exocomet detection in Sect. 6. We finally summarize the study in Sect. 7

2 Mission parameters

This feasibility study concerns imaging of exocomets with LIFE. The design of LIFE is under active evaluation (e.g., Hansen et al. 2022), but here we will use the same reference case design of LIFE as in Quanz et al. (2022). The design involves four telescope units, which as a baseline in this study will be assumed to each be 2 m in diameter. The unit telescopes relay their received signals to a beam combiner unit for interferometric purposes. LIFE will offer a normal interferometric imaging mode, but it will also offer a nulling interferometry mode, which allows to effectively cancel out the light from a bright star through destructive interference, whilst conserving light signals arising off-axis from the star, such as from planets, or indeed, comets. The nulling interferometry mode offers the best contrast at small separations from bright stars, which means that it is the natural mode to use for exocomet observations. The telescope spacecraft are foreseen to move in a rotation-like manner to optimize UV-plane coverage.

The reference case design includes a low-resolution spectroscopic mode (R ~ 20) spanning a wavelength range of 4–18.5 μm. It therefore captures a large fraction of the thermal radiation emitted from a temperate object such as an exocomet visiting the inner parts of a stellar system (e.g., Shanklin 2000), and also incorporates potentially important spectral features such as the silicate feature around 10 μm. The interferometric baseline will be adjustable in any design, but is foreseen to mainly operate in the range of tens to hundreds of meters. The nulling mode of LIFE is foreseen to operate in a rectangular configuration (see Fig. 1), where the shorter axis is called the ‘nulling baseline’ with length b, and the longer axis is called the ‘imaging baseline’ with length qb where, following Dannert et al. (2022), q = 6. In this study we will simply assume a nulling baseline of b = 25 m (unless stated otherwise). This gives an effective size of a spatial resolution element of a few tens of mas within the LIFE operating wavelength range. The field of view (FOV) for interferometric observations is set by the beam size of the individual telescope units, λ/D. Thus, if D = 2 m and λ = 10 μm, the FOV is 1 arcsec wide. For S/N (signal-to-noise ratio) and integration time calculations, we use the LIFEsim software (Dannert et al. 2022), which simulates the nulling interferometry mode of LIFE to produce realistic performance estimates.

3 Observing exocomets

All observational evidence for exocomets that exist to date originate from indirect measurements. The first evidence arose from the detection of temporally varying narrow spectral features superimposed on the broad spectral lines of the primary star in the β Pic system (Ferlet et al. 1987; Lagrange-Henri et al. 1988). This was interpreted as being due to exocomets on eccentric orbits, passing so close to the star that they evaporate – so called Falling Evaporating Bodies (FEBs; see Beust et al. 1990). Subsequently, the same type of features have been identified in the 49 Cet (Montgomery & Welsh 2012) and HD 172555 (Kiefer et al. 2014a) systems. A list of a few tens of more tentative FEB detections can be found in e.g. Strøm et al. (2020). Further concrete evidence for comets around β Pic came with the discovery of a cometary transit in the system (Zieba et al. 2019). Similar transit events had previously been identified in the KIC 3542116 system (Rappaport et al. 2018). While the exact morphology of a cometary transit depends on orientation (Lecavelier Des Etangs et al. 1999), the shape can generally be described as asymmetric, usually with a relatively steep flux drop initially, followed by a smoother rise back toward the baseline. Such a transit is not likely to be mixed up with any alternative known mechanisms, as long as the S/N is reasonably high. Subsequently, analysis of TESS (Ricker et al. 2015) data has revealed several tens of additional cometary transits in the β Pic system, which has even made it possible to derive a size distribution for the underlying comet population (Lecavelier des Etangs et al. 2022).

In this paper, we focus on observability of comets through direct imaging, which has so far not been accomplished. In the following, we will discuss the general factors of importance for exocomet imaging, which will be applied particularly to the case of β Pic in Sect. 5.1.

A famous feature of comets in the Solar System is the long cometary tail, and for this reason it may seem natural to consider exocomets as having extended flux distributions. However, the tail is typically very faint relative to the head of the comet, and despite the fact that LIFE would have an excellent spatial resolution at infrared wavelengths compared to most existing facilities, the resolution is still very coarse relative to the sizes of bodies in a stellar/planetary system. The density of the tail downstream from the comet nucleus can be described as an exponential drop-off ρ = ce^−λx (adapted from Zieba et al. 2019), where c is a normalization constant and x is a spatial coordinate, such that 1/λ is the scale length of the tail. For the deepest and most well-modeled transiting exocomet, assuming a characteristic star-comet separation of 1 au, Zieba et al. (2019) find a scale length of ~2 × 10⁸ m, or approximately 10% of the stellar diameter². This in turn corresponds to ~0.07 mas at the distance of β Pic, and even for the closest possible system of α Cen, it is only 1 mas. This is more than an order of magnitude smaller than the spatial resolution for the reference case LIFE design, as outlined in Sect. 2. For all practical purposes in this study, we can therefore consider all exocomets as essentially perfect point sources.

The spectral energy distributions of exocomets cannot be determined from transit observations, so for this aspect we need to rely on information from comets in the Solar System (e.g., Gehrz & Ney 1992). For example, the flux distribution of Hale-Bopp (Shanklin 2000) during its active cometary phase could to first order be formulated as a blackbody plus a strong silicate emission feature around 10 μm (Williams et al. 1997). The silicate emission is however not a ubiquitous feature in cometary spectra (e.g., Hanner et al. 1994, 1996). Hence, except where otherwise stated, here we will simply model the exocomet spectral energy distribution (SED) as a blackbody. This should be a conservative estimation, since silicate emission would make a comet brighter in the wavelength range of maximum sensitivity for LIFE. The SED then depends strictly on the temperature and total effective optical surface of the dust in the coma. Following the modeling in Zieba et al. (2019), we assume that the dust is optically thin. The total optical surface area can then be determined simply from the depth of the transit light curve. We assume that the dust thermalizes quickly enough that the temperature of the grains can be determined from the equilibrium temperature T_eq at the instantaneous star-comet distance r at any given time: $$T_{\text{eq}}={\left(\frac{L_{\ast }\left(1-A\right)}{16\pi {\sigma }_{\text{SB}}{r}^2}\right)}^{1/4}$$(1)

where σ_SB is the Stefan-Bolzmann constant and L_* is the luminosity of the star. In order to calculate T_eq we need an estimate for the albedo A of the dust. Here we again need to rely on Solar System comets. Kolokolova et al. (2004) compile a sample of measured cometary dust albedos, from which we can deduce (for nonextreme phase angles), a range between ~0.1 and ~0.4, with a characteristic mean albedo of 0.25 (e.g., Gehrz & Ney 1992; Mason et al. 1998, 2001). This is similar to the albedo of Earth at ~0.3. Putting this information together, we can note that the median transit depth for the distribution of observable transiting comets in the β Pic system is 174 ppm (Lecavelier des Etangs et al. 2022). For reference, given the radius of β Pic of 1.5 R_⊙ (Zwintz et al. 2019), an Earth-sized planet would exhibit a transit depth of approximately 40 ppm in the system. Hence, a typical comet in the observable transiting population around β Pic has a total optical surface area much larger than that of an Earth-sized planets. Furthermore, at equal separations from the star (characteristically ~1 au for the comet population at time of transit), the comets would also have a similar equilibrium temperature, or perhaps slightly higher (due to the slightly lower albedo), than an Earth-like planet. Since LIFE is scoped to be able to image Earth-like planets in and around the habitable zones of stars out to ~20 pc, it naturally follows that it is also well scoped to image large exocomets such as those observed around β Pic at a distance of 19.7 pc (Gaia Collaboration 2016).

Fig. 1 Schematic representation of the LIFE nulling configuration, where the four unit telescopes U1 through U4 are plotted as circles. The configuration is rectangular with a short axis b and a long axis qb. During observations, the configuration rotates in the plane of the targeted sky patch.

4 Potential target systems

4.1 β Pictoris

β Pictoris is an A6V-type Gray et al. (2006) star with an estimated mass of 1.83 M_sun (e.g., Brandt et al. 2021) and luminosity of 8.5 L_sun (Zwintz et al. 2019). β Pic is one of the most well-studied young nearby stars, and is the defining member of the β Pictoris moving group (BPMG, see Zuckerman et al. 2001), which has a collectively determined age of approximately 24±3 Myr (Bell et al. 2015). Due partly to its youth and proximity, β Pic hosts a particularly prominent debris disk. It was the first circumstellar disk to ever be spatially resolved (Smith & Terrile 1984), after having been identified through the infrared excess of its host star (Aumann 1985). The disk extends out to at least 2000 au (Janson et al. 2021b) and contains a range of structures and asymmetries (e.g., Kalas & Jewitt 1995), including a prominent warp (e.g., Heap et al. 2000; Golimowski et al. 2006) due to disk-planet interactions. To date, two giant planets have been discovered in the system: β Pic b (Lagrange et al. 2010) at 10.3 au and β Pic c (Lagrange et al. 2019; Nowak et al. 2020) at 2.7 au; both with masses in the range of ~9 M_jup (Brandt et al. 2021; Lacour et al. 2021). The architecture of the system relative to the approximate field of view of LIFE is shown in Fig. 2.

Given the large amount of planetesimals available in the β Pic system as evidence by the bright debris disk, and the presence of at least two massive planets to dynamically stir the disk, there are ample opportunities for cometary activity in the system. Indeed, as we have seen in Sect. 3, β Pic is the system where exocometary activity is most concretely proven, and best characterized. The presence of transiting comet data, in particular in this system is what enables us for the first time to make informed estimations about the prospects for directly imaging exocomets. This is the reason for why we use β Pic as the baseline case in Sect. 5.1, which is then the basis for a broader discussion of exocomets in nearby systems in Sect. 6.

Fig. 2 Extent of the β Pic planetary system relative to the ~1 arcsec FOV of LIFE at 10 μm (black dashed circle), as seen from Earth. The orbits of the planets (Nowak et al. 2020) are plotted in blue for the outer planet ‘b’ and red for the inner planet ‘c’. The inner green dashed circle marks the region inside of 1.6 au which is of primary interest for exocomet studies, since coma activity may be very limited outside this range. The edge-on disk in the system extends 200 times beyond the LIFE FOV.

4.2 ϵ Eridani

At a distance of only 3.220±0.001 pc (Gaia Collaboration 2016), ϵ Eri is one of the nearest stars to us, especially in the context of relatively Sun-like stars. Its spectral type is K2V (Keenan & McNeil 1989), and it has an estimated mass of 0.78 M_sun (Carrión-González et al. 2021) and an age in the range of ~400–800 Myr (Janson et al. 2015). It hosts a prominent debris disk (Greaves et al. 1998), which appears to be arranged in several distinct belts (Backman et al. 2009). While ϵ Eri has a high chromospheric activity, which has made precise radial velocity (RV) studies challenging, a planet candidate named ϵ Eri b has been reported in RV observations of the system (Hatzes et al. 2000). The exact orbital parameters have been highly uncertain, with early results implying a highly eccentric orbit, but long-baseline campaigns are converging toward a 3.5 au planet with a mass of ~0.7 M_jup and a low eccentricity (e.g., Llop-Sayson et al. 2021). A (multi)ringed structure such as the one seen in the ϵ Eri disk can imply the presence of additional planets (e.g., Quillen 2006) on wider separations than ϵ Eri b, although direct imaging studies have not yet detected any such planets, setting upper limits of ~1 M_jup from 3 au and outwards (e.g., Janson et al. 2015; Mawet et al. 2019).

Much like β Pic, ϵ Eri therefore hosts a combination of a bright debris disk, indicating a large supply of planetesimals, and one or several wide giant planets, indicating the possibility for substantial dynamical stirring of the planetesimal population. It thus may be one of the most promising targets for exocomet searches. Unlike β Pic, it does not yet feature any direct evidence of such comets, but on the other hand, its extreme proximity offers many advantages for the purpose of detecting exocomets, if such comets do exist there.

4.3 Fomalhaut

Fomalhaut hosts one of the most prominent debris disks in the Solar neighborhood (Kalas et al. 2008). It is an A4V-type (Gray et al. 2006) star with a mass of 1.92 M_sun (Mamajek 2012) and is located at a distance of 7.7 pc (van Leeuwen 2007). A point source known as Fomalhaut b has been repeatedly seen in several epochs of optical imaging with the Hubble telescope (e.g., Kalas et al. 2008, 2013). The visible-light point source has no infrared counterpart, which excludes the possibility that it originates from a planetary surface (Janson et al. 2012), but instead, it must be scattered light from a local concentration of dust. There could still be a planet at the center of the dust cloud (Pearce et al. 2021). However, detailed measurements of the source properties in Hubble data has shown that the point source appears to in fact be mildly extended, and expands with time (Gaspar & Rieke 2020). It also appears to be possibly accelerating outwards as a result of radiation pressure. This implies that the dust cloud is transient, and caused by the destruction of a planetesimal-sized body (of order ~ 100 km) through a catastrophic collision with another planetesimal (Lawler et al. 2015) or through tidal disruption due to gravitational interaction with one of the suspected planets in the system (Janson et al. 2020).

Regardless of the underlying mechanism, the disruption of a planetesimal on an eccentric orbit, and the fact that it was detected serendipitously, indicates a large frequency of cometary-type bodies in the Fomalhaut system (Lawler et al. 2015). The favorable proximity of the system, combined with this high expected cometary activity, makes Fomalhaut a compelling target for exocomet searches.

4.4 AU Mic

AU Mic is an M1V-type star (Keenan & McNeil 1989) with a mass of ~0.5 M_sun (Plavchan et al. 2020), and is possibly the youngest star within 10 pc, given its distance of 9.7 pc (Gaia Collaboration 2016). AU Mic is established as a member of the BPMG (e.g., Malo et al. 2014), just like β Pic, which results in an estimated age of 24±3 Myr (Bell et al. 2015). Also much like β Pic, AU Mic hosts a very prominent edge-on debris disk (Kalas et al. 2004), with interesting substructures. In particular, several dust clumps have been identified in the disk that move outwards from the star (Boccaletti et al. 2015, 2018), possibly implying recurring dynamical phenomena in the inner parts of the system. Additionally, there are two known transiting planets in the system, AU Mic b and c (e.g., Plavchan et al. 2020; Szabó et al. 2022), at semi-major axes of 0.06 and 0.11 au and masses of 12 and 22 M_Earth, respectively (Zicher et al. 2022), further emphasizing opportunities for dynamical stirring and exocometary activity in the system. AU Mic therefore constitutes another concrete example of the a priori most interesting systems to study in the context of direct imaging of exocomets.

4.5 Other systems

Beyond the individual examples that are highlighted above, there are a large number of systems in the Solar neighborhood where exocomet searches might be expected to be particularly fruitful. Debris disk targets constitute a promising object class in this context, since the dust observed in such disks is produced continuously through planetesimal disruptions – i.e., a debris disk is evidence for the presence of a large number of planetesimal bodies in the system. Rebollido et al. (2020) do not find a correlation with debris disks in general for targets with FEB-type comet signatures, but do find a correlation with near-infrared excess.

Within 20 pc, there are 20 systems with known debris disks listed in the Circumstellar Disks Catalog³. Aside from the already highlighted systems, noteworthy examples include GJ 581, which hosts a debris disk along with several RV-detected planets (e.g., Bonfils et al. 2005; Robertson et al. 2014); η Crv, which hosts a particularly bright warm debris disk, implying planetesimal activity around the habitable zone region (e.g., Defrère et al. 2015); and Vega, which hosts a debris ring with a wide gap inside of it (Hughes et al. 2012), similar in structure to systems like Fomalhaut or ϵ Eri.

Another path toward candidate exocomet imaging systems is if there are other indicators for exocomets, such as spectral lines that can be associated with FEBs, or exocomet transits. As discussed in Sect. 3, some such targets are already known, although they are not always suitable for observations with LIFE. For example, the exocomet transit target KIC 3542116 (Rappaport et al. 2018) is much too distant at ~260 pc to be a feasible target for imaging around the ~1 au range, and the FEB target HD 172555 (Kiefer et al. 2014a) is in principle feasible, but not ideal at 28.8 pc (Gaia Collaboration 2016). Both TESS (Ricker et al. 2015) and CHEOPS (Benz et al. 2021), and in the future PLATO (Rauer et al. 2014) are sensitive to exocomet transits around nearby stars, and may thus provide additional proven exocometary systems well ahead of the LIFE mission.

Beyond targeted exocomet searches, exocomets could conceivably be detected serendipitously as part of the larger LIFE program. LIFE will likely spend a large fraction of its mission time studying relatively Sun-like main-sequence stars, with a sufficient contrast to detect both terrestrial planets at various stages of their evolution and exocomets. It is not possible to make an accurate estimation of how common observable exocomets will be in such systems, not least because the Solar System may not necessarily be representative of a typical main-sequence star in this regard. However, as a rough estimation, we can consider the fact that in the Solar System, so-called ‘great comets’ with comparable sizes to what has been inferred from the exocomet distribution in Lecavelier des Etangs et al. (2022) occur approximately once per decade. Since the LIFE mission will extend for several years, and perhaps up towards a decade, it could be possible for individual exocomet events to be discovered in the main sample through sheer serendipity, but most likely not at any high rate.

5 Simulations

Here we present the outcomes of a series of LIFE observation simulations for the specific case of β Pic, as a baseline case for considering exocomet observability.

5.1 Frequency and detectability of comets

As mentioned previously, 30 exocomet transits have already been identified in the β Pic system. This allows us to both assess their individual prospects for being directly imaged, as outlined in Sect. 3, but also to estimate the rate at which they might be expected to occur. 156 days of TESS observations revealed 30 individual cometary transit events, with absorption depths between 107 ppm and 1963 ppm. For interpreting and simulating the corresponding size of the effective optical surface, it is useful to formulate this surface in terms of the radius (in units of Earth radii, R_E) that a circular disk with the same effective surface area as the cometary dust coma would have. For the observed β Pic exocomet population, this gives effective radii between 1.6 R_E and 7.0 R_E.

While exocometary transits can only be observed at a narrow range of orbital inclinations, direct imaging can reach exocomets at any inclination. Thus, to estimate the frequency of exocomets around β Pic amenable to direct imaging, we need to make some assumption about the inclination distribution of the comets. Since the β Pic circumstellar disk and planets are all encompassed by a narrow inclination distribution which is close to edge-on (~85–88 deg in different parts of the warped disk, e.g. Kalas & Jewitt 1995; Ahmic et al. 2009; Janson et al. 2021b), it might be tempting to assume that the exocomets follow a similar distribution. In this case, a relatively large fraction of the exocometary population may already be captured in the transit monitoring. However, it is important to note that assuming such a distribution would give drastically incorrect results for comets in the Solar system. We have compiled a list of all object classified as comets with semi-major axes <100 au in the Small-body database⁴, and show the distribution of inclinations in Fig. 3. The observed transiting comets are expected to typically reside at ~1 au from the star when the transit occurs (e.g., Zieba et al. 2019). At this separation, a transit occurs if the comet has an inclination within ± 0.48 deg from perfectly edge-on. Hence, the observable inclination band is 0.96 deg. We therefore use 0.96 deg as the sampling size for the distribution in Fig. 3, and normalize by the total number of objects, in such a way that the y-axis reflects the fraction ϕ of Solar systems comets that would be detectable in transit if viewed from a particular direction (as represented by the x-axis). Although there are many more comets close to the ecliptic plane than far from it, the peak of the distribution occurs quite far (~ 10–20 deg) from the ecliptic, and the distribution is so broad than no particular ~0.96 deg inclination interval captures any substantial fraction of the comets. The few most favorable directions capture up to ~4–5% of the comets, and the least favorable directions much less than 1%. If we assume that β Pic has a similar cometary inclination distribution as in the Solar system, it seems reasonable to assume ϕ = 0.01 as a characteristic value for the fraction of comets that exhibit transits. We therefore adopt this value as a baseline, but at the end of the section we will also note the implications of adopting other values of ϕ.

In order to evaluate the detectability of any individual exocomet, we use the LIFEsim (Dannert et al. 2022) software, which simulates the achievable signal-to-noise ratio S/N for a given integration time t_int, and a range of stellar, observatory, and “planetary” (in our case, cometary) parameters. The required stellar system parameters are the effective temperature T_*, radius R_*, distance d_*, ecliptic latitude l_*, and exozodiacal level z_* in units of Solar system zodiacal brightness (zodi). We use T_* = 8100 K and R_* = 1.5 R_sun from Zwintz et al. (2019), and d_* = 19.7 pc, and l_* = 4.98 rad from the NASA exoplanet archive based on Gaia measurements (Gaia Collaboration 2016). The exozodiacal level in the β Pic system is unknown, so we run simulations with z_* = 1, z_* = 10, z_* = 100 and z_* = 1000 zodi. Since β Pic has such a bright debris disk in general, it is likely that the zodiacal disk is also quite bright, so we use z_* = 100 zodi as a baseline value.

The observatory parameters refer to the spectral resolution, the sizes of the unit telescopes, and the interferometric baseline between the telescopes. For spectral resolution, we selected the default setting, which is R = 20. For telescope aperture sizes, a frequently used reference value for the mission is 2 m, which we choose as a characteristic size also in this study. For interferometric baseline, we experimented with different lengths, and settled for 25 m as a main option, since this offers good sensitivity across the most interesting range for detecting exocomets (projected separations of ~0.5–1 au).

The relevant planetary/cometary parameters, assuming a blackbody flux distribution, are the projected separation, the effective temperature, and the effective radius. For studying a range of possible parameter values, we implemented a 2-dimensional grid system, where effective radii ranged from 2 R_E to 7 R_E in steps of 1 R_E, and projected separations ranged from 20 to 80 mas in steps of 20 mas. To calculate an effective temperature for each grid point, we assumed a template comet with nearly identical orbital parameters as Halley’s comet (as listed in the Small-body database), except that the orbital period was scaled according to the β Pic system mass based on a fixed semi-major axis of 17.834 au, and the inclination was set to 60 deg in order to reflect a median inclination for a highly scattered (effectively uniform) hypothetical exocometary distribution. We simulated the orbit based on these parameters, and with a sampling frequency of 1 day, we calculated the effective temperature (based on the instantaneous star-comet separation, see Sect. 3) and projected separation at each sampled time-step. Choosing (arbitrarily) the post-periastron section of the orbit, we then have an estimate for the effective temperature as function of projected separation for comets on Halley-like orbits, which we interpolate to sample at the desired grid points of 20–80 mas, with a step size of 20 mas.

The S/N values calculated with LIFEsim across the separation-radius grid with our default parameters are shown in Fig. 4. It is immediately clear that a wide range of cometary bodies would be detectable at high significance in the β Pic system in 10 hours of integration time. Our baseline setting has z_* = 100 zodi, but we also plot calculations for other exozodiacal levels in the same figure. As expected, lower exozodiacal levels give even higher S/N values in the same amount of time, but even with substantially higher levels (z_* = 1000 zodi), detection remains possible. In Fig. 5, we show alternative cases with default stellar and cometary properties, but varying the telescope unit aperture size, to account for the boundary scenarios of possible mission sizes in Quanz et al. (2022). The S/N scales approximately proportionally to unit telescope size, as expected. We also checked the detectability impact of the interferometric baseline length. The overall impact changing the interferometric baseline is minor, but the sensitivity increases at smaller separations and decreases at larger separations for longer interferometric baselines, and vice versa.

With S/N values across the full size-separation grid in hand, and prior knowledge of the frequency and effective size distribution of the β Pic comet population from previous works, we can now evaluate how commonly observable exocomets can be expected to occur. We do this by once again assuming the same Halley-like template orbit as before, considering the half orbital arc that occurs after periastron for simplicity. We then generate 1000 simulated exocomets where each comet is randomly assigned a size from the sample of known exocomet effective sizes generated from the absorption depth distribution in Lecavelier des Etangs et al. (2022). Out of the 30 exocomets listed in Lecavelier des Etangs et al. (2022), We only use the 18 exocomets with effective sizes >2 R_E for representing the size distribution, since smaller comets are never observable with S/N > 5 in 10 hours of observations for β Pic under the baseline assumptions. For each comet, we interpolate the S/N grid at the locations corresponding to the comet’s size and the projected separations given by the simulated orbit (which is sampled with a 1-day cadence). We calculate for how many days each exocomet is observable, where S/N > 5 is the threshold for what is considered observable. The average duration of observability is 83.5 days across the simulated orbital half-arc, and therefore approximately two times that number for a full arc. The average duration of observability (in days) multiplied by the number of relevant exocomets that enter the observability zone per day gives the number of exocomets that are observable at any given time. The number of entering comets can be calculated by considering the fact that the 18 selected transiting cometary events occurred over a TESS observing window of 156 days in total. Hence, 0.115 sufficiently large transiting comets enter the window each day, so if the fraction of comets that transit is ϕ = 0.01 as discussed above, 11.5 potentially imageable comets enter the same window on a daily basis. In total, there is then an average of ~960 exocomets that are bright enough to be observable with LIFE at any random point in time.

This remarkably high number would in fact mean that the observations would be heavily confusion limited, both since the number of exocomets would far exceed the number of resolution elements inside of ~80 mas, and since LIFE has a limited effective UV-plane coverage, making it difficult to accurately reconstruct a complex mix of point sources, compared to reconstructing a field containing just one or a few point sources. We will discuss this in more detail in Sect. 5.2 and outline the implications in a broader context in Sect. 6, but as a basic first conclusion, it appears clear that individually detectable exocomets are plentiful in the β Pic system, which can probably be seen as constituting an upper envelope for the number density of exocomets around nearby stars.

The above calculations depend heavily on the value for the parameter ϕ, which as we have discussed is essentially unknown for β Pic, but is around ~1% if the solar system distribution of comets is representative. Since β Pic is younger and its disk and planets more massive, the orbital distribution of comets could also be quite different. For example, the close-in FEB population has been suggested to arise from resonances with the planets in the system (e.g., Beust & Morbidelli 2000). This might force subpopulations of the exocomets into relatively narrow orbital families. In the FEB population, two subfamilies have been identified, which have scatters in the argument of periapsis of ±8 deg and ±25 deg respectively (Kiefer et al. 2014b). The corresponding scatter in inclination distributions has not been measured from FEB data, and the FEB population is by no means necessarily representative of the exocomet population in total or the TESS transiting subpopulation. Hence, accurate clustering properties in the potentially imageable population cannot be predicted, but as a general principle, if the β Pic comets happen to be both narrowly distributed and fortuitously aligned toward the line of sight, ϕ could be higher and consequently, the number of imageable exocomets could be lower than calculated above. Still, since there are 960 sufficiently bright comets in the observing window if ϕ = 0.01, it follows that there are 9.6 comets in the same window even in the extreme case where ϕ = 1, such that every single comet in the system transits. Since all orbits would be edge-on in this circumstance, projection effects would bring those 9.6 comets closer to the effective inner working angle than what would be the case for the template orbit, but the majority of them would still be detectable. Hence, regardless of the specific inclination distribution of comets in the β Pic system, the simulations indicate that at least several exocomets would be formally imageable with LIFE at any given time.

Fig. 3 Histogram of the inclinations of <100 au comets in the Solar system. While low inclinations are generally statistically favored, the spread of the distribution is quite wide. Thus, if the comets in a system like β Pic share a similar orbital distribution, then only a small fraction of comets will transit as seen from any given orientation.

Fig. 4 S/N values (color scale) as function of temperature and effective coma radius for simulated exocomets in the β Pic system, in 10 h of observing time. Upper left: a case with low zodiacal contamination (1 zodi). Upper right: a higher level of 10 zodi. Lower left: the baseline value of 100 zodi. Lower right: a very high level of 1000 zodi.

Fig. 5 Same as Fig. 4, but varying the sizes of the unit telescopes from the baseline value (2 m diameter). The observing time in each case is 10 h. Left: 3.5 m apertures. Right: 1.0 m apertures.

5.2 Cometary crowding

The configuration and signal combination planned to be used for LIFE is described in detail in Dannert et al. (2022). In short, the signal combination from the four unit telescopes nulls out the on-axis light from the parent star, and for the rest of the field, it effectively creates a transmission pattern on the sky with alternating positive and negative regions. As the telescope configuration rotates, so does the transmission pattern, so a planetary or cometary off-axis companion moves in and out of positive and negative areas as seen by the telescope array. An off-axis companion at a given separation and position angle will therefore give rise to a unique modulation of the interferometric signal, as illustrated in Fig. 6. Thus, the astrometric and spectrophotometric properties of a single isolated off-axis source can be easily backed out from the nulling signal, even with the very limited uυ coverage provided by a four-element nulling array.

If several off-axis sources exist in the same field, their contributions blend into a mutual interferometric signal, which makes the signal interpretation much more challenging. How to most efficiently back out information about the individual sources in this context is work in progress, and beyond the scope of this study. However, since we can expect some exocometary fields to be crowded, at least in the case of β Pic as outlined above, we will discuss an example case of a moderately crowded field here in order to illustrate the impact of such crowding on the interferometric signal.

A simulation with 18 exocomets at a wavelength of 10 μm is shown in Fig. 7. The comets have the same size distribution as the 18 largest comets observed in the β Pic system (see Sect. 5.1), and the mass and distance of the star is also set to the same as β Pic. In order to make the illustration as clear as possible, we consider a case where all comets share the same orbital plane, and this plane is viewed face-on by LIFE. This means that the projected separation is always the same as the physical separation, such that comets that lie closest to the star in image space are also physically the closest to the star, and therefore hotter than the more distant comets. We also assume that all comets have exactly the same orbital parameters, except for the argument of periapsis ω which is randomly distributed between 0 and 360 deg. The orbits are Halley-like with semi-major axes of 17.834 au and eccentricities of 0.967.

Three different dates are simulated: an original epoch of observation called Day 0, a follow-up observation one week later called Day 7, and another observation 1 month from the original epoch called Day 30. In all three epochs, the total signal is dominated by the largest comet in the field (plotted in blue). This is especially true at Day 30 when the largest comet is also one of the closest to the parent star, and therefore one of the hottest in the field. This increases the total interferometric signal, and also makes it cleaner in the sense that it is dominated by a single source whose properties could probably be quite accurately backed out. On Day 0 and Day 7, the signal would be considerably more degenerate and prone to larger uncertainties. Crowding can therefore impose difficulties in the interferometric signal interpretation, and in the case of the β Pic system this is of course emphasized by the fact that, as we have seen in Sect. 5.1, the crowding there could be up to 50 times larger than in this simulated example. This would mean that even many individual resolution elements are crowded.

In a full signal extraction, it is possible> to mitigate the ambiguities resulting from a sparse sample of the uυ plane by accounting for all the spectroscopically observed wavelengths simultaneously, instead of just the 10 μm signal as in this illustrative example. It would also be possible to extend the coverage further by making additional observations with different baseline lengths. Nonetheless, we expect that the extreme crowding expected in the particular case of the β Pic system would make it a challenging target for LIFE. In the meantime, in the following two sections, we will disregard any possible issues associated with overcrowding of comets, and investigate what information can be acquired from an isolated exocomet, if observed over a handful of epochs.

Fig. 6 Transmission for the nulling interferometric mode used in this study. Upper: transmission pattern in polar coordinates at 10 μm of the nulling interferometric configuration used in this study. A comet at a separation of 40 mas would project along the black dashed line as the telescope array rotates. Lower: a cross section of the transmission pattern for a separation of 40 mas (i.e., along the dashed line in the upper panel). A comet residing at a position angle of 0 deg would map the transmission pattern of the black solid line. A comet at a different PA would result in a cyclically phase-shifted pattern of otherwise the same shape. See Dannert et al. (2022) for further details.

5.3 Astrometric constraints

Given that exocomets appear to be potentially detectable under certain circumstances, it is relevant to examine to which extent their properties can be constrained, and to which extent they can be distinguished from other possible point sources in the system (such as habitable planets, for example). Here, we focus on constraints related to astrometric measurements of exocomets.

In Dannert et al. (2022), it is found that the S/N-limited astrometric precision of LIFE for a point source with S/N = 9.7 is 1.5 mas. Beyond a limited S/N, astrometric precision is typically subject to calibration uncertainties related to the physical scaling and orientation of the astrometric grid. While such calibration details have not been examined in detail for LIFE yet, experience with astrometry from other facilities shows that providing calibration uncertainties at the ~1 mas level or far below is a demonstrably manageable task (e.g., Maire et al. 2016; GRAVITY Collaboration 2021). We therefore assume that the measurements would be S/N-limited. Under this circumstance, the astrometric uncertainty δ typically scales as δ ∝ (S/N)⁻¹ (e.g., Lindegren 1978).

In order to assess astrometric performance in a reasonably realistic context, we imagine a scenario in which a point source, suspected as a possible exocomet, is detected in observations acquired at date T₁. Follow-up observations are then acquired at two additional epochs (T₂ and T₃), taken 20 and 40 days later. Each observation contains 10 h of integration time. Astrometric data points are simulated by setting the “true” orbital parameters of the comet as corresponding to a Halley-type orbit with a semimajor axis a_T of 17.834 au, an eccentricity e_T of 0.967, and an inclination i_T of 60 deg just like in Sect. 5.1. The period for a system mass of 1.75 M_sun comes out to P_T = 56.9 years. The time of periastron T_T was set to an MJD of 58240.7 days, which in turn was set to three days before the first observation, T₁ = T_T + 3. The argument of periastron ω_T was set to 111.33 deg, and the ascending node Ω_T to 238.42 deg. The effective size of the exocomet coma was set to 5 R_Earth.

For each simulated observation date, we extracted the astrometry relative to the central star as ∆x₁ = 7.5 mas (eastward) and ∆y₁ = 19.6 mas (northward) at T₁, ∆x₂ = 38.2 mas and ∆y₂ = 22.9 mas at T₂, and ∆x₃ = 57.6 mas and ∆y₃ = 17.9 mas at T₃. Based on the physical and projected separations at each date, we then calculated S/N values using LIFEsim, acquiring values of 37.2, 29.6 and 15.0, leading to estimated astrometric uncertainties of 0.4 mas, 0.5 mas and 1.0 mas, respectively. Based on these uncertainties, we generated random errors with the corresponding σ, assuming Gaussian distributed errors, and added the errors of such a random set to the true astrometric positions. This resulted in simulated measurements $\text{Δ}{\tilde{x}}_1=7.8\pm 0.4$ mas and $\text{Δ}{\tilde{y}}_1=19.4\pm 0.4$ mas, $\text{Δ}{\tilde{x}}_2=37.9\pm 0.5$ mas and $\text{Δ}{\tilde{y}}_2=22.1\pm 0.5$ mas, and $\text{Δ}{\tilde{x}}_3=57.2\pm 1.0$ mas and $\text{Δ}{\tilde{y}}_3=18.3\pm 1.0$ mas.

It is clear from basic considerations that it would be impossible to determine detailed orbital properties of the comet based on the astrometric information given in this simulated example: with only three astrometric epochs, accurate orbital fitting is challenging under any circumstance, but even more so in this context, where the observational baseline only covers a tiny fraction (~0.2%) of the full cometary period. Nonetheless, the observations can be expected to provide some useful constraints on the orbit. In particular, the speed and acceleration of a highly eccentric object near periastron are significantly different from (e.g.) a smaller orbit of low eccentricity at the same instantaneous separation. The instantaneous orbital velocity υ_orb is given by: $${\upsilon }_{\text{orb}}=\sqrt{G{M}_{\ast }\left(\frac{2}{r}-\frac{1}{a}\right)}$$(2)

where G is the gravitational constant, M_* is the stellar mass, r the instantaneous separation and a the semi-major axis. Thus, for example, an object in a circular 1 au orbit around β Pic has an orbital velocity of 40 km s⁻¹, while a high-eccentricity object with a semi-major axis of 100 au has an orbital velocity of 57 km s⁻¹ when it crosses 1 au. Meanwhile, for large a (small 1/a), 2/r becomes the dominant term and υ_orb becomes very weakly dependent on a – as a result, the difference in velocity between a 10 au modestly eccentric orbit and a 100 au highly eccentric orbit when crossing 1 au is only ~2%. As a consequence, while we can expect to be able to distinguish between low and modest-to-high eccentricities in general, we cannot expect to distinguish between modestly high and very high eccentricities, with the limited astrometric information provided in the simulated example.

To quantify this, we ran orbital fitting with the OFTI code within the orbitize package⁵ (Blunt et al. 2020) on the simulated data points. We use a standard set of uninformed priors for the fitting, and run the simulations for 10⁴ orbits. The results of the fitting are shown in Fig. 8, with example orbits in Fig. 9. As expected, the true orbital parameters are poorly recovered. Of particular importance is the disagreement between the true eccentricity of 0.967 and the recovered eccentricity of ${0.71}_{-0.11}^{+0.10}$.

While it is straightforwardly expected that the best-fit eccentricity fails to provide a close match to the true value as argued above, it is perhaps more intuitively surprising that the true value is not recovered within the estimated error bars of the posterior (1.7σ deviation). We believe that the primary reason for this is the combination of a very high eccentricity with a phase very close to periastron. This is a priori a highly unlikely configuration, which means it cannot acquire sufficient probabilistic weight to contribute substantially to the posterior. In fact, the only way to observe a Halley-like exocomet is when it is extremely close to periastron, so the reason that we expect Halley-like exocomets to be observable in the first place is closely connected to this strong selection effect, which the fitting algorithm has no good way to account for. In a theoretical setting, we can of course adjust the priors in the fitting procedure to account for this selection effect – however, this solution is not applicable to the real world, where we have no prior information on the orbital distribution of comets. For example, modestly eccentric bright comets may be much more common in the young and dynamically active β Pic system than in the Solar system, and skewing the priors toward extreme eccentricities may yield erroneous results for such a population.

As a result, we find that it would be impossible to set any stringent upper limit on the eccentricity based on the astrometric data in our simulated example, due to unavoidable biases in the procedure. However, there is no corresponding bias for the lower limit on the eccentricity – that part of the parameter space is disfavored in the posterior simply on the basis that such orbits are not consistent with the existing simulated data. We thus conclude that the lower limit on the eccentricity is reliable, and that it is therefore possible to at least distinguish orbits of low eccentricity from orbits of intermediate or higher eccentricity, to a good degree of accuracy. Fitted parameters other than the eccentricity generally cannot be relied on in this context – since the procedure misses the “true” solution in terms of the eccentricity, it is likely to falsely estimate the global optima for other orbital parameters as well.

Fig. 7 Simulations of a crowded exocometary field. The three rows correspond to three different dates. Upper: day 0. Middle: day 7. Lower: day 30. The left column contains maps of the spatial distribution of the exocomets. Each circle is an exocomet, where the radius of the circle is directly proportional to the effective radius of the cometary coma. Each comet has its own color, and the dashed lines of each color corresponds to the orbital track of that particular comet. The right column shows the interferometric signal as a black thick line. For illustrative purposes, the contribution of each individual comet is plotted as a thin line that is color-matched with the corresponding symbol in the left column. The individual colored components cannot be measured directly; they can only be deduced from the total signal and its evolution in time, wavelength, and (if relevant) along different probed interferometric baselines. See the text for discussion.

Fig. 8 Corner plot for the most important orbital parameters from the orbit fitting test described in the text, with a 95% cut on the axes for the 104 orbital solutions. The true orbital parameters are not well recovered, due in large part to the short baseline probed with the simulated observations.

Fig. 9 20 example fitted orbits (colored lines), randomly sampled from the full population of fitted orbits. The orbits are generally significantly smaller and less eccentric than the input orbit, underlining the fact that the probed baseline is too short to determine exact parameters, although meaningful limits can be obtained (see text for details). Left: the orbital fits in two dimensions. The red star symbol marks the location of the star, and the yellow points mark the locations of the comet at three measured epochs. Upper right: separation as function of time. Upper left: position angle as function of time. Purple dots mark the location of the comet at the respective epochs.

5.4 Spectroscopic constraints

As mentioned previously, some Solar system comet spectra are featureless blackbody spectra, while others exhibit strong silicate emission. So far in this paper, we have only considered the pure blackbody case. Even in this case, useful information can be extracted from the spectra – the 4–18.5 μm spectral range of LIFE includes the blackbody peak for temperate exocomets, and the peak wavelength scales inversely with the temperature of the comet. Hence, the temperature of the dust can be derived at each measured epoch. By monitoring the temperature as function of separation over several epochs, it may in turn be possible to derive an estimate for the albedo of the dust in many cases.

However, a perhaps even more compelling set of information could be extracted if the comet shows silicate emission, similar to e.g. comet Hale-Bopp (Williams et al. 1997) in the Solar system. This would allow one to place constraints on the mineralogy of the dust, Mineralogical studies have been made of dust in a wide range of disks (e.g., Bouwman et al. 2008), including the disk of β Pic itself (Lu et al. 2022). Such studies provide constraints on the average composition of dust in the disk, which is a mixture of dust originating from many different bodies. Through spectroscopic studies of exocomets, dust from an individual body can be studied, providing for a potentially very clean analysis of its composition. Here, we have performed a test case to evaluate to which extent an exocomet exhibiting silicate emission features might be spectroscopically characterized.

For the purpose of this test, we assume a comet with the same effective coma size as in Sect. 5.3 (5 R_E) at an instantaneous star-comet separation of 0.6 au. We then model the spectral output of such a comet, assuming that it exhibits silicate emission, based on different compositions of the silicates. Following Lu et al. (2022), we source models of Enstatite from Chihara et al. (2002), Forsterite from Zeidler et al. (2015), and olivine and pyroxene from Jaeger et al. (1994) and Dorschner et al. (1995). We sequentially model one species at a time and note its impact on the output spectrum. This is obviously a highly idealized procedure, since a comet may in reality contain a mix of different silicate species; but the point of the analysis is to probe whether mineralogical variations can be confidently measurable in the first place, rather than attempt to realistically model the actual composition of an exocomet, which we have no way of knowing prior to any actual measurements.

The spectrum of the hypothesized exocomet is a combination of blackbody emission and silicate emission, so we start from a blackbody spectrum in the same way as for the featureless comets in Sect. 5.1. We then redistribute half of this flux into silicate emission – in other words, we halve the blackbody flux, normalize the pure silicate emission spectrum (of an individual species) to the same total flux (within the LIFE spectral range), and add the two together. The equal flux redistribution is arbitrary, but it represents a case which has less fractional silicate emission than, e.g., Hale-Bopp, but more fractional silicate emission than a featureless spectrum like that of Halley’s comet. It is therefore a mix encompassed by objects known to exist in nature. We then use the composite models in LIFEsim, using the same observational parameters as in the preceding sections. In addition to the standard 10 h simulated observations as before, we also ran 100 h observations to illustrate the effect of increased observing time on the achievable S/N.

The output S/N as function of wavelength for the simulations of each separate silicate species is shown in Fig. 10. In order to more clearly show how these S/N values affect the measured spectra, we selected Enstatite and Forsterite as two example cases and added random Gaussian noise to the respective spectral models with noise levels based on the simulated 10 h S/N values at each sampled wavelength. The results are shown in Fig. 11. It can be seen that at least for this relatively large comet, a single 10 h observation is sufficient to start distinguishing pure Enstatite from pure Forsterite, with e.g. Forsterite being distinguishable from an Enstatite model by more than 5σ in some spectral channels. Hence, some basic mineralogy can be performed in this circumstance, and distinguishing the presence of silicate emission from pure blackbody emission is easier still.

Fig. 10 Simulated S/N values per wavelength bin using LIFEsim, for four cases of possible cometary emission with different silicate species: Enstatite (blue), Forsterite (red), Olivine (green) and Pyroxene (purple). Solid lines denote 10 h observations, while dashed lines represent 100 h observations.

6 Discussion

As we saw in Sects. 5.1 and 5.2, large comets may potentially be so abundant in the β Pic system that exocomet observations with LIFE could become confusion limited – i.e., each resolution element close to the central star may contain several comets simultaneously. Such a scenario would essentially yield a diffuse cloud around β Pic, similar to an exozodiacal disk. Since the beam combination nulling method used for detecting point sources with LIFE cancels out symmetric features, the majority of this cloud would simply not get picked up by the observations. Since the expected brightness distribution of the exocomets is highly skewed, with a large number of small comets and a small number of large ones, the brightness distribution observed by LIFE could generally be expected to be dominated by a small number of discrete point sources, corresponding to the brightness peak of the exocomet distribution. Since the flux distribution observed by LIFE has to be synthesized from a limited number of interferometric baselines, mapping such a distribution might still be a challenge. Future work will be required to evaluate to which extent a large simultaneous number of point sources can be distinguished and mapped out by LIFE.

The above discussion is based on a scenario in which the β Pic cometary distribution is relatively uniformly spread in terms of orbital inclinations (ϕ = 0.01). If the distribution is much more concentrated closely to the line of sight ϕ ~ 1, overcrowding may be much less of an issue. In either case, the high implied rate of bright comets around β Pic clearly indicates that exocomet detection is a feasible task, at least in young systems with a high dynamical activity. If the β Pic system is in fact over-crowded with bright exocomets, this speaks in favor of exocomet detection in other, slightly less active systems. In this context, ϵ Eri and Fomalhaut may be particularly promising targets, for several reasons: both have bright debris disks but are an order of magnitude older than β Pic, which means that their cometary activity is probably intermediate between β Pic and regular main-sequence stars. Thus, they may feature a favorable balance between under- and overcrowding of exocomets in the field of view. Furthermore, both are at much smaller distances from us than β Pic, such that the achievable S/N (for equal-size comets) is much higher, and the inner solar system is much more finely resolved. Hence, while β Pic is the template object for this study (simply due to the fact that it is the only system for which we have detailed information about the cometary distribution), it appears probable that a system like ϵ Eri is probably an ideal target for exocomet detection and characterization purposes.

The detectability of exocomets have two separate but similarly important implications: on one hand, it opens up the possibility to study a new class of astrophysical objects directly, potentially allowing for new insights into the building blocks of planetary systems. On the other hand, the overlap in observational properties between exocomets and terrestrial exoplanets means that we also need to consider that they may potentially be mistaken for each other. For example, a survey primarily directed toward habitable exoplanet detections may need to consider exocomets as one possible source of false positive signal. This should generally not be an issue for old single stars like the Sun, since observable exocomet events are very rare in such cases, as we noted in Sect. 4.5. However, in bright debris disk systems in particular, exocometary interloper events could be quite common, as we have seen, so among such targets it is an important consideration. It may also be relevant for certain types of multiple systems.

For example, it has been hypothesized that the orbit of Proxima Cen around α Cen A and B could cause an enhanced cometary activity in the system, at least during certain epochs (Wertheimer & Laughlin 2006). In this context, we note that a particularly deep infrared imaging survey has been conducted around α Cen A and B with VLT/VISIR (Kasper et al. 2019). The proximity of the system and depth of the observations allow for sensitivity to much smaller and more temperate planets than would otherwise be reachable with a ground-based 8m-class facility. A possible candidate for such a planet in the system has been proposed in Wagner et al. (2021), in the form of an extended feature close to the coronagraphic edge, which could possibly be interpreted as the trace of a planet undergoing orbital motion during the 33 days over which the observations spanned, at some point disappearing behind the coronagraph edge. The planet would be temperate (semi-major axis of ~ 1.1 au), and have a size in the range of ~3.3–7 R_E. However, these ranges also overlap observationally with the exocomet population around β Pic. Hence, if any similar population exists in the α Cen system, even if much less numerous than the β Pic population, the observations could be equally well explained by a large exocomet. In this case, we do not consider either the planet or the comet scenario to be probable, since the reported candidate in Wagner et al. (2021) is explicitly tentative, and could also be explained simply as a residual noise feature. Nonetheless, the example highlights that we are currently at a phase in high-contrast imaging and interferometry where we need to start seriously consider the observability of exocomets, both as scientific targets and as potential contaminants for other science cases.

The two most common classes of astrophysical false positives for directly imaged exoplanets are chance alignments of background stars, and unresolved background galaxies. In both of these cases, the most common method to distinguish between true and false positives is to observe the target system at least twice, with a time baseline of months or years between the epochs (e.g., Janson et al. 2011). If the candidate point source does not share a common proper motion with the target star, it can be established as a background contaminant. For exocomets, this method could not be used in quite the same way for distinguishing them from planets, since the exocomet is bound to the star and therefore does share a common proper motion with it, just like a planet. However, as we saw in Sect. 5.3, astrometry at more than one epoch can provide a lower limit on the orbital eccentricity. While planets can have comet-like eccentricities (e.g. >0.7), classically habitable planets cannot, since they would spend a large fraction of their orbits outside of the habitable zone. It is therefore possible to at least distinguish habitable planets from exocomets in this manner. Furthermore, if followed up on a timescale longer than a few months, the exocomet would exit its near-periastron phase and turn quiescent, thereby disappearing entirely from view. A habitable planet can also temporarily disappear from view if on a highly inclined orbit, but would return back into view again over timescales of a few months, unlike a comet.

If an exocomet exhibits strong silicate emission (see Sect. 5.4), then it can be very concretely distinguished from planets, which would not show such features. However, many Solar system comets do not exhibit clear silicate emission, but rather show blackbody-like spectra. Much like in the astrometric case, in this scenario we would not necessarily be able to distinguish exocomets from planets in general, but we would be able to distinguish them from classically habitable planets, since they would have atmospheres and therefore detectable atmospheric features (e.g., Quanz et al. 2022). Distinguishing silicate emission features in exocomets would open a novel line of inquiry into the composition and differentiation processes of rocky and icy planetesimals, complementary to the results from polluted white dwarf, which indicate a fraction of exo-planetesimals to be geophysically processed (Bonsor et al. 2020, 2022). Constraining the composition and internal processing of wide-orbit planetesimals and exocomets can thus give a handle on late accretion phases and bombardment epochs of rocky exoplanets, altering their secondary atmospheres (Lichtenberg & Krijt 2021; Lichtenberg & Clement 2022).

From a probabilistic viewpoint, we do not expect exocomets as potential false positives for terrestrial planets to be a big problem for studies around mature Sun-like or low-mass stars. However, infrared emission from magma ocean atmospheres in the aftermath of giant impacts among rocky exoplanets in young systems promise to reveal insights into the climate state of young, Hadean-like exoplanets (Lupu et al. 2014), with distinct emission features in the optical (Boukrouche et al. 2021) and infrared (Bonati et al. 2019). Therefore, it is an issue that should be kept in mind when planning observational strategies for a mission such as LIFE.

While the focus of this study is infrared detection and characterization of exocomets, particularly with LIFE, there are of course also prospects for detecting exocomets through scattered light in the visible-light regime, if a sufficient contrast and sensitivity can be reached. As we have seen, comets of the type seen in transit around β Pic have similar effective optical areas as rocky planets, and based on Solar system comets, their albedos are somewhat lower than, but broadly similar to, Earth-like planets. Thus, much like for LIFE in the infrared, a mission scoped to image Earth-like planets at visible wavelengths (such as LUVOIR, see e.g. The LUVOIR Team 2019), would probably have good prospects for being able to also image exocomets. Due to the potentially low comet albedo and fewer distinctive potential spectral features than in the infrared, visible light detection or characterization of exocomets might be more challenging than in the infrared case for a mission scoped primarily towards Earth-like exoplanets, but the yield may nonetheless be substantial in suitable debris disk systems. Any detailed discussion on the potential exocometary yields of a LUVOIR-like mission is beyond the scope of this study, but we note that if exocomets could be discovered at both visible and infrared wavelengths, the broader wavelength coverage would further enhance the information that can be extracted about their characteristics.

Fig. 11 Example of exocomet spectroscopy with LIFE. Top: simulated LIFE spectra for two hypothetical exocomet cases, with noise levels set based on LIFEsim calculations for 10 h observations. The two cases correspond to pure Enstatite (blue) and pure Forsterite (red) compositions, respectively. The clearest distinction between the two cases is the peaking at shorter wavelengths for the Enstatite case than in the Forsterite case. Bottom: residuals of the same two spectra after subtraction of a pure Enstatite model, leaving significantly smaller residuals in the true Enstatite case than the Forsterite case, relative to the error bars.

7 Conclusions

As high-contrast direct imaging of exoplanets steadily improves with enhanced capabilities and ambitious plans for future missions, it is becoming increasingly important to consider the observability of exocomets in a direct imaging context – both as a scientific opportunity, and as a possible source of planetary false positives. In this study, we have examined several aspects of exocomet direct imaging and spectroscopy with LIFE. We have found that in systems with extreme cometary activities, such as β Pic, there may be as much as hundreds of detectable exocomets at any given time, which would render the observations confusion-limited, although the resulting signal would be largely dominated by the few brightest comets in the distribution. Somewhat older, more nearby, and slightly less active debris disk systems such as ϵ Eri or Fomalhaut may be the ideal targets for maximizing actual exocomet yield and S/N.

Orbital characterization of detected exocomets is difficult due to the small fraction of orbital arc that is covered during its active phase at periastron passage, but setting a lower limit on the eccentricity is feasible, which is useful, not least for avoiding false identification of detected point sources in the target system. Some exocomets might also be expected to exhibit strong silicate emission features, which would be detectable within reasonable amounts of integration time, and could allow for the first mineralogical studies of individual exocometary bodies.

Acknowledgements

M.J. gratefully acknowledges funding from the Knut and Alice Wallenberg Foundation. The authors thank Laetitia Rodet for useful discussions and helping to verify the orbital analysis. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work has made use of the Small Body Database managed by the Solar System Dynamics group at the Jet Propulsion Laboratory: https://ssd.jpl.nasa.gov/. The study has also made use of the CDS, NASA ADS, and NASA exoplanet archive services, as well as packages matplotlib (Hunter 2007) and numpy (Harris et al. 2020).

References

All Figures

	Fig. 1 Schematic representation of the LIFE nulling configuration, where the four unit telescopes U1 through U4 are plotted as circles. The configuration is rectangular with a short axis b and a long axis qb. During observations, the configuration rotates in the plane of the targeted sky patch.
In the text

Fig. 2 Extent of the β Pic planetary system relative to the ~1 arcsec FOV of LIFE at 10 μm (black dashed circle), as seen from Earth. The orbits of the planets (Nowak et al. 2020) are plotted in blue for the outer planet ‘b’ and red for the inner planet ‘c’. The inner green dashed circle marks the region inside of 1.6 au which is of primary interest for exocomet studies, since coma activity may be very limited outside this range. The edge-on disk in the system extends 200 times beyond the LIFE FOV.

In the text

	Fig. 3 Histogram of the inclinations of <100 au comets in the Solar system. While low inclinations are generally statistically favored, the spread of the distribution is quite wide. Thus, if the comets in a system like β Pic share a similar orbital distribution, then only a small fraction of comets will transit as seen from any given orientation.
In the text

	Fig. 4 S/N values (color scale) as function of temperature and effective coma radius for simulated exocomets in the β Pic system, in 10 h of observing time. Upper left: a case with low zodiacal contamination (1 zodi). Upper right: a higher level of 10 zodi. Lower left: the baseline value of 100 zodi. Lower right: a very high level of 1000 zodi.
In the text

	Fig. 5 Same as Fig. 4, but varying the sizes of the unit telescopes from the baseline value (2 m diameter). The observing time in each case is 10 h. Left: 3.5 m apertures. Right: 1.0 m apertures.
In the text

Fig. 6 Transmission for the nulling interferometric mode used in this study. Upper: transmission pattern in polar coordinates at 10 μm of the nulling interferometric configuration used in this study. A comet at a separation of 40 mas would project along the black dashed line as the telescope array rotates. Lower: a cross section of the transmission pattern for a separation of 40 mas (i.e., along the dashed line in the upper panel). A comet residing at a position angle of 0 deg would map the transmission pattern of the black solid line. A comet at a different PA would result in a cyclically phase-shifted pattern of otherwise the same shape. See Dannert et al. (2022) for further details.

In the text

Fig. 7 Simulations of a crowded exocometary field. The three rows correspond to three different dates. Upper: day 0. Middle: day 7. Lower: day 30. The left column contains maps of the spatial distribution of the exocomets. Each circle is an exocomet, where the radius of the circle is directly proportional to the effective radius of the cometary coma. Each comet has its own color, and the dashed lines of each color corresponds to the orbital track of that particular comet. The right column shows the interferometric signal as a black thick line. For illustrative purposes, the contribution of each individual comet is plotted as a thin line that is color-matched with the corresponding symbol in the left column. The individual colored components cannot be measured directly; they can only be deduced from the total signal and its evolution in time, wavelength, and (if relevant) along different probed interferometric baselines. See the text for discussion.

In the text

	Fig. 8 Corner plot for the most important orbital parameters from the orbit fitting test described in the text, with a 95% cut on the axes for the 104 orbital solutions. The true orbital parameters are not well recovered, due in large part to the short baseline probed with the simulated observations.
In the text

Fig. 9 20 example fitted orbits (colored lines), randomly sampled from the full population of fitted orbits. The orbits are generally significantly smaller and less eccentric than the input orbit, underlining the fact that the probed baseline is too short to determine exact parameters, although meaningful limits can be obtained (see text for details). Left: the orbital fits in two dimensions. The red star symbol marks the location of the star, and the yellow points mark the locations of the comet at three measured epochs. Upper right: separation as function of time. Upper left: position angle as function of time. Purple dots mark the location of the comet at the respective epochs.

In the text

	Fig. 10 Simulated S/N values per wavelength bin using LIFEsim, for four cases of possible cometary emission with different silicate species: Enstatite (blue), Forsterite (red), Olivine (green) and Pyroxene (purple). Solid lines denote 10 h observations, while dashed lines represent 100 h observations.
In the text

Fig. 11 Example of exocomet spectroscopy with LIFE. Top: simulated LIFE spectra for two hypothetical exocomet cases, with noise levels set based on LIFEsim calculations for 10 h observations. The two cases correspond to pure Enstatite (blue) and pure Forsterite (red) compositions, respectively. The clearest distinction between the two cases is the peaking at shorter wavelengths for the Enstatite case than in the Forsterite case. Bottom: residuals of the same two spectra after subtraction of a pure Enstatite model, leaving significantly smaller residuals in the true Enstatite case than the Forsterite case, relative to the error bars.

In the text