# The polarisation properties of the HD 181327 debris ring Evidence for sub-micron particles from scattered light observations

J Milli, E Choquet, R Tazaki, F Ménard, J.-C Augereau, J Olofsson, P

Thébault, O Poch, Anny Chantal Levasseur-Regourd, J Lasue, et al.

► To cite this version:

J Milli, E Choquet, R Tazaki, F Ménard, J.-C Augereau, et al.. The polarisation properties of the HD 181327 debris ring Evidence for sub-micron particles from scattered light observations. Astronomy and Astrophysics - A&A, 2024, 683, pp.A22. 10.1051/0004-6361/202347933. insu-04333442v2

HAL Id: insu-04333442

https://insu.hal.science/insu-04333442v2

Submitted on 1 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 polarisation properties of the HD 181327 debris ring
## Evidence for sub-micron particles from scattered light observations*

J. Milli<fnref n="1" />, E. Choquet<fnref n="2" />, R. Tazaki<fnref n="1" />, F. Ménard<fnref n="1" />, J.-C. Augereau<fnref n="1" />, J. Olofsson<fnref n="4" />,<fnref n="5" />, P. Thébault<fnref n="6" />, O. Poch<fnref n="1" />,
A.-C. Levasseur-Regourd<fnref n="7" />,†, J. Lasue<fnref n="8" />, J. B. Renard<fnref n="9" />, E. Hadamcik<fnref n="10" />, C. Baruteau<fnref n="8" />, H. M. Schmid<fnref n="3" />, N. Engler<fnref n="3" />,
R. G. van Holstein<fnref n="11" />, E. Zubko<fnref n="12" />, A. M. Lagrange<fnref n="6" />,<fnref n="1" />, S. Marino<fnref n="13" />, C. Pinte<fnref n="1" />, C. Dominik<fnref n="14" />, A. Boccaletti<fnref n="6" />,
M. Langlois<fnref n="15" />, A. Zurlo<fnref n="16" />,<fnref n="17" />,<fnref n="18" />, C. Desgrange<fnref n="1" />,<fnref n="4" />, L. Gluck<fnref n="1" />, D. Mouillet<fnref n="1" />, A. Costille<fnref n="2" />, and J. F. Sauvage<fnref n="2" />,<fnref n="19" />

Received 11 September 2023 / Accepted 22 November 2023

ABSTRACT

**Context.** Polarisation is a powerful remote-sensing tool to study the nature of particles scattering the starlight. It is widely used to characterise interplanetary dust particles in the Solar System and increasingly employed to investigate extrasolar dust in debris discs' systems.

**Aims.** We aim to measure the scattering properties of the dust from the debris ring around HD 181327 at near-infrared wavelengths.

**Methods.** We obtained high-contrast polarimetric images of HD 181327 in the *H* band with the SPHERE/IRDIS instrument on the Very Large Telescope (ESO). We complemented them with archival data from HST/NICMOS in the F110W filter reprocessed in the context of the Archival Legacy Investigations of Circumstellar Environments (ALICE) project. We developed a combined forward-modelling framework to simultaneously retrieve the scattering phase function in polarisation and intensity.

**Results.** We detected the debris disc around HD 181327 in polarised light and total intensity. We measured the scattering phase function and the degree of linear polarisation of the dust at 1.6 µm in the birth ring. The maximum polarisation is 23.6% ± 2.6% and occurs between a scattering angle of 70° and 82°.

**Conclusions.** We show that compact spherical particles made of a highly refractive and relatively absorbing material in a differential power-law size distribution of exponent -3.5 can simultaneously reproduce the polarimetric and total intensity scattering properties of the dust. This type of material cannot be obtained with a mixture of silicates, amorphous carbon, water ice, and porosity, and requires a more refracting component such as iron-bearing minerals. We reveal a striking analogy between the near-infrared polarisation of comets and that of HD 181327. The methodology developed here combining VLT/SPHERE and HST/NICMOS may be applicable in the future to combine the polarimetric capabilities of SPHERE with the sensitivity of JWST.

**Key words.** polarization – scattering – instrumentation: high angular resolution – comets: general – Kuiper belt: general – planetary systems

<footnote n="1">¹ Université Grenoble-Alpes, CNRS, IPAG, 38000 Grenoble, France
e-mail: julien.milli@univ-grenoble-alpes.fr</footnote>

<footnote n="2">² Aix-Marseille Univ., CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="3">³ Institute for Particle Physics and Astrophysics, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland</footnote>

<footnote n="4">⁴ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="5">⁵ European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="6">⁶ LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris-Diderot, Sorbonne Paris-Cité,
5 place Jules Janssen, 92195 Meudon, France</footnote>

<footnote n="7">⁷ LATMOS, Sorbonne Université, CNRS, CNES, Paris, France</footnote>

<footnote n="8">⁸ IRAP, Université de Toulouse, CNES, CNRS, UPS, Toulouse, France</footnote>

<footnote n="9">⁹ LPC2E, Université d'Orléans, CNRS, Orléans, France</footnote>

<footnote n="10">¹⁰ LATMOS, CNRS, 11 bld d'Alembert, 78280 Guyancourt, France</footnote>

<footnote n="11">¹¹ European Southern Observatory (ESO), Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago, Chile</footnote>

<footnote n="12">¹² Planetary Atmospheres Group, Institute for Basic Science (IBS), 55 Expo-ro, Yuseong-gu, Daejeon 34126, South Korea</footnote>

<footnote n="13">¹³ Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter, EX4 4QL, UK</footnote>

<footnote n="14">¹⁴ Anton Pannekoek Institute for Astronomy, Science Park 904, 1098 XH Amsterdam, The Netherlands</footnote>

<footnote n="15">¹⁵ CRAL, UMR 5574, CNRS, Université de Lyon, Ecole Normale Supérieure de Lyon, 46 Allée d'Italie, 69364 Lyon Cedex 07,
France</footnote>

<footnote n="16">¹⁶ Millennium Nucleus on Young Exoplanets and their Moons (YEMS), Chile</footnote>

<footnote n="17">¹⁷ Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile</footnote>

<footnote n="18">¹⁸ Escuela de Ingeniería Industrial, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile</footnote>

<footnote n="19">¹⁹ DOTA, ONERA, 13661, Salon Cedex Air, France</footnote>

<footnote>* The reduced images (FITS files) presented in Fig. 2 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/683/A22</footnote>

<footnote>† Deceased.</footnote>
# 1. Introduction

Debris discs are massive analogues of the asteroids and Edgeworth-Kuiper belts in the Solar System (e.g. Hughes et al. 2018, for a review). They are made of a population of large kilometre-sized planetesimals, which produce smaller dust particles through mutual collisions. The dust’ thermal emission creates an infrared excess in the spectral energy distribution (SED) of the star detectable above the photospheric emission with space-based infrared telescopes. Large infrared surveys revealed that at least 15% to 30% of main-sequence stars host a debris disc (e.g. Matthews et al. 2014; Montesinos et al. 2016), and this detection rate can reach 75% for early-type stars in young moving groups (F stars in the β Pic moving group; Pawellek et al. 2021). This means that kilometre-sized planetary embryos are a common outcome of stellar formation. Due to limitations in sensitivity, the detected debris discs have a dust mass several orders of magnitude higher than that in our Solar System at the present time (Krivov & Wyatt 2021), and they mostly belong to younger systems of a few tens to a few hundreds of millions of years. As collisional activity decays with age (Wyatt 2006) and might peak after planet formation, those discs are compatible with our view of the young Solar System after a few hundred million years when it was collisionally very active (Booth et al. 2009). Dust is therefore an important ingredient in these planetary systems and represents a valuable observable that can shed light on the architecture of the underlying planetary system and on the nature of the building blocks of planets to constrain the system formation. Major advances in high-angular resolution and high-contrast imaging techniques now allow one to angularly resolve and isolate the discs from the glare of the central star, opening up new perspectives to characterise the properties of the dust particles through their scattered starlight. About fifty-five debris discs have been detected in scattered light so far, but only a dozen of them have a high enough surface brightness and a favourable geometry or inclination allowing one to extract the scattering properties such as the phase function dependence with the scattering angle. This is the case for Fomalhaut (Min et al. 2010), HD 181327 (Stark et al. 2014), HD 61005 (Olofsson et al. 2016), HR 4796 (Perrin et al. 2015; Milli et al. 2017, 2019; Olofsson et al. 2020; Chen et al. 2020; Arriaga et al. 2020), HD 35841 (Esposito et al. 2018), HD 191089 (Ren et al. 2019) or TWA 7 (Ren et al. 2021). They typically exhibit a prominent peak of forward-scattering (HD 61005, HR 4796, HD 35841, or HD 191089) with a width that varies significantly from one system to the other, and a mild back-scattering, starting either at relatively small scattering angles (from 50° for HR 4796 and maybe even less for Fomalhaut), or at angles larger than 90° (HD 191089 and HD 181327).

Even fewer debris discs have been analysed in polarised light, jointly to total intensity, to extract the degree of linear polarisation as a function of the scattering angle (see Table 1), mostly because obtaining high-fidelity and flux-calibrated images simultaneously in total intensity and linearly polarised light is difficult and possible so far only for the brightest systems. This is the case for the debris disc detected around the star HD 181327, a young (18.5 Myr; Miret-Roig et al. 2020) F5/6V star member of the β Pictoris moving group and located at a distance of 48.2 ± 0.2 pc (Gaia Collaboration 2018). Its infrared excess $L_{IR,disc}/L_*$ is estimated at 0.2% (Lebreton et al. 2012). Hubble Space Telescope (HST) images in the near-infrared from NICMOS (1.1 µm; Schneider et al. 2006) and in the optical from ACS (0.6 µm) and later from STIS (0.4–0.8 µm; Stark et al. 2014) revealed a wide ring inclined by 28.5°+2.1°−2.0°, with an inner edge at 76.6 ± 1.0 au and a maximum brightness at 84.2 ± 1.0 au

<table><caption>**Table 1.** Discs with the degree of linear polarisation extracted over a range of scattering angles.</caption><thead><tr><th>Host star</th><th>Range of scat. angles</th><th>Band</th><th>Reference</th></tr></thead><tbody><tr><td>HR 4796</td><td>20°–120°</td><td>K</td><td>Arriaga et al. (2020)</td></tr><tr><td>HD 35841</td><td>22°–125°</td><td>H</td><td>Esposito et al. (2018)</td></tr><tr><td>HD 191089</td><td>30°–130°</td><td>H</td><td>Ren et al. (2019)</td></tr><tr><td>HD 114082</td><td>25°–145°</td><td>H</td><td>Engler et al. (2023)</td></tr><tr><td>HD 181327</td><td>60°–120°</td><td>H</td><td>This work</td></tr></tbody></table>

(deprojected semi-major axis values from Stark et al. 2014, taking into account the revised distance to the star). The ring is about 25 au wide, with a cleared interior and an extended nebulosity detected in the optical up to ~400 au (Schneider et al. 2006). Stark et al. (2014) extracted the scattering phase function (hereafter SPF), showing changes as a function of the distance to the star, compatible with dust segregation in the system. The width of the forward-scattering peak increases with the distance to the star, which suggests that smaller particles dominate the scattered light at larger distances. This is an expected behaviour for a collisionally active ring beyond which should lie a halo of small particles, produced in the main ring and placed on high-eccentricity orbits by radiation pressure (Lecavelier Des Etangs et al. 1996; Strubbe & Chiang 2006; Thébault & Wu 2008). This hypothesis is further supported by millimetre observations with ALMA at 1.3 mm (Marino et al. 2016), showing a ring with a semi-major axis at maximum dust density that is 4.2 ± 1.1 au smaller compared to optical observations (after correcting for the latest distance estimate).

A detailed modelling of the SED of the dust constrained by the resolved images of the disc was performed in Lebreton et al. (2012). It suggests the disc contains micron-sized particles (minimum size 0.9 µm), made of porous amorphous silicates (~12% in volume) and carbonaceous material (~23%) surrounded by an important layer of ice (~65%) in a ~65% porous structure, possibly fluffy.

Schneider et al. (2006) and Stark et al. (2014) however revealed some tensions between the minimum grain size suggested by SED modelling of about 1 µm and the scattered light behaviour which favours smaller, sub-micron particles. In an attempt to explore this inconsistency, we present here the first polarimetric observations of this system and extract the degree of linear polarisation. We describe our observations in Sect. 2, analyse the morphology in Sect. 3.1, extract the phase function and degree of linear polarisation in Sect. 4 and discuss the scattering properties of the ring in Sect. 5 before concluding in Sect. 6.

# 2. Observations

## 2.1. VLT/SPHERE instrumental setup

The star HD 181327 was observed on the night of 15 May, 2017, with the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE; Beuzit et al. 2019), as part of the Guaranteed Time Observations of the instrument consortium<fnref n="1" />. SPHERE is a high-contrast imager fed by an extreme adaptive optics (AO) system (Sauvage et al. 2016) to correct for the atmospheric turbulence and static aberrations. We used the near-infrared arm of SPHERE in dual-beam polarimetric imaging (DPI; Langlois et al. 2010; van Holstein et al. 2020a; de Boer et al. 2020). The

<footnote n="1">¹ ESO programme 099.C-0147(B).</footnote>
<table><caption>**Table 2.** Log of the SPHERE observations of HD 181327 and the reference star HD 202917 obtained on the UT date of 16 May, 2017.</caption><thead><tr><th>UT time (HH:MM)</th><th>Target</th><th>DIT (s) × NDIT<sup>(a)</sup></th><th># HWP positions</th><th># HWP cycles</th><th>Seeing<sup>(b)</sup> (′′)</th><th>τ<sub>0</sub><sup>(c)</sup> (ms)</th><th>Wind<sup>(d)</sup> (m s<sup>-1</sup>)</th><th>Strehl<sup>(e)</sup> (%)</th><th>Par. angle<sup>(f)</sup> start/end (°)</th></tr></thead><tbody><tr><td>05:59–06:18</td><td>HD 181327</td><td>32 × 1</td><td>2</td><td>15</td><td>1.3</td><td>1.8</td><td>12.3</td><td>71</td><td>-58.1/-52.1</td></tr><tr><td>06:23–07:34</td><td>HD 181327</td><td>32 × 2</td><td>4</td><td>15</td><td>1.2</td><td>2.4</td><td>12.5</td><td>71</td><td>-50.6/-23.7</td></tr><tr><td>07:56–09:34</td><td>HD 202917</td><td>64 × 2</td><td>4</td><td>11</td><td>1.2</td><td>1.9</td><td>11.2</td><td>58</td><td>-59.5/-23.0</td></tr></tbody></table>

**Notes.** (a)DIT is the individual Detector Integration Time in seconds and NDIT is the number of DITs. (b)From the DIMM. (c)Coherence time measured by the MASS-DIMM. (d)Measured 30 m above the platform by the Automated Weather Station part of the Astronomical Site Monitoring. (e)Estimate from the Real Time Computer of SPHERE. (f)Parallactic angle.

IRDIS imager splits the incoming light in two channels, and in DPI, a set of polarisers with orthogonal transmission axes is inserted in the dual-filter wheel in order to simultaneously image the light linearly polarised in two orthogonal directions. IRDIS provides a 11"×11" field of view with a pixel scale of 12.25 mas (Maire et al. 2016). To benefit from a high-Strehl atmospheric correction in the near infrared, we used the broad-band *H* filter (λ = 1.625 µm, Δλ = 0.29 µm). Observations were carried out in field-stabilised mode, because at the time of observations, the pupil-stabilised mode for DPI was not yet offered (van Holstein et al. 2017). They used the apodised Lyot coronagraph of radius 92.5 milliarcsecond (mas) to block the stellar light.

The details of the coronagraphic data obtained are summarised in Table 2 along with the atmospheric conditions. One polarimetric cycle is made of images recorded at four half-wave plate (HWP) switch angles to measure the Stokes parameters *Q*⁺, *Q*⁻, *U*⁺, and *U*⁻. It is referred to as a HWP cycle hereafter and in Table 2. Unfortunately for the first 15 HWP cycles out of 41 in total, only two positions of the HWP for *Q*⁺ and *Q*⁻ were recorded.

Because the observations were field-stabilised, Angular Differential Imaging (ADI; Marois et al. 2006) is not applicable to remove the glare of the central star and reveal the disc in total intensity. In addition, ADI is not applicable to objects azimuthally extended such as discs seen under a low inclination, because of severe self-subtraction of the astrophysical signal (Milli et al. 2012). We used an alternative approach by using a reference star observed in the same mode (same filter and coronagraph) and with similar AO correction. This type of data reduction is referred to as Reference Differential Imaging (hereafter RDI; e.g. Ruane et al. 2019). It is commonly applied for space-based imaging (for instance on HST/NICMOS with the ALICE programme; Soumm er et al. 2014; Choquet et al. 2014; Hagan et al. 2018), where the point-spread function (hereafter PSF) is very stable. From the ground, it was also one of the first PSF subtraction technique applied to AO-assisted observations (Mouillet et al. 1997) but was later superseded by ADI once pupil-stabilisation was recognised as the most efficient observing strategy to calibrate the halo of speckles of the central star. In field-stabilised observations, the spiders holding the secondary mirror of the telescope are not aligned with the Lyot stop, therefore they diffract the stellar light and relevant reference frames should match the parallactic angles of the telescope in order to have a similar orientation of the spider diffraction pattern. Among all the stars observed in the night of 15 May, 2017, HD 202917 observed immediately after HD 181327 provided the best result when used as a reference star to subtract the PSF for the sequence of images of HD 181327. HD 202917 was observed along the same range of parallactic angles, under similar atmospheric conditions (see Table 2), and is only slightly fainter

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Raw coronagraphic frames of HD 181327 and HD 202917. Left: first (top) and last (bottom) frames of the science target HD 181327. Right: two examples of frames of the reference star HD 202917, observed immediately after HD 181327, and with parallactic angles matching that of the first (top) and last (bottom) raw frame of HD 181327. North is up and east is to the left. The colour scale is logarithmic and the flux is shown in the raw detector counts called analogue to digital units (ADU).</figcaption></figure>

(G magnitude of 8.5 vs. 6.9 for HD 181327). An example showing two raw frames of the science star and two raw frames of the reference star is shown in Fig. 1. This reference star is known to host a faint debris disc, with a low fractional luminosity of 2.5 × 10⁻⁴, only detected in scattered light from space (Schneider et al. 2016; Soumm er et al. 2014) with a very faint peak surface brightness of 0.2 mJy arcsec⁻² along the semi-minor axis and a semi-major axis of 1.46". It is undetectable in our images neither in the raw frames nor in the reduced polarised or total intensity data. It therefore does not impact our choice of this star as a reference star.

## 2.2. Polarimetric data reduction

For each HWP cycle, we obtained the Stokes images *I*, *Q*, and *U* and corrected for the instrumental polarization (IP)
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Final images of the intensity (Stokes *I*, left image), azimuthal Stokes *Q*<sub>φ</sub> (middle) and *U*<sub>φ</sub> (right), calibrated in mJy arcsec<sup>-2</sup> (linear colour scale). North is up and east to the left. The features in the south-west at ~2” separation are a cluster of bad pixels on the IRDIS detector. A numerical mask of inner radius 1.25” and outer radius and 2.39” was applied on the intensity image.</figcaption></figure>

effects of the complete optical system using the IRDAP pipeline (van Holstein et al. 2020a,b). The polarisation of the central stellar halo was estimated to 0.08% ± 0.07% at the 1σ level by the pipeline. Therefore the star can be considered to be unpolarised, but we prefer using the image subtracted from this tiny stellar polarisation to have a slightly enhanced data quality.

As we expect the disc polarisation signal to be purely tangential or radial in the case of single-scattering by an optically thin disc illuminated by a single central illumination source, we use the azimuthal Stokes parameter $Q_{\phi}$ and $U_{\phi}$ defined as $Q_{\phi} = -Q \cos(2\phi) - U \sin(2\phi)$ and $U_{\phi} = +Q \sin(2\phi) - U \cos(2\phi)$ (de Boer et al. 2020), where $\phi$ is the polar angle between the north and the point of interest, measured from the north over east (the position angle). $Q_{\phi} > 0$ is equivalent to an azimuthal polarisation component while $Q_{\phi} < 0$ indicates a radial polarisation. The component $U_{\phi}$ describes the polarisation in the directions ±45° with respect to the radial direction.

The Stokes $Q_{\phi}$ and $U_{\phi}$ are shown in Fig. 2 (middle and right image), after conversion to mJy arcsec$^{-2}$. For the conversion, we estimated the star flux as the total flux of the mean off-axis PSF, encircled in a circular aperture of radius 170 px (2.1”), and we took into account the transmission of the neutral density filter and the difference in DIT between the off-axis PSF and the deep coronagraphic images. We assumed a stellar flux density of 4.22 ± 0.17 Jy (from 2MASS Johnson H filter; Cutri et al. 2003) and a pixel surface area of 12.25 × 12.25 mas$^2$.

For the polarimetric data reduction, we did not use the incomplete HWP cycles (first 15 cycles), because the improvement in signal-to-noise (S/N) was only marginal and we were concerned that combining together *Q* and *U* images with different S/N levels might bias the surface brightness distribution of the disc. In total intensity, this issue is irrelevant and all the 150 frames could be combined to form a master cube.

## 2.3. Total intensity data reduction

As explained in Sect. 2.1, we applied RDI to reveal the total intensity of the disc. The raw frames of the science target HD 181327 and its reference star HD 202917 were flatfielded, sky-subtracted, bad-pixel corrected and re-centred using the pipeline provided by the SPHERE Data Centre

(Delorme et al. 2017). This generated a science data cube consisting of 150 frames for HD 181327. We further selected the best frames, with no strong speckles and diffraction features at the expected location of the disc, which reduced the data cube to 98 frames. We looked for good references with matching parallactic angles in the data cube of HD 202917 and selected a reference cube consisting of 88 frames. We divided each of the science and reference cubes in four smaller data cubes sharing a similar range of parallactic angles. We applied for each set of science and reference cubes an RDI algorithm based on Principal Component Analysis (PCA; Soummer et al. 2012; Amara & Quanz 2012), to remove the quasi-static pattern of the PSF. The residual science images were then stacked together to obtain the image shown in Fig. 2 (left). For maximum efficiency, the PCA was optimised in a ring between 1.25” and 2.39”. The disc is clearly detected at the same location as in the polarised light image (Fig. 2, middle), although at a lower S/N because of the difficulty to remove efficiently the stellar halo in RDI compared to PDI. The inner part of the image is not shown because of strong contamination by unsubtracted stellar residuals. The features that extend radially at position angles 0°, 180° and 270° are instrumental artefacts not well corrected by the RDI post-processing. They correspond to phase aberrations induced by the pitch of the actuators of the SPHERE deformable mirror at specific spatial frequencies (40 cycles/pupil corresponding to a separation of 1.6” in the *H* band). On the other hand, the spiders diffraction pattern was well subtracted.

## 2.4. HST/NICMOS archival data

The HST/NICMOS were originally presented in Schneider et al. (2006). These observations were originally taken as part of the imaging survey GO-10177 (PI: Schneider), a search for debris discs around 26 targets with strong infrared excess. We reprocessed the deepest coronagraphic sequence obtained with the coronagraphic imaging mode of the NIC2 camera (0.’07565 pixel$^{-1}$, focal plane mask radius 0.’3), in the F110W filter ($\lambda$ = 1.104 µm, $\Delta\lambda$ = 0.5915 µm), obtained at two field orientations in a single spacecraft orbit on 2005 May 2 UT. These data were reduced and combined (de-rotated, stacked) using an advanced version of the pipeline developed for the
ALICE programme (PI: R. Soummer), a consistent reanalysis of the HST/NICMOS coronagraphic archives with advanced starlight subtraction methods (Choquet et al. 2014; Hagan et al. 2018), which allowed the discovery of 13 other new debris discs in scattered light (Soummer et al. 2014; Choquet et al. 2016, 2017, 2018; Marshall et al. 2018, 2023). The final image is shown later in Fig. 8 (left).

## 3. Analysis of the morphology and optical depth

The analysis of the morphology and radial profile of the dust density is required to constrain some geometrical properties of the disc and guide us for the extraction of the scattering properties. We used the polarimetric image, free from artefacts inherent to the use of star subtraction techniques in intensity, to investigate those two aspects.

### 3.1. Morphology

The disc detected around HD 181327 appears as an elliptical ring in the sky-projected plane. To derive the morphological parameters of this ring, we used the polarised intensity image $Q_\phi$ where the ring is detected at a S/N higher than in the total intensity image. We determined the radial location of the disk peak flux and fitted an ellipse to determine elliptical parameters of the dust ring. A similar technique was already used by Stark et al. (2014). As the disc is found in this study to be non-eccentric, we considered that the location of the belt peak surface brightness reflects the location of the underlying peak dust density and did not implement the iterative approach to correct for the $1/r^2$ illumination factor. We extracted radial profiles passing through the star and crossing the ring every 2° in position angle. For each profile, we determined the radial location of the maximum brightness of the disc by fitting a two-component power law (Eq. (1) in Milli et al. 2017). To find the best ellipse passing through the maximum brightness locations, we implemented the non-linear geometric fitting approach described in Ray & Srivastava (2008) within a Markov chain Monte Carlo (MCMC) framework (Foreman-Mackey et al. 2013) using uniform priors. The result is shown in Fig. 3. The ring has a semi-major axis $a = 1709 \pm 10$ mas oriented at a position angle PA of 99.1° ± 1.6°, a semi-minor axis $b = 1478 \pm 9$ mas, and its centre is offset by $x_0 = -6.8 \pm 8.2$ mas and $y_0 = 0.5 \pm 7.4$ mas in right ascension and declination ($x_0 < 0$ means an east offset, $y_0 > 0$ means a north offset). Using the Kowalsky deprojection technique described in Smart (1930) for binary systems and also applied by Stark et al. (2014); Rodigas et al. (2015); Milli et al. (2017, 2019) on debris discs, we derived the parameters of the true ellipse described by the dust particles in the orbital plane: the true semi-major axis $a$, the eccentricity $e$, the inclination $i$, the argument of pericentre $\omega$ and the longitude of ascending node $\Omega$. The result is given in Table 3. The derived PA and inclination are in agreement with the geometry derived from ALMA (Marino et al. 2016; Pawellek et al. 2021). They are also compatible within $1\sigma$ to the parameters derived from the STIS optical images of the ring, except for the semi-major axis, which is found to be ~2% larger in the optical. This small discrepancy could be real: optical and near-infrared images probe different dust populations and Stark et al. (2014) confirmed a spatial segregation of dust particles with smaller particles detected at large separations and larger particles within the birth ring. Our analysis also constrains the eccentricity to be below 1.1%, hence a poorly constrained argument of pericentre $\omega$.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Marginal probability distribution of the parameters of the best ellipse fitting ring: the E-W and N-S offset the ellipse centre $x_0$ and $y_0$, the semi-major and semi-minor axis $a$ and $b$, and the position angle PA. The graph in the upper right inset shows the data points used as input for the fit in blue, and the best ellipse in red.</figcaption></figure>

<table><caption>**Table 3.** Deprojected ellipse parameters.</caption><thead><tr><th>Parameter</th><th>IRDIS</th><th>STIS<sup>(a)</sup></th></tr></thead><tbody><tr><td><i>a</i> (au)</td><td>82.4 ± 0.5</td><td>84.2 ± 1.0<sup>(b)</sup></td></tr><tr><td><i>e</i></td><td>0.004<sup>+0.007</sup><sub>-0.001</sub></td><td>0.02 ± 0.01</td></tr><tr><td><i>i</i> (°)</td><td>30.2 ± 1.0</td><td>28.5<sup>+2.1</sup><sub>-2.0</sub></td></tr><tr><td><i>ω</i>(°)</td><td>-15<sup>+75</sup><sub>-60</sub></td><td>-70<sup>+32</sup><sub>-33</sub></td></tr><tr><td>Ω (°)</td><td>99.1 ± 1.6</td><td>11.2 ± 4.6<sup>(c)</sup></td></tr></tbody></table>

**Notes.** The uncertainty is given at $1\sigma$ and contains only the statistical error from the fit and no systematic error from the true north or star registration. $^{(a)}$Stark et al. (2014). $^{(b)}$This value uses the revised Gaia distance to the star of 48.2 pc instead of 51.8 pc from Stark et al. (2014). $^{(c)}$The argument of the ascending node $\Omega$ from Stark et al. (2014) is 90° offset from the value derived in this work. This is likely a problem of definition of the origin, taken here as east of north, as the literature mentions various conventions for the deprojection of discs (Chen et al. 2020).

### 3.2. Optical depth profile

We followed a methodology similar to Stark et al. (2014) in order to extract the optical depth $\tau$ of the ring from the polarised intensity IRDIS data. We deprojected the ring and corrected for the $1/r^2$ stellar illumination factor. We measured the median radial profile of the resulting map to estimate the optical depth of the disk. Because the disc brightness varies azimuthally due to the anisotropy of scattering by the dust, we extracted the radial profile in sectors (four wedges centred respectively on the forward-scattering, backward-scattering and 90°-scattering sides, with an opening angle of 90° each), normalised to one each sector independently and show a weighted average in Fig. 4. We note that the optical depth retrieved this way corresponds to the real geometrical optical depth only if all grains have the same
<figure><img src="image_5.png" /><figcaption>Fig. 4. Deprojected normalised optical depth of the disk at mm wavelengths (orange line extracted from data published in Pawellek et al. 2021), in the optical with STIS (green line Stark et al. 2014, assuming a 2% uncertainty), and in the near-infrared with IRDIS (blue line). The top panel shows the data (plain lines) with the 1σ uncertainty (shade) and the double-power law fit (dotted line limited to the range of separation where the fit was performed). The bottom panel is a zoom in the region where the optical depth peaks, highlighting with vertical dotted line the location of the maxima.</figcaption></figure>

scattering properties, which is probably not the case if the system has a large population of very small grains (see Sect. 5.4), but we use this parameter to allow easy comparison with the STIS and ALMA results.

To compare with the distribution of mm-sized grains we derive the deprojected intensity profile using the python tool FRANKENSTEIN (Jennings et al. 2020)<fnref n="2" />, using the ALMA band 7 data (λ = 0.88 mm) and the best-fit disc orientation presented in Pawellek et al. (2021). The profile is reconstructed from the interferometric visibilities directly, such that the actual resolution is higher than the actual beam size of 0.2" stated in Pawellek et al. (2021). Subsequently, we multiply the intensity profile by $\sqrt{r}$ to estimate the optical depth at mm wavelengths since the intensity is proportional to the optical depth and $T(r) \propto 1/\sqrt{r}$ at this wavelength. The result is shown in Fig. 4 (orange line). The radial profile distribution of mm dust peaks at 81 au (1.69"), consistent with previous estimates for mm-sized grains (Marino et al. 2016; Pawellek et al. 2021). We note that Marino et al. (2016) detected extended emission at 1.3 mm out to 200 au (4.1"), which is not recovered significantly here due to the lower sensitivity to large-scale structure in these observations compared to those at 1.3 mm.

We then compared the normalised optical depth retrieved from STIS, IRDIS and ALMA by fitting double-power profiles parameterised by the equation

$$ \tau(r) = \left( \frac{2}{\left(\frac{r}{r_0}\right)^{-2\alpha_{\text{in}}} + \left(\frac{r}{r_0}\right)^{-2\alpha_{\text{out}}}} \right)^{1/2} \quad (1) $$

The best fit parameters along with their 1σ uncertainty are shown in Fig. 4 (top panel). The fit was performed between 1.47" and 2.70", as shown by the dotted lines. The interest of this global fit, rather than two independent power-laws, is to yield the location of the peak value $r_{\text{max}}$ (which is usually different from $r_0$ unless $\alpha_{\text{out}} = \alpha_{\text{in}}$). A zoom around the maximum value of the optical depth is shown in the bottom panel, with the label giving the location of the peak value $r_{\text{max}}$.

We note first that the steepness of the inner and outer profiles as well as the location of the maximum $r_{\text{max}}$ change significantly with the wavelength of observation. The inner edge profile is steeper in the near-infrared than in the optical, with $\alpha_{\text{in,STIS}} = 9.2 \pm 0.3$ (green dotted line, steeper than the value of $7.0^3$ reported in Stark et al. 2014) and $\alpha_{\text{in,IRDIS}} = 16.2 \pm 0.5$ in the near-infrared with IRDIS (blue dotted line). The STIS and IRDIS instruments do not have the same spatial resolution (~68 mas sampled by only 1.4 px, vs. ~45 mas for IRDIS sampled by 3.7 px) and this difference is likely an effect of the convolution with the instrumental PSF. This could also be physical, as we expect smaller particles to be more sensitive to the Poynting-Robertson drag. The comparison between the outer edge IRDIS and STIS profiles show that both profiles are compatible within error bars between 1.72" (maximum optical depth) and 2.2". Beyond that, the STIS profile appears steeper than that of IRDIS, but one should take this result with caution because the noise is higher in the near-infrared IRDIS image than in the optical STIS image and the background noise is amplified by the $r^2$ multiplication. This noise already starts to affect the profile beyond 2", which may make the outer IRDIS profile less steep, so that the exponent $\alpha_{\text{out,IRDIS}} = -2.8 \pm 0.1$ is likely too shallow. Indeed the difference with STIS ($\alpha_{\text{out,STIS}} = -4.0 \pm 0.1$) is surprising because, in the canonical ring+halo scenario, we would expect the outer halo to be mostly populated by grains on high-eccentricity orbits that are bigger than the blow-out size $s_{\text{blow}}$. For an F5/6V star, $s_{\text{blow}}$ is ~1 µm for compact particles and could reach ~5 µm for 75% porous agglomerates (Arnold et al. 2019), meaning that halo grains should be much larger than $\lambda/2\pi$ for both STIS and IRDIS observations, and we would expect the radial profiles at these two wavelengths to be relatively similar. This similarity might break down if the halo has a significant population of sub-micron grains (see Sect. 5.4), which are poor scatterers in the near IR but not in the optical, but this should lead to a radial profile that is flatter at shorter wavelengths, in contradiction with our results. Another possible reasons for the difference may be the fact that the STIS optical depth profile is estimated from the total intensity image while the IRDIS optical depth profile uses the polarised intensity image. Compact sub-micron particles are expected to scatter more polarised light at near-infrared wavelengths (Rayleigh regime) compared to the optical (Mie regime). In any case, both the STIS and IRDIS profiles are steeper than the canonical $\tau(r) \propto r^{-1.5}$ expected in halos made of $s > s_{\text{blow}}$ grains (Strubbe & Chiang 2006), and a significant population of unbound $s < s_{\text{blow}}$ grains would only flatten the profile instead of render it steeper (Thebault et al. 2023). This could be an indication that an additional process is shaping the outer regions of this system.

Interestingly, the mm dust peak density is offset from the micron-size dust peak density. It is located $48 \pm 4.1$ mas ($2.3 \pm 0.2$ au) closer to the star, suggesting that there is a radial grain

<footnote n="2">² We use an $\alpha$ and $w_{\text{smooth}}$ parameters of 1.03 and $10^{-4}$, respectively. The extracted profile and in particular the peak location were not very sensitive to these choices.</footnote>

<footnote n="3">³ Stark et al. (2014) fit a single-power law to the increasing side of the optical depth profile, explaining the difference with the double-power law fit from Eq. (1).</footnote>
size segregation in the parent belt of the star. While not as pronounced as in the halo, such a size segregation is expected within the parent belt, especially for belts having a significant width (Thebault et al. 2014), as in the case for HD181327⁴.

In addition, smaller dust particles could be subject to a radial migration because of the presence of gas in the disc, which was indeed detected via ¹²CO (2–1) observations with ALMA (Marino et al. 2016). Gas drag can affect dust grains in different ways. Large particles for which radiation forces are negligible ($s \ge 100$ µm) are typically unaffected by gas drag unless gas densities are large enough in which case they tend to migrate inwards. In contrast, small grains tend to have orbital velocities smaller than the gas and thus gas-drag increases their angular momentum inducing an outward migration (Takeuchi & Artymowicz 2001). The end result is that, under the influence of gas, the distribution of small grains can be shifted outwards and concentrated in the outer region of the gas disc, where the gradient in the gas density becomes larger. Unfortunately, the CO distribution and thus gas radial profile were not well constrained due to insufficient signal-to-noise and thus we cannot quantitatively assess this scenario, which would require significantly more gas from another yet undetected species (Olofsson et al. 2022).

# 4. Extraction and analysis of the scattering properties

## 4.1. Strategy

Two main families of techniques exist to extract the SPF from discs (Olofsson et al. 2020): direct extraction using aperture photometry (as done for instance in Milli et al. 2017; Engler et al. 2019; Ren et al. 2019), or forward-modelling with parameterised SPF (as described for instance in Chen et al. 2020; Mazoyer et al. 2020). Olofsson et al. (2020) also introduced a non-parametric iterative method in which the phase function is an output of the modelling, but only applicable in polarimetry. In our case, the IRDIS dataset allows us to constrain the polarimetric and total intensity SPF simultaneously from the same observations, reducing considerably the biases and difficulties arising from combining different observations with possibly different instruments, spatial resolutions, observing conditions, etc. However, the type of noise dominating each image is very different, with a polarimetric intensity image dominated by the disc photon noise in the brightest part of the ring and by detector noise in the rest of the image, while the total intensity image is dominated by instrumental diffraction artefacts and other low-spatial frequency noise that the RDI reduction could not eliminate. The SPF direct extraction is therefore possible for the polarimetric image (see Sect. 4.2) but too uncertain for the total intensity image. We therefore decided to implement a forward-modelling strategy where we used the NICMOS image of the disc to constrain the SPF in total intensity, the IRDIS polarimetric data set was used to constrain the morphology and polarised SPF, and the IRDIS intensity image was only used to scale the disc model implementing the IRDIS morphology and NICMOS SPF to the correct surface brightness (Sect. 4.3). We do not expect significant differences in the SPF between 1.6 µm (IRDIS) and 1.1 µm (NICMOS), hence the choice of this strategy technique.

## 4.2. Extraction of the polarised scattering phase function

To extract the polarised scattering phase function (hereafter pSPF), we started from the polarised intensity image $Q_\phi$ shown in Fig. 2 (middle). The $U_\phi$ image shows, as expected, no astrophysical signal, and this image was therefore used to estimate the noise. The noise was assumed to be azimuthally symmetric, and we estimated a $1-\sigma$ noise map (called $\sigma_{pI}$) and a radial noise profile by computing the root-mean-square (RMS) of the pixel values in 1-px-wide annuli. We then used two different approaches to extract the polarised SPF from the $Q_\phi$ image: 1) we built a disc model and parametrised the pSPF before iterating to find the best model matching the data, or 2) we extracted the pSPF directly from the data.

### 4.2.1. Forward modelling of the polarised intensity image

To find the best model reproducing the data, we used the python implementation of the GraTeR code (Augereau et al. 1999) available as part of the Vortex Imaging Pipeline (VIP<fnref n="5" />; Gomez Gonzalez et al. 2017). The complete description of the input parameters is given in Appendix B of Milli et al. (2017) or in the online VIP documentation<fnref n="5" />. We set here the inclination $i$ and the longitude of the ascending node $\Omega$ to the values described in Table 3, and set the eccentricity to zero to avoid degeneracy in $\omega$ and limit the number of free parameters. Regarding the vertical and radial dust density distribution, we set the reference scale height to $\xi_0 = 6$ au at the reference radius $r_{0,pI}$ (variable parameter close to the semi-major axis $a$ of ~82 au), we used a Gaussian vertical profile $\gamma = 2$, a linear flaring $\beta = 1$. A similar parametrisation was used in Schneider et al. (2006), and these authors constrained the scale height to lie between 4 au $\le \xi_0 \le 8$ au, confirmed later by Stark et al. (2014), hence our choice to set $\xi_0$ to 6 au.

In a preliminary exploration of the parameter space, we noticed the models favoured a very steep inner edge of the disc, which translates into a large inner slope of the radial dust density distribution $a_{in}$ close to or larger than 15, in agreement with the inner edge slope measured directly in the convolved image and shown in Fig. 4. Large values of $a_{in}$ are difficult to constrain due to the limited angular resolution of the image, and discs with $\alpha_{in}$ larger than 15 would appear the same after convolution with the instrumental PSF. We therefore set this parameter to 15.5.

We let the outer radial density slopes $a_{out,pI}$ as a variable parameter. We used a custom phase function that we optimised to best reproduce the data. As we want to keep the number of free parameters as low as possible, we observed that the extracted polarised SPF can be well approximated by a cubic interpolation between five points at 0°, 60°, 90°, 120° and 180° scattering angle (this assumption will be validated later in Sect. 4.2.2). Three of those points are kept as free parameters, for scattering angles of 60°, 90° and 120° corresponding to the range seen from the Earth. The pSPF at 180° and 0° are set to zero and twice the value at 60° respectively, to provide a smooth SPF between 60° and 120° but the value has no physical meaning because it is not probed from the Earth, and it has no impact on the result of the fit.

In total there are five free parameters for the fit: the outer slope $a_{out,pI}$, the reference radius $r_{0,pI}$, and the 3 interpolation points of the polarised SPF called here $s_{60,pI}$, $s_{90,pI}$, $s_{120,pI}$ that

<footnote n="4">⁴ Marino et al. (2016) found a width of 23 au in the mm, which corresponds to a relative width of ~25%.</footnote>

<footnote n="5">⁵ A complete description and tutorial of the disc forward-modelling capabilities provided by VIP is available in the online documentation at https://vip.readthedocs.io/en/latest/tutorials/06_fm_disk.html</footnote>
<table><caption>**Table 4.** Summary of the free parameters (first column) and the interval over which a uniform prior was assumed (second column) and their most likely value for each of the three fits performed, to minimise the residuals in polarised intensity with IRDIS (third column), in total intensity with NICMOS (fourth column) or the combination of polarised intensity with IRDIS, and total intensity with NICMOS and IRDIS (fifth column).</caption><thead><tr><th>Parameter</th><th>Range</th><th>Best fit χ<sup>2</sup><sub>pI</sub></th><th>Best fit χ<sup>2</sup><sub>I,NIC</sub></th><th>Best fit χ<sup>2</sup><sub>sum</sub></th></tr></thead><tbody><tr><td>s<sub>60,pI</sub></td><td>[0, 7]</td><td>3.42 ± 0.48</td><td></td><td>3.52 ± 0.50</td></tr><tr><td>s<sub>90,pI</sub></td><td>[0, 5]</td><td>2.30 ± 0.32</td><td></td><td>2.34 ± 0.30</td></tr><tr><td>s<sub>120,pI</sub></td><td>[0, 3]</td><td>1.33 ± 0.34</td><td></td><td>1.34 ± 0.29</td></tr><tr><td>a<sub>out,pI</sub></td><td>[-9, -1]</td><td>-4.32 ± 0.95</td><td></td><td>-4.59 ± 1.07</td></tr><tr><td>r<sub>0,pI</sub> (au)</td><td>[65, 100]</td><td>79.9 ± 1.4</td><td></td><td>80.02 ± 1.08</td></tr><tr><td>s<sub>60,I</sub></td><td>[0, 7]</td><td></td><td>4.66 ± 0.37</td><td>4.47 ± 0.33</td></tr><tr><td>s<sub>90,I</sub></td><td>[0, 5]</td><td></td><td>2.86 ± 0.22</td><td>3.00 ± 0.21</td></tr><tr><td>s<sub>120,I</sub></td><td>[0, 3]</td><td></td><td>2.40 ± 0.30</td><td>2.24 ± 0.25</td></tr><tr><td>a<sub>out,I</sub></td><td>[-9, -1]</td><td></td><td>-7.37 ± 0.65</td><td>-7.39 ± 0.58</td></tr><tr><td>r<sub>0,I</sub> (au)</td><td>[65, 100]</td><td></td><td>81.81 ± 0.83</td><td>81.80 ± 0.85</td></tr><tr><td>Λ</td><td>[0, 1.8]</td><td></td><td></td><td>0.84 ± 0.09</td></tr></tbody></table>

Notes. The uncertainty is stated at $5\sigma$.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Interpretation of the IRDIS polarised intensity. Left: polarised intensity $Q_\phi$. Middle: best scattered light model explaining the data (unconvolved). Right: residuals after subtraction of the best model.</figcaption></figure>

combine the polarised SPF interpolated at 60°, 90° and 120° respectively and the disc total scattering cross-section to avoid introducing an additional scaling factor to match the brightness of the disc. The subscript ‘pI’ refers to polarised intensity, and is used to distinguish with the total intensity parameters using the subscript ‘I’.

Rigorously, one should construct synthetic Stokes Q and U images from the ray-traced image of a disc model, assuming purely tangential linear polarisation, then convolve those Q and U images with the measured IRDIS PSF, and reconstruct a synthetic $Q_\phi$ image to be compared to the data. This guarantees the effect of the convolution is correctly accounted for (see Appendix A of Engler et al. 2018, for a detailed description of that procedure). In our case, the disc is extended and shows a cavity much larger than the PSF size, therefore one can directly compare the data with the convolved ray-traced image from the disc model. We therefore implemented this strategy, which is faster, after checking that the difference is negligible with the slower but more rigorous approach<fnref n="6" />. We iterated on the

five free parameters to derive the best model explaining the data, by minimising the chi squared $\chi^2_{pI}$. It is defined as

$$ \chi^2_{pI} = \sum_i \left[ \frac{Q_\phi - M_{pI}(s_{60,pI}, s_{90,pI}, s_{120,pI}, a_{out,pI}, r_{0,pI})}{\sigma_{pI}} \right]^2, \quad (2) $$

where $M_{pI}$ is the disc model in polarised intensity, $\sigma_{pI}$ is the noise map, and $i$ is the index of the pixels used to compute the sum. A mask was created to encompass only the pixels between 1.3" and 2.2" where disc signal is detected.

We used an MCMC with 100 walkers and 500 steps. The chains converged, we considered the first 20% of the steps as burn-in phase. The results from the optimisation are summarised in Table 4 (column Best fit $\chi^2_{pI}$). The full posterior probability distribution of the free parameters is available in Fig. A.1. The best disc polarised intensity model is shown in Fig. 5 (middle), the reduced $\chi^2$ is 0.76 indicating a good match to the data, as visible in the residual image after subtraction of the best model Fig. 5 (right). We reproduced in Fig. 6 the best scattering phase function, as well as 100 models randomly drawn from the posterior distribution of the parameters $s_{60,pI}$, $s_{90,pI}$ and $s_{120,pI}$.

<footnote n="6">⁶ The relative change between the $Q_\phi$ image computed with both methods is less than 0.6% in the area of the image where there is detectable disc signal. This is negligible with respect to uncertainty associated with other sources of noise, such as the variability of the PSF during the observations.</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Polarised SPF estimated using the direct extraction method (blue curve) and the parametric model (black curve). The 1-σ uncertainty is shown in blue and black shade respectively. Both curves have been normalised to the same value at a scattering angle of 90° to allow a meaningful comparison.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. SPF in intensity parametrising the dust scattered light measured by NICMOS the best. The 3 degrees of freedom of the SPF are shown as black dots, and the model performs a cubic interpolation (black curve).</figcaption></figure>

### 4.2.2. Direct extraction

To avoid depending on the parametrisation of a disc model, the SPF can also be estimated directly from aperture photometry after accounting for various effects. We implemented this alternative approach as a sanitary check and followed the methodology detailed in Milli et al. (2017). First, we regularly sampled the best ellipse (as defined in Table 3). The spacing between each point was set to one resolution element. We associated to each point in the plane of the sky a unique scattering phase angle assuming a thin disk with zero scale height. We performed aperture photometry in elliptical apertures to account for the inclination of the system. The relative effect of the convolution and the offset from the star has no impact within our error bars. The limb brightening effect described in Olofsson et al. (2020), that can create brightness enhancements in the ansae of a disc has a negligible impact here because of the low inclination of HD 181327, we therefore did not implement any correction for that. The flux encircled in the aperture is therefore directly proportional to the polarised SPF of the disc. It is overplot in Fig. 6 (blue line) after normalising it to the same value as the pSPF extracted from the model-based approach (black line). Both approaches agree within error bars, which validates our extraction procedures. The direct extraction yields higher uncertainty because it is a non-parametric approach.

The pSPF varies linearly with the scattering angle between 60° and 120°. This linear decrease is also visible in the pSPF of HR 4796 (Arriaga et al. 2020) in the H and K band with a similar slope.

### 4.3. Total intensity

A direct extraction of the SPF in total intensity is not possible from the IRDIS total intensity image alone. This image suffers from various artefacts (AO features and low-frequency noise left behind by the RDI reduction) varying azimuthally along the ring that bias any measurement. We illustrate that problem in Appendix B. The model-based SPF retrieval technique is also not applicable here, because the retrieved error bars are too large to allow a meaningful interpretation. This image however still contain some important information on the overall surface brightness of the disc in total intensity, that is key in order to calibrate the polarised intensity and extract the degree of linear polarisation.

We therefore used the NICMOS image at a slightly different wavelength (1.104 µm vs. 1.625 µm for IRDIS) to model the SPF. The model used is the same as that used in polarised intensity, with five free parameters: the outer slope $a_{out,I}$, the reference radius $r_{0,I}$, and the 3 interpolation points of the SPF called $s_{60,I}$,

$s_{90,I}$, $s_{120,I}$. We decided to allow the reference radius and the outer slope to be different from the best fit values found in polarised light $a_{out,PI}$ and $r_{0,PI}$, because we found that NICMOS favoured a steeper outer edge than IRDIS implying also a slightly larger reference radius ($a_{out}$ and $r_{0}$ are correlated), and our focus is here on the SPF that can be best constrained if the morphology is also well described.

The effect of the PCA data reduction was taken into account through the following steps: for each disc on-sky orientation in the raw images, the disc model is convolved with the instrumental PSF, projected onto the Karhunen-Loève modes, and the resulting images are de-rotated and combined to create a forward model called $FM_{I,NIC}$ that can be directly subtracted from the NICMOS image $I_{NIC}$ to test the quality of the model. We used a similar metric as for IRDIS to estimate the quality of a model: the chi squared $\chi^2_{I,NIC}$ computed in a circular annulus where the disc signal is detected:

$$ \chi^2_{I,NIC} = \sum_i \left[ \frac{I_{NIC} - FM_{I,NIC}(s_{60,I}, s_{90,I}, s_{120,I}, a_{out,I}, r_{0,I})}{\sigma_{I,NIC}} \right]^2, \quad (3) $$

where $\sigma_{I,NIC}$ is the NICMOS noise map. This noise map was estimated by processing the reference star images from the NICMOS PSF libraries with the same method and reduction parameters as the NICMOS science images of HD 181327, as explained in Choquet et al. (2018). The PSF-subtracted libraries were then partitioned into sets with the same number of frames as HD 181327, rotated with the target image orientations, and combined. The noise maps were computed from the pixel-wise standard deviation across these sets of processed reference star images.

We implemented again a maximum likelihood approach, using an MCMC framework. The results from the optimisation are summarised in Table 4 (column Best fit $\chi^2_{I,NIC}$). The full posterior probabilities are shown in Appendix A.2 and we show in Fig. 7 the NICMOS image, the best model and the residuals after subtraction of the forward model. The SPF parametrising the best model is reproduced in Fig. 7, as well as 100 models randomly drawn from the posterior distribution of the parameters $s_{60,I}$, $s_{90,I}$ and $s_{120,I}$. The SPF is monotonously decreasing with scattering angle and seems to flatten to large angles. No backward scattering is detected beyond 110°, as visible in the optical (Stark et al. 2014) for some orbital distances. However we note that the SPF is similar to the SPF measured in the optical at a separation of 105 au (Fig. 7, top panel in Stark et al. 2014, corresponding to a radius of 97.7 au with the revised Gaia distance). Inspecting our residual image (Fig. 8, right) and given the SPF uncertainty (Fig. 7), our NIR observations support the absence of backward-scattering at scattering angles less than 120°.
<figure><img src="image_9.png" /><figcaption>Fig. 8. Interpretation of the NICMOS total intensity data. From left to right: NICMOS image reduced with PCA RDI, best model explaining the data, same model after forward modelling, residuals left after forward modelling, and subtraction from the data.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Interpretation of the IRDIS total intensity data. Left: IRDIS total intensity image. Middle: best unconvolved model from the combined MCMC fit. Right: residuals after forward modelling and subtraction of the best model (right).</figcaption></figure>

## 4.4. Combining total intensity and polarimetry

As explained in Sect. 4.1, we now combine the constraints from the IRDIS $Q_\phi$ image and the NICMOS $I_{NIC}$ image in order to find the best model reproducing the IRDIS intensity image $I_{IRD}$. We do not expect strong chromatic variations of the SPF in the range of scattering angles $60^\circ-120^\circ$. Chen et al. (2020) found for instance no dependence of the near-infrared SPF between 1.1 µm and 2.2 µm on the system HR 4796. We therefore assumed that the SPF derived from NICMOS at 1.104 µm is still valid at the IRDIS wavelength of 1.625 µm, and used the parameters $s_{60,1}$, $s_{90,1}$, $s_{120,1}$ in the model reproducing the disc seen in IRDIS total intensity. We refer the reader to the Appendix C.2 for a discussion on the validity of this assumption using the dust model introduced later in Sect. 5.1. As these parameters also encode the overall brightness of the disc, we need to introduce a scaling factor $\Lambda$ to match the surface brightness of disc model in the $I_{IRD}$ image. We re-used the parameters $\alpha_{out,pi}$ and $r_{0,pi}$ because we do not expect the morphological parameters to vary between the polarised and total intensity scattered light of the disc. Similarly as for NICMOS, we estimated the quality of the model using a chi squared $\chi^2_{IRD}$ computed in a circular annulus where disc signal is detected with a specific mask avoiding the artefacts in the image:

$$ \chi^2_{\text{I,IRD}} = \sum_i \left[ \frac{I_{\text{IRD}} - F M_{\text{I,IRD}}(s_{60,1}, s_{90,1}, s_{120,1}, \alpha_{\text{out,pi}}, r_{0,pi}, \Lambda)}{\sigma_{\text{I,IRD}}} \right]^2, \quad (4) $$

where $\sigma_{\text{I,IRD}}$ is the noise map associated to the $I_{IRD}$ image and $FM_{\text{I,IRD}}$ is the forward model obtained after convolution and processing of a disc model by the RDI reduction algorithm.

In order to propagate correctly the error bars derived from each individual fit, we performed a global optimisation of the sum of the individual chi squared $\chi^2_{\text{sum}} = \chi^2_{\text{pi}} + \chi^2_{\text{NIC}} + \chi^2_{\text{IRD}}$. We maximised the log likelihood defined as $-\frac{1}{2}\chi^2_{\text{sum}}$ in an MCMC framework, using the same uniform priors as those described for the individual optimisation of the IRDIS polarised intensity model and the NICMOS total intensity model, on top of which we added a uniform prior on the scaling factor $\Lambda$. The results are summarised in the last column of Table 4. The results are compatible within error bar with those obtained on the individual minimisation. As this global optimisation is the only way to obtain robust error bars on our fitted parameters, including the polarised and total intensity phase function, we use these results later on in the paper.

We show in Fig. 9 the best IRDIS total intensity model (middle panel), and the residuals (right panel) after computation and subtraction of the forward model. Thanks to this combined fit, the pSPF and SPF can be shown on the same scale in Fig. 10. The disc in polarised light is about five times fainter than in total intensity.

## 4.5. Polarised fraction

The polarised fraction is the ratio of the pSPF by the SPF. It is shown in Fig. 11 after propagating the error bars from the
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** SPF and pSPF on the same vertical scale (in arbitrary units) as a function of the scattering angle.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Polarised fraction as a function of the scattering angle.</figcaption></figure>

MCMC. The polarised fraction is maximum at about 23.6% ± 2.6%. The position of the maximum is not well constrained and occurs between scattering angles of 70° and 82°. The polarisation curve is asymmetric: it is decreasing beyond 100° and flatter below 80°.

An additional result from this modelling is the total scattered intensity of the disc relative to the star flux, $F_{\text{disc},\lambda}/F_{*\lambda}$, expressed in the H band. Using the best model presented in Fig. 9 and the stellar flux estimation presented in Sect. 2.2, we obtain $F_{\text{sca},\lambda,\text{disc}}/F_{*\lambda} = 0.082\% \pm 0.004\%$. This value is twice smaller than that found in the optical with STIS (Schneider et al. 2014), which is consistent with the fact that the disc is seen as blue between STIS and NICMOS, with $\Delta\text{mag}(\text{STIS} - \text{F110W}) = -0.9$ (Ren et al. 2023), equivalent to a flux ratio of ~2.3. The same quantity can be computed for the total polarised scattered intensity of the disc, $(F_{\text{sca},\lambda,\text{disc}})_{\text{pol}}/F_*$, using the best model presented in Fig. 5. We derive $(F_{\text{sca},\lambda,\text{disc}})_{\text{pol}}/F_{*\lambda} = 0.019\% \pm 0.001\%$. The uncertainty assumes a 5% accuracy on the estimation of the stellar flux. These numbers can be compared to the infrared excess of the system of 0.2% (Lebreton et al. 2012), to estimate an effective scattering albedo $\omega_{\lambda,\text{disc}}$ defined as the relative contribution of scattering by the disc into the line of sight to the total amount of stellar flux attenuation by scattering and absorption (Engler et al. 2023). Using their Eq. (9), this yields

$$ \omega_{\lambda,\text{disc}} = \frac{F_{\text{sca},\lambda,\text{disc}}/F_{*\lambda}}{F_{\text{sca},\lambda,\text{disc}}/F_{*\lambda}+L_{\text{IR,disc}}/L_*} = 0.29 \pm 0.01. $$

# 5. Interpretation and discussion

We discuss in this section the interpretation of the derived scattering properties with various scattering theories and experimental measurements or observations.

## 5.1. Interpretation with spherical particles made of amorphous carbon, silicates and water ice

Lebreton et al. (2012) showed that the SED of the disc can be well reproduced assuming spherical dust particles made of a mixture of silicates (12%), carbonaceous material (23%), amorphous ice (65%), with 65% porosity. Their best fit model is found for particles larger than $a_{\text{min}} = 0.81 \ \mu\text{m} \pm 0.31 \ \mu\text{m}$, with a differential power-law size distribution of exponent $\nu = -3.41 \pm 0.09$. In a first step we therefore investigated the compatibility of this type of dust population with the scattered light properties extracted from our new observations. In Milli et al. (2015, 2019) or Singh et al. (2021), we computed the theoretical SPF and polarised fraction for a sample of 6500 dust compositions and sizes, using the Mie theory and the distribution of hollow spheres (DHS; Min et al. 2005) as provided in the radiative transfer code MCFOST (Pinte et al. 2006). We reused these models to investigate their compatibility with the observations of HD 181327. These models are based on a porous dust particle composed of a mixture of astronomical amorphous silicates (Draine & Lee 1984), carbonaceous refractory material (Rouleau & Martin 1991), and water ice (Li & Greenberg 1998) partially filling the holes created by porosity. However, in this work we find that none of these models could satisfactorily reproduce both the total intensity SPF and polarised fraction of the dust particles. We refer to Appendix C for the details on the modelling approach and results. We note that a tension between the scattering properties and the SED was already reported in Lebreton et al. (2012) and Schneider et al. (2006): the best SED model predicts an albedo 4.5 times larger than what can be inferred from the NICMOS scattered light image, and the total intensity SPF is only weakly forward-scattering while the models predict for particles about the size of the wavelength much higher degrees of forward-scattering. Here, we reveal a new tension when taking into account the polarised fraction. We therefore look for solutions by relaxing one of the two assumptions made initially: we kept a similar composition but assumed a different shape for the particles (Sect. 5.2) or we kept the spherical shape assumption but explored a wider range of optical indices (Sect. 5.3).

## 5.2. Interpretation with a model of agglomerated debris particles

In an attempt to find more compatible models for compact irregular particles, we explored the agglomerated debris particles developed in Zubko et al. (2009), and used to interpret the scattering properties of cometary or circumplanetary dust (e.g. Zubko et al. 2014, 2020b; Arnold et al. 2019). They are a model of cosmic dust particles that approximately reproduces the highly irregular morphology of dust in the Solar System. Their morphology is similar to what was detected by Rosetta in the micron-sized dust particles of comet 67P/Churyumov-Gerasimenko (Bentley et al. 2016). This model uses the discrete dipole approximation (DDA; Purcell & Pennypacker 1973) to compute the scattering properties.

A population of aggregated debris particles is generated with size parameters $x$ between 1 and 32, where $x = \frac{2\pi s}{\lambda}$ in a differential power-law size distribution of exponent $\nu$. At $\lambda = 1.6 \ \mu\text{m}$, the particle sizes range therefore between 0.25 $\ \mu\text{m}$ and 8 $\ \mu\text{m}$, but the exact minimum and maximum size does not impact the results for the exponents $\nu$ between -1.5 and -3 explored here (cf. the study in Zubko et al. 2020a, for details). The material is made of a mixture of silicates supposed to be represented by an optical index $(n, k) = (1.6, 0.005)$, where n and
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Comparison of the best model among the grid detailed in Table 5 (dashed line) with the IRDIS polarised fraction (left, plain red line) and the NICMOS total intensity SPF (right, plain red line). This model provides the smallest reduced $\chi^2_{\text{pf+SPF}}$ of 4.4. The STIS total intensity SPF from Stark et al. (2014) is also shown (right, plain blue line) and compared with the model prediction at the STIS wavelengths (dashed blue line). The two panels also show in green dotted lines one model of aggregated debris particle from Zubko et al. (2009), matching well the polarised fraction but slightly less the total intensity SPF. This model employs 40% of the transparent material 1.6 + 0.001*i* and 60% of more absorbing material 2.43 + 0.59*i* in a size distribution parameterised with $v = -2.3$, and provides a reduced $\chi^2_{\text{pf+SPF}}$ of 10.3.</figcaption></figure>

*k* represent the real and imaginary part, and amorphous carbon with (*n*,*k*) = (2.43, 0.59). The agglomerated debris particles have a highly irregular morphology (see Zubko et al. 2020a, for details and illustrations) resembling what has been found in cometary dust in situ. We manually optimised the volume fraction of silicates and amorphous carbon to find a model matching both the polarised fraction and total intensity SPF (Fig. 12).

We define a goodness of fit using a reduced chi squared to measure the distance between the model and our measurements, both for the polarised fraction ($\chi^2_{\text{pf}}$) and for the total intensity SPF ($\chi^2_{\text{SPF}}$), and overall ($\chi^2_{\text{pf+SPF}} = \chi^2_{\text{pf}} + \chi^2_{\text{SPF}}$). We note that while the polarised fraction is an absolute value that can be directly compared to that of a model, the model SPF is scaled to match the measured NICMOS SPF in the least square sense before computing the $\chi^2_{\text{SPF}}$. We found that the model presented in Fig. 12 achieves a good fit for the polarised fraction with $\chi^2_{\text{SPF}} = 2.7$ and is slightly worse for the total intensity SPF with $\chi^2_{\text{SPF}} = 7.6$, for an overall $\chi^2_{\text{pf+SPF}} = 10.3$. The peak polarisation of 23.5% at a scattering angle of (or below) 80° can be well reproduced by the model. The total intensity SPF is more forward-scattering than the actual observations. This difficulty was already mentioned in Stark et al. (2014) with a model based on purely spherical particles. We considered alternative highly absorbing component with (*n*,*k*) = (1.855, 0.45) (Jenniskens 1993) or more absorbing silicates with (*n*,*k*) = (1.6, 0.01) instead of (1.6, 0.005), meaning for instance a higher fraction of iron in the generally Mg-rich silicate. However none of the those alternative models helped

to improve the fit for the total intensity SPF. A similar tension was already observed for the disc around HR 4796 where no Mie or DHS models could simultaneously reproduce the SPF and polarised fraction (Arriaga et al. 2020), despite a thorough grid search through all possible optical indices independently of the underlying physical material making up the dust. Inspired by this study and their grid search, we explored the optical indices of the dust over a much larger range, and kept the assumption of spherical particles.

### 5.3. Interpretation with spherical particles with a parameterised optical index

We generated the scattering properties of spherical particles in a differential power-law size distribution of exponent *v* between *s*<sub>min</sub> and 1 mm, and with optical indices (*n*, *k*). We used the tool OpTool (Dominik et al. 2021), implementing both the Mie and DHS theory (Min et al. 2005). The range of parameters explored are summarised in Table 5.

While the Mie theory failed to give an appropriate model with an overall $\chi^2_{\text{pf+SPF}}$ below 15 (i.e. better than the aggregated debris particle shown before), the DHS theory yields four models below this threshold, the best one having $\chi^2_{\text{pf+SPF}} = 4.4$. It is shown in Fig. 12.

To provide an estimate of the range of acceptable models for each goodness of fit estimators, we carry out a Bayesian analysis (e.g. Pinte et al. 2006; Duchêne et al. 2010) and assign a probability that the data are drawn from the model parameters.
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Marginal probability distributions of the four free parameters of our models, based on the polarised fraction (orange bars), on the NICMOS SPF (blue bars) or on the combination of both (green bars). These distributions are derived from models created using the DHS theory.</figcaption></figure>

<table><caption>**Table 5.** Parameters for the 16 640 models generated.</caption><thead><tr><th>Parameter</th><th>Min. value</th><th>Max. value</th><th>N<sub>sample</sub></th><th>Sampling</th></tr></thead><tbody><tr><td>Theory</td><td>Mie / DHS</td><td>/</td><td>/</td><td>/</td></tr><tr><td>s<sub>min</sub> (μm)</td><td>0.1</td><td>100</td><td>13</td><td>log.</td></tr><tr><td>n</td><td>1</td><td>4</td><td>16</td><td>linear</td></tr><tr><td>k</td><td>1e-7</td><td>1e2</td><td>10</td><td>log.</td></tr><tr><td>ν</td><td>-2.5</td><td>-5.5</td><td>4</td><td>linear</td></tr></tbody></table>

We do not have any a priori information on these parameters, we therefore assumed a uniform prior, corresponding to a uniform sampling of our parameters by our grid (see Table 5). The probability that the data corresponds to a given parameter set is given by $\Psi = \Psi_0\exp(-\frac{\kappa^2}{2})$ where $\Psi_0$ is a normalisation constant introduced so that the sum of probabilities over all models is unity. The probability given here is only valid within the framework of our modelling and parameter space. Figure 13 shows the inferred probability distribution for each of our four free parameters, after marginalisation against all other three parameters. It is shown here using the DHS theory because the best models are obtained with this theory.

While the total intensity SPF taken alone is not very constraining, the polarised fraction is, and when combined together those two observables show that a power-law exponent of -3.5 is very much favoured by the data. This is consistent with the value predicted by the theoretical collisional cascade in the birth ring (Dohnanyi 1969), although theoretical models of size distributions show very complex distributions near the blow-out size (Thebault & Kral 2019). Sub-micron particles are favoured with $s_{min} < 0.2$ µm, together with a relatively large imaginary part $k$ in the range 1–10, meaning highly absorbing material.

## 5.4. Discussion on the size

Assuming a spherical shape for the particles, the 0.1–0.2 µm minimum particle size is a robust conclusion of our analysis. A sub-micron minimum size of 0.81 ± 0.31 µm was also favoured by the SED analysis of Lebreton et al. (2012), in a power-law distribution consistent with the theoretical collisional cascade. Stark et al. (2014) had already noticed that the total intensity SPF extracted from optical observations was compatible with the Mie theory assuming ~0.1 µm particles but noted that such small grains are below the blow-out size and would not survive long enough in the system to dominate the scattering cross section over the bound grains. Our conclusion here is even more robust, because we combine the SPF with the polarised fraction which also strongly favours sub-micron particles. It could also explain why the disc colour is blue between the optical and the near-infrared (Ren et al. 2023). The blow-out size is estimated to be ~5 µm, based on the best Mie model of Lebreton et al. (2012) assuming a 65% porous ice-rich composition, that is 50 times larger than the particles inferred here. Even if ones assumes compact particles with different compositions, the smallest blow-out size could be as low as 0.8 µm (Arnold et al. 2019) which is still five times larger than what is inferred here, and no realistic compositions yield a $\beta$ parameter lower than the blowout threshold of 0.5 for particle radii above 0.05 µm (see Fig. 6b in Arnold et al. 2019). Such small sub-micron grains should thus be blown away from the system on dynamical timescales and their presence might thus appear counter-intuitive.

The possibility that a population of sub-blowout particles leaves a detectable imprint in scattered light and in the thermal SED of the dust is a scenario that was already proposed in Lebreton et al. (2012) and Stark et al. (2014), or even earlier for debris discs such as HD 141569A (Augereau & Papaloizou 2004). This scenario gained credence when the numerical and theoretical study of Thebault & Kral (2019) showed that, for debris discs at collisional steady-state, there is in fact always
a non-negligible fraction of unbound grains present in the system. This is especially the case for bright debris discs with infrared excess higher than 0.1%, for which the population of sub-blowout particles always leaves a detectable signature in scattered light. The contribution of unbound particles with respect to the overall disc flux at 1.6 µm reaches 40% in their simulation of an A6V star with a very bright disc of infrared excess 0.5% (see their Fig. 4b). In the thermal regime, this contribution can even be higher, reaching 80–90% at wavelengths in the narrow 10 to 20 µm domain, which may very well explain the result of Lebreton et al. (2012). There is, however, a limitation to these simulations, which is that the scattering phase function was assumed to be isotropic. Particles larger than the wavelength are indeed more forward-scattering than sub-micron particles, and for low-inclination discs such as HD 181327, most of the peak of forward-scattering is not visible, which tends to limit the contribution of those larger particles in the flux scattered towards an observer (Mulders et al. 2013). In a recent study, Thebault et al. (2023) explored these issues and found that, for a belt+halo system such as HD 181327, realistic scattering phase functions generally tend to increase the contribution of unbound grains to the flux in scattered light. This contribution can reach 10–30% in the parent ring and up 50 to 70% in the extended halo (see Fig. 2 of that paper). In addition, the presence of gas may also increase the fraction of sub-blowout particles by maintaining them longer than expected in the main parent belt because gas friction might slow down, or even halt their outward motion (Takeuchi & Artymowicz 2001). This scenario was proposed for HD 32297 (Bhowmik et al. 2019).

Another possibility is to invoke a recent catastrophic collision, which would produce a population of particles outside collisional equilibrium, with a large fraction of unbound grains (Kral et al. 2015). Such a massive collisional event was proposed in Stark et al. (2014) to explain a density enhancement seen in the disc at optical wavelength, but we do not detect such an asymmetry in the near-infrared. We note also that the large quantities of unbound grains produced by such a transient catastrophic event would be blown out of the system on relatively short timescales (Jackson et al. 2014), thus making the probability to witness such an event very low.

We can, however, not rule out the fact that the preference for sub-micron particles might be an artefact from the DHS or Mie theory, assuming perfectly spherical particles, while real disc particles are likely irregular and rough at their surface. If those irregularities or this surface roughness is of the order of 0.1 µm, DHS may lead to a better fit with $s_{min} = 0.1$ µm particles, as discussed in Min et al. (2005).

## 5.5. Discussion on the composition and shape of the particles and comparison to Solar-System comets

Our exploration of the optical index of the particles favours a relatively absorbing material ($k \sim 1$) and highly refractive ($n \sim 3.4$). In Fig. 14, we show a probability map of the optical index ($n, k$) based on the combination of the polarised fraction and SPF. It was obtained by marginalising (i.e. integrating) our probability distribution over the minimum particle size $a_{min}$ and the slope $v$. The sum of the probability of this map is 1, and the map sharply peaks at 83% for ($n, k$) = (3.4, 1). We see that most standard material used to interpret dust in protoplanetary or debris discs are less refractive ($n \le 3$), meaning that even a mixture of those would not yield a resultant material refractive enough. Troilite (FeS) is an interesting exception as the real part $n$ of the optical index is as large as 4.2 (Pollack et al. 1994) or

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

**Fig. 14.** Marginal probability map of the optical index ($n, k$) of the dust based on the polarised fraction and NICMOS SPF. We overplotted optical constants of amorphous carbon (Zubko et al. 1996; BE = BEnzene burning, ACAR = arc discharge between Amorphous Carbon electrodes in an ARgon atmosphere), water ice, olivine organics, iron and troilite from Pollack et al. (1994) and astrosilicates from Draine & Lee (1984). The red (respectively purple) line shows a mixture of astrosilicates (resp. amorphous carbon) with troilite.

even 6.8 using the optical constants from Henning & Stognienko (1996), not represented in Fig. 14 as it is beyond the exploration range. There is solid evidence for the widespread occurrence of sulfide minerals, including troilite in primitive bodies of the Solar System. In Antarctic micro-meteorites and in the grains of the 81P/Wild 2 comet returned by the Stardust mission, the presence of Fe–Ni sulfides is abundant and dominated by troilite (Dobrica et al. 2009). Fe–Ni sulfides are also very good dark and opaque mineral candidate to explain the low reflectance of the comet 67P/Churyumov-Gerasimenko (Quirico et al. 2016). Its presence in a small proportion could explain the large value of $n$ favoured by our model. For a reference, the averaged atomic mass fraction of Fe is 7.5% in the comet 67P (Bardyn et al. 2017) and the main carrier of Fe is supposed to be in the mineral phase, in the form of anhydrous sulfides and Fe–Ni alloys (Quirico et al. 2016). Such a proportion of opaque minerals, including sulfides such as pyrrhotite ($Fe_{1-x}S$) and troilite, may possibly also be present in the HD 181327 dust particles.

Interestingly, we also notice on Fig. 14 that the best composition lies relatively close to the amorphous carbon optical constants from Zubko et al. (1996), especially the sample labelled ‘BE’ with ($n, k$) = (2.78, 1.097), which is widely used to interpret circumstellar dust properties (e.g. Tazaki et al. 2023). The best model of Lebreton et al. (2012) used the amorphous carbon optical constants labelled ACAR, slightly further away from the region favoured by the dust scattered light properties, with a volume fraction of 23% for this component. We may therefore conclude that the amorphous carbon is a promising candidate if combined with a more highly refractive material such as troilite or other type of dark absorbing mineral. In the ($n, k$) parameter space represented in Fig. 14, we added the track where a compact (non-porous) mixture of amorphous carbon and troilite would lie (purple line, each symbol along the line representing an increment of 10% in one of the end members). The Bruggeman mixing rule was used to compute the optical constant of the resulting material. The best material candidate for the HD 181327 dust particles correspond to a mass ratio of 60% amorphous carbon and 40% troilite. This is an illustration of one possible mixture leading to a real material compatible with the scattered light properties extracted in this work. This is not
enough to confirm the presence of troilite or estimate its mass
abundance for several reasons. First, there is an infinite num-
ber of solutions, using for instance, different end members, more
than two materials, or adding porosity to the particles. The
red track shows the mix of astrosilicates with troilite. In gen-
eral, when one assumes only two end members mixed together
without porosity, the mixture track is a smooth and curved line
joining the two pure end members. Second, if the particles are
not homogeneous and consist for instance of a core and an exter-
nal coating, the bulk mass abundance may be different from
relative abundance retrieved from the scattering properties. Last,
the shape of the particles may be significantly different from
perfect spheres, implying deviations from the Mie or DHS the-
ory not investigated in this study. Muñoz et al. (2021) showed
that the use of the Mie theory may lead to overestimation of the
refractive index *n* and *k* and underestimation of the particle size
when it is used to fit the polarised fraction of particles of the
order of the wavelength or larger. In this respect, the investiga-
tion of the agglomerated debris particles presented here show
that the an irregular geometry with only amorphous carbon and
silicates alone, already provides a relatively good match to the
data without requiring high optical indices.

The main constituent of the model proposed in Lebreton et al. (2012) is water ice, with a volume fraction of 65%. Here, we see in Fig. 14 that the optical constant of water ice lies relatively far from the likely parameter space. This was also a conclusion from Sect. 5.1 and Appendix C, where no model including water ice could well reproduce the data. The main impact of the large amount of water ice is to produce a very high albedo of 65% for the dust. Stark et al. (2014) showed that such a high albedo is compatible with the scattered light observations, if one assumes a realistic SPF to account for the flux scattered outside the line of sight towards the Earth. Our current work cannot rule out the presence of water ice, if it is mixed with other more absorbing materials that would dominate the scattering properties. Observations at longer wavelengths, for instance in the 3.0 µm water absorption band, could confirm the presence of water ice, as it is now possible thanks to the sensitivity of the *James Webb* Space Telescope observations (see e.g. McClure et al. 2023).

Comets in our Solar System represent an interesting point of
comparison because cold debris rings are considered as a reser-
voir of cometary material releasing smaller particles through
collisions (see Kral et al. 2018; Hughes et al. 2018, for a review).
The dust released in the tail of comets close to their perihelion
passage when they sublimate has been widely studied for a num-
ber of bright comets, in the optical or in the near-infrared (see
Kiselev et al. 2015, for a review). We show in Fig. 15 the
comparison between the polarisation fraction extracted from the
HD 181327 particles and that of comets in the *J*, *H* or *K* band.
For geometrical reasons, it is easier to observe comets at high
scattering angles larger than 120°, especially beyond 160° close
to what is referred to as the negative polarisation branch. But the
NIR positive branch is now sufficiently documented to allow a
meaningful comparison, with an average maximum value *P*<sub>max</sub>
in the 20–30%. The HD 181327 dust properties interestingly
show a similar behaviour, with *P*<sub>max</sub> = 23.6% ± 2.6%. This
is similar to the debris disc HD 35861 (*P*<sub>max</sub> ~ 25%–30%;
Esposito et al. 2018), but larger than the HD 114082 disc
(*P*<sub>max</sub> ~ 17%; Engler et al. 2023) and smaller than HR 4796
(*P*<sub>max</sub> ~ 50%; Arriaga et al. 2020). This analogy with comets
can be used to extrapolate extrasolar dust properties from
what we know from cometary dust based on in situ or remote
observations, but also to guide our modelling approach based
on scattering theories already proven successful at modelling

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Comparison of the HD 181327 polarisation fraction with that from comets in the near-infrared (*J*, *H* or *K* band, adapted from Kiselev et al. 2015).</figcaption></figure>

cometary dust properties. In situ observations of the comet 67P
by the Rosetta mission showed that the cometary particles have
a variety of morphologies, including compact single elongated
grains and larger porous aggregate particles probably formed by
the hierarchical agglomeration of these smaller compact grains
(Bentley et al. 2016).

Numerical simulations have been developed to reproduce
the polarised fraction of some extensively observed comets,
such as 1P/Halley, C/1995 O1 Hale-Bopp or 67P/Churyumov-
Gerasimenko. Satisfactory fits to polarised observations have
been obtained over a wide range of wavelengths with aggregates
of submicron-sized grains mixed with spheroidal particles. Both
are consisting of absorbing organic-type material and weakly
absorbing silicate-type material (Levasseur-Regourd et al. 2008;
Lasue et al. 2009). With the additional constraint from the
total intensity SPF, two other models have proven successful to
reproduce cometary dust. First, the rough spheroid model, rep-
resenting a polydisperse mixture of randomly oriented smooth
and rough spheroids of a variety of aspect ratio, could reproduce
the intensity and polarisation properties of most cometary dust
(Kolokolova et al. 2015). The Gaussian Random Sphere model
of Markkanen et al. (2018) could also reproduce the photopo-
larimetric behaviour of 67P. Those successful numerical models
seem to show that the non-spherical nature of the particles seem
important to properly model the properties of cometary dust, and
their application to circumstellar debris dust is highly relevant
now that detailed data have been extracted on a couple of bright
debris discs.

6. Conclusions

In this paper, we extract the morphology of the HD 181327 debris disc, together with its polarised and total intensity scattered light properties, using near-infrared observations with VLT/SPHERE in the *H* band. We show first that the dust particles scattering predominantly at this wavelength have a different radial distribution from the particles responsible for the thermal emission
at 0.88 mm. Their peak density is offset by 2.3 au ± 0.2 au away from the star. This size segregation can be natural in belts sufficiently wide like HD 181327, or can result from the low amount of gas detected in this disc in previous ALMA observations. The outer slope is also less steep, as expected from the radiation pressure of the central star bringing smaller particles more sensitive to the radiation pressure on eccentric orbits.

We use forward modelling to extract both the total intensity SPF and the polarised intensity SPF. Due to the challenge at obtaining a high S/N total intensity image to extract the SPF of the dust from the SPHERE data alone, we develop an innovative technique, combining HST/NICMOS to constrain the total intensity SPF with SPHERE to constrain the polarised intensity SPF and morphology, yielding a reliable estimate of the polarised fraction. Such an approach may be used in the future to extract the polarisation fraction of other disks, in order to benefit from the low level of stellar residuals from ground-based high-contrast polarised coronagraphic imaging, and the high sensitivity of space-based total intensity imaging. An interesting prospect is for instance to combine SPHERE/IRDIS polarimetry with JWST/NIRCAM imaging to retrieve the polarised fraction and total intensity SPF of either additional discs or at additional near-infrared wavelengths.

In the range of scattering angles 60°–120° accessible for this disc given the system inclination, the maximum polarisation degree is 23.6% ± 2.6%. The total intensity SPF is monotonic and smoothly forward scattering, without evidence for backward scattering beyond 90°, as detected at optical wavelengths. Models based on compact spherical particles made of silicates, amorphous carbon, water ice and porosity that could well reproduce the SED cannot explain simultaneously the polarised and total intensity scattered light properties. We show that sub-micron particles in a power-law size distribution of exponent –3.5, made of a highly refractive material and as absorbing as amorphous carbon is required to explain both the polarisation and total intensity SPF in the assumption of spherical particles. Such a material can be obtained, for instance, by mixing iron sulphide with amorphous carbon, but there exists an infinite number of possibilities. The dust degree of polarisation is strikingly similar to Solar System comets. We therefore speculate that the dust particles may share similar properties, in particular irregular shapes with a mixture of compact grains and aggregated particles, as shown by the in situ measurements of Rosetta on the comet 67P/Churyumov-Gerasimenko. Future work may therefore look into the application of advanced scattering models developed to interpret the irregular compact and aggregated cometary dust to circumstellar dust. That could provide more insights into the composition and shape especially if those models can be coupled to radiative transfer to simultaneously reproduce the dust thermal emission and disk SED.

Access to the optical wavelength range represents another avenue to test the validity of the interpretation in terms of dust composition and size. The optical polarimetric channel of SPHERE, called ZIMPOL (Schmid et al. 2018), can be used to complement NIR polarimetric observations of protoplanetary and bright debris discs (e.g. Milli et al. 2019; Ma et al. 2023). The upcoming *Nancy Grace Roman* Space Telescope Coronagraph instrument (Kasdin et al. 2020) is expected to provide a higher sensitivity to faint debris discs and an accurate calibration of the polarisation fraction better than 3% to allow high-fidelity retrieval of the dust properties in the optical (Anche et al. 2023).

for that work during the interpretation steps of these observations, as part of the EPOPEE project (Etude des POussières Planétaires Et Exoplanétaires). She contributed significantly to the EPOPEE team, a research group funded by the Programme National de Planétologie (PNP) of CNRS-INSU in France, from which we acknowledge financial support. A.C.L.R. was one of the most enthusiastic members, always willing to share her wide expertise with the younger generation, and she has vastly contributed to this scientific collaboration between the disc and Solar System community. She sadly passed away only a few weeks after attending our June 2022 EPOPEE working group meeting in Grenoble. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (The Netherlands), ONERA (France) and ASTRON (The Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (The 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 RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). 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). Finally this project is co-funded by the European Union (ERC) through the projects Dust2Planets no. 101053020, ESCAPE no. 101044152 and COBREX no. 885593. Views and opinions expressed are however those of the authors 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. The project is also supported by the French National Research Agency (ANR-21-CE31-0015-03) as part of the DDISK project (Disk Discoveries: Investigating morphology and duSt properties with breaKthrough data science). R.T. acknowledges financial support from the CNES fellowship. A.Z. acknowledges support from the ANID – Millennium Science Initiative programme – centre Code NCN2021_080. The work by E.Z. was supported by the Institute for Basic Science (IBS-R035-C1).

## References

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

* Anche, R. M., Douglas, E., Milani, K., et al. 2023, PASP, 135, 125001

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

* Arriaga, P., Fitzgerald, M. P., Duchêne, G., et al. 2020, AJ, 160, 79

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

* Augereau, J. C., Lagrange, A. M., Mouillet, D., Papaloizou, J. C. B., & Grorod, P. A. 1999, A&A, 348, 557

* Bardyn, A., Baklouti, D., Cottin, H., et al. 2017, MNRAS, 469, S712

* Bentley, M. S., Schmied, R., Mannel, T., et al. 2016, Nature, 537, 73

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

* Bhowmik, T., Boccaletti, A., Thébault, P., et al. 2019, A&A, 630, A85

* Booth, M., Wyatt, M. C., Morbidelli, A., Moro-Martín, A., & Levison, H. F. 2009, MNRAS, 399, 385

* Chen, C., Mazoyer, J., Poteet, C. A., et al. 2020, ApJ, 898, 55

* Choquet, E., Hagan, J. B., Pueyo, L., et al. 2014, in IAU Symp., eds. M. Booth, B. C. Matthews, & J. R. Graham, 299, 30

* Choquet, É., Perrin, M. D., Chen, C. H., et al. 2016, ApJ, 817, L2

* Choquet, É., Milli, J., Wahhaj, Z., et al. 2017, ApJ, 834, L12

* Choquet, É., Bryden, G., Perrin, M. D., et al. 2018, ApJ, 854, 53

* Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246

* de Boer, J., Langlois, M., van Holstein, R. G., et al. 2020, A&A, 633, A63

* Delorme, P., Schmidt, T., Bonnefoy, M., et al. 2017, A&A, 608, A79

* Dobrica, E., Engrand, C., Duprat, J., et al. 2009, Meteor. Planet. Sci., 44, 1643

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

* Dominik, C., Min, M., & Tazaki, R. 2021, Astrophysics Source Code Library [record ascl:2104.010]

* Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

* Duchêne, G., McCabe, C., Pinte, C., et al. 2010, ApJ, 712, 112

* Engler, N., Schmid, H. M., Quanz, S. P., Avenhaus, H., & Bazzon, A. 2018, A&A, 618, A151

* Engler, N., Boccaletti, A., Schmid, H. M., et al. 2019, A&A, 622, A192

* Engler, N., Milli, J., Gratton, R., et al. 2023, A&A, 672, A1

* Esposito, T. M., Duchêne, G., Kalas, P., et al. 2018, AJ, 156, 47

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

*Acknowledgements.* We dedicate this paper to the memory of Anny-Chantal Levasseur-Regourd (ACLR), who provided significant support and inspiration
Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, AJ, 154, 7

Hagan, J. B., Choquet, É., Soummer, R., & Vigan, A. 2018, AJ, 155, 179

Henning, T., & Stognienko, R. 1996, A&A, 311, 291

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

Jackson, A. P., Wyatt, M. C., Bonsor, A., & Veras, D. 2014, MNRAS, 440, 3757

Jennings, J., Booth, R. A., Tazzari, M., Rosotti, G. P., & Clarke, C. J. 2020, MNRAS, 495, 3209

Jenniskens, P. 1993, A&A, 274, 653

Kasdin, N. J., Bailey, V. P., Mennesson, B., et al. 2020, SPIE Conf. Ser., 11443, 114431U

Kiselev, N., Rosenbush, V., Levasseur-Regourd, A.-C., & Kolokolova, L. 2015, in Polarimetry of Stars and Planetary Systems (Cambridge University Press), 379

Kolokolova, L., Das, H. S., Dubovik, O., Lapyonok, T., & Yang, P. 2015, Planet. Space Sci., 116, 30

Kral, Q., Thébault, P., Augereau, J. C., Boccaletti, A., & Charnoz, S. 2015, A&A, 573, A39

Kral, Q., Clarke, C., & Wyatt, M. C. 2018, Circumstellar Discs: What Will Be Next?, eds. H. Deeg, & J. Belmonte (Springer, Cham), 165

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

Langlois, M., Dohlen, K., Augereau, J.-C., et al. 2010, SPIE Conf. Ser., 7735, 77352U

Lasue, J., Levasseur-Regourd, A. C., Hadamcik, E., & Alcouffe, G. 2009, Icarus, 199, 129

Lebreton, J., Augereau, J.-C., Thi, W.-F., et al. 2012, A&A, 539, A17

Lecavelier Des Etangs, A., Vidal-Madjar, A., & Ferlet, R. 1996, A&A, 307, 542

Levasseur-Regourd, A. C., Zolensky, M., & Lasue, J. 2008, Planet. Space Sci., 56, 1719

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

Ma, J., Schmid, H. M., & Tschudi, C. 2023, A&A, 676, A6

Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, SPIE Conf. Ser., 9908, 990834

Marino, S., Matrà, L., Stark, C., et al. 2016, MNRAS, 460, 2933

Markkanen, J., Agarwal, J., Väisänen, T., Penttilä, A., & Muinonen, K. 2018, ApJ, 868, L16

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

Marshall, J. P., Milli, J., Choquet, É., et al. 2018, ApJ, 869, A10

Marshall, J. P., Milli, J., Choquet, E., et al. 2023, MNRAS, 521, 5940

Matthews, B. C., Krivov, A. V., Wyatt, M. C., Bryden, G., & Eiroa, C. 2014, Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (University of Arizona Press, Tucson), 521

Mazoyer, J., Arriaga, P., Hom, J., et al. 2020, SPIE Conf. Ser., 11447, 1144759

McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, Nat. Astron., 7, 431

Milli, J., Mouillet, D., Lagrange, A.-M., et al. 2012, A&A, 545, A111

Milli, J., Mawet, D., Pinte, C., et al. 2015, A&A, 577, A57

Milli, J., Vigan, A., Mouillet, D., et al. 2017, A&A, 599, A108

Milli, J., Engler, N., Schmid, H. M., et al. 2019, A&A, 626, A54

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

Min, M., Kama, M., Dominik, C., & Waters, L. B. F. M. 2010, A&A, 509, A6

Miret-Roig, N., Galli, P. A. B., Brandner, W., et al. 2020, A&A, 642, A179

Montesinos, B., Eiroa, C., Krivov, A. V., et al. 2016, A&A, 593, A51

Mouillet, D., Lagrange, A.-M., Beuzit, J.-L., & Renaud, N. 1997, A&A, 324, 1083

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

Muñoz, O., Frattin, E., Jardiel, T., et al. 2021, ApJS, 256, 17

Olofsson, J., Samland, M., Avenhaus, H., et al. 2016, A&A, 591, A108

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

Olofsson, J., Thébault, P., Kral, Q., et al. 2022, MNRAS, 513, 713

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

Perrin, M. D., Duchene, G., Millar-Blanchaer, M., et al. 2015, ApJ, 799, 182

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Pollack, J. B., Hollenbach, D., Beckwith, S., et al. 1994, ApJ, 421, 615

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

Quirico, E., Moroz, L. V., Schmitt, B., et al. 2016, Icarus, 272, 32

Ray, A., & Srivastava, D. C. 2008, J. Struct. Geol., 30, 1593

Ren, B., Choquet, É., Perrin, M. D., et al. 2019, ApJ, 882, 64

Ren, B., Choquet, É., Perrin, M. D., et al. 2021, ApJ, 914, 95

Ren, B. B., Rebollido, I., Choquet, É., et al. 2023, A&A, 672, A114

Rodigas, T. J., Stark, C. C., Weinberger, A., et al. 2015, ApJ, 798, 96

Rouleau, F., & Martin, P. G. 1991, ApJ, 377, 526

Ruane, G., Ngo, H., Mawet, D., et al. 2019, AJ, 157, 118

Sauvage, J.-F., Fusco, T., Petit, C., et al. 2016, J. Astron. Telesc. Instrum. Syst., 2, 025003

Schmid, H. M., Bazzon, A., Roelfsema, R., et al. 2018, A&A, 619, A9

Schneider, G., Silverstone, M. D., Hines, D. C., et al. 2006, ApJ, 650, 414

Schneider, G., Grady, C. A., Hines, D. C., et al. 2014, AJ, 148, 59

Schneider, G., Grady, C. A., Stark, C. C., et al. 2016, AJ, 152, 64

Singh, G., Bhowmik, T., Boccaletti, A., et al. 2021, A&A, 653, A79

Smart, W. M. 1930, MNRAS, 90, 534

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

Soummer, R., Perrin, M. D., Pueyo, L., et al. 2014, ApJ, 786, L23

Stark, C. C., Schneider, G., Weinberger, A. J., et al. 2014, ApJ, 789, 58

Strubbe, L. E., & Chiang, E. I. 2006, ApJ, 648, 652

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

Tazaki, R., Ginski, C., & Dominik, C. 2023, ApJ, 944, L43

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

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

Thebault, P., Kral, Q., & Augereau, J. C. 2014, A&A, 561, A16

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

van Holstein, R. G., Snik, F., Girard, H. H., et al. 2017, SPIE Conf. Ser., 10400, 1040015

van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020a, A&A, 633, A64

van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020b, Astrophysics Source Code Library [record ascl:2004.015]

Wyatt, M. C. 2006, ApJ, 639, 1153

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

Zubko, E., Kimura, H., Shkuratov, Y., et al. 2009, J. Quant. Spec. Radiat. Transf., 110, 1741

Zubko, E., Muinonen, K., Videen, G., & Kiselev, N. N. 2014, MNRAS, 440, 2928

Zubko, E., Videen, G., Arnold, J. A., et al. 2020a, ApJ, 895, 110

Zubko, E., Zheltobryukhov, M., Chornaya, E., et al. 2020b, MNRAS, 497, 1536
# Appendix A: Model fitting

## A.1. Extraction of the SPF from IRDIS data

The posterior probability distribution of the five free parameters. parametrising the polarised intensity disc model constrained by the IRDIS image are shown in Fig. A.1.

<figure><img src="image_17.png" /><figcaption>**Fig. A.1.** Posterior distribution of the parameters describing the polarised intensity of the disc.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. A.2.** Posterior distribution of the parameters describing the total intensity of the disc with NICMOS.</figcaption></figure>

## A.2. Extraction of the SPF from NICMOS data

The posterior probability distribution of the five free parameters. parametrising the total intensity disc model constrained by the NICMOS image are shown in Fig. A.2.

## A.3. Combined extraction of the SPF and pSPF from NICMOS and IRDIS data

The posterior probability distribution of the 11 free parameters. parametrising the total intensity disc model and polarised intensity disc model constrained by the NICMOS and IRDIS image are shown in Fig. A.3.
<figure><img src="image_19.png" /><figcaption>**Fig. A.3.** Posterior distribution of the parameters used in the combined fit of the NICMOS and IRDIS total intensity of the disc and the IRDIS polarised intensity.</figcaption></figure>
**Appendix B: Extraction of the total intensity SPF from IRDIS data**

To illustrate the impossibility to retrieve the SPF from IRDIS
total intensity data, we show in Fig. B.1 the residual image after
subtracting from the data two models with different SPF. Both
residual images provide a similar chi squared $\chi_l^2$ compatible with
the data. Even though the sharp adaptive optics features cre-
ated by the deformable mirror are masked and not part of the
optimisation area, diffraction patterns from the spiders are only
partially removed by the RDI reduction techniques. RDI cre-
ates in some region over-subtraction and the forward-modelling
approach supposed to account for this effect is degenerate.

<figure><img src="image_20.png" /><figcaption>**Fig. B.1.** Degeneracy in the SPF extraction from IRDIS alone. The left and middle images show residuals from IRDIS after the subtraction from the data of two disc models using different SPF parametrisations (right graph).</figcaption></figure>
# Appendix C: Incompatibility of the silicate, amorphous carbon and water ice porous particle model

## C.1. Model description and results

We parameterised the composition of such a dust model by the porosity without ice *P*, a fraction of vacuum removed by the ice $p_{H2O}$, and a silicate over organic refractory volume fraction $q_{Sior}$. The size of the smallest particles is written $s_{min}$ and the size distribution follows a differential power-law of exponent $\nu$. The grid of these five free parameters is described in Table C.1.

<table><caption>**Table C.1.** Grid of parameters for the 13 000 models generated.</caption><thead><tr><th>Parameter</th><th>Min. value</th><th>Max. value</th><th>N<sub>sample</sub></th><th>Sampling</th></tr></thead><tbody><tr><td>Scattering theory</td><td></td><td></td><td></td><td></td></tr><tr><td>s<sub>min</sub> (µm)</td><td>0.1</td><td>100</td><td>13</td><td>log.</td></tr><tr><td>p<sub>H2O</sub> (%)</td><td>1</td><td>90</td><td>5</td><td>log.</td></tr><tr><td>P (%)</td><td>0</td><td>80</td><td>5</td><td>linear</td></tr><tr><td>q<sub>Sior</sub></td><td>0</td><td>∞</td><td>5</td><td>[0, 0.5, 1, 2, ∞]</td></tr><tr><td>ν</td><td>-2.5</td><td>-5.5</td><td>4</td><td>linear</td></tr></tbody></table>

For each model, we define a goodness of fit using a reduced chi squared to measure the distance between the model and our measurements, for the polarised fraction ($\chi^2_{pf}$), for the total intensity SPF ($\chi^2_{SPF}$), and overall ($\chi^2_{pf+SPF} = \chi^2_{pf} + \chi^2_{SPF}$). The parameters and goodness of fit of the best models are shown in Table C.2 and their scattering properties are compared to our measurements in Fig. C.2. As also demonstrated in Sect. 5.3, we note that the polarised fraction is much more constraining than the SPF alone, due to the required scaling of the SPF: there are only five model with $\chi^2_{pf} < 10$ while there are 1912 models with $\chi^2_{SPF} < 10$. There is no model able to correctly match both the polarised fraction and the SPF: all models compatible with the polarisation fraction have a SPF steeper than observed. This tension was already revealed in Lebreton et al. (2012) and Schneider et al. (2006), where the models compatible with the spectral energy distribution had an SPF steeper than observed.

<table><caption>**Table C.2.** Goodness of fit estimates and corresponding parameters for the best models with respect to the polarised fraction, total intensity SPF, or both.</caption><thead><tr><th></th><th>best polar. frac</th><th>best SPF</th><th>best combination</th></tr></thead><tbody><tr><td>Theory</td><td>Mie</td><td>Mie</td><td>Mie</td></tr><tr><td>ν</td><td>-3.5</td><td>-3.5</td><td>-4.5</td></tr><tr><td>q<sub>Sior</sub></td><td>0</td><td>0.5</td><td>∞</td></tr><tr><td>p<sub>H2O</sub></td><td>0<sup>(a)</sup></td><td>0<sup>(a)</sup></td><td>0<sup>(a)</sup></td></tr><tr><td>s<sub>min</sub> (µm)</td><td>0.18</td><td>0.56</td><td>3.16</td></tr><tr><td>P</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td>χ<sup>2</sup><sub>pf</sub></td><td><strong>7.1</strong></td><td>77.7</td><td>84.8</td></tr><tr><td>χ<sup>2</sup><sub>SPF</sub></td><td>393</td><td><strong>1.6</strong></td><td>395</td></tr><tr><td>χ<sup>2</sup><sub>pf+SPF</sub></td><td>394</td><td>442</td><td><strong>65.3</strong></td></tr></tbody></table>

**Notes.** (a) When the porosity without ice *P* is 0, there is no empty holes to be filled by water ice, hence no water ice in the composition and the fraction of vacuum removed by the ice $p_{H2O}$ is unconstrained.

To provide an estimate of the range of acceptable models for each goodness of fit estimators, we carry out a Bayesian analysis

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

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

<figure><img src="image_23.png" /><figcaption>**Fig. C.1.** Comparison of the best models with the observations for the polarisation fraction (left) and total intensity SPF (right). The top, middle and bottom row shows respectively the best dust model for the polarisation fraction, for the NICMOS SPF and for the two observables simultaneously.</figcaption></figure>

and assign a probability $\Psi$ that the data are drawn from the model parameters. We do not have any a priori information on these parameters, we therefore assumed a uniform prior, corresponding to a uniform sampling of our parameters by our grid (see Table C.1). The derivation of $\Psi$ follows the same methodology as described in Sect. 5.3. Fig. C.2 shows the inferred probability distribution for each of our five free parameters, after marginalisation against all other four parameters. It is shown here using the Mie theory because the best models are obtained with this theory. It highlights the absence of models compatible with both the polarised fraction and the SPF. The minimum particle size
<figure><img src="image_24.png" /></figure>

**Fig. C.2.** Marginal probability distributions of the five free parameters of our model, based on the polarised fraction (orange bars), on the NICMOS SPF (blue bars), or on the combination of both (green bars). These distributions were derived from models created using the Mie theory

and the silicate to organic ratios favour indeed very different values whether we consider the polarised fraction alone (with pure organic sub-micron particles) or whether we consider also the SPF (favouring pure silicate micron-sized particles in this case). This discrepancy can be due to the inability of the Mie or DHS theory to capture the scattering properties of irregular or porous particles, or to the inadequacy of the dust composition made of only three basic constituents : amorphous carbon, silicates, ice.

## C.2. Validation of the assumption on the achromatic SPF between the J and the H band

When extracting the polarisation fraction from the IRDIS and NICMOS data, we made the assumption that the shape of the SPF in total intensity is the same between the IRDIS H band centred at 1.625 µm and the NICMOS wavelength centred at 1.104 µm (J band). The model developed in Appendix C.1 allows to test this assumption.

Among all the 6500 Mie models, the root mean square error (hereafter RMSE) between the SPF at 1.104 µm and 1.625 µm is 0.85% on average. It is less than 1% in 81% of the 6500 models, and less than 10% in 99% of them, with a maximum at 50% for one specific model. For the best models matching either the polarised fraction alone, the SPF in total intensity alone, or both simultaneously (cf Fig. C.1, the RMSE amounts respectively to 0.12%, 0.38% and 0.57%. The comparison between those SPF is shown in Fig. C.3. We conclude that in most cases the differences between the SPF at J and H is smaller than 1% and does not dominate our uncertainty.

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

**Fig. C.3.** Comparison of the total intensity SPF in the H (plain lines) and J band (dashed lines). The comparison is shown for three dust models representing the best dust model for the polarisation fraction alone, for the NICMOS SPF alone or for the two observables simultaneously.