A PR drag origin for the Fomalhaut disc’s pervasive inner dust: constraints on collisional strengths, icy composition, and embedded planets

ABSTRACT

Recent JWST observations of the Fomalhaut debris disc have revealed a significant abundance of dust interior to the outer planetesimal belt, raising questions about its origin and maintenance. In this study, we apply an analytical model to the Fomalhaut system, which simulates the dust distribution interior to a planetesimal belt, as collisional fragments across a range of sizes are dragged inward under Poynting–Robertson (PR) drag. We generate spectral energy distributions and synthetic JWST/MIRI images of the model discs, and perform an extensive grid search for particle parameters – pertaining to composition and collisional strength – that best match the observations. We find that a sound fit can be found for particle properties that involve a substantial water ice component, around 50–80 per cent by total volume, and a catastrophic disruption threshold, |$Q_\mathrm{D}^\star$|⁠, at a particle size of |$D\!\approx \!30\, \mu {\rm m}$| of (2–4)|$\, \times 10^6\, {\rm erg\, g}^{-1}$|⁠. Based on the expected dynamical depletion of migrating dust by an intervening planet, we discount planets with masses |$>1\, M_{\rm Saturn}$| beyond |$\sim$|50 au in the extended disc, though a planet shepherding the inner edge of the outer belt of up to |$\sim 2\, M_{\rm Saturn}$| is reconcilable with the PR-drag-maintained disc scenario, contingent upon higher collisional strengths. These results indicate that PR drag transport from the outer belt alone can account for the high interior dust contents seen in the Fomalhaut system, which may thus constitute a common phenomenon in other belt-bearing systems. This establishes a framework for interpreting mid-planetary system dust around other stars, with our results for Fomalhaut providing a valuable calibration of the models.

1 INTRODUCTION

The Fomalhaut system is known as one of the archetype debris discs owing to its luminosity and proximity (7.7 pc). Following the first detection of the disc’s substantial infrared excess emission by the Infrared Astronomical Satellite (IRAS) (Gillett 1986), Fomalhaut’s disc garnered significant attention due to its apparent multicomponent structure. Observations across the electromagnetic spectrum – from scattered light captured by Hubble (Kalas, Graham & Clampin 2005; Kalas et al. 2013), to the infrared with Spitzer and Herschel (Stapelfeldt et al. 2004; Acke et al. 2012), and at submillimetre/millimetre wavelengths with the James Clerk Maxwell Telescope (JCMT) and the Atacama Large Millimeter Array (ALMA) (Holland et al. 2003; MacGregor et al. 2017), amongst others – indicated an intricate system, including an outer belt at |$\sim$|140 au, as well as an inner, warm (⁠|$\sim$|150 K) disc component localized at 8–15 au, as inferred from the spectral energy distribution (SED) of the total emission (Su et al. 2016). This configuration seemed to resemble the planetesimal reservoirs of the Solar system, with a cold Kuiper belt and a warmer, relatively narrow asteroid belt.

However, recently, the JWST changed this conception of the Fomalhaut system, which for the first time resolved the distribution of dust interior to the outer belt (Gáspár et al. 2023). These observations reveal a prominent drawn-out inner disc from around 10 to 90 au, in addition to an ‘intermediate belt’, an eccentric (⁠|$e\!\approx \!0.3$|⁠) ring-shaped brightness feature stretching from roughly 60 to 110 au, which together fill up much of the space interior to the outer belt.

Such a continuous distribution of dust permeating the space interior to the outer planetesimal belt raises questions about its origin, in particular, whether it can be maintained by dust produced in the outer belt and being dragged inward by Poynting–Robertson (PR) drag, or whether other scenarios, such as dust delivery by comets (e.g. Faramaz et al. 2017; Marino et al. 2017), are necessary. This also pertains to the intermediate belt, which could conceivably be caused by a planet trapping migrating dust in orbital resonances, or by in situ dust production in a spatially confined, collisionally grinding planetesimal reservoir.

Wyatt (2005) provided an idealized approach to gauge the abundance of PR-drag-supplied dust interior to a cold belt – assuming discs to be composed of only the smallest, barely bound grains – which indicated the mechanism’s ineffectiveness to produce signals detectable with the contemporary instrumentation. Then, in light of more sensitive interferometric observing methods becoming available, Kennedy & Piette (2015) concluded that PR drag transport could no longer be considered insignificant, using a similar approach. Recent advancements of these analytical methods to incorporate a distribution of particle sizes (Rigley & Wyatt 2020), now allow a more rigorous reassessment of the PR drag scenario. In this study, we adapt and apply the model developed by Rigley & Wyatt (2020) to the Fomalhaut disc, to evaluate whether a pure PR drag transport scenario can account for the newly observed characteristics of the extended disc. Specifically, we simulate discs with a broad range of particle material parameters – facilitated by the fast analytical model – seeking outcomes consistent with both the disc’s SED, and the radial brightness profiles extracted from JWST images. In doing so, we incorporate more definitive properties of the outer belt derived from ALMA imaging as fixed model inputs. Given the simplicity of the model, the aim of this study is to assess the general validity of a PR-drag-fed extended disc, and does not target the reproduction of more intricate features, such as the eccentricity of the outer belt or the intermediate ring, or the enigmatic near-infrared excesses detected in interferometric observations attributed to a population of close-in hot dust (Lebreton et al. 2013; Mennesson et al. 2013; Su et al. 2016).

Section 2 provides an overview of the observational data. Section 3 introduces the PR drag disc model and details our method to generate and evaluate synthetic observations, across a grid of model parameters. Section 4 analyses the parameter dependence of suitable model outcomes. Section 5 presents the best-fit model outcome, and Section 6 compares spectra simulated from it to actual Spitzer spectra. Section 7 discusses our parameter constraints and the implications of the results. Section 8 summarizes our findings.

2 OBSERVATIONS

In this section, we examine the observational data utilized to fit our analytical model of the Fomalhaut debris disc. The data include overall excess flux measurements at different wavelengths, representing the disc’s SED, which are summarized in Table 1 and Fig. 1, and the radial brightness distribution captured by JWST’s Mid-Infrared Instrument (MIRI). The disc’s SED encodes comprehensive information on the disc’s overall thermal emission, predominantly reflecting the characteristics of the larger, colder particles in the outer belt, which emit in the far-infrared and (sub-)millimetre regime and are less prone to PR drag migration. Conversely, radial brightness profiles in the mid-infrared provide detailed spatial information about the smaller and warmer dust particles, which have a distribution more significantly influenced by PR drag. Together, these data sets encapsulate the coupling of particle sizes and optical properties across the entire debris disc, which we aim to reproduce simultaneously. ALMA observations are additionally used to inform our fixed model parameters concerning the location of the dust-producing planetesimals, as discussed in Section 3.3.

We extract the mid-infrared radial surface brightness profiles from the high-resolution images of the Fomalhaut disc recently obtained by JWST/MIRI at three different wavelengths with the 15.5 and 23 |$\mu {\rm m}$| coronagraphic filters (F1550C and F2300C) as well as the 25.5 |$\mu$|m filter (F2550W) (Gáspár et al. 2023). These images are reproduced in Fig. 2 from the published data. Following Gáspár et al. (2023), we deproject the reduced images using a position angle (PA) of 336.28|$^\circ$| and an inclination of 67.52|$^\circ$|⁠, and then find the median surface brightness at bins of astrocentric distance. Due to the unavailability of detailed error data, we use the standard deviation of the azimuthal variation at each given radius as a proxy for the uncertainty. Given an apparent inconsistency of the F2300C reduced fluxes with previously taken spectra, Gáspár et al. (2023) suspect an issue with the F2300C calibration pipeline (see their supplementary section 1.2). Therefore, while we display the F2300C data for comparison with our models, we do not use this data for model fitting.

3 DISC MODELLING

3.1 PR drag disc model

We employ the analytical model developed by Rigley & Wyatt (2020) to simulate the two-dimensional dust distribution – across radial extent and particle size – in a disc formed by a collisionally grinding planetesimal belt. The model combines previous analytical approaches to determine (1) the size distribution arising within the radial confines of the planetesimal belt under collisions and PR drag loss (Wyatt, Clarke & Booth 2011), as well as (2) the radially dependent size distribution interior to the belt under further collisional and drag-induced evolution (Wyatt 2005). The model reasonably approximates the results of numerical kinetic models (e.g. van Lieshout et al. 2014) whilst being significantly faster, thus enabling exploration of large parameter spaces. For more details about the model, the reader is referred to Rigley & Wyatt (2020).

We compute the dust grains’ material-dependent optical properties, including their |$\beta$|-factors – the ratio of radiation pressure to stellar gravity that governs a particle’s PR drag evolution – and thermal emission in a manner similar to Wyatt & Dent (2002). Following Li & Greenberg (1997), the method treats dust as aggregates of core-mantle grains, applying Maxwell–Garnett effective medium theory to derive the optical constants for the composite material. We then determine the grains’ optical efficiency coefficients using Mie theory, Rayleigh–Gans theory, or geometric optics, in the respective wavelength regimes. Finally, the absorption coefficients define the dust grain temperatures, which are computed by balancing the absorbed stellar radiation with the emitted thermal radiation.

The PR drag disc model has previously been applied to target systems of the HOSTS survey to assess if sufficient quantities of warm inner dust – required to explain a number of exozodi detections by the LBTI nulling interferometer (Ertel et al. 2020) – can be sustained by material being dragged inwards from an outer belt, or if alternative dust transport scenarios are necessary (e.g. via comets). In five of nine systems with known outer belts, this ‘PR-drag-transport-only’ model could reasonably reproduce the detected exozodi levels (Rigley & Wyatt 2020). Similarly, here we apply this model to the Fomalhaut system, while leveraging the more extensive observational data recently gathered by JWST/MIRI.

Both components of this model, that is, our implementations of the PR drag disc model from Rigley & Wyatt (2020) and of the method to find the dust optical properties and temperatures, are publicly available as individual python packages (see Section Data Availability).

3.2 Synthesizing MIRI images

For an adequate comparison of our model outcomes with the MIRI data, we need to generate synthetic MIRI images from the simulated dust spatial distribution and optical properties. Whereas Rigley & Wyatt (2020) focussed on the total emission from the disc (which we also employ in form of the SED), here we obtain the spatially resolved emission by first determining the radial profile of the disc’s face-on surface brightness at a specific wavelength as

where |$Q_\mathrm{abs}(\lambda , D)$| are the grains’ absorption efficiencies, |$B_\nu$| their spectral radiance, |$T(D,r)$| their temperatures, and |$\tau (D,r) \mathrm{d}D$| is the geometrical depth in particles of size |$D\!\rightarrow \!D\!+\!\mathrm{d}D$| at radius r.

For a given face-on radial brightness profile, disc inclination, and scale height aspect ratio, we produce the corresponding astrophysical scene (i.e. how the inclined disc would appear to an ideal observer, free from observational distortions) using the rave package (Han, Wyatt & Matrà 2022). This method draws sample points that are vertically Gaussian, azimuthally uniform, and radially distributed according to the radial profile, before binning them into a pixel grid projected at the appropriate orientation in space. We adopt a pixel scale of 0.26 arcsec per pixel (corresponding to 2 au at the distance of Fomalhaut), slightly more than twice the native scale of MIRI (0.11 arcsec per pixel). This coarser sampling serves to limit the computational costs imposed by the point spread function (PSF) convolution step for the coronagraphic filters, where the PSF varies across the field of view.

To produce synthetic MIRI images, we convolve the astrophysical scenes with the PSF of the respective MIRI filter, obtained using the webbpsf package (v. 1.3.0; renamed STPSF as of v. 2.0, Perrin et al. 2012). For the coronagraphic filters, where the PSF varies across the image, we pre-generate a data cube storing the PSF for each individual pixel, requiring 1.3 and 2.0 GB of memory for F1550C and F2300C, respectively.

We consider only the emission at the central wavelength of each filter. Extending our pipeline to multiple wavelengths per filter, which would also require caching the spatially variant PSFs for each wavelength, would be computationally prohibitive for our purpose. Given the bandwidth of the filters, this simplification might introduce inaccuracies. However, we found that, at least for the unconvolved face-on radial brightness profiles, the difference between central-wavelength and band-integrated profiles (with wavelength-dependent throughputs for the latter obtained from Pandeia, the JWST exposure time calculator system; Pontoppidan et al. 2016) is negligible for our purpose, reaching at most a few per cent. To replicate the saturation at the stellar core in the observations, we mask the pixels within 1.2 arcsec in the F2550W image. No noise is added to the images. We also do not replicate the differential imaging in the observations. Only one synthetic image is generated for each filter with the detector’s vertical axis aligned with the disc’s minor axis, comparable with the observations.

Since we model only the axisymmetric radial distribution of dust, we do not compare the observed and synthetic images directly but instead the radial brightness profiles extracted from the deprojected images, thereby ignoring any azimuthal features observed in the actual disc, in particular the slight eccentricity of the outer planetesimal belt, as well as the eccentric ‘intermediate ring’ feature superimposed on the disc (Gáspár et al. 2023). For deprojection, we again use the disc PA and inclination found by Gáspár et al. (2023), so that we are comparing profiles extracted from the observations and simulated model images in exactly the same way.

3.3 Application to Fomalhaut

To apply the PR drag model to the Fomalhaut disc, it is necessary to set appropriate input parameters for the central star, the planetesimal belt, and the particle characteristics.

For the central star, we adopt a mass |$M_{\star }\!=\!1.92\, {\rm M}_{\odot }$| (Mamajek 2012), and a luminosity |$L_{\star }\!=16.6\, {\mathrm{ L}_{\odot }}$| and spectrum (shown in Fig. 1, |$T_\mathrm{eff}\!=\!8523\, {\rm K}$|⁠) modelled with PHOENIX (Allard et al. 2012), as found via the SED fitting method by Yelverton et al. (2019, 2020) (spectrum provided by Kennedy, personal communication).

For the planetesimal belt, the model requires inputs for its geometric characteristics – namely, its radial extent and opening angle – as well as its dust mass, which is defined as the total mass up to a particle size of 1 cm, the largest size included in the model. The radial dimensions of the belt are informed by millimetre imaging data from ALMA observations (MacGregor et al. 2017), which delineate the location of the planetesimals. These data suggest a moderately eccentric belt (⁠|$e\!=\!0.12 \pm 0.01$|⁠) with an inner edge at pericentre at |$119.9 \pm 0.8$| au and an outer edge at apocentre at |$152.6 \pm 1.0$| au. Given that our model is axisymmetric, we approximate the observed belt as a wider, circular ring with an inner and outer edge fixed at 119.9 and 152.6 au, respectively, thereby preserving the radial extent of the actual belt.

The opening angle of the disc may also be estimated from the ALMA imaging data. We employ a parametric belt model to fit the full ALMA map by MacGregor et al. (2017), ignoring the belt’s eccentricity and instead assuming an axisymmetric ring centred on the geometric centre of the projected disc emission. Further assuming a radially and vertically Gaussian ring with peak intensity, peak radius, radial width, inclination, and scale height aspect ratio as free parameters we find a vertical half-opening angle of |$1.4\pm 0.1^\circ$| under this simplified model. This is comparable to the half-opening angle of |$1.0 \pm 0.25\circ$| found by Boley et al. (2012) by fitting an axisymmetric ring to a partial ALMA map the disc, although a vertically exponential profile was assumed in their model rather than a Gaussian profile. On the other hand, Kennedy (2020) find a half-opening angle of 1.7|$^\circ$| derived from fitting particle orbit distributions to the full ALMA map.¹ For simplicity, here we adopt a value of 1.5|$^\circ$|⁠. We discuss how variations of this parameter affect our model outcomes in Section 7.1.

The original model assumed the collisional velocity to be dominated by the particles’ vertical motion, with |$v_\mathrm{rel}\!=\!v_\mathrm{k} I_\mathrm{max}$| (Rigley & Wyatt 2020, equation 5), where |$v_\mathrm{k}$| is the circular Keplerian velocity and |$I_\mathrm{max}$| is the particles’ maximum inclination in radian, corresponding to the semi-opening angle of the disc. To also account for an eccentricity distribution within the eccentric belt, we add a velocity component depending on the particles’ mean proper eccentricity; |$v_\mathrm{rel}\!=\!v_\mathrm{k} \sqrt{1.25 \, \smash{e_\mathrm{p}^{\scriptscriptstyle 2} \, +\, I_\mathrm{max}^{\scriptscriptstyle 2}}}$| (Lissauer & Stewart 1993). For the mean proper eccentricity, |$e_\mathrm{p}$|⁠, we adopt a value of 0.019, as found by Kennedy (2020).

The belt dust mass – unlike other parameters that may be varied freely or fixed from the start – is set individually for each run such that the resulting total flux from the belt at |$\lambda \!=\!1.3\, {\rm mm}$| matches the belt flux measured by ALMA at that wavelength. Specifically, we scale the mass to match the mean of the fluxes found by MacGregor et al. (2017) via fitting in the image plane (24.7 |$\pm$| 0.1 mJy), and by White et al. (2017) via fitting the visibility data (30.8 |$\pm$| 0.1 mJy). The mass is found by rerunning the model iteratively, until the difference between the modelled and measured flux is |$<$|1 per cent (usually achieved within three iterations).

We assume the PR-drag-maintained disc to extend down to 0.5 au. This is roughly comparable to the sublimation distance of refractory dust (0.2–0.3 au for silicate; Lebreton et al. 2013). The observables we simulate and compare in this study – mid-infrared up to sub-mm wavelengths – proved relatively insensitive to the exact location of the innermost edge of the disc, provided it is |$\lesssim \!1~\mathrm{ au}$|⁠.

3.4 Particle parameters

The particle parameters in our model pertain to the grain composition and structure, which determine their optical properties, as well as to the grains’ collisional strength, and are treated as free parameters. The dust grain model we employ (Li & Greenberg 1997) assumes that particles are porous aggregates of core-mantle grains with a silicate core and an organic refractory mantle, as well as water ice filling the voids. It is parametrized by three fractions: (1) The volume fraction of silicate in the combined volume of silicates and organic refractory material, |$q_\mathrm{sil}= v_\mathrm{Si}/(v_\mathrm{Si} + v_\mathrm{C})$|⁠. (2) The matrix porosity, p, defined as the volume fraction of voids within the total volume of the silicate/organic matrix. These voids may be empty or partially filled with water ice, as determined by (3) the |$q_\mathrm{ice}$| parameter, defined as the ratio of ice volume to total void volume. For each of these parameters, we test the full range (i.e. 0–1, or rather 0–0.99 in the case of porosity).

For the collisional strength, the model of Rigley & Wyatt (2020) adopts a power-law determining the critical specific energy |$Q_\mathrm{D}^\star$|⁠, that is, the kinetic energy threshold required for an impactor to catastrophically disrupt a target, following Benz & Asphaug (1999). This two-parameter prescription is of the form |$Q_\mathrm{D}^\star = Q_\mathrm{0}(D/\mathrm{cm}) ^ {b_\mathrm{s}}$|⁠, with a normalization parameter, |$Q_\mathrm{0}$|⁠, and a slope parameter, |$b_\mathrm{s}$|⁠. Conventionally assumed parameters for |$Q_\mathrm{D}^\star$| typically correspond to a critical specific energy for a 1 cm target (i.e. |$Q_\mathrm{0}$|⁠) within one order of magnitude from |$10^7$| erg g|$^{-1}$| [e.g. Krivov, Löhne & Sremčević (2006): |$3.02\times 10^6\, {\rm erg}\, {\rm g}^{-1}$|⁠; Grün et al. (1985): |$\sim 9\times 10^6\, {\rm erg\, g}^{-1}$|⁠;² Heng & Tremaine (2010): |$10^7\, {\rm erg\, g}^{-1}$|⁠; Löhne, Krivov & Rodmann (2008): |$2.45 \times 10^7\, {\rm erg\, g}^{-1}$|⁠). However, a recent study of the collisional lifetimes of Solar system meteoroids finds |$\sim \!\!3$| orders of magnitude higher collisional strengths than conventionally assumed are necessary to fit observed meteor dynamics (Pokorný et al. 2024). To incorporate this large span of estimated strengths, we choose a dynamic range of a factor of |$10^4$| for |$Q_\mathrm{0}$|⁠. For |$b_\mathrm{s}$|⁠, we explore a range from |$-0.9$| to 0.45, comfortably accommodating the span of typically assumed values of around |$-0.3\,\,{\rm to}\,\,0$| (e.g. Krivov et al. 2006; Löhne et al. 2008; Heng & Tremaine 2010). We distinguish two fitting scenarios; Scenario A, across the entire range of tested values for |$b_\mathrm{s}$|⁠, and Scenario B, where only the more conventional range of |$b_\mathrm{s}\!\le \!0$| is considered.

With 9–14 values per parameter tested (see Table 2), this amounts to 131 670 parameter combinations, including |$1463$| unique compositions after excluding redundant settings (i.e. variations of |$q_\mathrm{ice}$| for |$p\! =\!0$|⁠). We precompute the optical properties and radially dependent temperatures for the |$1463$| materials across 45 particle sizes logarithmically spaced between 4 |$\mu {\rm m}$| and 1 cm, taking around 1 h on a server with a 40-core processor. We then run the PR drag model, including the generation of the three synthetic MIRI images, with the full parameter grid, taking around 6 h on the same machine. The per-core computing time of each individual model is around 6 s, with most of this time dedicated to image generation and convolution.

4 MODEL FITTING

4.1 Goodness-of-fit metric

To quantify how well the model reproduces the observations under a given set of parameters, we construct goodness-of-fit functions |$\chi ^2$|⁠, where the best fit has the minimum |$\chi ^2$|⁠. To simultaneously match both our observational data sets reasonably well, we will introduce |$\widehat {\chi ^2}$| as a combination of |$\chi _\mathrm{SED}^2$| and |$\chi _\mathrm{rad}^2$|⁠, which are the respective goodness-of-fit metrics for the SED and for the MIRI-retrieved radial brightness profiles.

We start from the general form of the |$\chi ^2$| statistic used to quantify the goodness-of-fit between observed data and model predictions, which is given by

where |$O_i$| are the observed values, |$E_i$| are the expected (model) values, and |$\sigma _i$| are the uncertainties in the observed values. To appropriately handle the logarithmic nature of our data (both in values, i.e. flux, and in sampling points, i.e. wavelength or, respectively, radial distance from the star), we first transform the observed and model values, as well as the uncertainties, to logarithmic scales, such that we replace:

In addition, given that the observational data points may not be evenly spaced, we incorporate a weighting mechanism based on data point density in logarithmic space. We calculate weights as the inverse of the distance between preceding and following data points:

where |$x_i$| are the sampling points of the observed values. This ensures that closely spaced (in log-scale) data points do not excessively dominate the fit. We obtain

where |$\sigma _{\log _{10}(O_i)} = \sigma _i/[O_i \ln (10)]$|⁠.

In the case of the two radial profiles retrieved with the F1550C and F2550W filters (note that we omit the F2300C profile, see Section 2), we compute one |$\chi ^2$| value for each and add them up to form |$\chi _\mathrm{rad}^2$|⁠:

Lastly, to ensure that both data sets equitably contribute to the combined goodness-of-fit metric, we normalize the ‘unscaled’ |$\chi ^2$| metrics by dividing them by their respective best-fit values across all parameter sets:

This normalization is necessary for two reasons. First, to factor in the varying numbers of sampling points in both observational data sets – compounded by the fact that multiple radial profiles were retrieved – and, secondly, to account for the different magnitudes of the uncertainties accompanying the data sets. The stated uncertainties for the SED are roughly an order of magnitude smaller than the uncertainties of the radial profiles, for which we use the azimuthal variation at a given radius (see Section 2). This normalization essentially equalizes the uncertainty scales of the two data sets, while upholding the relative weighting effect the uncertainties have among data points within each data set.

Finally, we define the combined goodness-of-fit metric as the mean of the normalized |$\chi ^2$| values for each data set:

This metric approaches one for the best achievable outcome, although it can only reach one if the best fits for both data sets are achieved simultaneously. Due to the normalization, this metric does not signify absolute goodness of fit but rather represents a relative measure of model performance that allows us to easily compare different model configurations across the entire parameter grid.

4.2 Exploring the parameter space for good fits

To explore which combinations of parameters optimize the fit between our model predictions and the observational data, and to identify interactions among these parameters, we assess the |$\widehat {\chi ^2}$| distribution across our models using a corner plot (Fig. 3). The pseudocolour diagrams in the lower triangle of the plot illustrate the interplay of any two parameters by showing the minimum achieved |$\widehat {\chi ^2}$| for each pairing while marginalizing over the remaining parameters. This approach helps us identify regions within the parameter space where the model fit is most favourable. The plots on the diagonal highlight how |$\widehat {\chi ^2}$| varies as each parameter is adjusted independently, indicating the sensitivity of the model fit to changes in individual parameters. By utilizing a |$\widehat {\chi ^2}$| cutoff, we may establish confidence intervals for these parameters. Given the lack of an absolute goodness-of-fit measure, we empirically set a cutoff at |$\widehat {\chi ^2} <\!1.5$|⁠, which we found to yield acceptable fits. Table 3 presents the best-fit parameters (those achieving the overall minimum |$\widehat {\chi ^2}$| value) along with the corresponding uncertainties, defined by the range of parameter values that fall within our |$\widehat {\chi ^2}$| cutoff.

First looking into the dependence of the goodness of fit on individual parameters, we note that suitable parameters suggest fairly icy particles – characterized by high porosity and a high water ice filling factor – with a weakly constrained silicate-organics ratio. The conventional range for |$Q_\mathrm{0}$| fits well, while the higher collisional strengths found for Solar system grains by Pokorný et al. (2024) are rejected. The commonly assumed slightly negative values for |$b_\mathrm{s}$| do allow for a suitable fit, yet the tendency for |$b_\mathrm{s}$| to favour positive values is noteworthy. To assess how a more conventional |$Q_\mathrm{D}^\star$| power-law slope would change our parameter constraints, we can also explore Scenario B where we impose |$b_\mathrm{s}{} \le 0$|⁠. This constraint worsens the fit only marginally, yet tightens the requirements for a high |$q_\mathrm{ice}$| and low |$Q_\mathrm{0}$| (see Table 3).

4.2.1 Parameter correlations

The marginalized |$\widehat {\chi ^2}$| distributions across two parameters, seen in the lower triangle in Fig. 3, let us identify parameter correlations, which manifest as skewed, elongated regions of suitable fits. While most parameter combinations show no evidence for correlation, a pronounced correlation between the collisional strength parameters |$Q_\mathrm{0}$| and |$b_\mathrm{s}$| is evident – apparently, |$b_\mathrm{s}\!\propto \!\log {Q_\mathrm{0}}$|⁠. We find that the region of suitable fits is described by a narrow corridor along a straight line that is approximately described by

Along this line, |$Q_\mathrm{D}^\star$| can thus be expressed as depending solely on |$b_\mathrm{s}$|⁠:

Here, given the exponential growth of the terms |$10^{2.5 b_\mathrm{s}+ 6.625}$| and |$(D/{\rm cm})^{b_\mathrm{s}}$|⁠, one may suspect that there exists a specific value of D, where the growth rates of these terms effectively counterbalance each other, resulting in a constant product |$Q_\mathrm{D}^\star$| across different values of |$b_\mathrm{s}$|⁠. To determine D at which |$Q_\mathrm{D}^\star$| is invariant with respect to variations in |$b_\mathrm{s}$|⁠, we differentiate |$\log (Q_\mathrm{D}^\star)$| with respect to |$b_\mathrm{s}$| and set it to zero, yielding |$D\!=\!10^{-2.5}\, {\rm cm}\!\approx \!32\, \mu {\rm m}$|⁠. At this value, |$Q_\mathrm{D}^\star$| is invariant for any combination of |$Q_\mathrm{0}$| and |$b_\mathrm{s}$| that satisfies equation (9), which is where the most suitable fits are found. In other words, a favourable model outcome is less determined by |$Q_\mathrm{0}$| or |$b_\mathrm{s}$| per se, but rather solely by |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$|⁠, which is optimal at around |$\sim 4\times 10^6\, {\rm erg}\, {\rm g}^{-1}$|⁠. The marginalized distribution of |$\widehat {\chi ^2}$| over this derived parameter, |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$|⁠, is shown in Fig. 4 (top left). The similar distributions for Scenarios A and B confirm that, for good fits, |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$| is rather independent of |$b_\mathrm{s}$|⁠. The fact that the size regime at which we obtain a constraint on |$Q_\mathrm{D}^\star$| is just above the blow-out size around Fomalhaut [|$D_{\mathrm{bl}}\!:\!\beta (D_{\mathrm{bl}})\!=\!0.5, D_{\mathrm{bl}}\!\approx \!13\, \mu {\rm m}$|⁠, for the best-fit material] is consistent with the expectation that, in a collisional cascade, the geometrical depth is dominated by the smallest particles that can remain bound, which therefore also contribute most to the emission we are fitting.

Another noteworthy parameter correlation exists between the porosity p and the ice filling fraction |$q_\mathrm{ice}$|⁠. There appears to be a region, within which suitable fits may be found that is roughly characterized by |$p\!\propto \!1/q_\mathrm{ice}$|⁠. This indicates a constraint on the product of p and |$q_\mathrm{ice}$|⁠, which is the total ice volume fraction of the particle material, |$\phi _\mathrm{ice}= v_\mathrm{ice} \, /\, v_\mathrm{tot} = p\cdot q_\mathrm{ice}$|⁠. Plotting the marginalized |$\widehat {\chi ^2}$| distribution over |$\phi _\mathrm{ice}$|⁠, as given in Fig. 4 (top right), reveals an extended minimum spanning roughly |$0.5\!<\!\phi _\mathrm{ice}\!<\!0.9$|⁠. Considering only models with |$b_\mathrm{s}\!\le \!0$| (Scn. B), this constraint tightens to |$0.7<\!\phi _\mathrm{ice}\!<\!0.9$|⁠.

The constraints on these secondary model parameters, |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$| and |$\phi _\mathrm{ice}$|⁠, as well as the overall grain porosity |$\phi _\mathrm{vac}= v_\mathrm{vac} \, /\, v_\mathrm{tot} = p\cdot (1-q_\mathrm{ice})$|⁠, are also summarized in Table 3.

4.2.2 Constraints on lifetimes and zodi level

We may also analyse how certain derived quantities are constrained by the fitting process. The respective marginalized |$\widehat {\chi ^2}$| distributions are depicted in Fig. 4 (bottom).

For one, the constraint we obtained on |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$| suggests that the fitting actually constrains the collisional lifetimes of particles of that size – or more specifically, the ratio of their PR drag lifetimes (the time it takes particles to migrate from the belt to the star) to their collisional lifetimes (were they to remain in the belt), |$\eta _0$|⁠. This ratio, which depends on |$Q_\mathrm{D}^\star$| among other factors, is a crucial intermediate quantity within the model that influences the determination of the geometrical depth profile. We find that the |$\eta _{0,\, {32}{\mu m}}$| of good fits is on the order of several tens, with Scenario B favouring a slightly higher |$\eta _{0,\, {32}{\mu m}}$| than Scenario A.

Another interesting parameter we can derive from the models is the zodi level. This quantity is defined as the disc’s surface density at the Earth-equivalent insolation distance – calculated as |$r_\mathrm{\, EEID} = \sqrt{L_\star /L_\odot }\, \mathrm{au}$|⁠, approximately 4.1 au for Fomalhaut – relative to the surface density of the Solar system’s zodiacal cloud at 1 au (Kennedy et al. 2015).³ As seen in Fig. 4 (bottom-right), well-fitting models have a zodi level of around |$60$|⁠, meaning that, at the radius where a planet would receive the same energy as Earth does from the Sun, the surface brightness is approximately |$60$| times that of the zodiacal cloud at the orbit of Earth. There is no significant difference between Scenarios A and B concerning this parameter. However, we do not expect the disc at 4 au to strongly influence the observables we are fitting, making this prediction more of an extrapolative indicator of disc density rather than a direct observational constraint.

5 THE BEST-FITTING MODEL

In this section, we inspect the best-fit model outcome (Scenario A) more closely. We first turn to the model’s ability to reproduce the spatial distribution of emission as observed by MIRI. Fig. 5 shows the synthetic images created from the best-fit model, next to unconvolved images of the inclined disc, that is, the ‘pristine’ astrophysical scenes generated with rave from the modelled face-on radial surface brightness profiles, before the PSF of the respective filter is applied. In the F2300C and F2550W images, we can easily discern the prominent outer belt, while we note that a substantial amount of emission stems from the region within the belt. In the case of the F1550C filter, only the hotter, innermost region of the disc is visible. Comparing before/after PSF convolution, indicates a marked difference in the 15.5 |$\mu$|m band, where the original concentrated source is severely spread out to form a notable diffraction pattern, as previously described by Boccaletti et al. (2024).

For the comparison with the MIRI images we resort to the radial brightness profiles, shown in Fig. 6, which we obtain after deprojecting the synthetic images in the same way as we do with the MIRI images. We first note that the profiles derived from the synthetic images (dashed lines) differ substantially from the pristine modelled profiles (dotted line), which correspond to an ideal, face-on observer (i.e. without projection of the disc, PSF convolution, and subsequent deprojection). In particular the steeper 15.5 |$\mu$|m profile is considerably flattened, again illustrating the blurring effect of the PSF already noticed in Fig. 5. In the cases of the coronagraphs, this blurring is confounded by the drop-off in transmission when approaching their inner working angles (IWAs),⁴ which correspond to radial distances of roughly 4 and 17 au on the disc major axis, and 10 and 43 au on the disc minor axis for the F1550C and F2300C filters, respectively (⁠|$\mathrm{IWA}_{F1550C}\!=\!0.49\,\mathrm{ arcsec}$| and |$\mathrm{IWA}_{F2300C}\!=\!2.16\,\mathrm{ arcsec}$|⁠; Boccaletti et al. 2015). These distances are also indicated in Fig. 6. Note that, while its IWA might be smaller, the transmission drop-off of the 4QPM coronagraph, F1550C, extends relatively further out, owing to its tapered transmission profile (with |$\sim$|66 per cent transmission at |$2\mathrm{IWA}_{F1550C}$|⁠; Boccaletti et al. 2015, fig. 6 therein).

Comparing the synthetic image profiles (dashed lines) to those retrieved from actual MIRI images (solid lines), we find considerable agreement. Significant deviations only appear exterior to the planetesimal belt (⁠|$>$|168 au), where the model brightness drops off more steeply than was observed. This is explained by the fact that our simplified model does not include a representation of the disc halo, that is, the ensemble of small dust near and below the blow-out size that is pushed by radiation pressure beyond the outer belt. Nevertheless, due to the PSF blurring effect combined with the inclination of the disc, considerable brightness even beyond the belt – where no dust exists in the model – is present in our synthetic images and derived radial profiles, forming a faint pseudo-halo. Conversely, we can infer that part of the halo observed by MIRI is PSF-blurred emission from the planetesimal belt, rather than actual halo dust.

Note that the model was only fitted to the profiles retrieved from the F1550C and F2550W images, given the issues with the F2300C calibration pipeline suspected by Gáspár et al. (2023). Nevertheless, we have plotted here the radial profile retrieved from the F2300C image, however, only after reversing the additional, ‘corrective’ scaling factor of |$3.14$| that Gáspár et al. (2023) had applied in order to align F2300C and F2550W fluxes. The consistent offset of a factor |$\sim 3$| between the resolved F2300C and F2550W fluxes in their original reduction (apparent in the gap between the solid lines in Fig. 6) has led Gáspár et al. (2023) to surmise (quite reasonably) an issue in the calibration pipeline, given that Spitzer spectrographic observation of the inner disc had revealed a nearly flat spectrum at those wavelengths, with only a few percentage points increase from 23 to 25.5 |$\mu$|m. (We will discuss the Spitzer spectral data in more detail in Section 6.) However, we see here that, while the underlying surface brightness (dotted lines) at those wavelengths is indeed very similar at shorter radial distances, the corresponding brightness profiles retrieved from the synthetic images (dashed lines) differ consistently by a factor of 2–3, owing in particular to the aforementioned transmission drop-off at short working angles. We therefore conclude that there is no evidence for an inherent issue with the F2300C reduction pipeline.

On the whole, we conclude that the modelled brightness distribution matches the one observed by MIRI remarkably well. The underlying radial surface density distribution of the model disc is illustrated in Fig. 7.

Secondly, we inspect disc’s simulated SED, which is given in Fig. 8. Overall, we find that the model SED aligns closely with the observed excess data points, which the model was fitted to. Only in the regime of 100–300 |$\mu$|m there is a consistent underprediction of the observed fluxes, with deviations up to |$\sigma \!=\!3.3$|⁠, which, given the simplifications of the model, we deem acceptable. We thus note that the total disc emission is reasonably well reproduced by the best-fit model. Fig. 8 also shows the contribution to the total emission from different (radial) parts of the disc. Up to the wavelength of 18 |$\mu {\rm m}$|⁠, almost all emission comes from regions within 30 au. Beyond 18 |$\mu {\rm m}$|⁠, the emission from within 30 au stagnates, while contributions from the outer regions significantly rise. At 33 |$\mu {\rm m}$|⁠, the outer belt (regions |$>$|120 au) becomes the dominant contributor with 50 per cent of the total emission, rising to 70 per cent at 50 |$\mu {\rm m}$|⁠, where emissions from the inner regions level off, and 90 per cent at 100 |$\mu {\rm m}$|⁠. At 1.3 mm, the wavelength of the ALMA data points, regions within 120 au contribute only |$1.1$| per cent to the total flux. (Note that the ALMA data points represent only the direct flux from the belt, which due to the low contribution from the interior regions still align with the total emission SED shown in Fig. 8.) The face-on surface brightness generated by our model at 1.3 mm interior to the belt is minuscule, due to a dearth of large particles, and within 100 au relatively constant at around 0.01 MJy sr|$^{-1}$| (⁠|$\approx \!2\times 10^{-4}~{\rm mJy\, arsec}^{-2}$|⁠).

6 COMPARISON WITH IRS SPECTRUM

For a more detailed examination around mid-infrared wavelengths we additionally compare the best-fit model to the spectra taken by the Spitzer Infrared Spectrograph (IRS), namely the high-resolution spectra taken of the central region (PID 90, reduction by Su et al. 2013) and the low-resolution spectrum taken of the south-east ansa of the outer ring [PID 1074, reduction from the Spitzer Heritage Archive (SHA), see also Stapelfeldt et al. 2004]. Given the complexity of how the disc emission contributes to the spectrum (as outlined below), we chose not to include this data in the fitting and rather to perform a consistency check on the best-fit model.

To simulate IRS spectra taken of our model disc, we generate pseudo-images of the disc and central star at a range of wavelengths (at 1 |$\mu {\rm m}$| steps) and measure the total flux within the extraction field of the respective spectrum recording. The pseudo-images are synthesized by first generating astrophysical scenes of the inclined disc (with rave) with a superimposed star represented by a point source rescaled to the stellar flux at each wavelength. These scenes are then convolved with the respective PSFs for each IRS module at each wavelength, which are obtained from the STinyTim tool (Krist 2006). An example of such a pseudo-image with outlines of the corresponding extraction fields is given in Fig. 9. The different IRS modules are designated SH (shorter wavelengths, high resolution) and LH (longer wavelengths, high resolution) used for the on-star spectrum of PID 90, and LL (longer wavelengths, low resolution) used for the on-ansa spectrum of PID 1074.

We determine the relative positioning of the extraction fields using the pointing information from the corresponding observation files on the SHA, combined with Fomalhaut’s position, calculated at the time of the observations using coordinates and proper motions provided by the Hipparcos mission (van Leeuwen 2007, data obtained from SIMBAD data base).

To mimic the reduction of the reference spectra, we also apply the slit loss correction for a point source to the simulated IRS spectra, even though the underlying sources are not strictly point sources. The corresponding correction factors were calculated as the inverse of the fraction of light that falls within the slit from a point source at the slit centre (or at the nod reference locations at 1/3 or 2/3 of the slit length in the case of the SH module) convolved with the respective STinyTim PSF. For the high-resolution spectrum, we merge the partial spectra of the SH and LH modules, computed individually up to, and down to 19 |$\mu {\rm m}$|⁠, respectively, where their wavelength ranges overlap. Note that we did not mimic the wavelength-dependent variable width extraction of the LL module extraction, but simply measured the total flux within the tapered slit aperture (seen in Fig. 9) resulting from the two-nod position reduction pipeline, which may affect their comparability. On the other hand, the full-slit extraction of the SH and LH module reduction pipeline should be comparable to our approach.

The resulting simulated spectra are shown alongside the observed ones in Fig. 10. In the case of the on-star spectrum from PID 90, we compare the excesses (i.e. we subtract our model photosphere fluxes from both the simulated and the observed spectrum) whereas the on-ansa spectrum from PID 1074 combines the local disc emission with the PSF wings of the central source. First off, it should be noted that the on-star spectra (both modelled and observed) consistently diverge from the total emission of the model disc at longer wavelengths, as expected, due to the only partial coverage of the cold outer ring by the slit field of view. The observed spectra are, naturally, more appropriately compared to our simulated IRS spectra.

For the on-star spectrum, we find that fluxes are underpredicted by our model, although the spectrum’s overall shape appears to be reasonably reproduced. In particular the marginal slope at wavelengths of 20–30 |$\mu$|m, which markedly diverges from the trend in the total disc emission, is reproduced well. This flattening of the spectrum is caused by the size and orientation of the LH module’s extraction field of view. At these wavelengths, the cold outer belt rapidly gains in brightness, which is largely missed by the LH module field of view. Tests with slit orientations aligned with the disc major axis, which better capture the outer belt, resulted in a spectrum that is more aligned with the steeper curve of the model disc’s total emission.

Moreover, the discontinuity at 19 |$\mu {\rm m}$| in our simulated spectrum appears to resemble the sudden flux increase in the observed spectrum at around 18–19 |$\mu$|m. This discontinuity is due to our discrete switchover from the SH to the LH module, which, given their different fields of view, are differently affected by the not strictly applicable point source correction and the imperfect pointing. For one, the point source slit loss correction can overcompensate for the somewhat – and at longer wavelengths increasingly so – extended source. (Note that for an entirely flat source, slit losses are exactly compensated by gains from off-slit sources making no correction necessary.) Compared to the SH slit, the wider LH slit captures more of that extended-source flux, which is also amplified by the point source correction, allowing the simulated spectrum to even slightly exceed the actual total disc flux at around 19–20 |$\mu$|m. On the other hand, the offset of the slit centre (or the nod reference positions) from the stellar position can cause the slit loss correction to undercompensate. While the offset of the reference positions from the star is similar in the nodded SH and the LH observations (⁠|$\sim$|0.5–1 arcsec), this effect is more significant in the case of the smaller field of view of the SH module, and even causes the simulated excess spectrum to become negative towards smaller wavelengths. We find that at |$<$|15 |$\mu$|m, where the disc flux is on the order of 1 per cent of the stellar flux, slit position changes of even just a few 100 mas severely impact the resulting simulated IRS excess spectrum.

The overall offset between the observed and the simulated on-star excess spectra may be reconciled by considering uncertainties in the absolute calibration of the total fluxes. With the stellar photosphere flux being the dominant contributor over the IRS wavelength range, even minor differences in the calibrated total fluxes (or the level of photosphere subtracted) lead to significant differences in the excess fluxes. Su et al. (2013) estimate a spectrophotometric accuracy of 5 per cent of their reduction which includes a calibration with custom spectral response functions derived from spectra of other stars (which we do not replicate here). The observed excess spectrum can be brought into reasonable alignment with the simulated one by reducing its total fluxes by |$\sim 4$| per cent, as demonstrated in Fig. 10. In fact, a re-calibration of the spectrum of Su et al. (2013) by Gáspár et al. (2023) has brought about a flux reduction of that magnitude, which would improve the fit of our simulated spectrum.

In the case of the on-ansa spectrum, we find that our model reproduces the overall spectrum shape well, yet at shorter wavelengths, it overpredicts the fluxes by a factor of up to two (at around 20 |$\mu {\rm m}$|⁠). This is to some extent surprising. If anything, we would rather have expected an underprediction by our model, since the south-east ansa of the actual Fomalhaut disc appears somewhat brighter than the north-west ansa (Gáspár et al. 2023) – an asymmetry which by construction we do not intend to reproduce given our axisymmetric model. However, at these shorter wavelengths, the spectrum is dominated by the PSF wings of the stellar source, not by the cold outer belt, which also gives rise to the well-reproduced negative slope (as also noted by Stapelfeldt et al. 2004). This suggests that the discrepancy more likely stems from inaccuracies in our simplified reduction or in the calibration of the observed spectrum, rather than from differences in the underlying disc emission.

We may also compare the IRS fluxes to the total fluxes in the MIRI images, observed and synthetic ones, which are overplotted in Fig. 10. It is evident that the total fluxes in the synthetic coronagraphic images, F1550C and F2300C, align neither with the on-star IRS spectrum, nor with the actual total emission of the model disc. This is because much of the flux originating from within several 10s of au is lost due to the transmission drop-off at short working angles, as already noted in Section 5. Additional losses occur due to the limited field of view of these filters, particularly in the case of the Lyot coronagraph, F2300C, where emission from the outer belt is significant. The total fluxes in the coronagraphic images are thus only 36 per cent (F1550C) and 30 per cent (F2300C) of the model disc’s total emission. In contrast, 79 per cent of the flux is recovered in the non-coronagraphic image, F2550W, all losses (157 mJy) being due to our masking of pixels within 1.2 arcsec from the star to mimic the saturation in the observed image.

We find that, in the case of the coronagraphic filters, the total image fluxes in the synthetic images align well with those from Gáspár et al. (2023) (after reversing their additional scaling factor of |$3.14$| for the F2300C image, see Section 5). The underprediction of the F2550W image flux by about 200 mJy stems from our model’s aforementioned omission of the disc halo, which was captured by the wider field of view of this filter. Including such a halo in our model would increase the total disc flux and the total flux in the synthetic F2550W image accordingly, while leaving observables with more constrained fields of view, such as the coronagraphic images and the on-star IRS spectra, relatively unaffected.

Given that we fitted our model to the F1550C and F2550W radial profiles up to the outer edge of the belt, the agreement (F1550C) and underprediction (F2550W) of total image fluxes are anticipated. The consistently low F2300C total fluxes, on the other hand, are reassuring and underscore the challenge of comparing such coronagraphic photometry with the IRS on-star excess spectrum.

7 DISCUSSION

This study sought to determine whether Poynting–Robertson (PR) drag transport from the outer belt alone could explain the significant amount of dust observed interior to the outer planetesimal belt of the Fomalhaut debris disc. To do this, we applied an analytical PR drag model to the Fomalhaut system, which allowed the exploration of a vast particle parameter space, pertaining to the dust composition and collisional strength. In a parameter grid search, the model was simultaneously fitted to the system’s excess emission SED and to radial brightness profiles extracted from resolved JWST/MIRI images.

We found that the observed radial characteristics of the extended disc could be explained solely by PR-drag transport, particularly when considering instrumental effects, namely PSF blurring and the transmission drop-off of the coronagraphs. Good fits to the observations postulate a significant water ice component in the dust grain material, and place firm constraints on the collisional strength of particles in the 10s of microns size regime.

Further validation of the best-fit model is given by comparison with IRS spectra. Our simulated spectra reproduce the general features of the observations, in particular the flatness of the star-centred spectrum at 20–30 |$\mu$|m, which we find is largely due to the orientation of the extraction aperture of that particular IRS observation, not necessarily reflecting the trends in the disc’s total emission.

The degree to which our rather simplistic model can reproduce the observations may appear surprising. However, the model’s ability to consistently explain various types of data – that is, the SED, MIRI-retrieved radial brightness profiles, and IRS spectra – suggests robustness rather than mere adaptation to specific data set anomalies. In fact, the fair reproduction of the F2300C profile and the IRS spectra occurred without them being fitted to. This robustness is supported by the limited number of free parameters and the analytical nature of the PR drag model, being applied broadly without tailoring excessively to any single feature of the Fomalhaut system.

Significant simplifications inherent to the model include: the absence of planets affecting the disc structure (e.g. through ejection or resonant trapping; Bonsor et al. 2018; Sommer, Yano & Srama 2020); the omission of other dust sources, such as comets or another collisional planetesimal reservoir (e.g. in the region of the warm inner disc component, termed ‘asteroid belt analogue’ by Gáspár et al. 2023); and the assumption that a single material composition and collisional strength prescription apply throughout the disc. Nonetheless, given the level of detail we are aiming to emulate – which notably excludes azimuthal features like the intermediate ring – the Fomalhaut disc may not require a more complex representation. Specifically, it might lack phenomena we have not modelled, such as giant planets capable of carving distinct cavities or comets delivering substantial dust quantities. In any case, our modelling supports the notion that the dust distribution in the Fomalhaut disc is consistent with a pure PR drag transport scenario. We will explore the implications of some of these simplifications further in the following sections, where we discuss aspects of the model and our obtained parameter constraints in more detail.

7.1 Collisional strength

Our model fitting strongly selects for the dust particles’ collisional strength parameters, |$Q_\mathrm{0}$| and |$b_\mathrm{s}$|⁠, which determine the size dependent critical dispersion energy: |$Q_\mathrm{D}^\star = Q_\mathrm{0}(D/{\rm cm}) ^ {b_\mathrm{s}}$|⁠. The correlation between these parameters revealed that it is the resulting |$Q_\mathrm{D}^\star$| at a few tens of micrometrs, rather than the individual parameters, that permits a well-fitting model outcome, thus constrained to a range of (2–4)|$\, \times 10^6\, {\rm erg\, g}^{-1}$|⁠. We assess this specific constraint by comparing it with extrapolations from conventional |$Q_\mathrm{D}^\star$| prescriptions, as illustrated in Fig. 11. These prescriptions are based on numerical and experimental simulations at target sizes greater than 1 cm, depicted with solid lines, and are commonly extrapolated to smaller sizes, indicated by dashed lines.

For instance, we include models by Benz & Asphaug (1999) for both solid water ice and basalt, derived from smoothed particle hydrodynamics (SPH) simulations across a broad range of target sizes (few cm to 100s of km). These models exhibit the well-known features of |$Q_\mathrm{D}^\star (D)$|⁠, that is, a strength regime with a negative slope for targets sizes |$<$|100 m, where gravitational forces are negligible, and a gravity regime with a positive slope for larger sizes, where the energy required for dispersal rises steeply due to the need to overcome gravitational binding. Benz & Asphaug (1999) found similar magnitudes of strength for basalt and ice, though basalt targets appeared easier to disrupt at slower impact velocities, while the opposite was true for ice. However, the ice strengths reported by Benz & Asphaug (1999) are considered questionable, as laboratory experiments have indicated strengths up to two orders of magnitude lower (Arakawa 1999; Ryan, Davis & Giblin 1999; Giblin, Davis & Ryan 2004). This discrepancy has prompted the adoption of downscaled ‘weak ice’ prescriptions that are more in line with experimental outcomes (Wyatt & Dent 2002; Leinhardt & Stewart 2009).

In comparison, our constraint on |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$| sits slightly below the extrapolations for ‘weak ice’. It is important to note that the collisional velocity for obtaining this constraint (127 m s|$^{-1}$|⁠, following our assumptions for the planetesimal belt, see Section 3.3), is smaller than those informing the depicted weak ice strength laws. However, no definitive trends with velocity were identified across laboratory impact experiments with ice within the 100–1000 m s|$^{-1}$| range, potentially due to variations in target preparation (Leinhardt & Stewart 2009).

For rocky targets, the speed-dependence of |$Q_\mathrm{D}^\star$| is better understood. Based on experimental data covering impact velocities from around 100 m s|$^{-1}$| to 3 km s|$^{-1}$|⁠, Nakamura et al. (2015) derived |$Q_\mathrm{D}^\star$| scaling laws for porous sulfate (gypsum, porosity |$\sim$|67 per cent) and non-porous silicate (pyrophyllite) dependent on target diameter as well as impact velocity. Their results, consistent with those for basalt by Benz & Asphaug (1999), indicate that the disruption threshold scales with impact velocity. Notably, this velocity dependence is more pronounced for porous targets where |$Q_\mathrm{D}^\star$| increases by a factor of |$\sim 6$| when the velocity increases tenfold, compared to the non-porous targets, where |$Q_\mathrm{D}^\star$| merely doubles for the same velocity change. The scaling laws provided by Nakamura et al. (2015) offer a |$Q_\mathrm{D}^\star {}(D)$| curve for these materials at velocities akin to those in our model, as illustrated in Fig. 11. The strength model for porous sulfate aligns with our |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$| constraint, whereas the model for non-porous rock is about an order of magnitude higher.

While the order-of-magnitude agreement of our |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$| constraint with the |$Q_\mathrm{D}^\star$| power laws for ‘weak ice’ and porous rocky material is reassuring, the validity of extrapolating these power laws from their tested regime at centimetre sizes down by 2–3 orders of magnitude to tens of micrometrs remains questionable. The possibility of an inflection point and even a positive slope regime at those small sizes (which in fact provides the best fit of our model) cannot be ruled out.

Moreover, the strength constraints obtained from our model are influenced by several factors. Since the analytical model primarily constrains the particles’ collisional lifetime, which depends on collisional velocity, adjustments to model parameters can significantly alter outcomes. For instance, reducing the semi-opening angle |$\epsilon$| by a factor of 1/3rd, that is, from 1.5|$^\circ$| to 1.0|$^\circ$| (as found by Boley et al. 2012), requires nearly halving the |$Q_\mathrm{D}^\star$| normalization to maintain the same optical depth model outcome. Conversely, an increase in |$\epsilon$| by 1/3rd to 2.0|$^\circ$| necessitates nearly doubling |$Q_\mathrm{D}^\star$|⁠. An increase of the collisional velocity might also result from considering the belt eccentricity or a higher assumed mean proper eccentricity, |$e_\mathrm{p}$|⁠, which will similarly raise our constraint on |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$|⁠. Furthermore, cratering collisions, which our current model omits, could significantly reduce the collisional lifetimes of particles. If we assume a universal reduction in collisional lifetimes due to cratering by a factor of 5 – a value in line with theoretical considerations (Kobayashi & Tanaka 2010) – |$Q_\mathrm{D}^\star$| would need to be increased by a factor of |$\sim 7$| to maintain equivalent model outcomes. Such an adjustment would position our |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$| constraint in between the plotted ‘weak ice’ and pyrophyllite models.

Beyond the more conventional |$Q_\mathrm{D}^\star$| prescriptions derived from impact experiments, our constraint can also be compared to the extreme collisional strengths recently deduced by Pokorný et al. (2024), who fitted models of Jupiter-family comet (JFC) micrometeoroids to Solar system observables such as the orbital element distributions of meteors. Similar previous studies also suggest that cometary meteoroids possess collisional lifetimes |$1\!-\!2$| orders of magnitude longer than those allowed by conventional collisional strengths (Nesvorný et al. 2011; Pokorný et al. 2014; Soja et al. 2019). Given that our |$Q_\mathrm{D,\, {32}\, {\mu }m}^\star$| constraint aligns with the lower end of conventional prescriptions, it starkly contrasts with these extreme Solar system strengths.

For reference, Fig. 12 shows the modelled collisional and PR drag time-scales for particles in the Fomalhaut belt, as a function of particle size, for our best-fit outcomes, which are broadly consistent with the time scales derived by Su et al. (2024) for the Vega debris disc – a system with remarkably similar properties to Fomalhaut.

One potential explanation for the discrepancy between the conventional strengths and those inferred from Solar system observations is the different collisional velocity regimes being probed. Collisional velocities in the Fomalhaut belt are on the order of 0.1 km s|$^{-1}$|⁠, while grains on JFC orbits colliding in the inner Solar system may easily reach relative velocities of 10s of km s|$^{-1}$| (see typical velocity distributions of JFC meteoroids at 1 au, e.g. Nesvorný et al. 2010, 2011; Szalay et al. 2019). Higher impact velocities and material porosity tend to increase collisional strengths for rocky materials (Benz & Asphaug 1999; Jutzi et al. 2010; Nakamura et al. 2015). However, extrapolating the high-porosity |$Q_\mathrm{D}^\star$| models of Nakamura et al. (2015) (tested up to 3 km s|$^{-1}$|⁠) to 10 km s|$^{-1}$| results in |$Q_\mathrm{D}^\star$| values ranging from 10|$^8$| to 10|$^9$| erg g|$^{-1}$| – still nearly two orders of magnitude lower than those reported by Pokorný et al. (2024).

To speculate, at velocities of 10s of km s|$^{-1}$|⁠, which exceed the capabilities of impact experiments with light-gas guns, collisional strengths might be higher than anticipated, perhaps due to largely damage-free penetration of fluffy grains by highly energetic but small projectiles. Such a hypothesis would be consistent with the fact that extreme collisional strengths are in particular required to fit the meteor orbit distributions and the mass flux to Earth (as indicated in Pokorný et al. 2024, fig. 4B therein) – both arguably governed by larger micrometeoroids, which spend more time near their highly eccentric source orbits and are thus more exposed to higher collisional velocities. Similarly, in the studies by Nesvorný et al. (2011) and Soja et al. (2019), only particle sizes above the 0.1–1 mm range require lifetimes significantly extended beyond conventional estimates.

This reasoning might offer a way to reconcile the extreme strengths required for the fitting of Solar system observations with those derived from impact experiments, which, conversely, appear compatible with observations of debris discs.

7.2 Grain composition

The fitting of our model also constrains the parameters of the dust grain model we employed. The results suggest that the particles are rather compact composites, consisting of porous aggregates of silicate-core, carbonaceous-mantle grains, with their voids filled with water ice, comprising |$50$|–|$80$| per cent of the total particle volume. Such a composition is generally consistent with what is known about the Solar system’s Kuiper belt, where densities measured in (binary) smaller trans-Neptunian objects, typically fall below 1 g cm⁻³, implying considerable ice volume fractions (Barucci & Fornasier 2023). More specifically, we may consider the composition of JFCs, which are known to stem from the Kuiper belt (Duncan, Quinn & Tremaine 1988; Fraser et al. 2024). Our constraint of the ‘refractories-to-volatiles’ ratio, represented as |$(m_\mathrm{sil}\!+\!m_\mathrm{carb})\, /\, m_\mathrm{ice}$| in our simplified grain model, of |$0.43\!-\!1.6$| (Scn. A) or |$0.43\!-\!0.77$| (Scn. B), is compatible with that found in the mass loss of comet 67P/Churyumov–Gerasimenko, which was scrutinized by the Rosetta probe (⁠|$0.7\!-\!2.3$|⁠, Biver et al. 2019), and falls on the lower end of the ratios inferred for its nucleus, which show significant discrepancy (e.g. |$>3$|⁠, Fulle et al. 2019; |$0.5\!-\!1.7$|⁠, Marschall, Morbidelli & Marrocchi 2025). Substantial amounts of water ice have also been inferred in other debris discs (e.g. Chen, Fitzgerald & Smith 2008; Lebreton et al. 2012; Lisse et al. 2012) and protoplanetary discs (e.g. Pontoppidan et al. 2005; Honda et al. 2008).

The relative abundance of carbonaceous material is not strongly selected by our model, as long as it is not too dominant (⁠|$q_\mathrm{sil}> 0.2\Leftrightarrow q_\mathrm{carb} = 1-q_\mathrm{sil}< 0.8$|⁠). This is consistent with the expected cosmic abundance of carbonaceous material relative to silicates, at a ratio of |$m_\mathrm{carb}/m_\mathrm{sil}\!\approx \!0.7$| (Li & Lunine 2003), equivalent to |$q_\mathrm{sil}\!\approx \!0.42$|⁠.

The structural grain model we employ is generally compatible with cosmic dust particles collected from the stratosphere, likely originating from the Kuiper belt region (their origin inferred from cosmic ray exposure ages, Keller & Flynn 2022; though see also Lin & Poppe 2024, who explore additional complexities and unresolved questions in track accumulation modelling), that is, porous aggregates of mineral grains bound by carbonaceous material. These may have held ices within their cavities initially, which would have naturally been lost due to their collection not only within the Solar system’s sublimation zone but also after atmospheric entry heating.

This leads us to another significant simplification in our model: the assumption that a single particle composition applies universally throughout the entire disc. However, in a realistic disc, the ice content of particles may change with radial location due to sublimation and/or UV photosputtering. The sublimation of the water ice is expected to occur once particles reach a temperature of |$\sim$|110 K (Kobayashi et al. 2008), which in the case of Fomalhaut corresponds to a radial distance of |$\sim$|30 au, depending on particle material and size. Yet, we observe no abrupt change of brightness in the resolved disc at around that distance, which might be expected from changes in the optical properties of particles as ices sublimate from aggregate voids. In particular, we see no evidence for a sublimation pile-up and subsequent depletion as theorized by Kobayashi et al. (2008). A distinct sublimation pile-up requires near-circular particle orbits (⁠|$e\!\lesssim \!0.01$|⁠; Kobayashi et al. 2008, 2011), which may explain its absence, given the Fomalhaut disc’s eccentric intermediate ring feature (⁠|$e\!\approx \!0.3$|⁠; Gáspár et al. 2023). Moreover, an optical depth reduction due to ice loss might not manifest when more realistic particle properties are considered. While the Kobayashi model assumes sublimating particles to shrink down to the compact refractory remnants, a grain model akin to ours may instead leave behind highly porous aggregates after losing its ices, potentially maintaining the original cross-section. Discrete dipole approximation simulations suggest that such irregular porous aggregates are relatively less sensitive to radiation pressure blow-out (with carbon aggregates of |$\sim 75$| per cent porosity exhibiting only twice the blow-out size of compact carbon grains around A-type stars, Arnold et al. 2019), which may enable these remnants to remain bound.

Exposed ices may also be removed from particles by stellar light far-UV photosputtering, which around A-type stars such as Fomalhaut could be effective even beyond 100 au (Grigorieva et al. 2007). Such early ice removal might seem at odds with the optical properties constrained by our model (and could help explain the absence of sublimation features). Indeed, the best fits across our models that exclude water ice entirely exhibit deficiencies in total emission at around 60–150 |$\mu$|m and in the resolved mid-infrared emission of the interior disc compared to the outer belt, indicating that the presence of ices is necessary to adequately reproduce these observables with the pure PR drag transport scenario. However, far-UV photosputtering (⁠|$\lambda _\mathrm{photon}\lesssim 200\, {\rm nm}$|⁠) can remove ices only from the outer layers of rock-ice mixture particles that are 10s or 100s of |$\mu$|m in size, likely leaving their optical properties at mid-infrared and longer wavelengths largely unaffected. Thus, it is conceivable that a more detailed grain model, in which ices are removed only from the outer layers of particles, could achieve fits similar to our current model, though exploring this goes beyond the scope of this paper. Explanations along these lines may also be necessary to reconcile previous studies that likewise invoked icy particles in debris discs around A-type stars to simultaneously fit excess SEDs and resolved infrared emission (Acke et al. 2012; Morales et al. 2016; Adams et al. 2018) with the photosputtering effect. Future JWST (near-infrared) scattered light observations may provide more definitive insights into the water ice content of particles (Kim et al. 2019; Kim, Kennedy & Roccatagliata 2024).

For reference, Fig. 13 shows the absorption coefficients, |$Q_\mathrm{abs}$|⁠, as computed with our optical properties model for our best-fit parameters (Scn. A), as a function of wavelength and particle size.

7.3 Planets

Our model has successfully replicated the observed smooth radial brightness distribution in the Fomalhaut disc without considering the presence of planets. Massive planets, however, would be expected to introduce distinct radial features, such as cavities carved by the ejection of particles from the system in close encounters, or ring-like overdensities caused by the trapping of particles in exterior mean-motion resonances. While prominent resonant rings may arguably be suppressed by particle collisions in a disc as dense as Fomalhaut’s (Kuchner & Stark 2010; Thebault, Kral & Ertel 2012), the process of cavity formation through particle ejection should be unaffected. One may thus be drawn to interpret the absence of clear brightness drops within the extended disc as an absence of giant planets.

To explore this further, we apply the formalism developed by Bonsor et al. (2018), which calculates the fraction of migrating particles lost to ejection or accretion by a planet based on the planet’s mass, semimajor axis, and the particles’ |$\beta$|-factor. Integrating this model across all particle size bins in our best-fit model disc – weighted by each bin’s contribution to the optical depth at specific radii – enables us to estimate the total reduction in optical depth for various planetary configurations, as shown in Fig. 14. This in turn approximates the expected brightness reduction a given planet would introduce. Generally, the more massive the planet and the larger its semimajor axis, the higher the loss in optical depth.

In exploratory tests, we post hoc applied optical depths reductions to our best-fit disc according to the Bonsor et al. (2018) formalism and examined the resulting synthetic images and radial brightness profiles. Crudely, we estimate that optical depth reductions of at least 30–70 per cent, depending on radial location, would have been discernible, as indicated by the dash–dotted line in Fig. 14. Brightness drops are more readily discerned at larger disc radii due to the blurring effect of the PSFs, obscuring smaller-scale structures, as well as the coronagraphs’ transmission drop-off. Based on this preliminary assessment, we are inclined to rule out the presence of planets more massive than Jupiter beyond |$\sim$|10 au, and those more massive than Saturn beyond |$\sim$|50 au.

Fig. 14 also shows the regime recently probed by JWST/NIRCam coronagraphy in a search for planets around Fomalhaut, achieving a sensitivity capable of detecting planets down to approximately Jovian mass (Ygouf et al. 2023). Out of 10 detected objects, 9 were confirmed as background objects. Object ‘S7’, however, although detected with only one of the two used filters, could thus-far not be attributed to a background object. Follow-up NIRCam observations to confirm or reject this planet candidate, seemingly located at a separation of |$\sim$|30 au from Fomalhaut, are underway (Beichman et al. 2023). Based on our analysis, we would expect a Jovian-mass planet at 30 au to carve a cavity with an optical depth drop on the order of 80 per cent. This is hardly compatible with the observed brightness distribution, as similarly argued by Ygouf et al. (2023). Nevertheless, the follow-up observations, which aim for a detection threshold as low as 0.3 |$M_{\rm Jupiter}$| (i.e. around Saturn’s mass), will serve as a test of the tentative upper limits we have placed on planetary masses in the Fomalhaut disc.

Chiang et al. (2009) constrained the parameters (⁠|$M_\mathrm{pl}$| and |$a_\mathrm{pl}$|⁠) of a hypothetical planet that could have sculpted the inner edge of the planetesimal belt as observed in scattered light. These parameters are represented in Fig. 14 by an orange dotted line. The constraints are based on the assumption that the outer boundary of the arising chaotic zone – where no stable orbits can be maintained – extends just to the inner edge of the planetesimal belt. This is described by the equation:

where |$a_\mathrm{inner}=133~ \rm {au}$| is the semimajor axis of the inner edge of the planetesimal belt. Chiang et al. (2009) additionally suggested that the mass of such a ring-shepherding planet should be |$< 3\, M_\mathrm{Jupiter}$|⁠, to avoid perturbing the belt particles onto more eccentric orbits, which would be inconsistent with the belt’s observed optical depth profile.

A shepherding planet would also be expected to eject dust migrating inward from the belt, and thus cause a reduction in optical depth. Indeed, the MIRI images of the disc indicate a brightness dip occurring near the inner edge of the belt, which however our planet-less model reproduces, as only a fraction of particles can effectively migrate from the belt before being ground down to blow-out sizes. Nevertheless, the magnitude of this ‘collisional depletion’ depends on the dust’s collisional strength, introducing a degeneracy; the observed dip might not solely result from collisions but also from planet-induced dynamical depletion. In our model, adjusting the collisional strength to be up to 10 times greater than in our best fit allows us to mimic the observed dip while accounting for the dynamical depletion caused by a shepherding planet with a mass up to |$0.5\, M_\mathrm{Jupiter}$|⁠. Higher levels of dynamical depletion, however, are increasingly difficult to offset with enhanced collisional strength. This suggests that the pure PR drag scenario for the Fomalhaut disc remains valid with the inclusion of a belt-shepherding planet up to |$0.5\, M_\mathrm{Jupiter}$|⁠, albeit with collisional strengths somewhat higher than those favoured in the planet-less case. The inclusion of a shepherding planet of a mass of |$0.2\, M_\mathrm{Jupiter}$| as hypothesized by Janson et al. (2020) to explain the transient dust cloud ‘Fomalhaut b’ (Kalas et al. 2008; Galicher et al. 2013; Gáspár & Rieke 2020) via tidal disruption, necessitates a doubling of collisional strengths to recover an equivalent outcome. Shepherding planets below |$0.1\, M_\mathrm{Jupiter}$| cause only marginal dynamical depletion and do not require an adjustment of collisional strengths to maintain the fit.

It is important to note that the ejection model of Bonsor et al. (2018) assumes planets to be on circular orbits, which may not align with the observed eccentricity of the planetesimal belt and the ‘intermediate ring’ brightness feature. Interestingly, this ‘intermediate ring’ visible in the F2300C and F2550W images, does not appear to have a significant impact on the azimuthally averaged radial brightness profiles. This suggests that if a planet is producing this feature (e.g. through resonant trapping), it is not massive enough to eject a substantial fraction of the migrating dust.

7.4 Implications for other debris discs

Let us put our findings about the Fomalhaut disc, namely, the plausibility of a pure PR drag transport scenario, into context. Before JWST’s revelations about Fomalhaut’s extensive interior disc, inferences about its structure typically relied on SED analysis of the unresolved excess, to determine the underlying dust temperature and thus location. With that approach, Su et al. (2013, 2016) derived a warm disc component with a blackbody temperature of |$\sim$|170 K, thus located at |$\sim$|8–15 au and presumably maintained by in situ dust production by a narrow, collisional asteroid belt analogue. Therefore, even when considering the degeneracy between the dust location and the dust’s uncertain emission properties, the discovery of the extended interior disc by Gáspár et al. (2023) came as a surprise, ruling out radially confined in situ production as the dust’s sole source.

Several 10s of debris discs with unresolved warm components have been identified whose infrared excesses can be fitted with two-component, or rather two-temperature (i.e. cold and warm) disc models, suggesting radially confined, warm belts seemingly associated with snowline locations (e.g. Morales et al. 2011; Ballering, Rieke & Gáspár 2014; Kennedy & Wyatt 2014; Morales et al. 2016; Ballering et al. 2017), thus often interpreted as exo-asteroid belts. The case of Fomalhaut begs the question: Is an extended, and thus potentially PR-drag-caused interior disc component the more common explanation for those unresolved infrared excesses?

In any case, the effectiveness of PR drag transport in explaining the appearance of the Fomalhaut disc underscores its relevance in the JWST era. Our findings support the conclusions of previous studies, which, even before the advent of highly resolved imagery, posited that PR drag alone could sustain the inferred abundances of warm dust around specific stars (Reidemeister et al. 2011; Löhne et al. 2012; Mennesson et al. 2014; Schüppler et al. 2014; Rigley & Wyatt 2020). Moreover, results on another debris disc captured by JWST/MIRI, namely that around A-type star Vega, were recently published (Su et al. 2024). These have likewise revealed a smoothly filled-in disc interior to Vega’s planetesimal belt (though without azimuthal features), consistent with pure PR drag evolution. More JWST observations specifically investigating the delivery mechanism of warm dust in debris disc systems are underway.

As established by Kennedy & Piette (2015), given the inevitability of PR drag and the relative insensitivity of interior PR-drag-delivered dust levels to belt density, all systems with detectable outer planetesimal belts, unless hosting a giant dust-ejecting companion, will be significantly permeated with dust. This has profound implications for the study of these exoplanetary systems. Notably, dust in habitable zones poses a significant challenge in the search for exo-Earths (e.g. Roberge et al. 2012; Currie et al. 2023), and deviations from the nominal (i.e. PR-drag-delivered) warm dust distribution can provide insights into the presence of planets or comets (e.g. Bonsor et al. 2018; Marino et al. 2018; Sommer et al. 2020). Recognizing these deviations and understanding their origins will require precise modelling of the PR drag scenario. Therefore, assuming our model adequately captures the nature of the Fomalhaut disc, the particle parameters we have constrained will serve as an essential calibration point for future efforts aimed at estimating or interpreting the distribution of warm dust in similar systems.

8 CONCLUSION

Our main conclusions are the following:

• The observed radial brightness distribution in the Fomalhaut disc is consistent with a pure PR drag scenario, that is, no other transport mechanisms appear to be necessary to generate the underlying dust distribution.

• We find that the faithful replication of instrumental effects and reduction methods is critical when trying to compare model outcomes with disc observables. This pertains for instance to the MIRI coronagraphy, where the transmission drop-off at small working angles and PSF blurring considerably skew the derived radial brightness profile of the inclined disc, or to Spitzer IRS spectrography, where the exact positioning of the extraction field in combination with the applied point-source correction significantly influences the obtained spectrum. Taking these effects into account, we find that our model naturally reproduces a significant offset between the resolved fluxes in the MIRI F2300C and F2550W images, previously noted by Gáspár et al. (2023, their supplementary section 1.2), despite more comparable fluxes at these wavelengths in the underlying source and synthesized spectra. This suggests that these datasets can be reconciled without necessarily invoking an issue with the F2300C reduction pipeline.

• We constrain the collisional strength of the particles in this pure PR drag scenario, specifically the catastrophic disruption threshold |$Q_\mathrm{D}^\star$| at a particle size of |$D\!\approx \!32\, \mu {\rm m}$| to a range of (2–4)|$\, \times 10^6\, {\rm erg\, g}^{-1}$|⁠. This roughly aligns with conventional strength scaling laws derived from impact experiments for solid ice (Leinhardt & Stewart 2009) as well as porous rock (Nakamura et al. 2015). The extreme strengths recently deduced for Solar system cometary grains (Pokorný et al. 2024) are contradicted by our findings, though this might be due to the different velocity regimes being probed. Accounting for cratering collisions, which our model omits, may increase our |$Q_\mathrm{D}^\star$| constraint by up to an order of magnitude.

• We constrain the disc particles to be rather compact composites, consisting of porous aggregates of silicate-core, carbonaceous-mantle grains, with their voids filled with water ice, comprising |$50$|–|$80$| per cent of the total particle volume. Notably, our simplified grain model, which does not account for changes in composition, effectively reproduces the observed radial brightness distribution, which in fact lacks the hypothesized sublimation depletion at the snowline of icy dust discs (Kobayashi et al. 2008). While we do not suggest that ices survive within the snowline (⁠|$\sim$|30 au), the absence of a discernible sublimation-induced depletion in the disc may indicate that particle behaviours – such as sublimation from the voids of a porous refractory matrix that largely preserves optical cross-section – differ from previous theoretical expectations. Further study is necessary to explore the impact of grain structure on sublimation dynamics and disc appearance.

• The radial (i.e. azimuthally averaged) brightness distribution of the Fomalhaut disc shows no clear indications of disturbances by massive planets. Given the expected cavity formations that such bodies would cause in a disc primarily sculpted by PR drag, our findings suggest the absence of Jovian-mass planets beyond 10 au and Saturnian-mass planets beyond 50 au within the extended inner disc. Nevertheless, an outer-belt-shepherding planet of up to |$\sim 0.5\, M_\mathrm{Jupiter}$| could still be compatible with a PR-drag-maintained interior disc, though requiring higher collisional strengths to offset the dynamical depletion. An upcoming JWST/NIRCam planet search with |$\sim$|Saturn-mass sensitivity will not only test these upper limits, but also, both the detection or non-detection of a shepherding planet will help refine the constraints we can place on collisional strengths.

• Given the inevitability of PR drag, all systems with detectable outer belts – unless obstructed by a dust-ejecting giant planet – are likely permeated with dust. Recently, high levels of pervasive inner dust have also been found around Vega (Su et al. 2024), supporting this picture. More JWST observation programmes focussing on warm dust in debris discs are underway (Han et al. 2024; Matra et al. 2024), which will test the generalizability of our results for Fomalhaut. The particle parameters we have constrained, particularly the collisional strength at the relevant sizes, which thus far relies solely on extrapolations, establish a calibration point for modelling PR-drag-delivered dust around other stars, which is essential for predicting dust levels in habitable zones and for interpreting system architectures from resolved images.

ACKNOWLEDGEMENTS

We are grateful for the support by United Kingdom Reasearch and Innovation (UKRI) - Science and Technology Facilities Council (STFC) grant no. ST/W000997/1. We thank Grant Kennedy for his support and for providing the Fomalhaut photosphere model, as well as Amy Bonsor and Marija Janković for insightful discussions. We also thank Kate Su and the reviewer, Petr Pokorný, for their valuable comments and suggestions. This work is based in part on observations made with the NASA/ESA/CSA JWST. The JWST observations are associated with programme 1193.

DATA AVAILABILITY

The JWST/MIRI data underlying this article are provided by Gáspár et al. (2023) and are available at https://github.com/merope82/Fomalhaut/. The core components of our modelling framework, that is, our implementation of the PR drag disc model from Rigley & Wyatt (2020) – extended with the capability to generate astrophysical scenes – and our method to find dust optical properties and temperatures, are available as separate python packages, pyrdragdisk and astrodust_optprops, respectively, at https://github.com/rigilkent/pyrdragdisk and https://github.com/rigilkent/astrodust_optprops.

Footnotes

REFERENCES