# The debris disc of HD 131488: bringing together thermal emission and scattered light

Nicole Pawellek, Attila Moór, Florian Kirchschlager, Julien Milli, Ágnes Kóspál, Péter Ábrahám, Sebastian Marino, Mark Wyatt, Isabel Rebollido, a Meredith Hughes, et al.

► To cite this version:

Nicole Pawellek, Attila Moór, Florian Kirchschlager, Julien Milli, Ágnes Kóspál, et al. The debris disc of HD 131488: bringing together thermal emission and scattered light. Monthly Notices of the Royal Astronomical Society, 2024, 527 (2), pp.3559-3584. 10.1093/mnras/stad3455. hal-04505683

HAL Id: hal-04505683

https://hal.science/hal-04505683v1

Submitted on 16 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The debris disc of HD 131488: bringing together thermal emission and scattered light

Nicole Pawellek<sup>1,2*</sup> Attila Moór,<sup>2,3</sup> Florian Kirchschlager<sup>1,4,5</sup> Julien Milli,<sup>6</sup> Ágnes Kóspál<sup>1,2,3,7,8</sup> Péter Ábrahám,<sup>2,3,7</sup> Sebastian Marino<sup>1,9</sup>, Mark Wyatt<sup>1,10</sup>, Isabel Rebollido,<sup>11</sup> A. Meredith Hughes<sup>1,12</sup>, Faustine Cantalloube,<sup>13</sup> and Thomas Henning,<sup>8</sup>

Accepted 2023 November 6. Received 2023 November 6; in original form 2023 June 12

## ABSTRACT

We show the first SPHERE/IRDIS and IFS data of the CO-rich debris disc around HD 131488. We use N-body simulations to model both the scattered light images and the spectral energy distribution of the disc in a self-consistent way. We apply the Henyey–Greenstein approximation, Mie theory, and the Discrete Dipole Approximation to model the emission of individual dust grains. Our study shows that only when gas drag is taken into account can we find a model that is consistent with scattered light as well as thermal emission data of the disc. The models suggest a gas surface density of $2 \times 10^{-5} M_{\oplus} \text{ au}^{-2}$ which is in agreement with estimates from ALMA observations. Thus, our modelling procedure allows us to roughly constrain the expected amount of gas in a debris disc without actual gas measurements. We also show that the shallow size distribution of the dust leads to a significant contribution of large particles to the overall amount of scattered light. The scattering phase function indicates a dust porosity of ~0.2...0.6 which is in agreement with a pebble pile scenario for planetesimal growth.

**Key words:** circumstellar matter – stars: individual (HD 131488) – infrared: stars.

## 1 INTRODUCTION

Circumstellar debris discs are optically thin collections of solids ranging from planetesimal size bodies down to dust grains. All of the components are thought to be part of a collisional cascade in which larger objects are gradually ground to smaller particles through mutual destructive collisions (Wyatt 2008). Observations are only sensitive to the lowest mass end of the population: thermal emission of dust is detectable at infrared (IR) and millimetre wavelengths, while the stellar light scattered by the disc is mostly observable in the optical/near-IR regime. Besides the gravitational force exerted by the star and possible planets, the observed second generation grains are also subject to additional non-gravitational forces related to stellar radiation and wind (Krivov 2010). Depending inversely on their size

the stellar radiation pressure can push dust grains on more and more eccentric orbits forming an extended halo of barely bound particles outside the planetesimal belt. Below a certain size dust is blown out from the system by this force. By causing an inward migration of grains the Poynting–Robertson effect and stellar wind drag can affect the spatial distribution of dust as well.

The presence of gas can also influence the dynamics and the spatial distribution of dust particles. Recently, detections of far-IR O I, C II, and particularly millimeter CO lines revealed gas in some 20 debris discs (e.g. Dent et al. 2014; Lieman-Sifry et al. 2016; Marino et al. 2016; Moór et al. 2017; Matrà et al. 2019a; Schneiderman et al. 2021). In most of these systems the observed gas is likely secondary and released through collisions of large volatile-rich bodies (Kral et al. 2017, 2019; Marino et al. 2020). Remarkably, as observations of less abundant CO isotopologues implied, in a subset of this sample the mass of CO gas is on a par with that of less massive protoplanetary discs (Kóspál et al. 2013; Péricaud et al. 2017; Moór et al. 2019;

<footnote n="1"><sup>1</sup>Institut für Astrophysik, Universität Wien, Türkenschanzstraße 17, 1180 Vienna, Austria</footnote>

<footnote n="2"><sup>2</sup>Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network (ELKH), Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary</footnote>

<footnote n="3"><sup>3</sup>CSFK, MTA Centre of Excellence, Budapest, 1121 Konkoly-Thege Miklós út 15-17., Hungary</footnote>

<footnote n="4"><sup>4</sup>Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="5"><sup>5</sup>Sterrenkundig Observatorium, Ghent University, Krijgslaan 281-S9, B9000 Gent, Belgium</footnote>

<footnote n="6"><sup>6</sup>Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="7"><sup>7</sup>ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary</footnote>

<footnote n="8"><sup>8</sup>Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="9"><sup>9</sup>Department for Physics and Astronomy, University of Exeter, Stocker Road, EX4 4QL Exeter, UK</footnote>

<footnote n="10"><sup>10</sup>Institute of Astronomy, University of Cambridge, Madingley Road, CB3 0HA Cambridge, UK</footnote>

<footnote n="11"><sup>11</sup>Centro de Astrobiología (CAB, CSIC-INTA), Camino Bajo del Castillo s/n, Villanueva de la Cañada, E-28692 Madrid, Spain</footnote>

<footnote n="12"><sup>12</sup>Astronomy Department and Van Vleck Observatory, Wesleyan University, 96 Foss Hill Drive, Middletown, CT 06459, USA</footnote>

<footnote n="13"><sup>13</sup>Aix Marseille Univ, CNRS, CNES, LAM, Pôle de l'Étoile Site de Château-Gombert 28, rue Frédéric Joliot-Curie 13388 Marseille cedex 13, France</footnote>

<footnote>*E-mail: nicole.pawellek@univie.ac.at</footnote>
Rebolledo et al. 2022). All of these CO-rich debris discs surround young, 5–50 Myr old, A-type stars; their observed gas material is at least partly co-located with the cold dust in these systems. Though we can measure only a few constituents of the complete gas mixture it is probable that the total gas mass is at least comparable to that of dust measured at millimeter wavelengths (e.g. Moór et al. 2017).

High spatial resolution scattered light images of several CO-rich debris discs revealed complex structures in the distribution of those small dust grains that could be most affected by gas. Optical and near-IR observations of HD 141569A have revealed complex morphology with two rings at ~245 and 400 au as well as spiral features in the disc (Biller et al. 2015, and references therein). Millimeter interferometric CO line observations of the system showed that the inner ring is located just at the outer edge of the gas disc (Flaherty et al. 2016; Di Folco et al. 2020). Recent imaging with VLT/SPHERE showed additional concentric ringlets between 47 and 93 au spatially with the gas disc (Perrot et al. 2016). By observing the disc around HD 131835 with SPHERE, Feldt et al. (2017) also discovered concentric dust rings that are co-located with the circumstellar gas material.

Though some of these structures could be the result of perturbations by planetary or stellar companions (e.g. Augereau & Papaloizou 2004; Feldt et al. 2017), the presence of gas in these systems provides alternative explanations. In an optically thin gaseous debris disc, the combined effect of stellar radiation and gas drag induces radial drift of dust. Assuming gas pressure decreases with radius, small dust particles migrate outward and can form a narrow ring at the outer edge of the gas disc (Takeuchi & Artymowicz 2001), as in the case of HD 141569A (Flaherty et al. 2016). Considering heating of gas by photo-electrons from nearby dust grains Klahr & Lin (2005) and Besla & Wu (2007) found that this effect can lead to strong local dust enhancements via a clumping instability. Depending on the gas and dust surface density such photoelectric instability can result in sharp concentric rings providing a feasible explanation for such features in HD 131835 and HD 141569A systems (Richert, Lyra & Kuchner 2018).

The usage of scattered light data is not limited to structural analysis, multiwavelength measurements allow to investigate the grain properties as well. Based on VLT/SPHERE imaging of the gaseous debris disc around HD 32297, Bhowmik et al. (2019) reported the presence of copious amount of grains smaller than the blowout size in this system and proposed that their pile-up is related to gas drag and/or avalanche mechanisms. According to this scenario by slowing down the motion of small unbound grains – that otherwise would leave the system on the orbital time-scale (Meyer et al. 2007) – gas drag can result in an overabundance of such particles with respect to a gas free case. Interestingly, the colour of the gas-bearing debris discs around HD 36 546 (Lawson et al. 2021) and HD 141 569 (Singh et al. 2021) also suggests the presence of copious submicron-sized or highly porous grains.

To further explore gas-dust interactions in an optically thin environment, in this paper we present the first spatially resolved scattered light images of the gaseous debris disc around HD 131488 obtained with the SPHERE instrument. HD 131488 is an A1-type star at a distance of 154.0 ± 2.5 pc (Gaia Collaboration 2016; Bailler-Jones et al. 2018; Lindegren et al. 2018), that likely belongs to the ~16 Myr old Upper Centaurus Lupus subgroup of the Scorpius-Centaurus association (Melis et al. 2013; Pecaut & Mamajek 2016). The IR excess emission of the system was first identified by Melis et al. (2013). Based on its spectral energy distribution (SED), the disc has a high fractional luminosity, the dust material is likely distributed in two belts (Melis et al. 2013). Using the ALMA interferometer at 1.3 mm, recently the disc was successfully resolved in continuum

and in $J = 2 - 1$ rotational transitions of $^{12}$CO, $^{13}$CO, and $C^{18}$O lines (Moór et al. 2017). Also, the star shows a gaseous CaII absorption associated with its circumstellar environment (Rebolledo et al. 2018).

Analysis of the continuum observation implied that large cold dust grains are confined in a ring with a radius of ~0.57 arcsec (~88 au). HD 131488 has the highest $C^{18}$O line luminosity of any gas-bearing debris disc found to date, in fact its measured $L_{C^{18}O}$ is even ~1.5 × higher than that of the well known protoplanetary disc around the Herbig Ae star, HD 100453 (van der Plas et al. 2019) and ~2 × higher than that in TW Hya (Favre et al. 2013). The outstandingly high CO mass makes HD 131488 an ideal choice for a detailed study of the gas-dust interaction.

In Section 2 we will discuss the observations of the disc around HD 131488, and the data reduction of the scattered light data. Section 3 will give an overview of the theoretical background used to generate our disc models starting with orbital parameters, scattered light models, and grain composition up to generating the final model images. We will discuss the grain size distribution (Section 3.5) including the influence of gas present within the disc. In Section 4 we present the results of our modelling effort which is then followed by a discussion in Section 5.

## 2 OBSERVATIONS AND DATA REDUCTION

We observed the disc around HD 131488 in the programme 0101.C-0753(B) (PI: A. Moór) on the night of 2018 April 7 for 1 h with the SPHERE instrument of the VLT (Dohlen et al. 2008; Beuzit et al. 2019), which is fed with an extreme adaptive optics system to reach a high contrast close to the star. We used the IRDIFS observing mode combining the near-IR dual-band camera Infra-Red Dual-beam Imager and Spectrograph (IRDIS; Dohlen et al. 2008) with the Integral Field Spectrograph (IFS; Claudi et al. 2008). The IRDIS observations were carried out using the dual band H23 filter with central wavelengths of 1.593 µm for H2 and 1.667 µm for H3 and a width of 139 nm. The IFS observations dispersed the Y-J band into 39 spectral channels from 958 nm to 1.329 µm. Both observations used the coronagraph N_ALC_YJH_S (Martinez et al. 2009; Carbillet et al. 2011) with a diameter of 185 mas and were performed in pupil tracking mode to allow for angular differential imaging (ADI; Marois et al. 2006). The observing conditions were slightly worse than average for the VLT site, with an average DIMM seeing of 0.88 arcsec and an average coherence time as measured by the Paranal MASS-DIMM of 3.5 ms. For a star of magnitude $G = 8$, this resulted in an average Strehl in the H band of about 70 per cent, as estimated by the adaptive optics system, while the direct measurement performed on the average non-coronagraphic images obtained before and after the coronagraphic sequences indicate a value of 66 per cent. Despite this performance being lower than average for an instrument like SPHERE, the conditions were very stable, leading to a good data set with homogeneous quality.

The raw IRDIS and IFS data were pre-processed by the High Contrast Data Centre (HC-DC)<fnref n="1" /> (Delorme et al. 2017). This pre-processing consists of flat fielding, bad-pixel correction, background subtraction, frame registration, and the IFS wavelength calibration. It uses native recipes from the ESO Data Reduction and Handling software (Pavlov et al. 2008) complemented by additional recipes

<footnote n="1">¹The HC DC, previously known as the SPHERE DC, performs data reduction on request and also processes all SPHERE public data to make them available publicly. More information is available at https://sphere.osug.fr/spip.php?rubrique16</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top: Classical ADI-reduced scattered light image (left) and PCA-reduced scattered light image (right) of HD 131488, obtained with IRDIS at 1.6 µm (average of the two IRDIS spectral channels). The surface brightness is given in mJy arcsec⁻². North is up and East is to the left. Bottom: Respective SNR maps for classical ADI and PCA images.</figcaption></figure>

developed by the SPHERE Data Center. This pre-processing results in spectro-temporal master cubes of images. For IRDIS, this represents a sequence of 76 images in two spectral channels, spanning 28.9° of field rotation for 34 min effective integration time. For the IFS, this represents a sequence of 60 images in 39 spectral channels, spanning 26.6° of field rotation for 32 min effective integration time.

We then processed the data with a classical ADI (Marois et al. 2006) reduction technique, which consisted of building a model of the coronagraphic image from the median of all pupil-stabilized images, which was then subtracted from each frame before de-rotating and stacking the images. To improve upon this reduction, we also performed a slightly more aggressive data reduction, where the model of the coronagraphic image is constructed using a principal component analysis (PCA; Amara & Quanz 2012; Soummer, Pueyo & Larkin 2012) retaining two principal components, a value found to maximize the signal-to-noise (SNR) of the disc. The reduction was performed over the whole frame in a single area extending from 36 mas to 1.23 arcsec radially. In Fig. 1, we show the result of both reductions for IRDIS.

The image was normalized to mJy arcsec⁻² in the following way. On the non-coronagraphic image, we measured the flux density encircled within a circle of radius 0.1 arcsec, encompassing the point spread function (PSF) core, wings, and diffraction spikes from the spiders. Then this flux density is corrected by the transmission of the neutral density filter used to obtain the non-coronagraphic image, and by the ratio between the detector integration time of the coronagraphic and non-coronagraphic images, to obtain a reference conversion value. To convert the coronagraphic image from ADU to mJy arcsec⁻², the coronagraphic image is divided by the reference conversion value, multiplied by the stellar flux density of HD 131488 found to be 6.9 Jy at the central wavelength of the H band and divided by the pixel surface area in arcsec². The pixel scale of IRDIS is 0.01225 arcsec per pixel (Maire et al. 2016). The image in Fig. 1 was not corrected by the throughput of the algorithm, which requires a disc model.

For the IFS data (Fig. 2), the master cube consists of temporal 60 frames and 39 spectral channels. We binned the spectral channels in three broader channels centred around 1.04, 1.18, and 1.29 µm, with
<figure><img src="image_3.png" /><figcaption>**Figure 2.** SPHERE/IFS data of HD 131488 reduced with PCA and binned in the following three spectral channels (from left to right): λ = 1.04, 1.18, 1.29 µm. Top: Surface brightness maps. The surface brightness is given in mJy arcsec⁻². North is up and east is to the left. Bottom: SNR maps of the respective surface brightness maps.</figcaption></figure>

a width of 0.16, 0.12, and 0.09 µm, respectively. We reduced each spectral channel independently with a PCA algorithm. The disc is clearly detected in each of those three spectral channels.

## 2.1 Radial profiles

The ADI and PCA reduced H23-band images (Fig. 1) clearly show the detected debris disc of HD 131488 between a radial distance of 0.14 arcsec (22 au) and 0.58 arcsec (90 au) in both eastern and western directions. Both reduction methods lead to similar results. To derive the position angle (PA) and inclination of the disc we used forward modelling applying the Henyey–Greenstein (HG) approach (see Section 4.3.1 for details). We found the PA of the disc in H23-band to be (97 ± 2)° and the inclination to be (84°±1.5)°.

By fitting the measured ALMA visibilities using a Gaussian ring model Moór et al. (2017) obtained comparable parameter values in thermal emission: PA = (96 ± 1)° and i = (82 ± 3)°. That study finds the maximum of the surface brightness at a distance of (88 ± 3) au, and a total disc width of (46 ± 12) au. The location of the peak of surface brightness is similar to the result of our scattered light observations from SPHERE (~90 au).

We extracted the radial profile of the surface brightness using the same method as described in Choquet et al. (2017) where H-band data of 49 Cet are analysed. In this method we produce slices along the semimajor axis with a length of 3 pixels above and below that axis and a width of 2 pixels. The length was found to provide the best SNR while covering the complete vertical disc extent. Then we calculate the mean value of the flux density for each slice. We estimate the noise level of the images by generating similar slices as for the radial profile itself, but along a line perpendicular to the disc’s semimajor axis. Thus, the slices are located outside of the disc. Then we calculate the standard deviation of each slice. The result is shown in Fig. 3.

Our observations reach an average disc SNR level of 4 using an ADI reduction between 22 and 80 au. This is a stronger detection than using a PCA reduction with a SNR of 3 for the same region. This is caused by a more aggressive reduction process of PCA leading to more over-subtraction of the disc and a lower SNR. Within a radius of 22 au the noise level is of the same order of magnitude as the disc signal. Thus, we will exclude the inner region from further analyses. The right panel of Fig. 3 suggests a possible detection beyond 90 au, especially in the western direction (~130 au). However, in this region the SNR is low so that the actual extent beyond 90 au remains uncertain.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left panel: Surface brightness as a function of radial distance to the star for the SPHERE observations in PCA reduction. The blue solid line shows the western part of the disc, the red dashed line the eastern part. Blue and red shaded areas show the 1σ noise level. The grey filled area shows the region where signal and noise are similar. Vertical black dashed lines give the location of the planetesimal belt at 88 au. On the western side a tentative detection of scattered light can be found up to ~130 au. Right panel: Same as left panel, but multiplied by radial distance squared.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Contrast curve for observations of HD131488 inferred from the ANDROMEDA code.</figcaption></figure>

## 2.2 Presence of planets

In order to specifically look for point-like sources, such as exoplanets, we made use of the ANgular Differential OptiMal Exoplanet Detection Algorithm (ANDROMEDA; Cantalloube et al. 2015), as implemented in the High-Contrast Data Centre (Delorme et al. 2017), which utilizes ADI and an inverse problem approach based on a maximum-likelihood estimator. It performs a pair-wise subtraction of frames with different rotation angles, models the expected signature that a planetary signal would leave in the residual image (using the off-axis PSF taken before and/or after the observing sequence) and tracks this signal within the pairs of residual images. We used a minimum rotation angle $\delta_{min}$ of $1\lambda/D$ between frames within a pair to limit self-subtraction, as recommended in Cantalloube et al. (2015). Some bright disc signal remain in the final ANDROMEDA SNR map (maximum SNR of about 9 for the front part of the ring) because the ring is narrow and its signal may appear point-like. Besides those, there is no point source above the $5\sigma$ contrast threshold shown in Fig. 4.

## 3 THEORETICAL BACKGROUND

In order to analyse the scattered light and thermal emission data of HD 131488, we make use of the MODERATO code (Wyatt et al. 1999; Lee & Chiang 2016; Olofsson et al. 2019; Pawellek et al. 2019b) calculating the orbits of dust particles influenced by

stellar gravity, radiation pressure, and collisional evolution. From their position within the disc the code infers the grains' flux density and generates disc images which can be compared to the actual observational images.

In Section 3.1 we explain the theoretical approach to calculate the particle orbits. In Sections 3.2–3.4 we describe the optical models and parameters, including dust compositions, used to infer the flux density of the dust. Section 3.5 focuses on the size distribution of the dust that is influenced by collisional forces and transport processes such as radiation pressure or gas drag. Finally, in Section 3.6 we show the resulting images.

### 3.1 Orbit parameters

The orbits of dust particles are altered by a number of processes, such as collisions or Poynting–Robertson drag. One of the strongest mechanisms for early-type stars is stellar radiation pressure shaping the overall dust grain distribution of a debris disc hosted by such a star. It is characterized by the radiation pressure parameter, $\beta$, which is defined as the ratio between the radiation force and stellar gravity (Burns, Lamy & Soter 1979). For arbitrary particles (including monomers and agglomerates), $\beta$ is given by

$$ \beta = \frac{|\vec{F}_{\text{rad}}|}{|\vec{F}_{\text{g}}|} = \frac{1}{4\pi G c} \times \frac{L_{\text{star}}}{M_{\text{star}}} \times \frac{\sigma_{\text{grain}} \bar{Q}_{\text{pr}}}{m_{\text{grain}}} \quad (1) $$

where $L_{\text{star}}$ and $M_{\text{star}}$ are the stellar luminosity and mass, $G$ the gravitational constant, $c$ the speed of light, $\bar{Q}_{\text{pr}}$ the radiation pressure efficiency averaged over the stellar spectrum, and $\sigma_{\text{grain}}$ and $m_{\text{grain}}$ the particle cross-section and mass. In this study, we will focus on spherical particles (including porous grains) so that $\beta$ can be calculated by

$$ \beta = \frac{3}{16\pi G c} \times \frac{L_{\text{star}}}{M_{\text{star}}} \times \frac{\bar{Q}_{\text{pr}}}{\varrho s (1-P)}, \quad (2) $$

where $\varrho$ is the bulk density of the dust material, $s$ is the grain radius (referred to as size), and $P$ is the porosity of the dust material (Kirchschlager & Wolf 2013). The parameter $\bar{Q}_{\text{pr}}$ is given as

$$ \bar{Q}_{\text{pr}} = \frac{\int \bar{Q}_{\text{pr}}(s, \lambda) F_{\lambda} \, d\lambda}{\int F_{\lambda} \, d\lambda}, \quad (3) $$
where $Q_{pr}$ describes the radiation pressure efficiency depending on $s$, wavelength, $\lambda$, and the stellar flux density, $F_{\lambda}$. For each grain and wavelength, $Q_{pr}$ depends on the absorption and scattering efficiencies, $Q_{abs}$ and $Q_{sca}$, respectively as well as the asymmetry parameter, $\langle \cos(\vartheta) \rangle$ (also called $g$), and is calculated by

$$Q_{pr}(s, \lambda) = Q_{abs}(s, \lambda) + Q_{sca}(s, \lambda) [1 - \langle \cos(\vartheta) \rangle(s, \lambda)]. \quad (4)$$

The asymmetry parameter depends on the scattering angle, $\vartheta$ and is calculated following Bohren & Huffman (1983),

$$\langle \cos(\vartheta) \rangle = g = \int_{4\pi} p \cos(\vartheta) d\Omega, \quad (5)$$

with $p$ being the phase function and $\Omega$ the solid angle. In total intensity, the phase function is given by $S_{11}$ of the Müller matrix (Bohren & Huffman 1983).

Knowing $\beta$, the orbital parameters of the dust grains can be inferred using the equations from Wyatt et al. (1999). Assuming that a dust particle is released from a planetesimal which possesses the orbital parameters semimajor axis, $a_p$, eccentricity, $e_p$, and true anomaly $f_p$, the orbit parameters of the dust grain (semimajor axis $a_d$ and eccentricity $e_d$) can be calculated by

$$a_d = \frac{a_p(1-\beta)(1-e_p^2)}{1-e_p^2-2\beta(1+e_p\cos(f_p))} \quad (6)$$

$$e_d^2 = \frac{\beta^2 + e_p^2 + 2\beta e_p \cos(f_p)}{(1-\beta)^2}. \quad (7)$$

From equation (7) we see that the particle's eccentricity reaches a value of larger than one when

$$\beta \ge \frac{1+e_p}{2[1+e_p \cos(f_p)]}.$$

Assuming that $e_p$ equals zero, i.e. the planetesimals possess circular orbits, this means that particles with $\beta \ge 1/2$ are expelled from the stellar system on either parabolic ($\beta = 1/2$) or hyperbolic orbits. The particle size where $\beta = 1/2$ is then called *blowout limit*. In the special case of $\beta \ge 1$, the trajectories of the particles become anomalous hyperbolas for which $e_d \le -1$ (Wyatt et al. 1999; Krivov, Löhne & Sremčević 2006).

We emphasize that the particles' orbit parameters (semimajor axis and eccentricity) are determined not only by the bulk density of the material, but also by the optical properties of the grains (e.g. absorption efficiency). There are different methods which can be used to infer those optical properties. The three of them used in this study will be discussed in the following section.

## 3.2 Scattered light models

The most common approach to calculate isotropic thermal emission is Mie theory (Mie 1908; Bohren & Huffman 1983) where the particles are assumed to be compact spheres. Due to the isotropy of thermal emission, the particles' shape is of no significant importance and disc models are usually in good agreement with observational data (e.g. Matrà et al. 2018; Moór et al. 2020; Pawellek et al. 2021). This looks different for scattered light data where we have to take into account the shape of the dust grains. Here the approach with Mie theory often leads to poor modelling results for debris discs, likely because the grains do not possess spherical shape (Pawellek et al. 2019b; see Section 3.2.2 for details).

Thus, alternative models are applied. A common approach is the HG approximation (Henyey & Greenstein 1941) which does not include any information on the shape of the grains, but there are

other methods as well such as the Discrete Dipole Approximation (DDA; Purcell & Pennypacker 1973). We will now introduce HG (Section 3.2.1), Mie (Section 3.2.2), and DDA (Section 3.2.3) as approaches to model scattered light observations.

### 3.2.1 HG approximation

The HG approach is used to calculate the scattering phase function, $p$, of the dust material assuming a simple analytical equation,

$$p(\vartheta) = \frac{1}{4\pi} \frac{1 - \langle \cos(\vartheta) \rangle^2}{[1 + \langle \cos(\vartheta) \rangle^2 - 2 \langle \cos(\vartheta) \rangle \cos(\vartheta)]^{3/2}}, \quad (8)$$

where the asymmetry parameter, $\langle \cos(\vartheta) \rangle$, is fixed to a certain value between -1 (back scattering) and 1 (forward scattering). Isotropic scattering implies $\langle \cos(\vartheta) \rangle = 0$. Applying this model to scattered light observations of debris discs, some studies infer the best-fitting asymmetry parameter to derive the general scattering properties of the dust material (e.g. Schneider et al. 2006; Millar-Blanchaer et al. 2015; Olofsson et al. 2016, 2020; Engler et al. 2017) which helps identify possible dust compositions. Other studies assume isotropic scattering to fit larger samples of discs (e.g. Esposito et al. 2020), in order to infer general scattering properties.

The HG approximation usually considers the bulk scattering properties of the dust, rather than considering the behaviour of different grain sizes in the disc, and is usually connected to simple geometric brightness profiles. To improve the HG method, some studies combined grain size distributions with HG properties to model debris discs more realistically (e.g. Esposito et al. 2016; Lee & Chiang 2016; Olofsson et al. 2016). The disadvantage of this approach is its inconsistency. By fixing the asymmetry parameter, $\langle \cos(\vartheta) \rangle$, the $\beta$ parameter (equation 4) is altered. This is because $\beta$ depends on the radiation pressure efficiency, $Q_{pr}$, which depends on $\langle \cos(\vartheta) \rangle$. Thus, a fixed HG parameter leads to a change in the spatial distribution of the dust (see appendix in Pawellek et al. 2019b). However, while this mixed approach does not provide reliable grain size information, it does allow the spatial dust distribution to be readily extracted from scattered light images.

Attempts have been made to scale $\beta$ correctly without taking into account any size information or optical properties of the particles (e.g. Adam et al. 2021; Olofsson et al. 2022). The advantage is that there is only a small number of free parameters to model, but any information on possible dust compositions remain unused. So far, a self-consistent calculation of optical (scattered light) and dynamical parameters (particle distribution, equations (6) and (7) is not possible with the HG approach. We note that this method is a good approach to analyse the material phase function when not focusing on individual particles though.

### 3.2.2 Mie theory

A solution to overcome the difficulties of the simple HG approximation is to apply a scattering model which includes different particle shapes e.g. Mie theory assuming compact spherical grains, or the hollow spheres model (Min, Hovenier & de Koter 2005).

While Mie theory is easy to implement into a code, it has the disadvantage of overestimating the forward scattering observed especially for large (tens of micron-sized) grains (e.g. Schuerman et al. 1981; Bohren & Huffman 1983; Weiss-Wrana 1983; Mugnai & Wiscombe 1986; McGuire & Hapke 1995). To circumvent this, it is possible to exclude large grains from the models as done in Pawellek et al. (2019b) where the maximum size included in the scattered
light model was fixed to 10 µm. However, it is possible that those
grains still contribute to the overall flux density of the debris disc (see
Section 3.5 for details) and thus, alter the results of the modelling.
This leads to an optimization problem where we need to find the
best maximum size so that all contributing particle sizes are taken
into account, and at the same time the overestimation of the forward
scattering is minimized.

A possibility to lower the forward scattering is to change the dust
material. Here we can apply the Effective Medium Theory (EMT)
using Bruggeman's mixing rule (Bruggeman 1935, 1936) to generate
mixtures of different sorts of dust. While Mie theory assumes
compact spheres as particles, we can simulate porous material with
EMT by generating inclusions of vacuum within the matrix of dust
grains. Then the usual Mie calculations can be applied.

3.2.3 DDA

Another way is to use a more complex model e.g. the DDA where
the optical properties of the grains are calculated by assuming that
a particle can be described by a spatial distribution of *N* discrete
polarisable dipoles (Purcell & Pennypacker 1973; Draine 1988). With
this method the particle shape is not limited to that of a simple sphere
but can represent nearly any arbitrary structure including porous
agglomerates or fluffy particles.

The DDA method is highly flexible and can accommodate a
huge variety of particle shapes. Thus, a number of free parame-
ters needs to be introduced, e.g. the dust composition, the grade
of porosity, particle shape, etc. In general, these free parameters
are barely constrained, however, in combination with information
from comets in our own Solar system it is possible to make
reasonable assumptions on those parameters (Section 3.3). While
the advantage of using DDA to model debris discs in scattered
light is evident – creating a self-consistent model of dynami-
cal and optical properties – the main disadvantage is its limited
applicability to particles of large grain size to wavelength ratio
(Draine & Flatau 2010). This limit defines a maximum grain size
of ≲10 µm for a wavelength around ~1 µm (Kirchschlager & Wolf
2013). Another caveat for DDA are highly conducting materials
that we will not take into account in this study (Michel et al.
1996).

So far, DDA is rarely used to model scattered light of debris
discs. Kirchschlager & Wolf (2013) investigated the influence of
grain porosity on the particles’ optical properties using DDA. They
found that the blowout size significantly increases for porous particles
compared to compact grains. The study presented by Brunngräber
et al. (2017) showed that the minimum grain size and the slope
of the grain size distribution are significantly overestimated when
modelling debris discs composed of porous dust with a disc model
assuming spherical, compact grains. Arnold et al. (2019) proved that
the blowout sizes of agglomerated particles and spherical grains
significantly differ but that the dust composition also plays an
important role (see Section 3.3 for details).

A few theoretical studies investigate the influence of porosity and
irregularity on particles' optical properties (e.g. Blum & Wurm 2008;
Kirchschlager & Wolf 2013, 2014; Ysard et al. 2018), while analysis
of protoplanetary discs assume non-spherical particles (e.g. Pinte
et al. 2008; Birnstiel et al. 2010; Ricci et al. 2012; Min et al. 2016).
There are also studies analysing the scattered light coming from
comets in our Solar System using DDA (e.g. Zubko 2013). However,
so far there is no direct application of DDA to actual debris disc
observations which we want to address in this study.

3.3 Dust composition

While the HG approach does not take into account individual dust
compositions such as astrosilicate, water ice, or carbon, the applica-
tion of Mie theory or DDA makes it possible to choose appropriate
compositions freely. Indeed, there is a whole zoo of possible materials
(e.g. Zubko et al. 1996; Henning & Mutschke 1997; Li & Greenberg
1998; Draine 2003; Jäger et al. 2008; Mutschke & Mohr 2019) which
makes it necessary to make assumption for the composition based
on Solar system data or spectra of debris discs.

Considering spectra, outer planetesimal belts observed with in-
struments like Spitzer/IRS typically do not reveal any solid state
features that would allow the dust composition to be constrained
(Chen et al. 2006). The reason is that the dust grains are usually
either too cold or too large to generate visible spectral features.
Of course, there are exceptions for debris discs with hot dust
components like HD 172555 (Chen et al. 2006), HD 36546 (Lisse
et al. 2017), or HD 145263 (Lisse et al. 2020) where silica and
carbon-rich material was found. On the other hand, studies of
comets and asteroids in the Solar system showed that the solid
material in our own debris disc often possesses high porosities
with values of *P* ~ 50 per cent (e.g. Fulle et al. 2015; Sakatani
et al. 2021). We note that it remains debatable whether the dust
composition in the inner region (element abundances and porosity)
is similar to that in the outer region due to different material
processing.

In terms of our study we want to combine dynamical and optical
properties in one model to infer grain sizes. The dust composition
is not the main focus here since no spectra or polarimetric data are
available for HD 131488 that would allow conclusions on the solid
material. Thus, reliable size information is more important. Arnold
et al. (2019) found that for absorptive particles like pure amorphous
carbon grains, porous spheres produce much larger blowout sizes
than dust agglomerates, while for weakly absorbing, pure silicate
grains, porous spheres produce slightly smaller blowout sizes than
agglomerates.

To explore the potential of DDA models in comparison to Mie
theory, we will use the porous grain model of Kirchschlager &
Wolf (2013) and consider particles of a basic spherical shape with
inclusions of vacuum to reflect grain porosity. We assume the
diameter of the vacuum inclusions ('voids') to be as large as 1/100
of the grain diameter, apply astronomical silicate (Draine 2003) and
vary the porosity of the dust grains. This approach is comparable to
that of Arnold et al. (2019) which also uses a basic spherical shape
for the particles. However, Arnold et al. (2019) uses different void
sizes so that the particle structure becomes more complex. With their
irregular shape the orientation of the grains within the debris disc also
becomes more important. In Appendix A1 we show the influence of
different void sizes and spatial distributions on the scattering phase
function.

We decided against a variation of the voids for our scattered light
models of the debris disc around HD 131488 as this would open
several dimensions in parameter space (size and spatial distribution of
the voids, three directions for the orientation of each particle) making
the modelling of the disc complicated and expensive in computational
time. Furthermore, from a statistical point of view we would assume
that while the orientation of irregular particles is important, averaging
over all of them would lead to optical properties similar to a (nearly)
spherical particle rendering the computational effort moot. However,
we note that based on Appendix A1 a dust composition using a
different void size might lead to different results than inferred in this
study.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Radiation pressure parameter, $\beta$, as a function of grain size, $s$, for different scattering models. The horizontal dotted line shows the blowout limit at $\beta = 0.5$. The solid lines show the results assuming EMT and Mie theory, the dashed lines assuming DDA.</figcaption></figure>

## 3.4 Comparing HG, Mie and DDA

### 3.4.1 Radiation pressure parameter

In Fig. 5 we compare $\beta$ (equation 2) for the case of HD 131 488 using Mie theory and DDA assuming a stellar luminosity of 13.9 $L_{\odot}$ and mass of 1.8 $M_{\odot}$. We apply the same stellar spectrum as used in our modelling of HD 131 488 (see Section 4.1 for details). We assume particles with a basic spherical shape and varying porosity. In the case of Mie theory, we mix the refractive indices of astrosilicate with those of vacuum to generate porous material (Section 3.2.2). In the case of DDA, we generate spherical particles of astrosilicate with voids of vacuum and calculate the optical properties directly without using any mixing rules. Hence, applying a porosity of 0 leads to comparable results for DDA and Mie (black lines in Fig. 5).

With increasing porosity the blowout size increases as well for both DDA and Mie theory. Considering HD 131 488 and assuming Mie theory to calculate the absorption and scattering efficiencies, the blowout sizes vary between 2.9 µm for compact grains and ~ 11 µm for particles with a porosity of 0.8. The increase of blowout size with increasing porosity is in agreement with results from other studies

(e.g. Kirschslager & Wolf 2013; Pawellek & Krivov 2015; Arnold et al. 2019). The blowout sizes of compact grains inferred from DDA and Mie show small differences which are due to different calculation methods and set-ups as well as averaging $Q_{pr}$ over a limited number of wavelengths.

The differences between Mie theory and DDA get more pronounced for sub-blowout grains and for larger porosities. DDA predicts smaller $\beta$-values than Mie theory which will influence the amount and orbits of sub-blowout particles present in the disc. We showed in Appendix A1 that different sizes and spatial distributions of voids can influence the scattering phase function of particles. This also leads to changes in the $\beta$ parameter as it depends on the optical properties of the material (see equations 2 and 4).

### 3.4.2 Scattering properties

Another important aspect are the scattering properties of grains. Fig. 6 shows the phase function $p(\theta)$, also $S_{11}$ from the Müller matrix, as function of scattering angle, $\theta$, for different particle sizes assuming a porosity of zero (left panel) and 0.4 (right panel), and for HG following equation (8) with $(\cos(\theta)) = 0.5$.

We see that in both panels the peak of $S_{11}$ at $\theta \sim 0.0$ increases towards larger sizes while the HG approach does not show this behaviour (it is grain size independent). The peak is the aforementioned strong forward scattering for big particles. Comparing EMT (solid lines) and DDA (dashed lines), the phase functions look similar in the case of compact spheres. This is expected since Mie theory can be viewed as a limiting case for both EMT and DDA when assuming compact particles rather than porous grains.

In the case of $P = 0.4$ the phase functions also look similar for (sub-)micron-sized particles indicating that EMT leads to similar scattering properties as DDA when assuming basic spherical particles with small void sizes. However, for larger grains (~10 µm) and scattering angles of $\theta \gtrsim 70^{\circ}$ the deviations of DDA and EMT become more pronounced.

## 3.5 Grain size distribution

For N-body dust models both in scattered light and thermal emission we need to define a size distribution including a minimum and a maximum size of particles that are present in the debris disc. A

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Phase function as a function of scattering angle for different grain sizes and porosities: $P = 0.0$ (left), $P = 0.4$ (right). Solid lines show results for EMT, dashed lines for DDA. The black dash-dotted line represents HG with $(\cos(\theta)) = 0.5$. Grains with $s = 0.1$ µm are shown in red, $s = 1.0$ µm in green, and $s = 10$ µm in blue.</figcaption></figure>
typical size distribution follows a power law,

$$N(s) \mathrm{d}s = N_0 \left(\frac{s}{s_0}\right)^{-q} \mathrm{d}s \quad (9)$$

for grains on bound orbits where $N_0$ and $s_0$ are normalization constants, and $q$ the size distribution index (see Section 4.4) usually set to 3.5 (Dohnanyi 1969). However, when taking into account collisional evolution we find that there is an overabundance of bound grains close to the blowout limit (e.g. Strubbe & Chiang 2006; Thébault & Wu 2008). This can be explained by the fact that smaller grains become unbound and leave the system. Thus, they cannot act as projectiles to destroy the larger bound grains. We take this into account and apply the correction factor introduced by Strubbe & Chiang (2006), $f(e_d) \propto (1 - e_d)^{-3/2}$.

### 3.5.1 Minimum size

While protoplanetary discs are often modelled applying (sub-)micron-sized dust grains that are coupled to the gas (e.g. Szulágyi et al. 2019; Vorobyov et al. 2022), the situation is different in typically gas-depleted debris discs where radiation pressure strongly affects the smaller dust particles that are no longer coupled to the gas. As shown in Section 3.1, grains with $\beta \ge 1/2$ are unbound and expelled from the stellar system on very short time-scales. However, models of scattered light data show that sub-blowout grains are often necessary to fit the debris disc data (e.g. Thebault & Kral 2019).

Assuming that the debris disc is in a quasi steady state, i.e. the production and destruction rates of grains due to collisions are equal, we can estimate the abundance of sub-blowout grains applying the collisional model from Wyatt et al. (2007b). The idea is that only unbound (sub-blowout) grains can get lost from the disc so that the total mass loss rate of the dust equals the production rate of the sub-blowout particles. The mass loss rate is given by

$$\dot{M} = \frac{M_{\mathrm{dust}}}{\tau_{\mathrm{max}}}, \quad (10)$$

where $\tau_{\mathrm{max}}$ is the collisional lifetime of the largest grains considered in our model which also contains most of the mass in the distribution $M_{\mathrm{dust}}$ (for details see e.g. Wyatt et al. 2007b; Löhne, Krivov & Rodmann 2008). The mass of the unbound grains seen at the moment the disc image is taken has to be equal to the mass loss rate so that we can use a normalization constant C,

$$C = \frac{\dot{M}}{M_{\mathrm{dust}}} = \frac{1}{\tau_{\mathrm{max}}}. \quad (11)$$

So far, we assumed that the number of grains follows equation (9) without taking into account blowout limits or mass loss rates. Now, we get a corrected number of grains $C N_i(s)$ for sub-blowout grains of the ith size based on the production rate not following equation (9) anymore. Finally, using the orbital information on each grain (e.g. mean anomaly) we can infer the number of particles produced at each location in the disc.

We investigate the influence of radiation pressure on the minimum grain size by comparing the outcome of the radiation pressure model (RP) to a model that ignores effects of radiation pressure (non-RP). We assume a size distribution of grains between 0.1 and 1000 µm and a total dust mass of those grains of 0.1 M⊕ in both cases. The mass loss rate and thus, the production rate of the sub-blowout particles depends on the dynamical excitation of the disc, i.e. the proper eccentricity of the planetesimals. We note that the planetesimal belt as a whole can exhibit an eccentricity of zero while individual planetesimals can deviate from a circular orbit. Only when they do,

destructive collisions are possible. To estimate the mass loss rate in the RP model we assume a proper eccentricity $\langle e \rangle$ of the colliding planetesimals of 0.1 comparable to the classical Edgeworth–Kuiper belt (e.g. Elliot et al. 2005; Vitense, Krivov & Löhne 2010) and a dispersion of inclinations $\langle i \rangle$ of 0.1 for a central radius of 88 au (Section 2.1) following the approximation $\langle e \rangle \approx \langle i \rangle$ from Wyatt et al. (2007a).

Fig. 7 shows the influence of radiation pressure on the size distribution and the total flux density for a dust disc made of compact, spherical grains ($P = 0.0$) assuming astronomical silicate as dust composition and a dust mass of 0.5 M⊕. The total flux density seen by the observer is calculated by

$$F_{\mathrm{v}}(s) \mathrm{d}s = N(s) F_{\mathrm{v,star}} \left(\frac{R_{\mathrm{star}}}{d_{\mathrm{star}}}\right)^2 \left(\frac{s}{2r}\right)^2 p(\theta) Q_{\mathrm{sca}} \mathrm{d}s, \quad (12)$$

where $R_{\mathrm{star}}$ and $d_{\mathrm{star}}$ are the stellar radius and the distance to the observer, respectively, $F_{\mathrm{v,star}}$ is the flux density of the star at $R_{\mathrm{star}}$ and at the observational wavelength $\lambda$, $r$ is the distance of the particle from the star, and $\theta$ is the scattering angle.

In the left panel of Fig. 7 the total cross-section of particles is shown as a function of grain size. In the non-RP model (dashed line) the total cross-section increases towards smaller sizes following the power law from equation (9). Thus, the smallest grains possess the largest cross-section in this model (grey shaded area), and contribute the major part of the total flux density (also grey shaded area in the right panel). More than 95 per cent of the total flux density come from particles of the assumed sub-blowout size in the non-RP model. In the right panel, we also see that the contribution of the smallest grains ($s \le 0.5$ µm) to the total flux density decreases again due to the decreasing scattering efficiency of the particles for which $\lambda > s$. The same effect is visible in the RP model.

Here, the sub-blowout grains leave the system on (anomalous) hyperbolic trajectories, and are re-produced by destructive collisions of larger bodies. As a result, the number of those grains is much smaller than that of bound grains. We see this effect as a steep decrease of the cross-section in the left panel of Fig. 7 for grains between ~1 and ~3 µm for which $1/2 \le \beta \le 1$ (red shaded area). For smaller grains the total cross-section increases again due to the power law distribution given by equation (9). Since their total cross-section is much smaller, the contribution of the small particles to the total flux density is significantly smaller compared to the non-RP case. Compared to the 95 per cent of the flux density coming from sub-blowout grains in the non-RP case, their contribution is only ~10 per cent in the RP model (red shaded area). This is in agreement with the results from Thebault & Kral (2019) which investigated the influence of sub-blowout grains on disc modelling results. The fraction of particles with $\beta > 0.5$ is still large enough that we should not exclude them completely, and therefore we will not fix our size distribution to the blowout limit but to a size of 0.1 µm to account for the presence of those sub-blowout grains. As can be seen in Fig. 7, the contribution of grains smaller than 0.1 µm is negligible due to the small scattering efficiency of those particles at a wavelength of $\lambda = 1.6$ µm.

### 3.5.2 Maximum grain size and size distribution index

Knowing the maximum size is important since the DDA method limits us to grain sizes smaller than 10 µm due to the number of dipoles necessary to calculate the optical parameters (Kirschlager & Wolf 2013) at an observational wavelength of ~1 µm. Here, we already assume a sphere as simplified basic grain shape and add small
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Influence of radiation pressure and collisions on the grain size distribution and the total flux density assuming a size distribution with $q = 3.5$, a porosity of $P = 0.0$, and a dust mass of $0.5 M_\odot$. Left panel: Total particle cross-section as a function of grain size; right panel: total flux density as function of grain size at $\lambda = 1.6 \mu m$. Red dashed lines indicate grains with $\beta = 1/2$ and $\beta = 1$; black dashed line: model without radiation pressure; black solid line: model including radiation pressure. Grey and red shaded areas: area with grains < blowout limit for both models. Blue shaded area: contribution to total flux density < 1 per cent assuming the case of radiation pressure.</figcaption></figure>

vacuum inclusions in order to represent porous grains (Section 3.3). More complicated shapes would lead to much higher and unfeasible computational times (Arnold et al. 2019) due to a larger number of free parameters. Thus, we need to consider if we need to fill up our size distribution by adding grains > 10 µm that were calculated using Mie theory since they might still significantly contribute to the total flux density of the disc.

Modelling cometary dust using DDA, Zubko (2013) found that grains larger than $15 \times \lambda/(2\pi)$ do not significantly contribute to the estimates of back scattering or geometric albedo, but that those particles increase the computational time. For HD 131488 this would mean to exclude all grain sizes larger than ~4 µm (the blowout size for compact spherical particles lies at ~3 µm). However, in our study we are more interested in the particles' contribution to the total flux density including physical mechanisms such as radiation pressure, and thus, the estimate from Zubko (2013) might not be valid in our case. Furthermore, assuming that dust grains are produced in destructive collisions of bigger bodies we need to include larger dust sizes. Pawellek et al. (2019b) inferred the maximum size considered in the disc models of 49 Cet by estimating the width of the forward scattering peak of compact spherical grains (equation 7 therein). Following a similar approach for the disc around HD 131488 we would get a maximum size of ~9 µm which is already more than 2 times larger than the estimate from Zubko (2013).

However, while the scattering efficiency for grains smaller than the observational wavelength decreases, it stays nearly constant for big particles (see appendix of Pawellek et al. 2019b). Therefore, their contribution to the total flux density is determined by their size distribution rather than their scattering properties. We make a rough estimate and assume that the total flux density coming from a certain size of grains $s$ is given by $F_v \times N_0 s^{3-q}$ based on equations (9) and (12). For simplicity we also assume that $Q_{sca}$ and $p(\theta)$ are constant for large particles. Thus, the ratio of flux densities coming from two different sizes $s_1$ and $s_2$ is given by

$$ \frac{F_v^1}{F_v^2} \approx \left( \frac{s_1}{s_2} \right)^{3-q} . \qquad (13) $$

For example, grains with a size ratio of 10 and a size distribution index of 3.5 reach a flux density ratio of 3 i.e. the contribution of

grains of size $s_1$ is only 3 times larger than that of particles $s_2 = 10 s_1$. This is also visible in Fig. 7 where we assume $q = 3.5$. We find that 95 per cent of the total flux density stems from particles smaller than 25 µm, and 99 per cent from particles smaller than 200 µm (blue shaded area in Fig. 7). This means that we need to 'fill up' our size distribution with spherical grains between 10 and 25 µm, or 200 µm, respectively, to include the part of the size distribution that still contributes significantly to the total flux density. A similar approach was considered for protoplanetary discs (Min et al. 2016). In terms of the size parameter which is defined as $x = 2\pi s/\lambda$, and assuming an observational wavelength of $\lambda = 1.6 \mu m$ this means values of either ~100 (for $s = 25 \mu m$) or ~800 (for $s = 200 \mu m$) as upper limits for the size distribution when studying debris discs compared to a value of 15 suggested by Zubko (2013) when studying cometary tails.

As shown by equation (13), the contribution of large particles to the total flux density depends on the size distribution index, and thus, the maximum size also depends on this parameter. In case of $q = 3.5$ (Fig. 7), we assume an ideal collisional cascade with constant impact velocities and material strength (Löhne 2020). However, SED modelling of debris discs showed that $q$ often differs from this value. It varies mostly between 3 and 4 (e.g. Pawellek et al. 2014, 2021; Löhne 2020). If $q = 3$, all grain sizes with $\beta < 0.5$ contribute the same total cross-section to the distribution (red dashed line in Fig. 8) which results in an equal contribution of flux density based on equations (12) and (13). In this case, the definition of a maximum grain size is rather difficult since larger grains still contribute significantly to the total flux density. For $q = 4$ each size contributes the same mass rather than cross-section, and thus, the total flux density is dominated by particles close to the blowout limit (blue dotted line in Fig. 8). Here, the maximum grain size could be even lower than the 200 µm suggested for the case $q = 3.5$.

To account for most of the cases we will assume a maximum size of $10^4 \mu m$ which also gives us the opportunity to model the disc's SED with the same size distribution at far-IR wavelengths (see Section 4.4 for details). This means in the case of DDA modelling, we will calculate grains with sizes $s \le 10 \mu m$ using DDA and grains with $s > 10 \mu m$ using Mie theory. We note that a similar approach was considered for protoplanetary discs (Min et al. 2016). As shown in Fig. 5 the blowout sizes for Mie theory and DDA are very close
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Flux density as a function of grain size for different size distribution indices. The dust mass is fixed to 0.5 M<sub>⊕</sub>.</figcaption></figure>

when assuming small void sizes as done in this study. Different void sizes and spatial distributions within the particles as used in Zubko (2013) or Arnold et al. (2019) would make it difficult to fill up the DDA size distribution with Mie grains and thus lead to inconsistencies in the model.

Considering the sub-blowout grains (Section 3.5.1), their contribution to the total flux density is determined by the mass loss rate (equation 10) which, in return, depends on the total dust mass and the collisional lifetime of the largest particles in the cascade. The collisional lifetime decreases with increasing q. Thus, the mass loss rate increases and the contribution of sub-blowout particles to the total flux density increases as seen in Fig. 8.

### 3.5.3 Influence of gas on the size distribution

The disc around HD 131488 was found to possess a high content of gas (e.g. Moór et al. 2017, 2019; Rebollido et al. 2022) that could have an impact on the dust distribution in both space and size. Depending on the surface density of the gas, particles up to a certain size can be dragged efficiently by the gas while larger grains remain unaffected. Thus, strong gas drag might lead to a much higher amount of sub-blowout grains in the disc that are usually expelled from the system. The (dimensionless) Stokes number gives the time-scale necessary to stop a grain from its relative motion towards the gas. It is given by

$$St = \frac{\pi s q}{2 \Sigma_{\text{gas}}}, \qquad (14)$$

where $\Sigma_{\text{gas}}$ is the surface density of the gas (e.g. Marino et al. 2020). If $St \ll 1$ the dust particle is stopped nearly instantly and follows the motion of the gas.

Fig. 9 shows the dust particle size as a function of the Stokes number for different gas surface densities. While the surface density of the gas is not well constrained, we can estimate a rough value for the CO gas surface density based on the observations presented by Moór et al. (2017). The study estimated a CO gas mass of $9 \times 10^{-2}$ M<sub>⊕</sub> and a reanalysis of the data constrains the radial extent of the gas to 30–130 au (see Appendix C). These translate to an average CO gas surface density of $2 \times 10^{-6}$ M<sub>⊕</sub> au<sup>-2</sup>.

However, the total gas mass is probably much larger than this. First, the above CO mass estimate should be treated as a lower limit, since it was derived assuming an ISM-like abundance of C<sup>18</sup>O,

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Stokes number as a function of grain size.</figcaption></figure>

which is likely to be an underestimate of the true value due to isotope-selective photo-dissociation (Moór et al. 2017). Moreover, this estimate considers only CO, which may not be the dominant species. If the gas has a residual primordial nature, then the gas composition is dominated by H<sub>2</sub> molecules, whose mass exceeds that of CO by orders of magnitude (Miotello et al. 2023, and references therein). According to current theories, however, it is more likely that the observed gas is of secondary origin and has been released from icy bodies, e.g. via collisions, sublimation, photo-desorption, and/or as an outcome of the thermal evolution of young large icy planetesimals (Kral et al. 2019; Marino et al. 2020; Bonsor et al. 2023). In our Solar system, H<sub>2</sub>O, CO, and CO<sub>2</sub> are the most abundant species in the cometary gas (Mumma & Charnley 2011). While self-shielding and shielding by C atoms can substantially increase the photo-dissociation lifetime of CO molecules, similar mechanisms are not available for CO<sub>2</sub> and H<sub>2</sub>O molecules, which are therefore rapidly dissociated due to UV photons. In order to determine the total gas mass of a CO-rich debris disc, we would therefore need to know not only the mass of CO, but also the amounts of the various photo-dissociation and photo-ionisation products (C, C<sup>+</sup>, O, H) of the main molecules. Although, thanks to ALMA, estimates of the C content of an increasing number of CO-rich discs are available (Cataldi et al. 2023), the amounts of O and H in such discs are not known.

Based on molecular abundances measured in cometary atmospheres in the Solar system (Mumma & Charnley 2011), the mass ratios of above photo-dissociation products to CO can vary over a wide range, with an upper bound of ~25 (taking into account that the rapid photo-dissociation of CO<sub>2</sub> results in CO gas). Although this is subject to a number of uncertainties when applied to the disc of HD 131488 – for example, not only is the composition of the ice bodies there unknown, but also the mechanisms that lead to the gas production, which can result in different gas mixtures for the same ice composition – it can be said that the average gas surface density can be as high as several times $10^{-5}$ M<sub>⊕</sub> au<sup>-2</sup>. Another aspect is that the gas distribution is probably not uniform. It is well possible that the density if gas in the planetesimal belt, which is the main production site for both gas and small dust particles, is higher than elsewhere.

In the case of compact grains, we find that all particles smaller than the blowout size possess a Stokes number significantly smaller than 1 only if $\Sigma_{\text{gas}} \gtrsim 10^{-4}$ M<sub>⊕</sub> au<sup>-2</sup>, while for coupling grains of size $s \sim 0.5 \mu\text{m}$ – that would be the brightest particles if radiation
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Scattered light models of a debris disc generated with the MODERATO code (Pawellek et al. 2019b). Left: Ignoring radiation pressure. Right: Including radiation pressure. The lobes visible in the radiation pressure model can be explained by the scattering phase functions of the dust grains.</figcaption></figure>

pressure was inefficient in expelling sub-blowout grains (Fig. 7) – needs gas surface densities higher than ~$10^{-5}$ M$_{\oplus}$ au$^{-2}$. These surface densities are much higher than the values estimated for CO only, but taking into account the other components, the total gas surface density may be quite similar to them.

## 3.6 Images of disc models

### 3.6.1 Disc appearance

We now analyse the influence of the aforementioned aspects (Sections 3.1–3.5) on the disc images. Fig. 10 shows the effect of radiation pressure for a grain size distribution of compact grains between 0.1 and $10^4$ µm and with a size distribution index of $q = 3.9$ indicating a high fraction of grains at the lower end of the size distribution.

In the left panel all particles possess $\beta = 0$ condition, implying that all grains stay close to their parent bodies on non-eccentric orbits. The close proximity of small grains to the planetesimals might be the case for a large surface density of gas (see Section 3.5.3). While the presence of gas does not change the $\beta$ values, the orbits of the grains will be altered so that the particles stay close to their parent body which can be roughly described with $\beta \approx 0$. In this case the scattered light is dominated by grains around ~0.5 µm (Fig. 7) that are not expelled from the system, i.e. radiation pressure is not efficient. These particles possess more or less isotropic scattering properties at $\lambda = 1.6$ µm (Fig. 6) so that the model disc also shows a more isotropic distribution of scattered light.

In the right panel of Fig. 10 radiation pressure is included, i.e. $\beta > 0$. Radiation pressure is effective when the surface density of gas is low. Grains smaller than the blowout limit are expelled from the system, and only contribute a minor fraction of the dust due to their reproduction by collisions (Fig. 7). Particles close to the blowout limit are moving on highly eccentric orbits forming a halo of roughly bound grains. Only larger grains for which $\beta \ll 0.5$ stay close to the parent belt and dominate the scattered light i.e. the surface brightness at the centre of the belt. The large particles possess strong forward scattering (Fig. 6) and thus, the model disc shows a peak close to the star in the centre of the image where the scattering angle is small. Azimuthal changes in brightness in the right panel are caused by

large particles. A finer grid of grain sizes will smooth out those variations.

### 3.6.2 Total flux density

Radiation pressure not only influences the appearance of the disc, but also affects the total flux density of the model disc which is connected to the disc mass. As shown in Fig. 7, the highest contribution to scattered light comes from small grains with sizes of ~0.5 µm when radiation pressure is not taken into account. These grains do not significantly contribute to the total dust mass and thus, only a small amount of material is needed to generate a high total flux density when such small grains are present within the disc.

If radiation pressure is included, these dominating sub-blowout grains are expelled from the system, and the majority of the scattered light comes from bound particles close to the blowout limit (Fig. 7). If we use the same total dust mass in both cases (RP and non-RP), the flux density of the bound grains is much smaller compared to that of the sub-micron-sized grains. Therefore, we would need a much larger dust mass to generate the same level of flux density when taking radiation pressure into account. This effect is even stronger for porous dust grains. The amount of scattered light coming from a compact grain is larger than the light coming from a porous particle of the same size (Section B).

# 4 MODELLING RESULTS

## 4.1 Stellar photosphere and dust composition

For all our approaches we apply an ATLAS9 model (Castelli & Kurucz 2004) as stellar photosphere to determine the influence of the host star HD 131488. Here, the stellar temperature, metallicity, and surface gravity provided by Rebollido et al. (2018) are taken into account to generate the synthetic spectrum. We assume a dust composition of pure astronomical silicate with a bulk density of 3.3 g cm$^{-3}$, and use porosities of $P = 0.0$ (for compact grains), 0.2, 0.4, 0.6, and 0.8 by using Bruggeman's mixing rule of EMT (Bruggeman 1935, 1936).
<figure><img src="image_12.png" /><figcaption>Figure 11. Best-fitting model assuming HG approximation. From left to right: PCA-reduced SPHERE/IRDIS observations; HG model; Residual image. The color scale is given in mJy arcsec⁻².</figcaption></figure>

## 4.2 Fitting approach

To find the best-fitting model, we use a $\chi^2$-minimization assuming that an ideal residual image should only contain white noise in each pixel. The $\chi^2$-parameter is then computed for each pixel by

$$ \chi^2 = \sum_{i=1}^{N_{\text{pixel}}} \left( \frac{F_{i, \text{residual}}}{F_{i, \text{noise}}} \right)^2 . \qquad (15) $$

The noise is estimated by computing a disc-free image with the same noise distribution as in the PCA-reduced image. This is done by derotating the IRDIS images in the opposite direction, compared to the correct reduction. The faint disc signal present in individual images is therefore diluted when the images are stacked, to produce a final disc-free image containing only residual noise.

We use the $\chi^2$ to estimate the uncertainties of our free parameters. Assuming a confidence level of 95 per cent, we infer the critical $\chi^2$-value for which we need to reject the hypothesis that our model represents the observations. In terms of a reduced $\chi^2$ this means a value of $\le 1.05$ in our case. Based on the best-fitting parameter values we change each free parameter individually until the final $\chi^2$ gets larger than the critical value.

## 4.3 Comparing scattered light models

### 4.3.1 HG model

The HG approach is useful to infer general scattering properties of the dust material by modelling the phase function. Fig. 11 shows the results of this procedure.

We use a simple geometric model with the free parameters PA, inclination, g, central radius $r_0$, and flux density. We apply a two-part power law as radial profile centred at $r_0$. The slopes of the power law were fixed to $\alpha_1 = 12$ and $\alpha_2 = -12$ so that the disc is narrow and $r_0$ corresponding to the peak density of the disc. The best-fitting value for $r_0$ is then found for $(110 \pm 25)$ au but remains rather uncertain. While Moór et al. (2017) found a best-fitting value of $(88 \pm 3)$ au it seems that in scattered light the disc peaks at a larger distance. We will discuss this issue in Section 5.2.

We included an *ad hoc* phase function and find a best-fitting value for the asymmetry parameter $g = \langle\cos(\theta)\rangle$ of $(0.67 \pm 0.07)$, a PA of $(97 \pm 2)^\circ$ and an inclination of $(84_{-2.0}^{+1.5})^\circ$ (see Section 2.1). We assume a zero disc eccentricity. The positive g-parameter indicates a material of forward scattering particles which is comparable to the results of other debris disc studies (e.g. Millar-Blanchaer et al. 2015; Olofsson et al. 2016, 2019; Engler et al. 2017, 2019). It also indicates that the particles might not resemble compact spherical bodies since we would expect an even stronger forward scattering around $g = 0.9$ when applying Mie theory.

### 4.3.2 Mie model

Now, we generate the semidynamical disc models using the MODERATO code. The code assumes surface density profiles for the parent belt following a Gaussian distribution. We use the following free parameters: disc width $\Delta r$, dust mass $M_d$, size distribution index q, and dust porosity P. We fix PA and i to the values inferred from the HG approach to keep the fitting process fast. Each planetesimal in the belt releases grains of different sizes following a power law (equation 9). Sub-blowout grains are produced following the collisional model of Wyatt et al. (2007b). Then the orbits of the individual dust particles are calculated and from their position the light scattered in the direction of the observer is inferred.

For our first model we consider grains calculated by Mie theory that range from compact spheres to high porosity ($0.0 \le P \le 0.8$ in steps of 0.2) and assume a size distribution index smaller than 4 (for $q = 4$ each size bin would contribute the same mass to the total dust mass). The particles are produced in a single axisymmetric ring, and we do only take into account the effect of radiation pressure on the particles' orbits.

We find a best-fitting porosity of 0.6. However, we cannot exclude porosities between 0.0 and 0.4 as they also lead to results within the confidence interval. Assuming a porosity of 0.8 did not lead to a well-fitting model. The best-fitting size distribution index is $q = 3.0 \pm 0.1$. For grains with sizes between $0.1 \mu m \le s \le 10^4 \mu m$ we
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Best-fitting model assuming Mie theory and radiation pressure. From left to right: PCA-reduced SPHERE/IRDIS observations; model generated by MODERATO assuming *P* = 0.6 and *q* = 3.0; forward model (PCA); residual image.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Radial profiles of the planetesimal belt. Dashed line: Profile inferred from ALMA data. Solid line: Profile used to model scattered light data. Both profiles use a Gaussian.</figcaption></figure>

need a dust mass of ~4.5 ± 0.8 $M_⊕$ to fit the scattering flux density observed. In Fig. 12 we show the results for the best-fitting model.

Based on ALMA observations Moór et al. (2017) used a Gaussian radial profile and found a best-fitting value for the central radius of (88 ± 3) au and a total disc width of (46 ± 12) au (after updating the distance of HD 131488 by Gaia data; Section 2.1). In contrast to that we find that a symmetric Gaussian radial profile does not lead to a fit within the expected $\chi^2$ confidence level of the scattered light data, i.e. $\chi^2_{reduced} > 1.05$. We found that in this case the dust at radial distances smaller than 88 au dominate the scattered light and lead to the disc appearing smaller than observed.

We tested different values as inner boundary and found that only when moving the inner boundary to (88 ± 5) au i.e. ignoring the inner part of the Gaussian and only taking the outer part into account, our models could reach a $\chi^2$-value within the appropriate confidence level ($\chi^2_{reduced} \le 1.05$). With this approach we find a best-fitting disc width of (30 ± 3) au (Fig. 13). We will discuss this discrepancy in disc width in more detail in Section 5.2.

### 4.3.3 DDA model

In a similar approach as described in Section 4.3.2 we now apply optical dust properties inferred from DDA. Again we use the free parameters: disc width $\Delta r$, dust mass $M_d$, size distribution index $q$, and dust porosity $P$. We consider grains ranging from compact spheres to high porosity (0.0 ≤ $P$ ≤ 0.8 in steps of 0.2) and assume a size distribution index smaller than 4. The results for the best-fitting model are shown in Fig. 14. Similar to the Mie model we find a best-fitting for $P = 0.6$. Again, we cannot exclude lower porosities ($P = 0.0, 0.2, 0.4$) as they also lead to models within the confidence interval, but with larger $\chi^2$. The higher porosity of 0.8 can be excluded as it leads to best fits outside the confidence interval. The size distribution index is found as $q = 3.0 \pm 0.2$, and the dust mass as $M_{dust} \approx 4.4 \pm 0.7 M_⊕$. Also the asymmetric radial profile with a disc width of 30 ± 3 au is similar to that found in Section 4.3.2 (see Section 5.2 for details).

Based on the scattering phase function and blow-out sizes (Figs 5 and 6) the differences between the Mie and DDA approach for the pure radiation pressure model were found to be minor when assuming grains of basic spherical shape and small sizes of the vacuum inclusions (Section 3.2.3). This is now confirmed by a similar quality of our Mie and DDA models ($\chi^2$ values are similar). This leaves us with the question whether the time consuming DDA approach is useful to model debris discs. We will discuss this question in Section 5.1.

## 4.4 SED

### 4.4.1 Model set-up

We now compare the scattering flux density inferred from the radiation pressure model with the results from modelling the thermal emission of the disc around HD 131488 at longer wavelengths. With the MODERATO code we generate thermal emission images at wavelengths smaller than $10^4 \mu$m, and then calculate the SED. The SED is calculated by the same approach as the SONATA code (Pawellek et al. 2014, 2021; Pawellek & Krivov 2015), but now includes the effect of radiation pressure on particles of different sizes. This approach guarantees that the models for thermal emission and scattered light are self-consistent. The photometric data of the
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Best-fitting model assuming radiation pressure. From left to right: PCA-reduced SPHERE/IRDIS observations; model generated by MODERATO assuming *P* = 0.6 and *q* = 3.0; forward-model (PCA); residual image. The differences between the residuals of the DDA and the Mie model (Fig. 12) are subtle. The Mie model leaves slightly more dark areas in the eastern part of the disc.</figcaption></figure>

<table><caption>**Table 1.** Continuum flux density.</caption><thead><tr><th>Wavelength [µm]</th><th>Flux density [mJy]</th><th>Instrument</th><th>Reference</th></tr></thead><tbody><tr><td>0.42</td><td>2229.66 ± 34.95</td><td>TYCHO B</td><td>1</td></tr><tr><td>0.43</td><td>2558.05 ± 23.47</td><td>APASS B</td><td>2</td></tr><tr><td>0.47</td><td>2671.94 ± 24.51</td><td>APASS G</td><td>2</td></tr><tr><td>0.51</td><td>2045.98 ± 19.39</td><td>Gaia BP</td><td>3</td></tr><tr><td>0.53</td><td>2335.20 ± 30.67</td><td>TYCHO V</td><td>1</td></tr><tr><td>0.54</td><td>2631.62 ± 24.14</td><td>APASS V</td><td>2</td></tr><tr><td>0.62</td><td>2165.84 ± 19.87</td><td>APASS R</td><td>2</td></tr><tr><td>0.64</td><td>2052.12 ± 18.87</td><td>Gaia G</td><td>3</td></tr><tr><td>0.78</td><td>1718.90 ± 16.76</td><td>Gaia RP</td><td>3</td></tr><tr><td>0.79</td><td>1527.60 ± 31.34</td><td>DENIS I</td><td>4</td></tr><tr><td>1.24</td><td>1123.65 ± 36.53</td><td>2MASS J</td><td>5</td></tr><tr><td>1.65</td><td>763.25 ± 33.65</td><td>2MASS H</td><td>5</td></tr><tr><td>2.16</td><td>500.87 ± 15.84</td><td>2MASS Ks</td><td>5</td></tr><tr><td>3.38</td><td>240.65 ± 8.31</td><td>WISE</td><td>6</td></tr><tr><td>4.63</td><td>163.45 ± 5.16</td><td>WISE</td><td>6</td></tr><tr><td>8.98</td><td>164.20 ± 7.16</td><td>AKARI</td><td>7</td></tr><tr><td>12.33</td><td>111.12 ± 4.80</td><td>WISE</td><td>6</td></tr><tr><td>22.25</td><td>153.15 ± 8.80</td><td>WISE</td><td>6</td></tr><tr><td>101.40</td><td>331.20 ± 19.84</td><td>Herschel/PACS</td><td>8</td></tr><tr><td>163.60</td><td>184.80 ± 25.54</td><td>Herschel/PACS</td><td>8</td></tr><tr><td>1322.42</td><td>2.91 ± 0.31</td><td>ALMA</td><td>9</td></tr><tr><td>1652.22</td><td>1.64 ± 0.17</td><td>ALMA</td><td>10</td></tr><tr><td>8750</td><td>0.0595 ± 0.0124</td><td>ATCA</td><td>11</td></tr></tbody></table>

References: [1] Høg et al. (2000); [2] Henden et al. (2016); [3] Gaia Collaboration (2018); [4] DENIS Consortium (2005); [5] Cutri et al. (2003); [6] Wright et al. (2010); [7] Ishihara et al. (2010); [8] Marton et al. (2017); [9] Moór et al. (2017); [10] This work; [11] Norfolk et al. (2021).

dust continuum used to constrain the thermal emission model are given in Table 1.

### 4.4.2 Results

In Fig. 15 we show the SED of HD 131488. In a first approach we use the SONATA code to fit a two component model to the observational data (red and blue dashed lines in Fig. 15), but we note that this code does not take into account radiation pressure. We find a best-fitting model for a dust mass of 1.0 ± 0.2 $M_⊕$ assuming the same porosity

<figure><img src="image_16.png" /><figcaption>**Figure 15.** SED of HD 131488. The purple circles show the measured flux density at different wavelengths. The lines show models using astronomical silicate with a bulk density of $1.3 \text{ g cm}^{-3}$ for $P = 0.6$. The blue and red dashed lines show the two component model done by the SONATA code ignoring radiation pressure. The green solid and the orange dash-dotted lines show the results from the MODERATO code where no second component was taken into account.</figcaption></figure>

(0.6) and size distribution index (3.0) found by our scattered light models.

In a second approach we use the MODERATO code which now includes the effects of radiation pressure. In the scattered light images we do not detect warm dust close to the host star since the coronagraph of the SPHERE instrument is blocking out the inner region. Hence, we do not model a possible Asteroid belt analogue with MODERATO but only the outer Kuiper belt analogue (green solid line and orange dash-dotted line in Fig. 15). When applying the best-fitting parameters inferred from the scattered light model ($P = 0.6, q = 3.0, M_⊕ = 5 M_⊕, r_0 = 88 \text{ au}, \Delta r = 30 \text{ au}$) we get the orange dash-dotted line seen in Fig. 15. This line is not consistent with the observational data at far-IR wavelengths. However, we find that the observational data can be fitted when using a dust mass of 1 $M_⊕$ instead of 5 $M_⊕$ (green solid line in Fig. 15) which is a comparable dust mass as found by the SONATA model.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Best-fitting model assuming radiation pressure and gas drag with a surface density of $\Sigma = 2 \times 10^{-5} M_{\oplus} \text{au}^{-2}$. From left to right: PCA-reduced SPHERE/IRDIS observations; model generated by MODERATO assuming $P = 0.6$ and $q = 3.0$; forward-model (PCA); residual image.</figcaption></figure>

The ATCA point at 8.7 mm shows a higher flux density than predicted by our models. One reason is that at a wavelength of ~1-cm particles with sizes larger than our applied maximum size of $10^4 \text{ µm}$ contribute to the emission. Furthermore a simple power law size distribution might not be valid at long wavelengths. A higher amount of large particles might be present in the disc. This seems to be a common occurrence based on observational data from other debris discs (Lestrade et al. 2020) and data from our own Solar system (e.g. Morbidelli et al. 2021).

As mentioned before, the dust mass necessary to model the thermal emission data is a factor 5 lower than the prediction made from scattered light models only. The dust mass is well constrained by the SED so that we need to find a scattered light model that can reproduce the total flux density with this mass value. This suggests that a higher amount of sub-blowout grains is retained than our pure radiation pressure model predicts. A possible explanation might be the presence of gas that we will investigate in Section 4.5.

## 4.5 Combining SED and scattered light results

As found for the pure radiation pressure model, a dust mass of $4.4 \pm 0.7 M_{\oplus}$ would be necessary to generate the amount of scattered light observed for HD 131488. However, SED models predict a dust mass of only $1.0 \pm 0.2 M_{\oplus}$. A possible explanation for the much higher flux density at short wavelengths might be the presence of gas within the disc that could increase the amount of sub-blowout grains dominating the scattered light. We showed in Section 3.5.3 that a gas surface density of $10^{-4} M_{\oplus} \text{au}^{-2}$ would be enough to strongly couple all the sub-blowout grains to the gas, and that a surface density of at least $\sim 10^{-5} M_{\oplus} \text{au}^{-2}$ seems realistic for HD 131488. We note that this is a rough estimate without uncertainties as we do not have enough data on gas species present in the disc. We now assume that sub-blowout grains up to a certain size are efficiently coupled to the gas and not expelled by radiation pressure. In a simple approach we vary the gas surface density and assume that particles for which the Stokes number is $St \le 1$ (equation 14) are retained within the disc by setting the $\beta$ parameter to zero. We note that in reality $\beta$ is not changed by the presence of gas, but the grains' orbits are.

In Fig. 16 we show the resulting best-fitting model for which the dust masses found in scattered light and thermal emission are equal (1 $M_{\oplus}$). We find that for a size distribution index of $q = 3.0$ we need

a surface density of $\Sigma = (2.0 \pm 0.1) \times 10^{-5} M_{\oplus} \text{au}^{-2}$ to fit both the SED and the scattered light data. As explained in Section 3.5.3 this value is in agreement with expectations/estimates from ALMA observations.

## 4.6 Deriving dust properties

### 4.6.1 Average phase function

In Section 3.4.2 and Appendix A2 we analyse the scattering phase function for different porosities and grain sizes independent of actual debris disc models. We find that for increasing porosity the phase function of specific sizes decreases for larger scattering angles, 9 (Fig. A2). At small angles ($9 < 10^{\circ}$) we see an increase of the phase function for grains with sizes of $s \ge 1 \text{ µm}$ which we mentioned before as *forward scattering*. In Fig. 10 the forward scattering is visible in the right panel where radiation pressure expels grains smaller than the blowout size leading to a dominance of grains with $s \sim 10 \text{ µm}$. In the left panel it is not visible as the size distribution is dominated by 0.1 µm-sized grains that do not show the forward scattering peak in the phase function (Fig. A2).

Now we combine those results on the scattering behaviour of the individual dust particles with the semidynamical disc model and infer an *average phase function* for the debris disc model. To do so we infer the scattering angle and flux density (equation 12) of each dust grain of size $s$ at position $r$ in the disc. Then we multiply the value by the number of grains within the same size bin and at the same location (equation 9) and sum over all particle sizes and distances to get a flux density that only depends on the scattering angle, 9. Finally, we normalize this flux density by the total flux density of the model to get the average phase function.

In Fig. 17 we show the resulting average phase function for our best-fitting models (pure radiation pressure, gas drag) and for different porosities. Similar to the results for individual grains the average phase function shows smaller values at larger scattering angles when the porosity increases. Also, at small angles the phase function increases. Comparing the grain model to the HG model (black-dash double-dotted line in Fig. 17) we find that a porosity of $P = 0.6$ leads to the closest match between Mie, DDA, and HG approach. For $P = 0.0$ the phase function shows a larger contribution at larger scattering angles while for $P = 0.8$ the contribution is
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Average phase function as function of scattering angle. The parameter *P* gives the porosity of the dust composition, HG the phase function assuming the HG approximation. Best-fitting models use *P* = 0.6. All models use *q* = 3.0.</figcaption></figure>

smaller. Thus, the *g* parameter gives a hint on the general porosity of the material being larger than 0.0. The very porous material of 0.8 could be excluded by our scattered light model already.

### 4.6.2 Reflectance

We do not only possess the SPHERE/IRDIS data for the disc around HD 131488 at $\lambda = 1.6$ µm that were used to generate the scattered light models in this paper, but we also have the data available from IFS at $\lambda = 1.04, 1.18, \text{and } 1.29$ µm (see Section 2). This gives us the opportunity to analyse the reflectance, i.e. the fraction of stellar light that is scattered by the dust of the disc at all these wavelengths. The reflectance gives us information on the colour of the debris disc in scattered light. We assume that the colour is determined by the grain size distribution, i.e. the amount of small grains compared to large grains.

To derive the reflectance we inferred the total flux density of the disc from the observational images at five different wavelengths (IFS: 1.04, 1.18, 1.29 µm; IRDIS: 1.593, 1.667 µm). For this we used the HG approach (Section 4.3.1) with *g* = 0.67 and optimized a scaling factor to minimize the residuals in the five spectral channels. For our DDA model predictions we only used four wavelengths (1.04, 1.18, 1.29, 1.6 µm) as the IRDIS wavelengths are very close to each other. We applied the best-fitting values (*q* and *P*) and calculated the expected flux density. Then we divided the total flux density by the respective stellar flux density at the respective wavelength. The stellar flux density of the observations was inferred directly by the PSF while for the models we interpolated the stellar spectrum used in the modelling process (Section 4.1). The result is shown in Fig. 18.

The observations (black squares) show a weak decrease of reflectance with increasing wavelength, but the uncertainties are large so that the colour of the disc – while suggested to be blue – remains uncertain. Recent studies on debris disc colours (e.g. Thebault & Kral 2019; Ren et al. 2023) show that most debris discs seem to be blue in scattered light at wavelengths longer than 1 µm. Ren et al. (2023) studied discs in the optical and near-IR and also found some discs to possess a red colour at wavelengths shorter than 1 µm. Based on the right panel of Fig. 7 we would expect this as sub-micron sized grains sensitively traced at visual wavelengths become less abundant than their still bound micron-sized counterparts.

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Reflectance for the different wavelength bands used by VLT/SPHERE. Black squares show the values inferred from observations including 1$\sigma$ errorbars, orange circles those predicted from the best-fitting model including gas drag. Additional lines show the predictions for models using radiation pressure only and different size distribution indices. These lines were scaled to a value of 0.1 per cent at 1.5 µm.</figcaption></figure>

To get an idea on how the reflectance would change for different size distribution indices we generated models that show the difference for flat ($q = 3.0$), intermediate ($q = 3.5$), and steep ($q = 4.0$) size distributions. We only assumed radiation pressure to influence the grains. Gas drag was ignored. The results are shown as red, green and blue lines in Fig. 18. The scaling for those models was done so that the lines cross 0.1 per cent at a wavelength of 1.5µm to keep the plot readable. It is not connected to the actual dust mass.

We find that for a flat size distribution (red line) there is no change in reflectance with wavelength. This is understandable as for $q = 3.0$ all particles independent of their size contribute the same total cross-section and thus, the same amount of light should be scattered at all wavelengths. For $q > 3.0$ we find the reflectance decreasing with increasing wavelength. The larger *q*, the steeper the decrease becomes as the ratio between small and large particles increases as well. For values of $q < 3.0$ the reflectance is expected to increase since here the total cross-section of big grains is larger than for small particles. Due to the large uncertainties of the slope we find that a constraint on the size distribution based on reflectance is not possible.

The orange line in Fig. 18 shows the result for our best-fitting model including our simple assumptions on gas drag. This model and the observations are scaled correctly. First, we see that our best-fitting model in general predicts a reflectance that is significantly higher than the observational reflectance even taking into account the large uncertainties. This can be explained by the differences in scattering phase function (Fig. 6). For angles smaller than 20° the phase function for the DDA model is larger than the HG phase function by more than a factor of two. The total flux density is dominated by grains with small scattering angles. However, for small angles observations are rather uncertain as for example the PCA reduction leads to strong over-subtraction. Therefore, we do not put too much emphasis on the different levels of reflectance, but are more interested in its slope.

We get a similar slope of radiation pressure model and gas drag model when assuming $q \sim 3.4$ for the former and $q = 3.0$ for the latter. This shows that even for a size distribution index of $q = 3.0$ we can get a decreasing slope and thus, a blue colour of a disc when including the effects of gas drag. However, due to the large uncertainties a comparison with the observations is not reliable. Based on the results
of our SED and scattered light models we assume that the reflectance indicates a slightly blue colour for the disc around HD 131488.

# 5 DISCUSSION

## 5.1 Mie versus DDA

Modelling scattered light observations of debris discs is often difficult as optical dust models such as Mie theory were found to give poor modelling results (e.g. Krist et al. 2010; Milli et al. 2017; Pawellek et al. 2019b). In this paper we analysed the benefit of using DDA to model scattered light observations of debris discs. This approach is not new as other studies investigated circumstellar discs applying DDA (e.g. Min et al. 2016; Arnold et al. 2019, 2022; Audu et al. 2023) already. However, with this study we present the first analysis using a semidynamical disc model and thus, combining optical properties with dust dynamics to generate scattered light models.

Instead of changing the dust material we varied the level of porosity to investigate the resulting disc models in both scattered light and thermal emission. Other studies showed that porosity influences the modelling outcome e.g. the ratio between dominant grain size and blowout size which hints at the disc's dynamical excitation (Pawellek & Krivov 2015; Brunngräber et al. 2017). In this study we stayed comparable to Mie theory when applying DDA. We assumed particles of basic spherical shape. For Mie grains we applied EMT to generate porous material. For DDA we generated small inclusions of vacuum with a size of 1/100 (Section 3.2.3). Since we cannot use DDA for grains with $s \ge 10$ µm, this approach allowed us to fill up the size distribution with Mie particles. To use more complex particle structures in DDA more work is needed to infer possible ways of filling up the size distribution.

We find that when assuming grains of spherical shape and small void sizes, DDA and Mie lead to similar results for pure radiation pressure models. Deviations between the models can be explained by different set-ups (e.g. different blowout sizes, phase functions, $\beta$, etc). This outcome is somewhat expected as Mie theory is a limiting case for DDA when assuming spherical shapes. For the special case of HD 131488 we find that particles of basic spherical shape and small vacuum inclusions can reproduce the observations very well. Thus, Mie theory seems a valid approach to model the scattered light data for this disc. As mentioned before, this seems not the case for many debris discs, although there are studies reaching a similar result (e.g. Ertel et al. 2011).

While the benefit of DDA is not particularly obvious for this study, we emphasize that we are now able of introducing arbitrarily shaped grains or dust aggregates similar to Zubko (2013) or Min et al. (2016) into the MODERATO code to model discs where Mie theory is not a good approximation for the dust particles. However, we note that a transition between DDA and Mie will be necessary to cover the whole grain size range.

## 5.2 Radial extent

### 5.2.1 ALMA versus SPHERE

The MODERATO code uses the location of the planetesimal belt as input to calculate images at different wavelengths. The largest dust grains traced by ALMA are barely affected by radiation pressure or other transport mechanisms (e.g. Pawellek et al. 2019a) so that we can assume the radial extent inferred from ALMA data to reflect the actual planetesimal belt location. Also assuming that the dust grains are produced in mutual collisions within the planetesimal

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Normalized flux density as function of radius inferred by the collisional code ACE for a planetesimal belt of different ages.</figcaption></figure>

belt, the radial extent of the disc at shorter wavelengths should be in agreement with those dust grains put on eccentric orbits due to radiation pressure.

While the radial profile inferred from ALMA favours a Gaussian with a central radius of (88 ± 3) au and a total disc width of (46 ± 12) au (Moór et al. 2017), all three modelling approaches for scattered light prefer a radial profile that deviates from the ALMA data. The HG approach (Section 4.3.1) used a narrow ring with a peak in surface brightness at 110 ± 25 au. Both DDA and Mie used an asymmetric Gaussian starting at the peak fixed to 88 au (similar to ALMA) and a disc width of 30 au (Fig. 13). Despite the differences in profile set-ups we found that HG, DDA, and Mie predict that within 88 au the amount of dust material has to be small in order to fit the scattered light data (see Sections 5.2.2 and 5.2.3 for details).

Due to the low spatial resolution of the ALMA data and the nearly edge-on orientation of the disc, the region within 88 au is not well constrained at long wavelengths. The uncertainty in disc width (46 ± 12 au) inferred from ALMA is one indicator for this. Ignoring the inner region, we find that the radial profile for the outer region ($r > 88$ au) seems to be consistent in ALMA and scattered light models [width for ALMA: (23 ± 6) au; width for SPHERE: (30 ± 5) au].

### 5.2.2 Inward transport – PR-drag

We found that all scattered light models predict that the amount of dust within 88 au is low. In this section we investigate how large the amount is that we would expect at these regions due to transport processes, and whether the dust should be visible in scattered light. Similar to Pawellek et al. (2019a) we analyse the effect of Poynting–Robertson drag on the radial distribution to estimate the amount of dust drifting inwards from the planetesimal belt. To do so we use the collisional code ACE (Löhne et al. 2017) for a planetesimal belt between 88 and 118 au and bodies up to 40 km in radius. We collisionally evolve the belt for several Myr and infer the surface density of the dust as function of radius. From this we then calculate the flux density applying equation (12).

In Fig. 19 we show the flux density as function of radius normalized to the maximum found at the location of the planetesimal belt. At the age of 16 Myr – the proposed age of HD 131488 – the amount of dust drifting inwards due to PR-drag leads to a flux density between 6 and 7 orders of magnitudes lower than that of the planetesimal
<figure><img src="image_21.png" /><figcaption>Figure 20. Disc models for a disc radius of 100 au and a disc width of 20 au applying different porosities P.</figcaption></figure>

belt which is well below the detection limit of VLT/SPHERE. Even for older systems the flux density within 88 au would be 5 orders of magnitude lower. The low amount of dust expected for distances close to the planetesimal belt is in agreement with the scattered light models suggesting that particles moving inwards due to PR-drag do not significantly contribute to the scattered light.

### 5.2.3 Inward transport – gas drag

In addition to PR-drag the gas present in the disc might cause the dust grains to migrate. For HD 131488 we found gas between 30 and 130 au (Section 3.5.3) so that it seems possible to have dust as close as 30 au. We use the classical approach of Takeuchi & Artymowicz (2001) to estimate the amount of dust within 88 au.

The ratio of the force supporting the gas against stellar gravity to the gravity force is given by the parameter $\eta$. The gas can add (remove) angular momentum to (from) the dust grains. As a result the particles migrate outwards (inwards) until they reach a certain stability distance, $r_s$. At this distance the gas pressure gradient and the stellar radiation pressure balance each other (Takeuchi & Artymowicz 2001) so that

$$ \beta(s) = \eta(s, r_s). \qquad (16) $$

We are interested in the particles that migrate inwards. For those $\beta < \eta$. In a similar approach to Krivov et al. (2009) we inferred the $\beta$-values for which this relation is fulfilled (see fig. 2 in Krivov et al. 2009). Note that $\eta$ is only a function of the gas surface density profile, and the gas temperature, $T$, (equation 10 in Krivov et al. 2009). For both we assume power laws with typical exponents: $\Sigma \propto r^{-3/2}$ and $T \propto r^{-1/2}$. We assume that $\eta \propto L_{\text{star}}^{-0.25}/M_{\text{star}}$. At a distance of 88 au it is expected that only grains with $\beta < 0.05$ are dragged inwards. All particles with $\beta > 0.05$ are expected to drift outwards. This means that there seems to be no inward-drift of (sub-)blowout particles for which $\beta \gtrsim 0.5$ (Fig. 5).

Additionally, we need to take into account the drag force of the gas which determines the migration time-scale of the large dust grains and thus, how many of them we would expect to drift inwards. Following Marino et al. (2020) we found that for particles with $\beta \le 0.05$ the migration time-scale would be longer than the collision time-scale. Thus, for these grains we do not expect inward-migration due to gas-drag.

We conclude that the effect of PR-drag is more dominant than gas drag when analysing the inner region of HD 131488 in scattered light. As shown in Section 5.2.2, the flux density coming from grains migrating inwards due to PR-drag is several orders of magnitude lower compared to that of the planetesimal belt. Based on our findings

we assume the effect of gas drag to be lower than that of PR-drag, and thus, we assume that the total amount of dust migrating inwards from the planetesimal belt is small and remains unseen in scattered light.

### 5.2.4 Projection effects

In Fig. 20 we show a disc model where $r_0 = 100$ au, $\Delta r = 20$ au, $i = 84^{\circ}$, and $M_d = 1 M_{\oplus}$ for different levels of porosity. We find that for an edge-on disc the radial extent of surface brightness seems to decrease with increasing porosity.

An explanation might be given by the scattering phase function. The lower the disc inclination the lower the range of scattering angles we can observe. For a face-on disc we only get particles with a scattering angle of $\theta = 90^{\circ}$. Thus, if the material becomes more porous, the total flux density of the disc decreases (Fig. 6; Appendix B; Samra, Helling & Birnstiel 2022). As a result, the sensitivity limit of our instrument is reached at smaller distances from the star already so that the radial extent might seem smaller.

If the disc is now close to edge-on, we nearly cover all scattering angles between 0 and $180^{\circ}$. For more porous material the phase function shows a steeper decrease at larger scattering angles compared to compact materials. Thus, with increasing porosity the particles at larger scattering angles contribute less to the total flux density which results in an apparently decreasing disc extent shown in Fig. 20.

## 5.3 Dust properties

### 5.3.1 Scattering phase function

Our first scattered light image was generated by using the HG approximation where we found $g = 0.67$ to give the best-fitting model (Section 4.3.1). This indicates a relatively high level of forward scattering compared to other debris discs that were modelled with HG in scattered light (e.g. Engler et al. 2020; Stark et al. 2023). A strong forward scattering is expected for spherical particles which is the reason that our Mie models were well suited to fit the scattered light data.

Assuming that the results from the HG model give the best approximation of the ‘real’ scattering phase function of HD 131488, we find that a porosity of $P = 0.6$ is the closest fit to the HG function for both DDA and Mie models (Fig. 17, Section A2), but that we cannot rule out porosities of 0.2 and 0.4. A porosity of 0.0 would lead to a higher fraction of backward scattering, a porosity of 0.8 would not exhibit enough backward scattering compared to the HG function.
### 5.3.2 SED

We modelled the SED of HD 131488 for the different porosities and found that for $P = 0.8$ the size distribution parameter is not well constrained as $q$ reaches a value of $-8.2 \pm 2.7$ (note the negative sign). The size distribution index is determined by the long-wavelength data (far-IR to mm). For higher porosities the decrease becomes steeper even if $q$ stays constant. A negative $q$-value shows the dominance of the largest particles in the size distribution ($10^4 \mu m$) in order to fit the long-wavelength data. This is not consistent with collisional evolution but rather a pure outcome of the fitting procedure.

Since lower porosities result in reasonable fits of the photometric data, and are in agreement with collisional evolution, we conclude that a very high porosity of $P = 0.8$ seems unrealistic for the material in the disc of HD 131488 based on the SED. This is confirmed by the scattered light models.

### 5.3.3 Conclusion

Combining the results from the scattering phase function, and SED modelling, we conclude that the material in the disc around HD 131488 probably possesses a porosity between 0.2 and 0.6. While we cannot rule out smaller porosities, we find that porosities larger than 0.6 seem unlikely. Porosities of 0.6 are consistent with results from our Solar system where porosities of ~50 per cent were found for ‘rubble-pile’ asteroids (e.g. Weidling et al. 2009; Walsh 2018; Omura & Nakamura 2021; Sakatani et al. 2021).

In comparison a study of the debris disc around AU Mic found a significantly higher porosity of 76 per cent (Arnold et al. 2022). The difference in porosity between the late and early-type stars might indicate a direct influence of dynamical excitation on the dust material. Pawellek & Krivov (2015) found that debris discs around earlier-type stars such as HD 131488 possess a higher excitation i.e. a higher collision velocity than discs around late-type stars such as AU Mic. Collisions between planetesimals lead to compaction of the material and thus might decrease the overall porosity down to 40 per cent (e.g. Housen, Sweet & Holsapple 2018; Walsh 2018). This might be the case for HD 131488. In contrast to that, AU Mic might be less dynamically excited so that the material might not be compacted to the same degree as HD 131488. So far, we are lacking a study investigating the porosity of debris discs as function of stellar luminosity to make any conclusive remarks on relations between disc excitation and porosity.

## 5.4 Size distribution index

Our scattered light and SED model predict a size distribution index of $q = 3.0$ which suggests that all particle size bins contribute the same cross-section so that even large particles contribute to the near-IR scattered light image (Fig. 8). Assuming a collisional cascade we would expect $q$ to lie between 3 and 4 (Section 3.5) with an ideal collisional cascade at $q = 3.5$ (Dohnanyi 1969). This puts the size distribution index of our debris disc to the lower boundary for such a cascade. While not very common, such low $q$-values can also be found for other debris discs such as HD 32297 and HD 131835 (Löhne 2020; Norfolk et al. 2021). Both of these discs contain a significant amount of gas (Moór et al. 2019) so that we cannot rule out a link between gas content and size distribution index. On the other hand, other CO-rich debris discs did not show such a low $q$-value (e.g. HD 9672, HD 21997; Pawellek et al. 2014) so that a more thorough study is needed to draw any conclusions on a possible

<figure><img src="image_22.png" /><figcaption>**Figure 21.** Flux density as function of grain size. Blue lines indicate different proper eccentricities $\langle e \rangle$, red lines different levels of gas surface density $\Sigma$.</figcaption></figure>

link between those parameters. Studying a sample of 22 debris discs including both gaseous and gas-poor discs, Norfolk et al. (2021) found that for gaseous discs the $q$-value tends to be lower compared to that of gas-poor discs. However, this is based on a small number of gaseous discs so that an actual trend is still debatable.

Considering the results on the reflectance (Section 4.6.2), we were not able to constrain the size distribution index due to the large uncertainties.

SED models of debris discs that included mm-data suggest that a single power law as size distribution might not be realistic (e.g. this work; Lestrade et al. 2020). Furthermore, from our own Solar system we know that the size distributions of the Asteroid and Edgeworth-Kuiper belt change with size and also take values of $q < 3$ (e.g. Yoshida & Nakamura 2007; Morbidelli et al. 2021). Thus, it is likely that HD 131488 also possesses a more complex size distribution.

## 5.5 Combining scattered light and SED models

There are several studies that analysed debris discs at several wavelengths (e.g. Ertel et al. 2011; MacGregor et al. 2015; Ballering et al. 2016; Thebault & Kral 2019; Pawellek et al. 2019b; Esposito et al. 2020; Thebault, Olofsson & Kral 2023). The study of Schneider et al. (2006) investigated the debris disc around HD 181327 and tried to combine results from thermal emission and scattered light. In their fig. 14 the study shows that there is nearly no overlap between SED and scattered light models. With our study we were able to generate a self-consistent model fitting both thermal emission and scattered light for the first time with a semidynamical disc model.

While we were able to fit the scattered light data with a pure radiation pressure model, we found that we would need an amount of dust five times higher than required by thermal emission data when assuming a proper eccentricity of the planetesimals of $\langle e \rangle = 0.1$ (see Section 3.5.1).

We investigated the influence of dynamical excitation on dust mass. In Fig. 21 we show the total flux density per size bin for different levels of proper eccentricity $\langle e \rangle$ in blue. All models assume a porosity of 0.6. We see that a higher dynamical excitation increases the flux density per size bin as a higher amount of sub-blowout particles is produced by collisions. We also find that for a pure radiation pressure model with dynamical excitation of $\langle e \rangle = 0.9$ we only need a total dust mass of $1.3 M_\oplus$ to reproduce the scattered light observations. This is close to the mass inferred by the SED
(1.0 ± 0.2 M<sub>⊕</sub>), but first, the χ<sup>2</sup>-value for this model is larger than the critical value and thus, does not lead to a good fit. And secondly, a proper eccentricity of 0.9 leads to a very short collisional time-scale for the planetesimals. For example for a km-sized body at a radius of 88 au the lifetime is ~30 Myr assuming ⟨e⟩ = 0.1. Assuming ⟨e⟩ = 0.9 instead this shortens to 0.5 Myr for the same body (Wyatt et al. 2007b; Löhne, Krivov & Rodmann 2008). Thus, for 16 Myr-old HD 131488, it would be likely that the debris disc was collisionally depleted if ⟨e⟩ = 0.9. In that case, the disc would not be detectable for our instruments anymore. As a result, we need a model with a lower dynamical excitation but with a mechanism retaining the sub-blowout grains (in our case gas drag) to explain the observational data of HD 131488.

In Fig. 21 we see that a higher gas surface density increases the size of the grains coupled to the gas i.e. the number of size bins, but not the flux density per size bin. Thus, gas drag and dynamical excitation are degenerate. Compared to our best-fitting model with ⟨e⟩ = 0.1 and Σ = 2 × 10<sup>-5</sup> M<sub>⊕</sub> au<sup>-2</sup> we find that for a very high dynamical excitation of ⟨e⟩ = 0.5 we still need a gas surface density of Σ = 1 × 10<sup>-5</sup> M<sub>⊕</sub> au<sup>-2</sup> to reproduce the observations. The model also gives a χ<sup>2</sup>-value below the critical value and thus, a good fit.

While this model can reproduce the data, there are several arguments against it. First, the dynamical excitation is related to the size distribution index *q* which, in an ideal case, is given as 3.5. The larger *q*, the higher the amount of micron-sized particles compared to mm-sized ones (Fig. 8). This implies that for larger *q* the disc might possess a higher dynamical excitation as more small grains are produced during collisions of the large counterparts. However, the SED model of HD 131488 and its scattered light data both led to a best-fitting of *q* = 3.0 ± 0.1. This is a low value indicating a small dynamical excitation.

Secondly, the high amount of gas in the disc leads to damping of the particles' eccentricities which we use in our gas drag model. And thirdly, other dynamically excited ('self-stirred') discs were found to possess proper eccentricities of ⟨e⟩ ≲ 0.2 rather than 0.5 (e.g. Krivov, Löhne & Sremčević 2006; Thébault & Augereau 2007; Löhne et al. 2012; Pawellek & Krivov 2015; Schüppler et al. 2015; Daley et al. 2019; Geiler et al. 2019; Matrà et al. 2019b). We therefore prefer the model with ⟨e⟩ = 0.1 and Σ = 2 × 10<sup>-5</sup> M<sub>⊕</sub> au<sup>-2</sup>. Anyhow, we cannot rule out smaller proper eccentricities.

## 5.6 Influence of gas

As mentioned in Section 5.5, a pure radiation pressure model is not able to reproduce the observational data of HD 131488. Only a retaining mechanism such as gas drag led to consistent scattered light and thermal emission models.

We found a gas surface density for CO of 2 × 10<sup>-6</sup> M<sub>⊕</sub> au<sup>-2</sup> but assume that other gas species might add to this. Smirnov-Pinchukov et al. (2022) tried to find additional molecules in CO-rich debris discs but did not detect any. It is likely that the species targeted by this study are not shielded from the stellar UV radiation and thus, dissociate very quickly (e.g. Matrà et al. 2018). Additionally, the selected molecules are not thought to be dominant components in the gas mixture so that they would hardly contribute to the total gas mass. We would need detections of C, O, and H to constrain the total gas mass reliably assuming a secondary origin for the gas. In the case of a primordial origin, a detection of H<sub>2</sub> would help constraining the gas mass.

A total surface density of 10<sup>-5</sup> M<sub>⊕</sub> au<sup>-2</sup> for our debris disc seems possible (Section 3.5.3), but we cannot offer a strong constraint on

<table><caption>**Table 2.** Summary of best-fitting results for HD 131488 and the sections of the paper where they were discussed.</caption><thead><tr><th>Parameter</th><th>Best fit</th><th>Reference section</th></tr></thead><tbody><tr><td>R<sub>central</sub> [au]</td><td>88 ± 5</td><td>Section 4.3.2</td></tr><tr><td>ΔR [au]</td><td>30 ± 3</td><td>Section 4.3.2</td></tr><tr><td>s<sub>min</sub> [μm]</td><td>10<sup>-1</sup></td><td>Section 3.5.1</td></tr><tr><td>s<sub>max</sub> [μm]</td><td>10<sup>4</sup></td><td>Section 3.5.2</td></tr><tr><td>q</td><td>3.0 ± 0.1</td><td>Sections 4.3.2, 4.3.3, 4.4</td></tr><tr><td>M<sub>dust</sub> [M<sub>⊕</sub>]</td><td>1.0 ± 0.2</td><td>Section 4.4</td></tr><tr><td>P</td><td>0.2...0.6</td><td>Section 5.3</td></tr><tr><td>Σ<sub>gas</sub> [M<sub>⊕</sub> au<sup>-2</sup>]</td><td>(2.0 ± 0.1) × 10<sup>-5</sup></td><td>Section 4.5</td></tr></tbody></table>

this number. However, with our simple assumptions on gas drag we found that a gas surface density of 2 × 10<sup>-5</sup> M<sub>⊕</sub> au<sup>-2</sup> is sufficient to couple a high amount of sub-blowout grains so that we can fit the SED and scattered light data with the same amount of dust. While this is an first order agreement with the rough estimate from ALMA, a more detailed work is needed to get reliable constraints on the total gas mass.

What we can say is that, when ignoring a mechanism that can retain small dust particles (in our case gas) it seems not possible to fit both thermal emission and scattered light of the disc around HD 131488 at the same time. Matching the amount of dust needed to fit the SED and the scattered light data opened a way to roughly estimate the surface density of the gas.

## 6 SUMMARY

In this study we analysed the scattered light and thermal emission data of the debris disc around HD 131488 applying a semidynamical disc model in combination with HG, DDA, and Mie theory. The SPHERE/IRDIS and IFS data of HD 131488 were presented for the first time. The modelling results are summarized in Table 2.

The high amount of CO-gas found in the disc (Moór et al. 2017) is capable of retaining a large fraction of sub-blowout grains. Only if we take into account these particles, we are able to generate a model that can fit all data available for this debris disc (thermal emission and scattered light). This opens a way of roughly estimating the amount of gas necessary to fit all data.

The radial profile of the planetesimal belt preferred by the scattered light models seems to deviate from the one inferred by ALMA observations in the way that in scattered light we do not expect a significant amount of dust within the central radius of 88 au. This is in agreement with expectations from PR and gas drag models. However, the deviation of the profiles might be attributed to the low spatial resolution of the ALMA data.

The disc possesses a flat size distribution (q = 3) and moderate level of porosity of ~20...60 per cent which is in agreement with a collisional cascade and results from Asteroid observations suggesting a pebble pile scenario for planetesimal growth within this system. Compared to the disc around AU Mic the material might have been more compacted by collisions. While the reflectance might indicate a slightly blue colour of the debris disc, the uncertainties of the observations are too large to draw any conclusions. However, the modelling results including gas drag indicate a blue colour as well and thus, might be in agreement with the observational results.

The modelling approach of DDA and Mie led to similar results when assuming particles of basic spherical shape and small sizes for vacuum inclusions. For HD 131488 Mie theory leads to well-fitting models indicating that the dust particles possess a scattering behaviour similar to spheres. To study the influence of more complex
particle shapes more work is needed as DDA is limited to small sizes
making Mie grains necessary to fill-up the size distribution.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their constructive criticism and
help improving the presentation of this study.

NP is grateful to Torsten Löhne, Philippe Thébault, Jürgen Blum,
and Peter Woitke for many useful discussions. NP also thanks
Kevin Wagner and Benoit Pairet for comments on data analysis.
AMH is supported by a Cottrell Scholar Award from the Research
Corporation for Science Advancement. SM is funded by the Royal
Society through a University Research Fellowship.

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.

SPHERE is an instrument designed and built by a consortium
consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Ger-
many), LAM (Marseille, France), LESIA (Paris, France), Laboratoire
Lagrange (Nice, France), INAF-Osservatorio di Padova (Italy),
Observatoire de Genève (Switzerland), ETH Zurich (Switzerland),
NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands)
in collaboration with ESO. SPHERE was funded by ESO, with
additional contributions from CNRS (France), MPIA (Germany),
INAF (Italy), FINES (Switzerland), and NOVA (Netherlands).
SPHERE also received funding from the European Commission
Sixth and Seventh Framework Programmes as part of the Optical
Infrared Coordination Network for Astronomy (OPTICON) under
grant number RI13-Ct-2004-001566 for FP6 (2004–2008), grant
number 226604 for FP7 (2009–2012) and grant number 312430
for FP7 (2013–2016). We also acknowledge financial support from
the Programme National de Planétologie (PNP) and the Programme
National de Physique Stellaire (PNPS) of CNRS-INSU in France.
This work has also been supported by a grant from the French
Labex OSUG@2020 (Investissements d'avenir – ANR10 LABX56).
The project is supported by CNRS, by the Agence Nationale de la
Recherche (ANR-14-CE33-0018). It has also been carried out within
the frame of the National Centre for Competence in Research PlanetS
supported by the Swiss National Science Foundation (SNSF). MRM,
HMS, and SD are pleased to acknowledge this financial support of
the SNSF. Finally, this work has made use of the the SPHERE
Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTH-
EAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice) and Obser-
vatoire de Paris/LESIA (Paris) and is supported by a grant from Labex
OSUG@2020 (Investissements d'avenir-ANR10 LABX56). This
work has made use of the High Contrast Data Centre, jointly operated
by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille),
OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and
Observatoire de Lyon/CRAL, and supported by a grant from Labex
OSUG@2020 (Investissements d'avenir – ANR10 LABX56). We
thank P Delorme and E Lagadec (High Contrast Data Centre) for
their efficient help during the data reduction process.

This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2015.1.01243.S. ALMA is a partnership of ESO
(representing its member states), NSF (USA), and NINS (Japan),
together with NRC (Canada), MOST and ASIAA (Taiwan), and
KASI (Republic of Korea), in cooperation with the Republic of Chile.
The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and
NAJO.

Based on observations collected at the European Southern Observatory under ESO programme 0101.C-0753(B).

DATA AVAILABILITY

The data underlying this article will be shared on request to the
corresponding author. The ALMA and VLT/SPHERE data are
publicly available and can be queried and downloaded directly
from the ALMA archive: https://almascience.nrao.edu/asax/ and the
SPHERE archive: https://archive.eso.org/wdb/wdb/eso/sphere/.

REFERENCES

Adam C. et al., 2021, A&A, 653, A88

Amara A., Quanz S. P., 2012, MNRAS, 427, 948

Arnold J. A., Weinberger A. J., Videen G., Zubko E. S., 2019, AJ, 157, 157

Arnold J. A., Weinberger A. J., Videen G., Zubko E. S., 2022, ApJ, 930, 123

Audu S., Weinberger A., Videen G., Arnold J., Zubko E., 2023, in American Astronomical Society Meeting Abstracts. p. 269.05

Augereau J.-C., Papaloizou J. C. B., 2004, A&A, 414, 1153

Bailer-Jones C. A. L., Rybikz J., Fouesneau M., Mantelet G., Andrae R., 2018, AJ, 156, 58

Ballering N., Rieke G., Su K. Y. L., Gaspar A., 2016, in American Astronomical Society Meeting Abstracts #227. p. 228.03

Besla G., Wu Y., 2007, ApJ, 655, 528

Beuzit J. L. et al., 2019, A&A, 631, A155

Bhowmik T. et al., 2019, A&A, 630, A85

Biller B. A. et al., 2015, MNRAS, 450, 4446

Birnstiel T. et al., 2010, A&A, 516, L14

Blum J., Wurm G., 2008, ARA&A, 46, 21

Bohren C. F., Huffman D. R., 1983, Absorption and Scattering of Light by Small Particles. Wiley and Sons, New York, Chichester, Brisbane, Toronto, Singapore

Bonsor A., Wyatt M. C., Marino S., Davidsson B. J. R., Kral Q., 2023, MNRAS, 526, 3115

Bruggeman D. A. G., 1935, Ann. Phys., Lpz., 416, 636

Bruggeman D. A. G., 1936, Ann. Phys., Lpz., 417, 645

Brunngräber R., Wolf S., Kirchschlager F., Ertel S., 2017, MNRAS, 464, 4383

Burns J. A., Lamy P. L., Soter S., 1979, Icarus, 40, 1

Cantalloube F. et al., 2015, A&A, 582, A89

Carbillet M. et al., 2011, Exp. Astron., 30, 39

Castelli F., Kurucz R. L., 2004, preprint (arXiv:astro-ph/0405087)

Cataldi G. et al., 2023, ApJ, 951, 111

Chen C. H., Sargent B. A., Bohac C., et al., 2006, ApJS, 166, 351

Choquet É. et al., 2017, ApJ, 834, L12

Claudi R. U. et al., 2008, in McLean I. S., Casali M. M., eds, Proc. SPIE Conf. Ser. Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II. SPIE, Bellingham, p. 70143E

Cutri R. M. et al., 2003, 2MASS All Sky Catalog of Point Sources

DENIS Consortium, 2005, VizieR Online Data Catalog, p. II/263

Daley C. et al., 2019, ApJ, 875, 87

Delorme P. et al., 2017, in Reylé C., Di Matteo P., Herpin F., Lagadec E., Lançon A., Meliani Z., Royer F., eds, SF2A-2017: Proc. Annual meeting of the French Society of Astronomy and Astrophysics. p. 347

Dent W. R. F. et al., 2014, Science, 343, 1490

Di Folco E., Péricaud J., Dutrey A., Augereau J. C., Chapillon E., Guilloteau S., Piétu V., Boccaletti A., 2020, A&A, 635, A94

Dohlen K., Saisse M., Origne A., Moreaux G., Fabron C., Zamkotsian F., Lanzoni P., Lemarquis F., 2008, in Advanced Optical and Mechanical Technologies in Telescopes and Instrumentation. p. 701859

Dohnanyi J. S., 1969, J. Geophys. Res., 74, 2531

Draine B. T., 1988, ApJ, 333, 848

Draine B. T., 2003, ARA&A, 41, 241

Draine B. T., Flatau P. J., 2010, preprint (arXiv:1002.1505)

Elliot J. L. et al., 2005, AJ, 129, 1117
Engler N. et al., 2017, A&A, 607, A90

Engler N. et al., 2019, A&A, 622, A192

Engler N. et al., 2020, A&A, 635, A19

Ertel S., Wolf S., Metchev S., Schneider G., Carpenter J. M., Meyer M. R., Hillenbrand L. A., Silverstone M. D., 2011, A&A, 533, A132

Esposito T. M. et al., 2016, AJ, 152, 85

Esposito T. M. et al., 2020, AJ, 160, 24

Favre C., Cleeves L. I., Bergin E. A., Qi C., Blake G. A., 2013, ApJ, 776, L38

Feldt M. et al., 2017, A&A, 601, A7

Flaherty K. M. et al., 2016, ApJ, 818, 97

Fulle M. et al., 2015, ApJ, 802, L12

Gaia Collaboration, 2018, VizieR Online Data Catalog, p. I/345

Gaia Collaboration, 2016, A&A, 595, A1

Geiler F., Krivov A. V., Booth M., Löhne T., 2019, MNRAS, 483, 332

Henden A. A., Templeton M., Terrell D., Smith T. C., Levine S., Welch D., 2016, VizieR Online Data Catalog, p. II/336

Henning T., Mutschke H., 1997, A&A, 327, 743

Henyey L. G., Greenstein J. L., 1941, ApJ, 93, 70

Hög E. et al., 2000, A&A, 355, L27

Housen K. R., Sweet W. J., Holsapple K. A., 2018, Icarus, 300, 72

Ishihara D. et al., 2010, A&A, 514, A1

Jäger C., Mutschke H., Henning T., Huisken F., 2008, ApJ, 689, 249

Kirchschlager F., Wolf S., 2013, A&A, 552, A54

Kirchschlager F., Wolf S., 2014, A&A, 568, A103

Klahr H., Lin D. N. C., 2005, ApJ, 632, 1113

Kóspál Á. et al., 2013, ApJ, 776, 77

Kral Q., Matrà L., Wyatt M. C., Kennedy G. M., 2017, MNRAS, 469, 521

Kral Q., Marino S., Wyatt M. C., Kama M., Matrà L., 2019, MNRAS, 489, 3670

Krist J. E. et al., 2010, AJ, 140, 1051

Krivov A. V., 2010, Res. Astron. Astrophys., 10, 383

Krivov A. V., Löhne T., Sremčević M., 2006, A&A, 455, 509

Krivov A. V., Herrmann F., Brandeker A., Thébault P., 2009, A&A, 507, 1503

Lawson K. et al., 2021, AJ, 162, 293

Lee E. J., Chiang E., 2016, ApJ, 827, 125

Lestrade J. F. et al., 2020, in MayetF, CatalanoA., Macías-Pérez J. F., PerottoL., eds, mm Universe @ NIKAZ – Observing the mm Universe with the NIKAZ Camera, Vol. 228. EPJ Web of Conf., p. 00015

Li A., Greenberg J. M., 1998, A&A, 331, 291

Lieman-Sifry J., Hughes A. M., Carpenter J. M., Gorti U., Hales A., Flaherty K. M., 2016, ApJ, 828, 25

Lindegren L. et al., 2018, A&A, 616, A2

Lisse C. M., Sitko M. L., Russell R. W., Marengo M., Currie T., Melis C., Mittal T., Song I., 2017, ApJ, 840, L20

Lisse C. M. et al., 2020, ApJ, 894, 116

Löhne T., 2020, A&A, 641, A75

Löhne T., Krivov A. V., Rodmann J., 2008, ApJ, 673, 1123

Löhne T. et al., 2012, AN, 333, 441

Löhne T., Krivov A. V., Kirchschlager F., Sende J. A., Wolf S., 2017, A&A, 605, A7

MacGregor M. A., Wilner D. J., Andrews S. M., Hughes A. M., 2015, ApJ, 801, 59

McGuire A. F., Hapke B. W., 1995, Icarus, 113, 134

Maire A.-L. et al., 2016, in EvansC. J., SimardL., TakamiH., eds, Proc. SPIE, Conf. Ser. Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI. SPIE, Bellingham, p. 990834

Marino S. et al., 2016, MNRAS, 460, 2933

Marino S., Flock M., Henning T., Kral Q., Matrà L., Wyatt M. C., 2020, MNRAS, 492, 4409

Marois C., Lafrenière D., Doyon R., Macintosh B., Nadeau D., 2006, ApJ, 641, 556

Martinez P., Dorrer C., Aller Carpentier E., Kasper M., Boccaletti A., Dohlen K., Yaitskova N., 2009, A&A, 495, 363

Marton G. et al., 2017, preprint (arXiv:1705.05693)

Matrà L., Marino S., Kennedy G. M., Wyatt M. C., Öberg K. I., Wilner D. J., 2018, ApJ, 859, 72

Matrà L., Öberg K. I., Wilner D. J., Olofsson J., Bayo A., 2019a, AJ, 157, 117

Matrà L., Wyatt M. C., Wilner D. J., Dent W. R. F., Marino S., Kennedy G. M., Milli J., 2019b, AJ, 157, 135

Melis C., Zuckerman B., Rhee J. H., Song I., Murphy S. J., Bessell M. S., 2013, ApJ, 778, 12

Meyer M. R., Backman D. E., Weinberger A. J., Wyatt M. C., 2007, Protostars and Planets V, p. 573

Michel B., Henning T., Stognienko R., Rouleau F., 1996, ApJ, 468, 834

Mie G., 1908, Ann. Phys., Lpz., 330, 377

Millar-Blanchaer M. A. et al., 2015, ApJ, 811, 18

Milli J. et al., 2017, A&A, 599, A108

Min M., Hovenier J. W., de Koter A., 2005, A&A, 432, 909

Min M., Rab C., Woitke P., Dominik C., Ménard F., 2016, A&A, 585, A13

Miotello A., Kamp I., Birnstiel T., Cleeves L. C., Kataoka A., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, ASP Conf. Ser. Vol. 534, Protostars and Planets VII. Astron. Soc. Pac., San Francisco. p. 501

Moór A. et al., 2017, ApJ, 849, 123

Moór A. et al., 2019, ApJ, 884, 108

Moór A. et al., 2020, AJ, 159, 288

Morbidelli A., Nesvorny D., Bottke W. F., Marchi S., 2021, Icarus, 356, 114256

Mugnai A., Wiscombe W. J., 1986, Appl. Opt., 25, 1235

Mumma M. J., Charnley S. B., 2011, ARA&A, 49, 471

Mutschke H., Mohr P., 2019, A&A, 625, A61

Norfolk B. J. et al., 2021, MNRAS, 507, 3139

Olofsson J. et al., 2016, A&A, 591, A108

Olofsson J. et al., 2019, A&A, 630, A142

Olofsson J., Milli J., Bayo A., Henning T., Engler N., 2020, A&A, 640, A12

Olofsson J., Thébault P., Kennedy G. M., Bayo A., 2022, A&A, 664, A122

Omura T., Nakamura A. M., 2021, PSJ, 2, 41

Pavlov A., Möller-Nilsson O., Feldt M., Henning T., Beuzit J.-L., Mouillet D., 2008, in Bridger A., Radziwill N. M., eds, Proc. SPIE Conf. Ser. Vol. 7019, Advanced Software and Control for Astronomy II. SPIE, Bellingham. p. 701939

Pawellek N., Krivov A. V., 2015, MNRAS, 454, 3207

Pawellek N., Krivov A. V., Marshall J. P., Montesinos B., Ábrahám P., Moór A., Bryden G., Eiroa C., 2014, ApJ, 792, 65

Pawellek N., Moór A., Pascucci I., Krivov A. V., 2019a, MNRAS, 487, 5874

Pawellek N. et al., 2019b, MNRAS, 488, 3507

Pawellek N., Wyatt M., Matrà L., Kennedy G., Yelverton B., 2021, MNRAS, 502, 5390

Pecaut M. J., Mamajek E. E., 2016, MNRAS, 461, 794

Péricaud J., Di Folco E., Dutrey A., Guilloteau S., Piétu V., 2017, A&A, 600, A62

Perrot C. et al., 2016, A&A, 590, L7

Pinte C. et al., 2008, A&A, 489, 633

Purcell E. M., Pennypacker C. R., 1973, ApJ, 186, 705

Rebollido I. et al., 2018, A&A, 614, A3

Rebollido I. et al., 2022, MNRAS, 509, 693

Ren B. B. et al., 2023, A&A, 672, A114

Ricci L., Testi L., Maddison S. T., Wilner D. J., 2012, A&A, 539, L6

Richert A. J. W., Lyra W., Kuchner M. J., 2018, ApJ, 856, 41

Sakatani N. et al., 2021, Nat. Astron., 5, 766

Samra D., Helling C., Birnstiel T., 2022, A&A, 663, A47

Schneider G. et al., 2006, ApJ, 650, 414

Schneiderman T. et al., 2021, Nature, 598, 425

Schuerman D. W., Wang R. T., Gustafson B. Å. S., Schaefer R. W., 1981, Appl. Opt., 20, 4039

Schüppler C. et al., 2015, A&A, 581, A97

Singh G. et al., 2021, A&A, 653, A79

Smirnov-Pinchukov G. V., Moór A., Semenov D. A., Ábrahám P., Henning T., Kóspál Á., Hughes A. M., di Folco E., 2022, MNRAS, 510, 1148

Soummer R., Pueyo L., Larkin J., 2012, ApJ, 755, L28

Stark C. C., Ren B., MacGregor M. A., Howard W. S., Hurt S. A., Weinberger A. J., Schneider G., Choquet E., 2023, ApJ, 945, 131

Strubbe L. E., Chiang E. I., 2006, ApJ, 648, 652
Szulágyi J., Dullemond C. P., Pohl A., Quanz S. P., 2019, MNRAS, 487, 1248

Takeuchi T., Artymowicz P., 2001, ApJ, 557, 990

Thébault P., Augereau J.-C., 2007, A&A, 472, 169

Thebault P., Kral Q., 2019, A&A, 626, A24

Thébault P., Wu Y., 2008, A&A, 481, 713

Thebault P., Olofsson J., Kral Q., 2023, A&A, 674, A51

van der Plas G. et al., 2019, A&A, 624, A33

Vitense C., Krivov A. V., Löhne T., 2010, A&A, 520, A32

Vorobyov E. et al., 2022, A&A, 658, A191

Walsh K. J., 2018, ARA&A, 56, 593

Weidling R., Güttler C., Blum J., Brauer F., 2009, ApJ, 696, 2036

Weiss-Wrana K., 1983, A&A, 126, 240

Wright E. L. et al., 2010, AJ, 140, 1868

Wyatt M. C., 2008, ARA&A, 46, 339

Wyatt M. C., Dermott S. F., Telesco C. M., Fisher R. S., Grogan K., Holmes E. K., Piña R. K., 1999, ApJ, 527, 918

Wyatt M. C., Smith R., Greaves J. S., Beichman C. A., Bryden G., Lisse C. M., 2007a, ApJ, 658, 569

Wyatt M. C., Smith R., Su K. Y. L., Rieke G. H., Greaves J. S., Beichman C. A., Bryden G., 2007b, ApJ, 663, 365

Yoshida F., Nakamura T., 2007, Planet. Space Sci., 55, 1113

Ysard N., Jones A. P., Demyk K., Boutéraon T., Koehler M., 2018, A&A, 617, A124

Zubko E., 2013, Earth Planets Space, 65, 139

Zubko V. G., Mennella V., Colangeli L., Bussoletti E., 1996, MNRAS, 282, 1321

## APPENDIX A: SCATTERING PHASE FUNCTIONS

### A1 Size and spatial distribution of voids

In Fig. A1 we show the scattering phase function as function of scattering angle analysing the influence of different void sizes (inclusions of vacuum) and spatial distributions of those voids. The void size of 1/100 was applied in all DDA scattered light models of this study (red solid line, top panel). We see that for small sizes (1/100 and 2/100) the phase function does not change significantly, but that the changes become more pronounced with larger sizes (5/100 and 20/100) which is in agreement with results from studies using more complex particle structures (e.g. Arnold et al. 2019).

In a similar fashion we kept the void size constant (20/100, bottom panel of Fig. A1) and analysed the influence of the spatial distribution of the vacuum inclusions. We see that for a large size even their spatial distribution can change the phase function significantly.

### A2 Different porosities

In Fig. A2 we show the scattering phase function for different porosities as inferred from our models: for $s \le 10 \mu m$ we use DDA, for $s > 10 \mu m$ we use Mie theory.

<figure><img src="image_23.png" /><figcaption>**Figure A1.** Scattering phase function as function of scattering angle for a particle of radius 1 µm, assuming a porosity of $P = 0.4$ at a wavelength of 1.6 µm. Top panel: Change of the void size but keeping the spatial distribution constant. Bottom panel: Change of the spatial distribution of voids with constant sizes of 20/100 voids per particle diameter. The different random distributions are called a, b, c, and d.</figcaption></figure>

We see that for sub-micron-sized dust the scattering is more or less isotropic and the behaviour independent of the level of porosity. This changes for micron-sized and larger grains. Up to 10 µm-sized grains the phase functions become more complex in structure.

For increasing porosity the values for back-scattering ($\theta \sim 180^\circ$) decreases for all grains. However, we note that for $s = 10 \mu m$ we cannot differentiate between the different phase functions. Interestingly, the forward-scattering ($\theta < 5^\circ$) does not change with porosity.
<figure><img src="image_24.png" /><figcaption>**Figure A2.** Scattering phase function as function of scattering angle for different porosities at a wavelength of 1.6 µm. The different panels show the results for different grain sizes.</figcaption></figure>

## APPENDIX B: FLUX DENSITIES

In Fig. B1 we show the contribution to the total flux density per size bin for different porosities. As expected from Fig. 5 we see

an increase of the blowout size with increasing porosity. For bound grains the level of flux density is comparable. The total flux density is decreasing with increasing porosity.
<figure><img src="image_25.png" /><figcaption>**Figure B1.** Flux density as function of grain size for different porosities.</figcaption></figure>

## APPENDIX C: $^{12}$CO POSITION-VELOCITY DIAGRAM

Here we reanalyse the $^{12}$CO J = 2–1 emission reported by Moór et al. (2017) to constrain the extent of the gas. Fig. C1 shows a positional-velocity diagram of $^{12}$CO obtained assuming the inclination and PA derived from the scattered light images, and a stellar mass of 1.8 $M_{\odot}$ (Matrà et al. 2018). We can constrain the radial distribution of CO by overlaying two diagonal lines representing the line-of-sight velocities as a function of separation along the major axis, at two fixed orbital radii, and assuming Keplerian rotation. The curves in white dotted lines connecting the two diagonal lines show the maximum line-of-sight velocities as a function of projected separation. By varying the two orbital radii such that the white wedges enclose most of the CO emission, we find that the CO gas is mostly contained between 30 and 130 au. The significant emission just outside the white wedges is due to the large beam size that smooths the radial extent of CO.

<figure><img src="image_26.png" /><figcaption>**Figure C1.** Position-velocity diagram of $^{12}$CO J = 2–1 emission. The grey contours represent emission at $3\sigma$. The diagonal white solid lines represent the line-of-sight velocity of gas at a fixed orbital radius and in Keplerian rotation as a function of projected separation. The white dotted lines show the maximum velocity along the line of sight for a Keplerian rotational profile and as a function of separation. The horizontal white line at the bottom left represents the beam FWHM of 0.51 arcsec in the direction of the disc PA.</figcaption></figure>

This paper has been typeset from a TeX/Latex file prepared by the author.