# ALMA millimetre-wavelength imaging of HD 138965: new constraints on the debris dust composition and presence of planetary companions

J P Marshall, S Hengst, A Trejo-Cruz, C del Burgo, J Milli, M Booth, J C

Augereau, E Choquet, F Y Morales, P Thébault, et al.

► To cite this version:

J P Marshall, S Hengst, A Trejo-Cruz, C del Burgo, J Milli, et al. ALMA millimetre-wavelength imaging of HD 138965: new constraints on the debris dust composition and presence of planetary companions. Monthly Notices of the Royal Astronomical Society, 2025, 541 (1), pp.71 - 84. 10.1093/mn-ras/staf984. hal-05469994

HAL Id: hal-05469994

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

Submitted on 21 Jan 2026

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 CC BY 4.0 - Attribution - International License
# ALMA millimetre-wavelength imaging of HD 138965: new constraints on the debris dust composition and presence of planetary companions

J. P. Marshall<sup>1,★</sup> S. Hengst,<sup>2</sup> A. Trejo-Cruz,<sup>3</sup> C. del Burgo,<sup>4,5</sup> J. Milli,<sup>6</sup> M. Booth,<sup>7</sup> J.C. Augereau,<sup>8</sup> E. Choquet,<sup>9</sup> F. Y. Morales,<sup>10,11</sup> P. Thébault,<sup>12</sup> F. Kemper<sup>10</sup>,<sup>13,14,15</sup> V. Faramaz-Gorka<sup>16</sup> and G. Bryden<sup>10</sup>

Accepted 2025 June 6. Received 2025 June 5; in original form 2025 January 5

## ABSTRACT

HD 138965 is a young A type star and member of the nearby young Argus association. This star is surrounded by a broad, bright debris disc with two temperature components that was spatially resolved at far-infrared wavelengths by Herschel. Here, we present Atacama Large Millimeter/Submillimeter Array (ALMA) millimetre-wavelength imaging of the cool outer belt. These reveal its radial extent to be $150_{-7}^{+10}$ au with a width ($\sigma$) of $49_{-6}^{+7}$ au ($\Delta R/R = 0.77$), at a moderate inclination of $49.9_{-3.7}^{+3.3}$. Due to the limited angular resolution, signal-to-noise and inclination we have no constraint on the disc's vertical scale height. We modelled the disc emission with both gravitational and radiation forces acting on the dust grains. As the inner belt has not been spatially resolved, we fixed its radius and width prior to modelling the outer belt. We find astronomical silicate is the best fit for the dust composition. However, we could not reject possible scenarios where there are at least 10 per cent water-ice inclusions. Combining the spatially resolved imaging by ALMA with non-detection at optical wavelengths by HST, we obtain a limit on the scattering albedo $\omega \le 0.09$ for the debris dust in the outer belt. Analysis of the outer belt's architecture in conjunction with simple stirring models places a mass limit of $2.3 \pm 0.4 M_{Jup}$ on a companion interior to the belt ($a \le 78$ au), a factor of two improvement over constraints from high contrast imaging.

**Key words:** planet–disc interactions – circumstellar matter – stars: individual: HD 138965 – radio continuum: planetary systems.

## 1 INTRODUCTION

Main-sequence stars play host to both planets and belts of planetesimals. Dynamical interaction between the components can excite the motions of the planetesimals to such a degree that a collisional cascade is triggered generating detectable amounts of dust. The dust produced in the collisional cascade is most commonly detected at optical or near-infrared wavelengths as scattered light (in total

intensity or polarization, e.g. Esposito et al. 2020; Ren et al. 2023; Marshall et al. 2023b; Crotts et al. 2024) or at infrared to millimetre wavelengths as thermal emission (e.g. Eiroa et al. 2013; Kennedy & Wyatt 2013; Montesinos et al. 2016; Holland et al. 2017; Matrà et al. 2025). Due to their origins in destructive processes, these dusty structures are collectively referred to as ‘debris discs’ (Wyatt 2008; Hughes, Duchêne & Matthews 2018).

These debris discs are larger and more massive analogues to the Solar system’s Asteroid or Edgeworth–Kuiper belts (Horner et al. 2020). There is some dependence of disc detection on the stellar luminosity, with the highest detection rate found around A

<footnote n="1"><sup>1</sup>Institute of Astronomy and Astrophysics, Academia Sinica, 11F of AS/NTU Astronomy-Mathematics Building, No.1, Section 4, Roosevelt Rd, Taipei 106319, Taiwan</footnote>

<footnote n="2"><sup>2</sup>School of Mathematics, Physics, and Computing, University of Southern Queensland, West Street, Toowoomba, QLD 4350, Australia</footnote>

<footnote n="3"><sup>3</sup>Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Antigua Carretera a Pátzcuaro # 8701, Ex-Hda. San José de la Huerta, Morelia, Michoacán, C.P. 58089, México</footnote>

<footnote n="4"><sup>4</sup>Instituto de Astrofísica de Canarias, Vía Láctea S/N, La Laguna 38200, Tenerife, Spain</footnote>

<footnote n="5"><sup>5</sup>Departamento de Astrofísica, Universidad de la Laguna, La Laguna 38200, Tenerife, Spain</footnote>

<footnote n="6"><sup>6</sup>Institut de Planétologie et d'Astrophysique de Grenoble, Université Grenoble Alpes, 621 Avenue Centrale, 38400 Saint-Martin-d'Hères, France</footnote>

<footnote n="7"><sup>7</sup>UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="8"><sup>8</sup>Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="9"><sup>9</sup>Aix-Marseille Université, CNRS, CNES, LAM, 13013 Marseille, France</footnote>

<footnote n="10"><sup>10</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="11"><sup>11</sup>STEM Division, Los Angeles Mission College, 13356 Eldridge Avenue, Sylmar, CA 91342, USA</footnote>

<footnote n="12"><sup>12</sup>LESIA-Observatoire de Paris, UPMC Univ. Paris 06, Univ. Paris-Diderot, 75014, France</footnote>

<footnote n="13"><sup>13</sup>Institut de Ciències de l'Espai (ICE, CSIC), Can Magrans, s/n, E-08193 Cerdanyola del Vallès, Barcelona, Spain</footnote>

<footnote n="14"><sup>14</sup>ICREA, Pg. Lluís Companys 23, E-08010 Barcelona, Spain</footnote>

<footnote n="15"><sup>15</sup>Institut d'Estudis Espacials de Catalunya (IEEC), E-08860 Castelldefels, Barcelona, Spain</footnote>

<footnote n="16"><sup>16</sup>Steward Observatory, Department of Astronomy, University of Arizona, 933 N. Cherry Ave, Tucson, AZ 85721, USA</footnote>

<footnote>*E-mail: jmarshall@asiaa.sinica.edu.tw</footnote>
stars (Thureau et al. 2014), an intermediate fraction around FGK stars (Sibthorpe et al. 2018), and the lowest detection rate around M dwarfs (Lestrade et al. 2009; Kennedy et al. 2018). However, these detections are strongly sensitivity limited at levels of fractional luminosity around ten times greater than that predicted for dust in the Solar system (≈10⁻⁷, Vitense et al. 2012). As such, the supposed lack of cold debris discs around M dwarf stars could be a result of that limitation (Luppe et al. 2020; Cronin-Coltsmann et al. 2023). Likewise, there is a dependence on disc detection with age. All else being equal, younger stars have brighter debris discs (Wyatt et al. 2007; Kains, Wyatt & Greaves 2011). This is consistent with debris discs originating from the steady collisional grinding of parent bodies into dust over the lifetime of the system (Najita, Kenyon & Bromley 2022). That said, the interstellar medium can influence the morphology of debris discs throughout their life (Heras et al. 2025). As such, young stars with reliable ages (e.g. in moving groups) make favourable targets for the detection and characterization of both debris discs and planetary companions.

HD 138 965 is a young A star lying at a distance of 78.60 ± 0.17 pc (Gaia Collaboration 2016; Babusiaux et al. 2023; Gaia Collaboration 2023). It has been identified as a member of the 40 to 50 Ma Argus association (Zuckerman 2019), but an older age of 350±³⁷/₅⁶ Ma based on fitting *Gaia* observations has been also proposed (David & Hillenbrand 2015). The first evidence of cool infrared excess was obtained from the *InfraRed Astronomical Satellite* (IRAS, Patten & Willson 1991; Rhee et al. 2007), with subsequent detections by *Akari* Far Infrared Surveyor (Liu, Wang & Jiang 2014) and *Spitzer* Multi-band Imaging Photometer for SIRTF (MIPS, Morales et al. 2009; Chen et al. 2014a). The *Spitzer* InfraRed Spectrograph (IRS) spectrum is featureless, but its shape is indicative that a substantial warm inner component to the debris disc is present within this system (Chen et al. 2014a). The cool, outer disc was spatially resolved in *Herschel* Photodetector Array Camera and Spectrotraph (PACS) observations, with a measured radius of 187 ± 6 au (Morales et al. 2016).

The disc has yet to be imaged in scattered light, and remains one of the few relatively bright debris discs ($L_d/L_* > 10^{-4}$) around nearby stars without such observations. Combining multiwavelength, spatially resolved imaging of debris disc systems with radiative transfer modelling is necessary to weaken degeneracies inherent between the dust properties and architecture (e.g. Marshall et al. 2014), thereby providing the greatest precision possible from the determination of the system's properties.

The architecture of HD 138965's debris disc, with two distinct temperature components, might imply the presence of two distinct planetesimal belts separated by (a) gap(s) produced by (a) planet(s), keeping the intervening region free of dust. Planetary companions influence the dynamics of planetesimal belts, stirring their motions and initiating the required collisional cascade and produce the observed debris disc (Mustill & Wyatt 2009). Alternatively, the planetesimal belt may be stirred only by small bodies within the belt (Krivov & Booth 2018), or by some combination of these two factors (e.g. Muñoz-Gutiérrez, Marshall & Peimbert 2023).

The radial and vertical structure of debris belts can potentially be used to constrain the mass of planetary companions stirring or sculpting them (e.g. Mustill & Wyatt 2009; Krivov & Booth 2018; Daley et al. 2019; Marino 2021). Similarly, limits on the mass(es) of any potential companion(s) can be inferred from the overall architecture of the debris disc system (Pearce et al. 2022) or from the analysis of the gas kinematics. The power of constraints based on the disc architecture are limited by how well the belt is resolved and the nature of the interaction between companion and belt (i.e. stirring or

sculpting). Existing imaging observations have not spatially resolved the width of the outer belt of HD 138965 (Marshall et al. 2021), so the constraints derived from such calculations remain weak. Theoretical models suggest that the extent of HD 138965's debris disc is consistent with self-stirring by planetesimals, independent of the assumed age of the system (Krivov & Booth 2018). The current best limits come from high-contrast imaging observations of HD 138965 which rule out the presence of a 10 $M_J$ companion beyond 70 au from the star, assuming a stellar age of $348_{-54}^{+39}$ Ma (Matthews et al. 2018).

Here, we present Atacama Large Millimeter/Submillimeter Array (ALMA) Band 6 (1.3 mm) observations of HD 138965's debris disc, for the first time spatially resolving the system at millimetre wavelengths. We interpret these observations in conjunction with the scattered light non-detection to place upper limits on the dust albedo through radiative transfer modelling. We further infer mass limits on planetary companions in the vicinity of the disc based on the belt architecture and stellar age. The remainder of the paper proceeds as follows. In Section 2 we summarize the observations of HD 138965 that underpin the subsequent analysis. Then, in Section 3 we present the modelling approach, for both the host star and the disc architecture, before presenting our results. Our findings are then discussed in context in Section 4. Finally, we summarize our results and present our conclusions in Section 5.

## 2 OBSERVATIONS

Here we describe the various data sets obtained and combined to model HD 138965 and its accompanying debris disc. Most of these observations are obtained from publicly available catalogues or archive as science-ready, reduced products.

### 2.1 ALMA

Our target source HD 138965 (ICRS: RA = 15h 40m 11s.395, DEC = −70d 13m 41.47s) was observed by ALMA on 2019 December 18, during Cycle 6. The observations, under project code 2019.1.01220.S (PI: J.P. Marshall), were carried out using the 12-m array with a total of 42 antennas and employing a single pointing (no mosaic).

A total of four spectral windows (SPW) were used, including one centred at ~230.1 GHz to target the CO (2-1) molecular line transition. The remaining three SPWs were centred at the topocentric frequencies of 227.2, 243.1 and 245.1 GHz and were employed to observe continuum emission. All SPWs fall in the ALMA Band 6 frequency coverage. Observations are sensitive to emission with a largest angular scale of ~12.5 arcsec, while the field of view is approximately a factor of two larger. The instrument setup chosen used a bandwidth of 2 and 1.875 GHz, for the continuum and line SPWs, respectively. The corresponding channel spacings are ~15625 and 976.6 kHz. Calibration and imaging of the single data set were performed following standard procedures from the ALMA Observatory; during the observation, quasars J1427−4206 and J1558−6432 were used as bandpass/flux and phase/amplitude calibrators, respectively. The only data set obtained for our ALMA programme suffered from bad weather, resulting in amplitude drops for the flux calibrator during its observing scan. Affected data was removed, and a difference of at least 5 per cent in flux level was found before and after data flagging. This indicates a larger than expected flux uncertainty for the calibration, and therefore imaging, of HD 138965.

A total on-source time of 40 min was achieved for HD 138965. All imaging was done using a robust value of 0.5 in the Briggs
weighting scheme. The final continuum image has a sensitivity of 3 × RMS ~70 µJy, with a restored (synthesized) beam of 1.39 arcsec × 1.17 arcsec, at a position angle $\phi = 10^\circ 5$. CO (2–1) line emission cubes were also inspected, with a channel width of ~1.3 km s$^{-1}$ (RMS = 1.5 mJy beam$^{-1}$) and binned to 10 km s$^{-1}$ (RMS = 0.67 mJy beam$^{-1}$). Detailed data analysis and properties of the continuum and CO line emission are described in Section 3.2 and 3.3.

## 2.2 Ancillary data

We compile a comprehensive set of photometric and spectroscopic data spanning optical to millimetre wavelengths. This includes optical and near-infrared photometry from Hipparcos (Høg et al. 2000) and 2-Micron All Sky Survey (2MASS, Skrutskie et al. 2006), near- and mid-infrared photometry from Akari (Ishihara et al. 2010) and Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010), Spitzer IRS spectroscopy and MIPS photometry (Chen et al. 2014a), IRAS and Herschel PACS far-infrared photometry (Morales et al. 2016), and the newly acquired ALMA millimetre photometry. Photometric measurements were colour corrected before modelling following the relevant instrument handbooks. A summary of the compiled data are presented in Table 3. The Spitzer IRS spectrum was taken from the CASSIS archive (Lebouteiller et al. 2011) and scaled to the stellar photosphere model assuming a negligible infrared excess at 5 µm. Spectrophotometry using PYPHOT (Fouesneau 2022) was carried out to check for consistency between the rescaled spectrum and the WISE and MIPS photometric bands. The rescaled Spitzer IRS spectrum is consistent with the photometry.

These data complement the non-detection in scattered light from Very Large Telescope Spectro-Polarimetric High-contrast Exoplanet REsearch observations and the spatially resolved imaging observations in thermal emission at far-infrared (Herschel) and millimetre (ALMA) wavelengths. They will be used to infer the stellar properties, in conjunction with the Gaia-derived stellar parallactic distance, and thereafter modelled in conjunction with the images to determine the dust grain minimum size and size distribution exponent through radiative transfer modelling.

## 3 MODELLING AND RESULTS

Here, we detail the data analysis and its outcomes. We first determine the fundamental parameters of HD 138965 to identify an appropriate stellar atmosphere model. The orientation and extent of the debris disc are then calculated from the ALMA imaging observations with constraints from the previous, lower angular resolution Herschel data. The stellar model in conjunction with the dust spatial distribution are used as inputs to a 1-D radiative transfer code to infer the dust grain properties ($s_{min}$, $M_{dust}$, $q$) consistent with the observed thermal emission.

### 3.1 Fundamental stellar parameters

In order to derive the fundamental stellar parameters of HD 138965 (Gaia DR3 5819811291755502976; 2MASS J15401155–7013403; IRAS 15351–7003), we applied the Bayesian inference code of del Burgo & Allende Prieto (2016, 2018) on a grid of stellar evolution models constructed from the PARSEC v1.2S tracks (Bressan et al. 2012; Tang et al. 2014; Chen et al. 2014b, 2015). This choice is based on the good statistical agreement for detached eclipsing binaries, especially for stars on the main-sequence, where their measured dynamical masses and the corresponding predictions are consistent

<table><caption>Table 1. Fundamental stellar parameters of HD 138965.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Ref.</th></tr></thead><tbody><tr><td>Π (mas)</td><td>12.722 ± 0.028</td><td>Gaia DR3</td></tr><tr><td>M<sub>G</sub> (mag)</td><td>1.949 ± 0.006</td><td>Gaia DR3</td></tr><tr><td>G<sub>BP</sub> – G<sub>RP</sub> (mag)</td><td>0.087 ± 0.005</td><td>Gaia DR3</td></tr><tr><td>[Fe/H]</td><td>–0.3 ± 0.3</td><td>Zhang, Green & Rix (2023)</td></tr><tr><td>Effective temperature (K)</td><td>8881 ± 31</td><td>this work</td></tr><tr><td>Radius (R<sub>☉</sub>)</td><td>1.557 ± 0.009</td><td>this work</td></tr><tr><td>Mass (M<sub>☉</sub>)</td><td>1.90 ± 0.04</td><td>this work</td></tr><tr><td>Surface gravity (log g, cgs)</td><td>4.332 ± 0.012</td><td>this work</td></tr><tr><td>Luminosity (L/L<sub>☉</sub>)</td><td>13.59 ± 0.09</td><td>this work</td></tr><tr><td>Bolometric magnitude (mag)</td><td>1.907 ± 0.007</td><td>this work</td></tr><tr><td>Age (Ma)</td><td>31 ± 11</td><td>this work</td></tr></tbody></table>

on average to within 4 per cent (del Burgo & Allende Prieto 2018).

The grid of PARSEC v1.2S models arranged for this analysis comprises ages ranging from 2 to 13800 Myr and steps of 5 per cent, and [M/H] from –2.18 to 0.51 with steps of 0.02 dex, adopting the photometric passband calibration of Riello et al. (2021), with the zero points of the VEGAMAG system. We used it to infer the stellar parameters of HD 138965 through the code of del Burgo & Allende Prieto (2016, 2018), fed by the three following input parameters: the absolute G magnitude $M_G$, the colour $G_{BP} – G_{RP}$ and the iron to hydrogen [Fe/H]. Table 1 lists them. We reasonably assumed null extinction when deriving the colour $G_{BP} – G_{RP}$ and $M_G$ from the Gaia DR3 photometry and astrometry (Gaia Collaboration 2023). $M_G$ was determined from the apparent G magnitude by subtracting the distance modulus, which was calculated from the distance of 78.60 ± 0.17 pc corresponding to the Gaia DR3 trigonometric parallax. We adopted [Fe/H] = –0.30 ± 0.30 dex, derived by Zhang et al. (2023) from Gaia XP spectra.

The Galactic space coordinates in units of parsec of HD 138965 are (X, Y, Z) = (55.66 ± 0.11, –53.09 ± 0.11, –16.195 ± 0.032), which lies on the perimeter of the nearby, 40–50 Ma Argus association (see Zuckerman 2019). In addition, the Galactic space velocities in km s$^{-1}$ of HD 138965 are (U, V, W) = (–24.69 ± 0.24, –10.58 ± 0.23, –4.84 ± 0.07), where U, V, and W are positive toward the Galactic centre, the Galactic rotation, and the north Galactic pole, respectively. We determined the spatial and velocity parameters of HD 138965 from its Gaia DR3 Equatorial Coordinates (ICRS) at Ep = 2016.0, proper motion, radial velocity (–9.38 ± 0.30 km s$^{-1}$), and parallax. Zuckerman (2019) used Gaia DR2 data instead, deriving a mean Galactic space velocity in km s$^{-1}$ of (–22.5 ± 1.2, –14.5 ± 2.1, –5.0 ± 1.6), excluding HD 188728, CD-52 9381, and HD 192640 because they may not be members of Argus, as well as HD 129496, HD 145689, and even HD 138965 because of the relative large uncertainties in UVW, and adding HD 102647. We calculated a mean Galactic space velocity, in km s$^{-1}$, from Gaia DR3 data without adding or removing any star, arriving at (–23.1 ± 2.0, –14.5 ± 2.7, –4.9 ± 1.6). We note that Zuckerman (2019) adopted a radial velocity of –2.0 ± 4.3 km s$^{-1}$ for HD 138965, which is significantly lower and much less precise than the most recent value from Gaia DR3. A full discussion of the membership of all stars collected by Zuckerman (2019) is out of the scope of this paper, but we find convincing evidence in the Gaia DR3-based spatial and velocity distributions that HD 138965 is particularly a member of the Argus association. This is supported by the Banyan Sigma analysis which assigns a 99.4 per cent probability of its membership (Gagné et al. 2018). Imposing the prior that it is on
the pre-main sequence and using the aforementioned inputs to feed the Bayesian inference code, we inferred the fundamental parameters of HD 138965 presented in Table 1.

Our effective temperature ($T_{\text{eff}} = 8881 \pm 31$ K) and surface gravity (log $g = 4.332 \pm 0.012$) agree within a 2-$\sigma$ level with those of Zhang et al. (2023) ($T_{\text{eff}} = 9814 \pm 442$ K, log $g = 4.69 \pm 0.20$). Our age estimation of 31 $\pm 11$ Ma is consistent with that of Rhee et al. (2007) (20 Ma), who considered the UVW analysis of Zuckerman & Song (2004) as well as the location of HD 138965 on the Hertzsprung-Russell diagram.

As a validation test, we fed our Bayesian code with the same three input parameters but without imposing the prior on the evolution stage, finding that the age significantly increases to 805 $\pm 468$ Ma. This value is within $1-\sigma$ uncertainty of the value of $350_{-56}^{+37}$ Ma from David & Hillenbrand (2015). The similarity is even greater if we adopt Solar metallicity, since we arrive at an age of 338 $\pm 185$ Ma. However, the rotational broadening $v \sin i = 94.8 \pm 1.7$ km s$^{-1}$ from Gaia DR3 and the radius of HD 138965 lead to a maximum rotational period of 0.831 $\pm 0.016$ d that supports the youth of this host star.

## 3.2 Millimetre continuum emission

We model the disc architecture as a single Gaussian belt. The model is defined by a flux density $f_{\text{disc}}$, radius $R_{\text{peak}}$, width (standard deviation) $\sigma_R$, scale height $h = z/R$, inclination $i$, and position angle $\phi$. Additionally, we include parameters for the stellar photospheric contribution $f_*$, and position offsets $\Delta RA$, $\Delta Dec$ between the star (assumed to be at the phase centre of the observations) and the disc centre. Previous marginally resolved Herschel observations of the disc at far-infrared wavelengths provide good constraints on the radius, inclination, and position angle for the priors for these parameters (Morales et al. 2016; Marshall et al. 2021). These priors are normally distributed with a standard deviation 3 per cent of the initial value and bounds informed by the likely range of each parameter.

We explore the parameter space of the model in a Bayesian framework using EMCEE (Foreman-Mackey et al. 2013), following the approach presented in previous, similar works (e.g. Marshall et al. 2018, 2023a). Walkers are initialized at the maximum likelihood values for the disc parameters determined from modelling the Herschel observations. Each set of parameters is used to generate a disc model using the radiative transfer code RADMC-3D (Dullemond et al. 2012) which is then converted into synthetic visibilities using GALARIO (Tazzari, Beaujean & Testi 2018). These model visibilities are binned and compared to the ALMA observations through least-squares fitting. In total, we created 90 000 realizations of the model (90 walkers, 1000 steps). The final posterior probability distribution used the last 10 000 realizations, from which the maximum likelihood parameters and their uncertainties were obtained. The posteriors are all monomodal and approximately normally distributed. We therefore take the 50th percentile of the posterior of each parameter as its model value and calculate the uncertainties from the 16th and 84th percentiles of each parameter distribution. We present the marginalized posteriors in Appendix B. The results of the model fitting are summarized in Table 2, and we show the observations and residuals (after subtraction of the maximum likelihood model) in Fig. 1.

We find that HD 138965's disc has a total flux density of 1.46 $\pm 0.23$ mJy at a wavelength of 1.269 mm. The cold belt is spatially resolved, with a radius of $150_{-7}^{+10}$ au and a width of $49_{-6}^{+7}$ au

<table><caption>**Table 2.** Results from visibility fits to the ALMA observation. We note that the relative scale height $h = (z/R)$ is unconstrained from the modelling and the stellar flux density $f_*$ is an upper limit.</caption><thead><tr><th>Parameter</th><th>Range</th><th>Value</th></tr></thead><tbody><tr><td>R<sub>peak</sub> (au)</td><td>50–300</td><td>150<sup>+10</sup><sub>-7</sub></td></tr><tr><td>σ<sub>R</sub> (au)</td><td>10–100</td><td>49<sup>+7</sup><sub>-6</sub></td></tr><tr><td>h</td><td>0.01–0.30</td><td>0.15<sup>+0.07</sup><sub>-0.08</sub></td></tr><tr><td>i (°)</td><td>0–90</td><td>49.9<sup>+3.3</sup><sub>-3.7</sub></td></tr><tr><td>φ (°)</td><td>0–180</td><td>173.3<sup>+3.7</sup><sub>-4.1</sub></td></tr><tr><td>f<sub>dust</sub> (mJy)</td><td>0.1–10.0</td><td>1.460 ± 0.230</td></tr><tr><td>f<sub>star</sub> (µJy)</td><td>0–300</td><td>&lt; 70</td></tr><tr><td>ΔRA (′′)</td><td>-1–1</td><td>0.11<sup>+0.5</sup><sub>-0.6</sub></td></tr><tr><td>ΔDec (′′)</td><td>-1–1</td><td>0.16<sup>+0.7</sup><sub>-0.7</sub></td></tr></tbody></table>

(FWHM = $115_{-14}^{+16}$ au). The radius is consistent with the Herschel-derived extent, and the width of the debris belt is marginally resolved with a fractional width $\Delta R/R = 0.77$. This places HD 138965 amongst the more radially extended discs that have been observed today spanning 0.1 to $>1.0$ (Matrà et al. 2025), and in comparison to the Solar system's Edgeworth-Kuiper belt ($\Delta R/R \simeq 0.3$ Kavelaars et al. 2009; Petit et al. 2011). A $2-\sigma$ brightness asymmetry between the NW and SE ansae of the disc is seen in the imaging data. Given the large uncertainties and calibration of these observations it is not necessarily a real feature of the disc, but should motivate future observations. The large relative vertical scale height, $h = 0.15_{-0.08}^{+0.07}$, might be taken as a further indication of the disc being stirred. The associated uncertainty, in conjunction with the moderate inclination and modest signal-to-noise, instead leads us to the interpretation that a broad range of scale heights tested by the model will adequately fit the observations. This motivates the need for higher spatial resolution, higher signal-to-noise observations in future to constrain the internal dynamics of this disc. Its inclination, $i = 49.9_{-3.7}^{+3.3}$, and position angle, $\phi = 173.3_{-4.1}^{+3.7}$, are the most precise measurements of the disc orientation, both consistent with previous estimates. As expected we do not detect either emission from the warm belt or the star and obtain a $3-\sigma$ upper limit of 90 µJy on the stellar photospheric emission. Furthermore, we find no evidence for an offset between the star (phase centre) and disc centre from these observations although the constraints on such an offset are low due to the limited spatial resolution of the observations.

## 3.3 Millimetre CO line emission

We use the spatially resolved extent and orientation of the disc from the continuum image to identify the region of the spectral window covering the CO (2-1) line at 230.538 GHz to search for gas emission from the disc. A spatial mask is applied to the spectral cube, with pixels lying within the region from $R_{\text{peak}}-\text{HWHM}$ to $R_{\text{peak}}+\text{HWHM}$ and with a flux density $\ge 3-\sigma$ being included in the summation. For each pixel included in the mask we shift the spectrum's frequency according to the projected Keplerian velocity, interpolate these individual spectra to a set of common velocities and sum these to produce a final CO spectrum. The resulting spectrum centred on the CO (2-1) line at 230.54 GHz is presented in Fig. 2, with a shaded region denoting the frequency space included in the search for CO line emission centred at the stellar radial velocity of $-9.38$ km s$^{-1}$ taken from Gaia DR3 (Gaia Collaboration 2023). We find no evidence for any CO emission associated with the disc, and
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Left: ALMA Band 6 continuum image of HD 138965. The image has been cleaned and reconstructed with a Briggs weight of 0.5. The instrument beam (1.39 arcsec × 1.17 arcsec, φ = 10.5°) is denoted by the white ellipse in the bottom left corner. Contours are in steps of 2−σ from +2−σ. Orientation is north up, east left. Right: Residuals after subtraction of the maximum likelihood model from the observations. We see there two features of interest in the residual map. A 2−σ bump is located at the position of the northern ansa of the ring and is responsible for the observed ‘asymmetry’, and another 2−σ bump is located south of the ring exterior to the disc location.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** ALMA CO spectrum centred on a rest frequency of 230.538 GHz. The shaded region denotes the 10 km s⁻¹ wide channel used to determine the presence, or rather absence, of significant CO emission, centred on the stellar radial velocity from Gaia DR3 (i.e. −9.38 km s⁻¹). Uncertainties in the unbinned spectrum channels (1.3 km s⁻¹ wide) are denoted by the vertical lines. The horizontal dashed lines denote the uncertainty in the unbinned spectrum from +2σ in steps of 1−σ.</figcaption></figure>

obtain a 5−σ upper limit of 3.35 mJy beam⁻¹ in a 10 km s⁻¹ wide channel, equivalent to 2.7 × 10⁻²³ W m⁻², which is around a factor of three above the predicted CO emission level of 1.1 × 10⁻²³ W m⁻² (MCO = 8.6×10⁻⁷ M⊕) of Kral et al. (2017).

## 3.4 Spectral energy distribution

We model HD 138965’s spectral energy distribution (SED) as a star plus two dust components, using two separate parametric models to represent the warm excess present at mid-infrared wavelengths,

<table><caption>**Table 3.** Photometry used in the radiative transfer modelling.</caption><thead><tr><th>Wavelength (µm)</th><th>Flux density (mJy)</th><th>Telescope/ Instrument/Filter</th><th>Reference</th></tr></thead><tbody><tr><td>0.44</td><td>10 700 ± 200</td><td>Johnson B</td><td>1</td></tr><tr><td>0.55</td><td>9838 ± 186</td><td>Johnson V</td><td>1</td></tr><tr><td>0.58</td><td>8680 ± 20</td><td>Gaia G</td><td>2</td></tr><tr><td>1.24</td><td>4654 ± 172</td><td>2MASS J</td><td>3</td></tr><tr><td>1.65</td><td>3052 ± 190</td><td>2MASS H</td><td>3</td></tr><tr><td>2.16</td><td>2059 ± 58</td><td>2MASS K<sub>s</sub></td><td>3</td></tr><tr><td>3.4</td><td>921 ± 88</td><td>WISE W1</td><td>4</td></tr><tr><td>4.6</td><td>551 ± 19</td><td>WISE W2</td><td>4</td></tr><tr><td>9</td><td>206 ± 10</td><td>Akari IRC9</td><td>5</td></tr><tr><td>12</td><td>103 ± 4</td><td>WISE W3</td><td>4</td></tr><tr><td>13</td><td>134 ± 7</td><td>Spitzer IRS</td><td>6</td></tr><tr><td>22</td><td>80 ± 5</td><td>WISE W4</td><td>4</td></tr><tr><td>24</td><td>81 ± 2</td><td>Spitzer MIPS</td><td>6</td></tr><tr><td>31</td><td>111 ± 2</td><td>Spitzer IRS</td><td>6</td></tr><tr><td>60</td><td>482 ± 63</td><td>IRAS 60</td><td>7</td></tr><tr><td>70</td><td>561 ± 29</td><td>Spitzer MIPS</td><td>6</td></tr><tr><td>100</td><td>527 ± 17</td><td>Herschel PACS</td><td>8</td></tr><tr><td>160</td><td>299 ± 13</td><td>Herschel PACS</td><td>8</td></tr><tr><td>1270</td><td>1.46 ± 0.23</td><td>ALMA Band 6</td><td>9</td></tr></tbody></table>

Note. References: 1. Høg et al. (2000) 2. Gaia Collaboration (2023) 3. Skrutskie et al. (2006) 4. Wright et al. (2010) 5. Ishihara et al. (2010) 6. Chen et al. (2014a) 7. Beichman et al. (1988) 8. Morales et al. (2016) 9. This work.

and the cold excess at mid-infrared to millimetre wavelengths. A summary of the photometry used in the radiative transfer modelling is provided in Table 3.

We obtained the best BT-NextGen stellar atmosphere model (AGSS2009, Allard, Homeier & Freytag 2012) that fits the HD
<table><caption>**Table 4.** Assumed inner belt parameters for HD 138965.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>s<sub>min</sub> (μm)</td><td>12</td></tr><tr><td>T<sub>dust</sub> (K)</td><td>168</td></tr><tr><td>q</td><td>3.5</td></tr><tr><td>M<sub>s</sub> (×10<sup>-3</sup> M<sub>⊕</sub>)</td><td>0.154</td></tr><tr><td>r<sub>m</sub> (au)</td><td>13.7</td></tr><tr><td>σ<sub>r</sub> (au)</td><td>1.75</td></tr><tr><td>Composition</td><td>Astronomical silicate</td></tr></tbody></table>

138965’s photometry from a classical least-squares approach. The fundamental stellar parameters derived are the effective temperature $T_{\text{eff}} = 8800$ K, metallicity [Fe/H] = -0.5, and the surface gravity of log g = 4.5, consistent with the values inferred in Section 3.1. This model was scaled to the available optical and near-infrared photometric data up to 5 µm.

To model the thermal emission of the grains we applied STARDUST, a model that takes into account the dust composition (determined by optical constants) and a combination of forces that influence the spatial redistribution of dust grains. These forces include the conservative forces of gravity competing with the outward radiation pressure (e.g. Burns, Lamy & Soter 1979), and the non-conservative force of the Poynting-Robertson effect (e.g. Poynting 1903; Robertson 1937; Klačka et al. 2007) that causes grains to spiral inwards. STARDUST includes modelling a simple collisional lifetime of the grains that depends on the optical depth of the debris disc (e.g. Kuchner & Stark 2010). We used the fractional luminosity of HD 138965’s disc as a proxy for the optical depth. We summarise the *Stardust* model in Appendix A.

For our purposes, we assumed the radial planetesimal surface density to follow a Gaussian profile and differential power-law grain size distribution (dn/ds ∝ s⁻q). In all considered scenarios, we assume the maximum size of the dust grains to be 3000 µm and the inner belt to have a grain composition of astronomical silicate (Draine 2003). The outer, cooler belt was assumed to have a fixed peak surface density at 150 au with a standard deviation of 49 au (taken directly from the ALMA observations, see Table 2).

We note that, since the inner belt is not spatially resolved, there will be a degeneracy between location and the temperature of the grains. We performed an initial analysis using STARDUST to estimate the inner belt characteristics. As such we had the following components to be free parameters for both the inner and outer belts: belt radius ($r_m$) and width ($\sigma_r$), minimum grain size ($s_{min}$), dust mass ($M_s$), and the size distribution component ($q$) to initially vary between 3 and 4 as is the case for collisional active discs (e.g. MacGregor et al. 2016; Marshall et al. 2017; Norfolk et al. 2021). Here we found that the inner belt had a minimum grain size of approximately 12 µm, however, as the value of $q$ did not converge.

To pin down the radial component of the inner belt, we ran STARDUST with the following assumptions for the inner belt: minimum grain size of 12 µm, size distribution component to be 3.5 (steady-state collisional cascade; Dohnanyi 1969), and the disc fractional width ($\Delta r/r$) to be 0.3. The mean radial distance of the inner belt and the mass were kept as free parameters. This scenario still did not converge to a solution, but it did produce a reasonable goodness-of-fit to the data where our reduced minimum $\chi^2$ test was less than one. The posterior probability distribution for the corresponding EMCEE run is presented in Fig. B2 (Appendix B). A summary of the fixed parameters that have been assumed for the inner belt component of HD 138965 are given in Table 4.

Now that we have our fixed assumptions for the inner belt, we then modelled the outer belt with different compositions and porosities. For all scenarios, we set the minimum grain size in the outer belt to be the blowout limit of the corresponding grain composition. For the outer belt component ($q$), we allowed to vary between 2 and 5 to account for potential size-dependent velocity distributions for all bodies (e.g. Pan & Schlichting 2012). The different compositions considered were a mixture of astronomical silicate with either crystalline water-ice (Warren & Brandt 2008), amorphous water-ice (Hudgins et al. 1993), amorphous carbon (Preibisch et al. 1993), or vacuum (as porosity). For each mixture (or inclusion matrix particles; IMPs), we assumed the grains to be astronomical silicate combined with an inclusion material with a volume fraction of 0.5, 0.25 and 0.1 for each material.

To determine the maximum likelihood values for each model we again adopt a Bayesian approach using the Markov Chain Monte Carlo (MCMC) code EMCEE (Foreman-Mackey et al. 2013) to explore the parameter space and determine the posterior probability distributions for each parameter from which the maximum amplitude probability and uncertainties were determined. We ran a total of 320 000 realizations for each model with 80 walkers and 4000 steps. We used the first 500 steps as burn-in for the MCMC chains and calculated probability distributions from the final 280 000 realizations. We used the 50th value from the posterior of each parameter to determine the reduced $\chi^2$ of a given architecture and dust composition for comparison between dust composition models. For each of these scenarios, all of them converged to a solution before reaching the 4000th step.

To determine which model to be the best fit, we applied the Bayesian Information Criterion (BIC; Schwarz 1978):

$$ \text{BIC} = \chi^2 + k \ln N, \quad (1) $$

where ‘k’ is the number of free parameters and ‘N’ is the number of data points used in the fit. In our scenarios, we have $k = 3$ (for outer belt composed of only astronomical silicate) or $k = 4$ (for outer belt composed of mixed grains), and $N = 15$. The lower the BIC value the better the fit. We then calculate the change of BIC ($\Delta$BIC; Burnham & Anderson 2004):

$$ \Delta \text{BIC} = \text{BIC}_{\text{mc}} - \text{BIC}_{\text{best}}, \quad (2) $$

where ‘mc’ is the model being compared with the ‘best’ model. In our case, the best model as determined by ‘BIC’ is the scenario where the outer belt is composed of pure astronomical silicate dust grains. If $\Delta$BIC < 2, the models are indistinguishable; 2 < $\Delta$BIC < 6 there is some support for the ‘best’ model but not conclusive; 6 < $\Delta$BIC < 10 there is stronger evidence to support the best model; ands $\Delta$BIC > 10 indicates that the ‘best’ model is highly preferred. See Table 5 for a summary of the model outcomes.

On this basis, the model corresponding to $P = 0.25$ is a reasonable fit to the observations, while 90:10 astronomical silicate:amorphous water-ice, and the remaining porosity scenarios are all plausible dust compositions. See Fig. 3 for the corresponding SED to the scenario of 90:10 astronomical silicate:amorphous water-ice. The 90:10 and 75:25 of astronomical silicate:crystalline water-ice were not categorically ruled out based on the $\Delta$BIC values. We can however rule out the 50:50 compositions for both amorphous and crystalline water ice, and all three of the carbon:astronomical silicate compositions tested here. It should be noted that the emission from the inner belt is fixed following the initial modelling where the outer belt grains are assumed to be solely astronomical silicate. This may introduce an inherent bias toward astronomical silicate being the best solution for the outer belt. Ultimately, we would need to either
<table><caption>**Table 5.** Input parameters and the corresponding goodness-of-fit values from the radiative transfer modelling. S: astronomical silicate; WC: crystalline water ice; WA: amorphous water ice; C: carbon; P: porosity; with the volume fraction of inclusion material mixed with the astronomical silicate are 0.10, 0.25 and 0.50, respectively. NB: Dust masses are for grains up to 3 mm in size.</caption><thead><tr><th>Parameter</th><th>S100</th><th>S90WC10</th><th>S75WC25</th><th>S50WC50</th><th>S90WA10</th><th>S75WA25</th><th>S50WA50</th><th>S90C10</th><th>S75C25</th><th>S50C50</th><th>P = 0.1</th><th>P = 0.25</th><th>P = 0.5</th></tr></thead><tbody><tr><td>s<sub>blow</sub> (μm)</td><td>2.70</td><td>2.86</td><td>3.16</td><td>4.04</td><td>2.86</td><td>3.16</td><td>3.87</td><td>2.86</td><td>3.04</td><td>3.46</td><td>2.93</td><td>3.41</td><td>4.80</td></tr><tr><td>q [2.0–5.0]</td><td>3.73<sup>+0.07</sup><sub>−0.07</sub></td><td>3.82<sup>+0.09</sup><sub>−0.08</sub></td><td>3.86<sup>+0.11</sup><sub>−0.10</sub></td><td>4.01<sup>+0.20</sup><sub>−0.17</sub></td><td>3.78<sup>+0.08</sup><sub>−0.08</sub></td><td>3.83<sup>+0.10</sup><sub>−0.09</sub></td><td>4.00<sup>+0.18</sup><sub>−0.15</sub></td><td>4.15<sup>+0.10</sup><sub>−0.09</sub></td><td>4.41<sup>+0.10</sup><sub>−0.09</sub></td><td>4.89<sup>+0.08</sup><sub>−0.16</sub></td><td>3.63<sup>+0.07</sup><sub>−0.07</sub></td><td>3.49<sup>+0.07</sup><sub>−0.06</sub></td><td>3.30<sup>+0.07</sup><sub>−0.07</sub></td></tr><tr><td>M<sub>s</sub> [1–20]</td><td>7.74<sup>+0.49</sup><sub>−0.48</sub></td><td>7.75<sup>+0.55</sup><sub>−0.54</sub></td><td>7.18<sup>+0.58</sup><sub>−0.56</sub></td><td>5.86<sup>+0.70</sup><sub>−0.62</sub></td><td>7.44<sup>+0.51</sup><sub>−0.50</sub></td><td>7.07<sup>+0.55</sup><sub>−0.54</sub></td><td>6.12<sup>+0.65</sup><sub>−0.61</sub></td><td>7.05<sup>+0.43</sup><sub>−0.42</sub></td><td>9.20<sup>+0.76</sup><sub>−0.74</sub></td><td>8.89<sup>+0.60</sup><sub>−0.63</sub></td><td>6.67<sup>+0.43</sup><sub>−0.43</sub></td><td>5.23<sup>+0.35</sup><sub>−0.34</sub></td><td>3.18<sup>+0.21</sup><sub>−0.21</sub></td></tr><tr><td>(×10<sup>−3</sup> M<sub>⊕</sub>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>χ<sup>2</sup></td><td>9.3</td><td>14</td><td>17</td><td>150</td><td>11</td><td>28</td><td>100</td><td>19</td><td>20</td><td>120</td><td>9.4</td><td>8.3</td><td>10</td></tr><tr><td>χ<sup>2</sup><sub>red</sub></td><td>0.66</td><td>0.90</td><td>1.1</td><td>9.8</td><td>0.78</td><td>1.9</td><td>6.9</td><td>1.3</td><td>1.3</td><td>7.7</td><td>0.62</td><td>0.55</td><td>0.67</td></tr><tr><td>ΔBIC</td><td>0.0</td><td>7.0</td><td>10</td><td>140</td><td>4.4</td><td>21</td><td>96</td><td>12</td><td>14</td><td>110</td><td>2.8</td><td>1.7</td><td>3.4</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Spectral energy distribution of HD 138965 for the scenario where grains subjected to forces of gravity, radiation pressure, and drag induced by the Poynting–Robertson effect. The inner and outer belts were composed of 90:10 astronomical silicate:amorphous water-ice. Data points are literature photometry, except for the ALMA data point. The *Spitzer/IRS* spectrum is shown from 5 to 33 µm. The stellar photosphere model is denoted by the dotted line, the dashed-dot lines denote the cold and warm disc components, and the total emission by the solid line.</figcaption></figure>

independently determine the architecture of the inner belt through spatially resolved imaging or obtain spectroscopic data for the dust revealing features associated with water ice or other species to better infer the composition of the debris dust.

Through our initial analysis using MCMC, we found the value of the inner belt's mean radial location to be at $15^{+3}_{-2}$ au, and minimum grain size $s_{min}$ of 12 µm, broadly consistent with the values reported in Morales et al. (2016). The $s_{min}$ adopted for the inner belt is much larger than that of the outer belt. We would expect $s_{min}$ to be a few times the blowout size and for a 13.6 $L_{\odot}$ star, that ratio would be $s_{min} \simeq 2.5s_{blow}$ or 7 µm (Pawellek et al. 2014). The effective temperature of 12 µm astronomical silicate dust grains is 168 K at 15 au from HD 138965, consistent with the sublimation temperature of water ice in vacuum (e.g. Kossacki 2021). This supports our assumption of a dry inner belt; if the minimum grain size were smaller the dust would be hotter such that the major implication would be the inner belt would need to be placed further out to match the observed mid-infrared emission. Our study of the outer belt composition implemented a refined model compared to Morales et al. (2016), where they assumed a single composition of 50:50 astronomical silicate:ice for the dust grains and a fixed density distribution. Our best-fit results correspond to the porosity scenarios and have a strong candidate with 90:10 astronomical silicate:crystalline water-ice mixture for the outer belt.

## 3.5 Constraints from scattered light

The large extent and broad width of HD 138965's debris disc combined with its intermediate inclination ($i = 49.9$) result in a low-surface brightness for the disc, despite its relatively high fractional luminosity ($L_d/L_* \simeq 4 \times 10^{-4}$). This makes it a challenging target for detection in scattered light. Previous attempts have been made with *HST/ACS* at optical wavelengths and with *VLT/SPHERE* in the near-infrared. These have both failed to detect any emission from the debris disc. Using eq. (2) from Marshall et al. (2018) in conjunction with the resolved disc architecture and known stellar brightness we can calculate a limit on the dust albedo, $\omega$. A point of caution must be made here in that the fractional luminosity ($L_d/L_*$) which is used to infer the albedo in Marshall et al. (2018) may be a poor approximation of the absorption efficiency at the observing wavelength of scattered light imaging. The dust albedo calculated in that work is a bulk measure across all dust grain sizes, not the albedo for a specific wavelength or grain size.

HD 138965 was observed by VLT/SPHERE IRDIS (Beuzit et al. 2019) using the BB_H filter for the SPHERE High-Angular Resolution Debris Disks Survey<fnref n="1" /> (SHARDDS; Wahhaj et al. 2016; Choquet et al. 2017; Marshall et al. 2018; Cronin-Coltsmann et al. 2021). These observations are deeper than those presented in Matthews et al. (2018), but the disc is still not detected. From the SHARDDS observation we obtain a $5-\sigma$ upper limit of 83 µJy arcsec$^{-2}$ at a separation of 2 arcsec from the star. To estimate the equivalent integrated disc brightness required to be detectable, we inject model discs with an architecture consistent with the ALMA observation into the data cube. The cube was then reduced and stacked using the classical ADI algorithm. The model brightness was varied until the integrated stack exhibited residuals in excess of the $5-\sigma$ limit. We find a minimum brightness of around 10 mJy is required for the disc to be detected in the SHARDDS data. Based on this calculation, we would need to go at least a factor of four deeper in imaging HD 138965 to detect the disc at near-infrared wavelengths. This assumes reducing the limit of the albedo to $\le 0.2$ would start to probe values consistent with other discs at near-infrared wavelengths (e.g. Ren et al. 2023).

We further estimate the contrast for point sources around HD 138965 at radial separations out to $\simeq 1$ arcsec. This contrast curve was calculated using the Regime Switching Model algorithm as presented in Dahlqvist et al. (2022). We convert the observed contrast to mass in order to determine the limits on planetary companions in the vicinity of the disc, as shown in Fig. 4. The mass limit as a function of separation was calculated for ages of 30 ± 10 Ma (our inferred age), 45 ± 5 Ma (Argus association membership), and $348^{+39}_{-54}$ Ma (from

<footnote n="1">¹ESO programs 096.C-0388 and 097.C-0394</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Mass limits from VLT/SPHERE SHARDDS observation of HD 138965 based on AMES-COND evolutionary models (Allard et al. 2001). Contrast as a function of separation was calculated using the RSM algorithm as presented in Dahlqvist et al. (2022). The dot, dash, and solid lines represent different assumed ages for the star, with the coloured shaded regions being the associated uncertainties. The shaded regions denote the expected extent of the two debris belts. We find limits at 80 au of 4 $M_{\text{Jup}}$ for an age of 30 Ma, or 5.5 $M_{\text{Jup}}$ for an age of 45 Ma.</figcaption></figure>

Gaia). The masses were converted from the measured contrast using the AMES-COND grid of models (Allard et al. 2001). In the vicinity of the inner edge of the cold debris belt (≈80 au) we obtain upper limits to the mass of any companion of 4 $M_{\text{Jup}}$ ($t_{\star} = 31$ Ma), 5 $M_{\text{Jup}}$ ($t_{\star} = 45$ Ma), and 15 $M_{\text{Jup}}$ ($t_{\star} = 348$ Ma).

A surface brightness limit of 20 mag arcsec⁻² at 2 arcsec (≡ 33 µJy arcsec⁻²) was obtained in the HST/ACS F606W filter (Krist et al. 2010). The HST/ACS measurement is more constraining than the VLT/SPHERE measurement, resulting in a limit of $\omega \le 0.09$ at optical wavelengths. This is consistent with the minimum grain size in the outer belt derived from modelling the SED, assuming an astronomical silicate composition for the dust. The calculated albedo limit for HD 138965’s disc is intermediate for those discs detected in scattered light (0.02 to ≥0.2; Choquet et al. 2018; Ren et al. 2023), from which we infer deeper observations would stand a good chance of detecting the disc.

We can obtain an alternative measure of the dust scattering albedo using the maximum likelihood dust model determined in Section 3.4. Here we take the $Q_{\text{abs}}$ and $Q_{\text{sca}}$ values derived from the dirty ice dust composition for grains with $s = 3.6$ µm, being the minimum grain size in our distribution and therefore the largest contribution to the disc surface area. We infer a scattering albedo of $\omega = 0.56$ over all angles at a wavelength of 1.6 µm, using the relation $\omega = Q_{\text{sca}}/(Q_{\text{abs}} + Q_{\text{sca}})$ (Bohren & Huffman 1998). Considering instead only a restricted range of scattering angles close to 90°, based on the orientation of the disc, we instead infer an effective dust albedo ≈0.1 for astronomical silicate dust grains around 10 µm in size (following e.g. Mulders et al. 2013; Choquet et al. 2018). This would result in non-detection of the disc and be marginally consistent with upper limits from HST.

## 4 DISCUSSION

The outer belt of HD 138965’s debris disc is found to be very broad with $\Delta R/R = 0.77$. This is far greater than the Solar system’s Edgeworth–Kuiper belt ($\Delta R/R \simeq 0.3$; Kavelaars et al. 2009; Petit et al. 2011) and comparable to HR 8799’s outer belt (Faramaz et al. 2021). A handful of radially broad discs followed-up at sufficient angular resolution have been found to contain substructures associated with disc–planet interaction believed to represent sculpting or scattering of planetesimals (e.g. Marino 2021). The pursuit of further observations of HD 138965 at higher angular resolution to determine the architecture of the outer belt are therefore recommended.

Despite the modest angular resolution and signal-to-noise of the ALMA observations the outer belt is spatially resolved in both radius and width. We can therefore place constraints on the presence of a planetary companion stirring or sculpting the belt following the analysis presented in Pearce et al. (2022). For these calculations we assume the disc inner edge to lie at 101 au (i.e. $R - \sigma_R$) and the stellar parameters from Table 1. For the single planet sculpting scenario, we calculate a planetary mass $\ge 2.3 \pm 0.4$ $M_J$ at a semimajor axis $\le 74^{+8}_{-6}$ au. This is a factor of five improvement over the previous limit from Matthews et al. (2018), and a factor of two better than obtained from the deeper SHARDDS image presented in Section 3.5, but still relatively weak. At the other extreme, assuming the disc is self-stirred by massive planetesimals, we can constrain the total mass of the debris disc to be at least 400 $\pm 100$ $M_\oplus$ in bodies up to 100 km in size. This mass is consistent with the maximum allowable debris disc masses as a function of disc properties (e.g. radius, fractional luminosity) determined by Krivov & Wyatt (2021).

Our analysis of the dust dynamics has been modelled by STARDUST, where we include P–R effected grains (i.e. inward spiralling of dust grains) to eventually result in collisional destruction on its way to the host star. Due to the inner belt not being spatially resolved we attempted to mitigate the degeneracy between grain temperature and location, we performed an initial analysis to estimate its characteristics under the assumption that the outer belt grains are composed of pure astronomical silicate. The mean values taken from the posterior probability distributions for this scenario, we estimated the location of the inner belt to be mean distance of $15^{+3}_{-2}$ au, this result is in line with the study conducted by Morales et al. (2016). With all subsequent scenarios modelled, we fixed the minimum grain size of the outer belt to be the blowout limit of the grain corresponding to the relative dust composition. We note that under the assumptions of inner belt characteristics there may be an inherent bias towards the best model where pure astronomical silicate grains are in the outer belt. On this basis, we cannot rule out the scenarios where there is at least 10 per cent crystalline or amorphous water-ice inclusions in the astronomical silicate grains in the outer belt. However, we can eliminate grains for the scenarios with 50 per cent water-ice inclusions and any grain that has amorphous carbon mixed in for HD 138965’s outer belt. To be sure of the inner belt radial location, higher spatial resolution images would be required. For example, at a radius of ≈20 au the inner belt would have an angular separation of 0.25 arcsec which is potentially resolvable by JWST/MIRI in imaging mode at 10 µm (Wright et al. 2023), or in future with ATLAST at 350 µm (Booth et al. 2024; Klaassen et al. 2024).

## 5 CONCLUSIONS

We have spatially resolved for the first time the radius and width of the outer cool belt of HD 138965’s debris disc at millimetre wavelengths with ALMA. We find the peak emission from the belt at 150 au
with a standard deviation of 49 au (FWHM 115 au, $ΔR/R \simeq 0.77$). The disc radius, inclination (49°9), and position angle (173°3) of the belt are in good agreement with previous marginally resolved Herschel observations. The disc is very broad, comparable in width to several systems that have exhibited indirect evidence of disc-planet interactions. However, the low signal-to-noise and relatively low angular resolution of the ALMA observations preclude any more detailed insights regarding the belt structure or its stirring.

Using the newly acquired ALMA imaging of the disc we modelled the disc spectral energy distribution using a novel code to redistribute dust grains based on their size and the gravitational and radiation forces acting upon them. As the inner belt is not spatially resolved, we assumed that both the inner and outer belt grains were composed of pure astronomical silicate to constrain the system's inner belt by determining its location to be $15_{-2}^{+3}$ au to help mitigate the degeneracy between the temperature and location for grains. We found, whilst the pure astronomical silicate grains were considered the best scenario for the composition of the dust in the outer belt, we could not categorically disregard scenarios with a 10 per cent crystalline or amorphous water-ice inclusions in the silicate grains.

We then undertook a study of archival scattered light observations from VLT/SPHERE to obtain some constraint on the dust optical properties, but the disc remains undetected. Using the non-detection by HST at optical wavelengths in conjunction with the spatially resolved architecture of the disc from ALMA we calculate an upper limit to the dust albedo of 0.09. That limit from observations is inconsistent with the same limit determined from the dust grain size and composition inferred from the SED. These two values can be somewhat reconciled by considering either a restricted range of scattering angles to reduce the inferred albedo from HST imaging, or instead invoking the broad belt width and self-subtraction inherent to the data reduction process for its non-detection in SPHERE imaging.

Despite non-detection of the disc in scattered light, analysis of these data provided deeper constraints on the mass of companions to the system. Interpretation of the contrast curve derived from the SHARDDS observation of HD 138965 places a mass limit of $4 M_{Jup}$ at 80 au from the star, assuming an age of 30 Ma. Using the architecture of the disc we instead obtain a mass limit of $2.3 \pm 0.4 M_{Jup}$ for a single planet sculpting the disc using the estimates of Pearce et al. (2022). However, those predictions may underestimate the planetary mass required to stir a disc once the mass of the belt itself is taken into account as demonstrated by Sefilian (2024). This comparison illustrates the predictive power of debris disc architectures for the inference of planetary companions, albeit with caution regarding the depth that can be achieved.

Future avenues for the characterisation of the HD 138965 system with current and planned facilities are promising. For example, in the infrared JWST is able to resolve both the inner and outer belts of the debris disc and simultaneously provide deep limits on the presence of planetary companions. At millimetre wavelengths, deeper observation with ALMA to resolve the outer belt's structure is possible, albeit demanding, in terms of the time required with its current sensitivity; future bandwidth upgrades will make such high-angular resolution observations more feasible.

## ACKNOWLEDGEMENTS

The authors thank the anonymous referee for their constructive and helpful comments.

This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). This research has also made use of NASA's Astrophysics Data

System. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.01220.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC 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 NAOJ.

JPM acknowledges support by the National Science and Technology Council of Taiwan under grant NSTC 112-2112-M-001-032-MY3.

CdB acknowledges support from the Agencia Estatal de Investigación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) under grant WEAVE: EXPLORING THE COSMIC ORIGINAL SYMPHONY, FROM STARS TO GALAXY CLUSTERS and the European Regional Development Fund (ERDF) with reference PID2023-153342NB-I00/10.13039/501100011033, as well as from a Beatriz Galindo Senior Fellowship (BG22/00166) from the MICIU. The Universidad de La Laguna (ULL) and the Consejería de Economía, Conocimiento y Empleo of the Gobierno de Canarias are also gratefully acknowledged for the support provided to CdB (2024/347).

FK acknowledges support from the Spanish Ministry of Science, Innovation and Universities, under grant number PID2023-149918NB-I00. This work was also partly supported by the Spanish program Unidad de Excelencia María de Maeztu CEX2020-001058-M, financed by MCIN/AEI/10.13039/501100011033.

EC acknowledges funding from the European Union (ERC, ESCAPE, project no. 101044152). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

VF-G acknowledges funding from the National Aeronautics and Space Administration through the Exoplanet Research Program under Grant Nos. 80NSSC21K0394 and 80NSSC23K1473 (PI: S. Ertel), and Grant No. 80NSSC23K0288 (PI: V. Faramaz).

**Facilities:** ALMA

**Software:** ASTROPy (Astropy Collaboration 2013, 2018, 2022), CORNER (Foreman-Mackey 2016), EMCEE (Foreman-Mackey et al. 2013), GALARIO (Tazzari et al. 2018), MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020).

## DATA AVAILABILITY

The data presented in this work are available upon reasonable request to the authors.

## REFERENCES

* Allard F., Hauschildt P. H., Alexander D. R., Tamanai A., Schweitzer A., 2001, ApJ, 556, 357

* Allard F., Homeier D., Freytag B., 2012, Phil. Trans. R. Soc. A, 370, 2765

* Astropy Collaboration, 2013, A&A, 558, A33

* Astropy Collaboration, 2018, AJ, 156, 123

* Astropy Collaboration, 2022, ApJ, 935, 167

* Augereau J. C., Beust H., 2006, A&A, 455, 987

* Babusiaux C. et al., 2023, A&A, 674, A32

* Beichman C. A., Neugebauer G., Habing H. J., Clegg P. E., Chester T. J., eds, 1988, Infrared Astronomical Satellite (IRAS) Catalogs and Atlases. Volume 1: Explanatory Supplement. 1

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

* Bohren C. F., Huffman D. R., 1998, Absorption and Scattering of Light by Small Particles. John Wiley and Sons
Booth M. et al. 2024, in Zmuidzinas J., Gao J.-R., eds, SPIE Conf. Ser. Vol. 13102, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy XII. SPIE, Bellingham, p. 1310206

Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127

Burnham K. P., Anderson D. R., 2004, Soc. Methods Res., 33, 261

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

Chen C. H., Mittal T., Kuchner M., Forrest W. J., Lisse C. M., Manoj P., Sargent B. A., Watson D. M., 2014a, ApJS, 211, 25

Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, MNRAS, 452, 1068

Chen Y., Girardi L., Bressan A., Marigo P., Barbieri M., Kong X., 2014b, MNRAS, 444, 2525

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

Choquet É. et al., 2018, ApJ, 854, 53

Cronin-Coltsmann P. F. et al., 2021, MNRAS, 504, 4497

Cronin-Coltsmann P. F., Kennedy G. M., Kral Q., Lestrade J.-F., Marino S., Matrà L., Wyatt M. C., 2023, MNRAS, 526, 5401

Crotts K. A. et al., 2024, ApJ, 961, 245

Dahlqvist C. H. et al., 2022, A&A, 666, A33

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

David T. J., Hillenbrand L. A., 2015, ApJ, 804, 146

del Burgo C., Allende Prieto C., 2016, MNRAS, 463, 1400

del Burgo C., Allende Prieto C., 2018, MNRAS, 479, 1953

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

Draine B. T., 2003, ApJ, 598, 1017

Dullemond C. P., Juhasz A., Pohl A., Sereshti F., Shetty R., Peters T., Commercon B., Flock M., 2012, Astrophysics Source Code Library, preprint(arXiv:1202.015)

Eiroa C. et al., 2013, A&A, 555, A11

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

Faramaz V. et al., 2021, AJ, 161, 271

Foreman-Mackey D., 2016, J. Open Source Softw., 24

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Fouesneau M., 2022, 'pyphot', https://github.com/mfouesneau/pyphot, (accessed 2024 August 8)

Gagné J. et al., 2018, ApJ, 856, 23

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

Gaia Collaboration, 2023, A&A, 674, A1

Harris C. R. et al., 2020, Nature, 585, 357

Heras A. M., Eiroa C., del Burgo C., Marshall J. P., Montesinos B., 2025, A&A, 694, A325

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

Holland W. S. et al., 2017, MNRAS, 470, 3606

Horner J. et al., 2020, PASP, 132, 102001

Hudgins D. M., Sandford S. A., Allamandola L. J., Tielens A. G. G. M., 1993, ApJS, 86, 713

Hughes A. M., Duchêne G., Matthews B. C., 2018, ARA&A, 56, 541

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

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

Kains N., Wyatt M. C., Greaves J. S., 2011, MNRAS, 414, 2486

Kavelaars J. J. et al., 2009, AJ, 137, 4917

Kennedy G. M. et al., 2018, MNRAS, 476, 4584

Kennedy G. M., Wyatt M. C., 2013, MNRAS, 433, 2334

Klaassen P. et al., 2024, Open Research Europe, 4, 112

Klačka J., Kocifaj M., Pástor P., Petržala J., 2007, A&A, 464, 127

Kossacki K. J., 2021, Icarus, 368, 114613

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

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

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

Krivov A. V., Booth M., 2018, MNRAS, 479, 3300

Krivov A. V., Müller S., Löhne T., Mutschke H., 2008, ApJ, 687, 608

Krivov A. V., Wyatt M. C., 2021, MNRAS, 500, 718

Kuchner M. J., Stark C. C., 2010, AJ, 140, 1007

Lebouteiller V., Barry D. J., Spoon H. W. W., Bernard-Salas J., Sloan G. C.,

Houck J. R., Weedman D. W., 2011, ApJS, 196, 8

Lestrade J. F., Wyatt M. C., Bertoldi F., Menten K. M., Labaigt G., 2009, A&A, 506, 1455

Liu Q., Wang T., Jiang P., 2014, AJ, 148, 3

Luppe P., Krivov A. V., Booth M., Lestrade J.-F., 2020, MNRAS, 499, 3932

MacGregor M. A. et al., 2016, ApJ, 823, 79

Marino S., 2021, MNRAS, 503, 5100

Marshall J. P. et al., 2014, A&A, 570, A114

Marshall J. P. et al., 2023a, MNRAS, 521, 5940

Marshall J. P., Cotton D. V., Bott K., Bailey J., Kedziora-Chudczer L., Brown E. L., 2023b, MNRAS, 522, 2777

Marshall J. P., Maddison S. T., Thilliez E., Matthews B. C., Wilner D. J., Greaves J. S., Holland W. S., 2017, MNRAS, 468, 2719

Marshall J. P., Milli J., Choquet É., del Burgo C., Kennedy G. M., Matrà L., Ertel S., Boccaletti A., 2018, ApJ, 869, 10

Marshall J. P., Wang L., Kennedy G. M., Zeegers S. T., Scicluna P., 2021, MNRAS, 501, 6168

Matrà L. et al., 2025, A&A, 693, A151

Matthews E. et al., 2018, MNRAS, 480, 2757

Montesinos B. et al., 2016, A&A, 593, A51

Morales F. Y. et al., 2009, ApJ, 699, 1067

Morales F. Y., Bryden G., Werner M. W., Stapelfeldt K. R., 2016, ApJ, 831, 97

Mulders G. D., Min M., Dominik C., Debes J. H., Schneider G., 2013, A&A, 549, A112

Müller S., Löhne T., Krivov A. V., 2010, ApJ, 708, 1728

Muñoz-Gutiérrez M. A., Marshall J. P., Peimbert A., 2023, MNRAS, 520, 3218

Mustill A. J., Wyatt M. C., 2009, MNRAS, 399, 1403

Najita J. R., Kenyon S. J., Bromley B. C., 2022, ApJ, 925, 45

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

Pan M., Schlichting H. E., 2012, ApJ, 747, 113

Patten B. M., Willson L. A., 1991, AJ, 102, 323

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

Pearce T. D. et al., 2022, A&A, 659, A135

Petit J. M. et al., 2011, AJ, 142, 131

Poynting J. H., 1903, MNRAS, 64, 1

Preibisch T., Ossenkopf V., Yorke H. W., Henning T., 1993, A&A, 279, 577

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

Rhee J. H., Song I., Zuckerman B., McElwain M., 2007, ApJ, 660, 1556

Riello M. et al., 2021, A&A, 649, A3

Robertson H. P., 1937, MNRAS, 97, 423

Schwarz G., 1978, Ann. Stat., 6, 461

Sefilian Antranik A., 2024. The Astrophysical Journal, 966, 140

Sibthorpe C., Kennedy G. M., Wyatt M. C., Lestrade J. F., Greaves J. S., Matthews B. C., Duchêne G., 2018, MNRAS, 475, 3046

Skrutskie M. F. et al., 2006, AJ, 131, 1163

Tang J., Bressan A., Rosenfield P., Slemer A., Marigo P., Girardi L., Bianchi L., 2014, MNRAS, 445, 4287

Tazzari M., Beaujean F., Testi L., 2018, MNRAS, 476, 4527

Thureau N. D. et al., 2014, MNRAS, 445, 2558

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Vitense C., Krivov A. V., Kobayashi H., Löhne T., 2012, A&A, 540, A30

Wahhaj Z. et al., 2016, A&A, 596, L4

Warren S. G., Brandt R. E., 2008, J. Geophys. Res. (Atmospheres), 113, D14220

Wenger M. et al., 2000, A&AS, 143, 9

Wolf S., Hillenbrand L. A., 2003, ApJ, 596, 603

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

Wright G. S. et al., 2023, PASP, 135, 048003

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

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

Zhang X., Green G. M., Rix H.-W., 2023, MNRAS, 524, 1855

Zuckerman B., 2019, ApJ, 870, 27

Zuckerman B., Song I., 2004, ARA&A, 42, 685
# APPENDIX A: DUSTY GRAIN MODEL

A quasi-static dynamics and radiative transfer code for dust grains in debris discs ('STARDUST') was developed to:

(i) predict the radial location of the dust grains depending on the forces acting on the grains due to the host star according to size (dynamics),

(ii) predict the thermal emission of the dust grains from the grain's location, size, and composition (radiative transfer).

The grains in the Stardust model are assumed to be spherical objects that absorb and emit radiation according to Mie theory (Bohren & Huffman 1998). In addition to gravity, these grains are subject to different forces that can move them far away from the point of production (e.g. Krivov et al. 2008). The model can switch 'on' and 'off' the force of stellar radiation pressure (e.g. Burns et al. 1979) and the inward drag due to the Poynting-Robertson (P-R) effect (e.g. Poynting 1903; Robertson 1937; Klačka et al. 2007) in comparison to the mean collisional time of the grain (e.g. Kuchner & Stark 2010).

The model assumes an optically thin disc with no gas present. When the grains are traversing along elliptical orbits due to radiation pressure, it is assumed that no collisions occur. However, when the model activates the P-R effect, i.e. when the grain loses angular momentum and spirals inwards towards the star, a simple collision assumption is implemented based on the optical depth of the disc. Additional force components, such as stellar wind, will be considered with future iterations of this model.

## A1 Radiation pressure affected grains

Grains in the model are initially on Keplerian orbits, where grains are under the influence of gravity. When the conservative radiation pressure force is switched on, this pushes grains out along elliptical orbits according to grain size. To determine the resultant orbit is by the ratio ($\beta$) of radiation pressure to gravitational forces (e.g. Krivov 2010):

$$ \beta = 0.574 \frac{L_{\star} Q_{\text{pr}}(s)}{M_{\star} \rho s}, \qquad (A1) $$

where, $L_{\star}$ and $M_{\star}$ are the solar luminosity and mass of the host star, grain density ($\rho$), and size ($s$) are expressed in units of 'g cm$^{-3}$' and 'µm' respectively, and the dimensionless quantity $Q_{\text{pr}}(s)$ is the radiation pressure efficiency dependent on grain size. The $\beta$ values in the model are calculated for each corresponding grain size.

The eccentricity ($e$) of the orbit can be derived from the $\beta$ (e.g. Augereau & Beust 2006):

$$ e = \frac{\beta}{1 - \beta} \qquad (A2) $$

and the semimajor distance ($a$) of the ellipse, with pericentre ($r_p$), calculated using:

$$ a = \frac{r_p}{1 - e}. \qquad (A3) $$

Our aim is to determine the probability of a grain at a particular point in time along the orbit. If we simply plot stellar distance as a function of True Anomaly ($\theta$), time is not considered (see left plot of Fig. A1 for a grain corresponding to $\beta = 0.3$). Therefore, to find a better representation we start with Kepler's Second Law in the following form:

$$ r^2 d\theta = mabdt \qquad (A4) $$

where $m = 2\pi/P$ is referred to as the mean motion, and $P$ is the orbital period in years for a given semi-major axis in au and host star mass in Solar masses ($P = a^{3/2}M_{\star}^{-1/2}$). The semimajor axis ($a$) and the semiminor axis ($b$) are in units of au, and the period ($P$) is in years. Eq. (A4) cannot be solved analytically, thus, we have adopted the following iterative case:

$$ d\theta = \Delta\theta = \theta_{n+1} - \theta_n \qquad (A5) $$

and the time-step defined by:

$$ dt = \Delta t (= t_{n+1} - t_n) \qquad (A6) $$

So, substituting eqs (A5) and (A6) into eq. (A4), yields the following relationship:

$$ \theta_{n+1} = \theta_n + \frac{abm}{r_n^2} \Delta t \qquad (A7) $$

eq. (A7) allows iterative solution to the stellar distance of the grain as a function of both True Anomaly and time. An appropriate time-step has been chosen to correspond to the $\beta$-value of the largest grain size. See middle plot of Fig. A1, for the grain orbit corresponding to $\beta = 0.3$. The next step is to group the radial data determined by solving eq. (A7) with a bin size of 1 au (although in principle this can be any size desired), resulting in the relative probability distribution corresponding to $\beta$ in the right-hand plot of Fig. A1. Inspecting this figure shows that the grain will spend most of its time on closest and farthest approaches of the star.

<figure><img src="image_6.png" /><figcaption>**Figure A1.** Grain corresponding to $\beta_{\text{pr}} = 0.3$ being released at a stellar distance of 30 au. Left panel: stellar distance (au) of the grain as a function of True Anomaly (radians). Middle panel: stellar distance (au) of the grain as a function of time (years). Right panel: probability of finding the grain along its orbit as a function of stellar distance (au) binned in sizes of 1 au.</figcaption></figure>
This may be counter-intuitive for the nearest approach as the grain has the fastest speed. However, these extreme points (pericentre and apocentre) are turning points for the grain and the corresponding stellar distances (r) at these locales are more frequent, as opposed to locales where the grain is travelling towards and away from these extremities.

## A2 Grain size distribution

The distribution of grain sizes typically favour smaller grains in debris discs. Dohnanyi (1969) presented a collisional model that showed the number density of particles per unit volume (dn/dV) in the mass range m to m + dm is:

$$ \frac{dn}{dV} \propto m^{-\alpha} dm \qquad (A8) $$

Where the exponent, $\alpha$, determines the rate of collisions. For example, for $\alpha$ equal to 11/6 corresponds to a steady-state collisional cascade. However, in the literature (Dohnanyi 1969), the steady-state collisional cascade model is often cited with the transform that shows the relative number of grains according to size:

$$ dN \propto s^{-q} ds \qquad (A9) $$

Where N is the number of grains corresponding to size s, and q is the exponent that determines the type of collision. In this transform, a steady-state collisional cascade occurs when q = 3.5. Generally, for debris discs, this exponent lies between 3.0 and 4.0 (Hughes et al. 2018).

## A3 P-R effect and collisions

The model can include grains that undergo orbital decay due to loss of energy from the incident electromagnetic radiation on a moving non-rotating spherical grain, i.e. the P-R Effect. The time taken for grain to spiral in from outer orbit of semimajor axis ($a_i$) to an inner orbit of reduced semimajor axis ($a_f$) due to the P-R Effect (Wyatt):

$$ t_{pr} = 400 M_{*}^{-1} (a_{i}^{2} - a_{f}^{2}) / \beta, \qquad (A10) $$

For simplicity, the P-R time is compared to the (estimated) mean collisional time ($t_{coll}$, Kuchner & Stark 2010):

$$ t_{coll} = \sqrt{(a_i^3/M)/4\pi\tau}, \qquad (A11) $$

Where $\tau$ is the effective optical depth of planetesimal belt.

## A4 Radiative transfer of grains

We have described how the model defines grain location, size, and composition including belt structure. We can now estimate the temperature and determine the direct thermal emission of the grains.

If thermal equilibrium is assumed, i.e. the energy that the grains absorbs is equal to its emission, then the grain temperature ($T_g$) at a specific stellar distance (r) is derived from the following expression (Wolf & Hillenbrand 2003):

$$ r = \frac{R_{*}}{2} \left[ \frac{\int_{0}^{\infty} Q_{\lambda}^{\text{abs}}(s) F_{\lambda}(T_{*}) d\lambda}{\int_{0}^{\infty} Q_{\lambda}^{\text{abs}}(s) B_{\lambda}(T_{g}) d\lambda} \right]^{1/2} \qquad (A12) $$

where $R_*$ and $T_*$ are star's radius and effective temperature respectively. $F_\lambda(T)$ is the photosphere stellar model for star, $B_\lambda(T_g)$ is the blackbody function for a given grain temperature ($T_g$), and $Q_\lambda^{\text{abs}}(s)$ is the absorption efficiency calculated for a grain radius (s) as a function of wavelength.

To determine the direct thermal emission, we will assume grains of a particular size (s) are within a rotationally symmetric disc heated by the host star observed from a distance D. These grains are within an annulus with mean stellar distance (r) with a thickness of one bin size. The direct thermal emission, measured in Janskys, is derived from a modified blackbody expression:

$$ F_{\lambda}^{\text{th}}(r, s) = 4\pi^2 N(r, s) \left(\frac{s}{D}\right)^2 Q_{\lambda}^{\text{abs}}(s) B_{\lambda}(T_g), \qquad (A13) $$

where the grain number density ($N(r, s)$) is defined as:

$$ N(r, s) = 2\pi n s r \qquad (A14) $$

corresponding to the number of grains (n) of a specific grain size (s) that are present within an annulus with a mean radius defined by the stellar distance r. It is noted that eq. (A13) is a modified form of the equation presented in both Krivov et al. (2008) and Müller, Löhne & Krivov (2010) where they calculate the flux for all grains at all distances across the entire disc.

The optical properties of the IMPs were determined using eq. (8.50) of Bohren & Huffman (1998) where f = 0.5, 0.75 and 0.9, respectively; note this method was used in Morales et al. (2016). These types of composition is commonly referred to as 'dirty ice'.

# APPENDIX B: POSTERIOR PROBABILITY DISTRIBUTIONS OF THE MODELLING

Here we present a corner plot (Fig. B1) showing the posterior probability distributions from the EMCEE runs (100 000 realizations of the model) and the values for the maximum likelihood model parameters presented in the text (Table 2, Section 3). We omit the posterior for the scale height parameter (h) from the corner plot as it is flat and uninformative. We also present the corner plot (Fig. B2) corresponding to the probability distribution to help pin down inner belt characteristics where the outer belt is composed of pure astronomical silicate.
<figure><img src="image_7.png" /><figcaption>**Figure B1.** Corner plot showing the posterior probability distributions for the disc architecture model parameters. The posterior for the scale height, *h*, has been omitted as it is flat and uninformative.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure B2.** Corner plot showing the posterior probability distributions from the SED model analysis. The 'cross-hairs' are the Maximum Likelihood Values from the analysis.</figcaption></figure>

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