Open Access
Issue
A&A
Volume 683, March 2024
Article Number A22
Number of page(s) 22
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202347933
Published online 01 March 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 Foma-lhaut (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 LIR,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 , with an inner edge at 76.6 ± 1.0 au and a maximum brightness at 84.2 ± 1.0 au (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.

Table 1

Discs with the degree of linear polarisation extracted over a range of scattering angles.

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 consortium1. 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 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; Soummer 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 (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−4, only detected in scattered light from space (Schneider et al. 2016; Soummer et al. 2014) with a very faint peak surface brightness of 0.2 mJy arcsec−2 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.

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.

thumbnail Fig. 1

Raw coronagraphic frames of HD181327 and HD202917. Left: first (top) and last (bottom) frames of the science target HD181327. Right: two examples of frames of the reference star HD202917, observed immediately after HD181327, and with parallactic angles matching that of the first (top) and last (bottom) raw frame of HD181327. 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).

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) 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ϕ and Uϕ defined as Qϕ = −Q cos(2ϕ) − U sin(2ϕ) and Uϕ = +Q sin(2ϕ) − U cos(2ϕ) (de Boer et al. 2020), where ϕ is the polar angle between the north and the point of interest, measured from the north over east (the position angle). Qϕ > 0 is equivalent to an azimuthal polarisation component while Qϕ < 0 indicates a radial polarisation. The component describes the polarisation in the directions ±45° with respect to the radial direction.

The Stokes Qϕ and Uϕ 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 mas2.

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.

thumbnail Fig. 2

Final images of the intensity (Stokes I, left image), azimuthal Stokes Qϕ (middle) and Uϕ (right), calibrated in mJy arcsec−2 (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.

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 flat-fielded, 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 (λ = 1.104 µm, Δλ = 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ϕ 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/r2 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 ± 10 mas oriented at a position angle PA of 99.1° ± 1.6°, a semi-minor axis b = 1478 ± 9 mas, and its centre is offset by x0 = −6.8 ± 8.2 mas and y0 = 0.5 ± 7.4 mas in right ascension and declination (x0 < 0 means an east offset, y0 > 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 ω and the longitude of ascending node Ω. 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σ 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 ω.

thumbnail Fig. 3

Marginal probability distribution of the parameters of the best ellipse fitting ring: the E–W and N–S offset the ellipse centre x0 and y0, 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.

Table 3

Deprojected ellipse parameters.

3.2 Optical depth profile

We followed a methodology similar to Stark et al. (2014) in order to extract the optical depth τ of the ring from the polarised intensity IRDIS data. We deprojected the ring and corrected for the 1/r2 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 centreed 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 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)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 to estimate the optical depth at mm wavelengths since the intensity is proportional to the optical depth and 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 (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 rmax (which is usually different from r0 unless αout = α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 rmax.

We note first that the steepness of the inner and outer profiles as well as the location of the maximum rmax change significantly with the wavelength of observation. The inner edge profile is steeper in the near-infrared than in the optical, with αin,STIS = 9.2 ± 0.3 (green dotted line, steeper than the value of 7.03 reported in Stark et al. 2014) and αin,IRDIS = 16.2 ± 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 r2 multiplication. This noise already starts to affect the profile beyond 2″, which may make the outer IRDIS profile less steep, so that the exponent αout,IRDIS = −2.8 ± 0.1 is likely too shallow. Indeed the difference with STIS (αout,STIS = −4.0 ± 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 sblow. For an F5/6V star, sblow 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 λ/2π 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 τ(r) ∝ r−1.5 expected in halos made of s > sblow grains (Strubbe & Chiang 2006), and a significant population of unbound s < sblow 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 ± 4.1 mas (2.3 ± 0.2 au) closer to the star, suggesting that there is a radial grain 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 HD1813274.

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 12CO (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 ≳ 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).

thumbnail 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.

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 an 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 NIC-MOS 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ϕ shown in Fig. 2 (middle). The Uϕ 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-σ noise map (called σ) 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ϕ 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 (VIP5; 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 documentation5. We set here the inclination i and the longitude of the ascending node Ω to the values described in Table 3, and set the eccentricity to zero to avoid degeneracy in ω and limit the number of free parameters. Regarding the vertical and radial dust density distribution, we set the reference scale height to ξ0 = 6 au at the reference radius r0,pI (variable parameter close to the semi-major axis a of ~82 au), we used a Gaussian vertical profile γ = 2, a linear flaring β = 1. A similar parametrisation was used in Schneider et al. (2006), and these authors constrained the scale height to lie between 4 au ≤ ξ0 8 au, confirmed later by Stark et al. (2014), hence our choice to set ξ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 ain 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 ain are difficult to constrain due to the limited angular resolution of the image, and discs with ain 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 aout,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 aout,pI, the reference radius r0,pI, and the 3 interpolation points of the polarised SPF called here S60,pI, S90,pI, S120,pI that 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 ‘’ 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ϕ 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 approach6. We iterated on the five free parameters to derive the best model explaining the data, by minimising the chi squared It is defined as (2)

where MpI is the disc model in polarised intensity, σ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 ). 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 χ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 S60,pI, S90,pI and S120,pI.

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).

thumbnail Fig. 5

Interpretation of the IRDIS polarised intensity. Left: polarised intensity Qϕ Middle: best scattered light model explaining the data (unconvolved). Right: residuals after subtraction of the best model.

thumbnail 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.

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 over-plot 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 Η and Κ 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 aout,I, the reference radius r0,I, and the 3 interpolation points of the SPF called S60,I, S90,I, S120,I We decided to allow the reference radius and the outer slope to be different from the best fit values found in polarised light aout,PI and r0,pI, because we found that NICMOS favoured a steeper outer edge than IRDIS implying also a slightly larger reference radius (aout and r0 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 FMi,nic that can be directly subtracted from the NICMOS image INIC 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 computed in a circular annulus where the disc signal is detected: (3)

where σ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 ). 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 S60,I, S90,I and S120,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°.

thumbnail 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).

thumbnail 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.

thumbnail 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).

4.4 Combining total intensity and polarimetry

As explained in Sect. 4.1, we now combine the constraints from the IRDIS Qϕ image and the NICMOS INIC image in order to find the best model reproducing the IRDIS intensity image IIRD. We do not expect strong chromatic variations of the SPF in the range of scattering angles 60º–120º. 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 s60,I, s90,I, s120,I 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 Λ to match the surface brightness of disc model in the IIRD image. We re-used the parameters αout,pI and r0,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 computed in a circular annulus where disc signal is detected with a specific mask avoiding the artefacts in the image: (4)

where σI,IRD is the noise map associated to the IIRD image and F MI,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 . We maximised the log likelihood defined as 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 Λ. 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 planel), and the residuals (right panel) after computation and subtraction of the forward model. Thanks to this combined fit, the pSFP 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 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, Fdisc,λ/F★,λ, 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 Fsca,λ,disc/F★,λ = 0.082% ± 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 ∆mag(STIS − F 110W) = −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, (Fsca,λ,disc)pol/F, using the best model presented in Fig. 5. We derive (Fsca,λ,disc)pol/F★,λ = 0.019% ± 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 ωλ,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 .

thumbnail Fig. 10

SPF and pSPF on the same vertical scale (in arbitrary units) as a function of the scattering angle.

thumbnail Fig. 11

Polarised fraction as a function of the scattering angle.

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 amin = 0.81 µm ± 0.31 µm, with a differential power-law size distribution of exponent v = −3.41 ± 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 in a differential power-law size distribution of exponent v. At λ = 1.6 µm, the particle sizes range therefore between 0.25 µm and 8 µm, but the exact minimum and maximum size does not impact the results for the exponents v 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 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 and for the total intensity SPF , and overall . 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 . We found that the model presented in Fig. 12 achieves a good fit for the polarised fraction with and is slightly worse for the total intensity SPF with , for an overall . 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.

thumbnail 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 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.001i and 60% of more absorbing material 2.43 + 0.59i in a size distribution parameterised with v = −2.3, and provides a reduced of 10.3.

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 smin 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 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 . 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.

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 where Ψ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 & Krai 2019). Sub-micron particles are favoured with smin < 0.2 µm, together with a relatively large imaginary part k in the range 1–10, meaning highly absorbing material.

thumbnail 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.

Table 5

Parameters for the 16 640 models generated.

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 β 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 & Krai (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 smin =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 ~ 1) and highly refractive (n ~ 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 amin 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 ≤ 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 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 81PA/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 (Fe1−xS) 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 number 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 general, 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 external 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 theory 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 investigation 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 reservoir 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 number 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 Pmax in the 20–30%. The HD 181327 dust properties interestingly show a similar behaviour, with Pmax = 23.6% ± 2.6%. This is similar to the debris disc HD 35861 (Pmax ~ 25%–30%; Esposito et al. 2018), but larger than the HD 114082 disc (Pmax ~ 17%; Engler et al. 2023) and smaller than HR 4796 (Pmax ~ 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 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, representing 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 photopolarimetric 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.

thumbnail 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.

thumbnail 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 etal. 2015).

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 Corona-graph 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).

Acknowledgements

We dedicate this paper to the memory of Anny-Chantal Levasseur-Regourd (ACLR), who provided significant support and inspiration 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 Osserva-torio 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 (Investisse-ments 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).

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.

thumbnail Fig. A.1

Posterior distribution of the parameters describing the polarised intensity of the disc.

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.

thumbnail Fig. A.2

Posterior distribution of the parameters describing the total intensity of the disc with NICMOS.

thumbnail 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.

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 compatible with the data. Even though the sharp adaptive optics features created 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 creates in some region over-subtraction and the forward-modelling approach supposed to account for this effect is degenerate.

thumbnail 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).

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 pH2O, and a silicate over organic refractory volume fraction qSior. The size of the smallest particles is written smin and the size distribution follows a differential power-law of exponent ν. The grid of these five free parameters is described in Table C.1.

Table C.1

Grid of parameters for the 13 000 models generated.

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 , for the total intensity , and overall . 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 while there are 1912 models with . 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 C.2

Goodness of fit estimates and corresponding parameters for the best models with respect to the polarised fraction, total intensity SPF, or both.

To provide an estimate of the range of acceptable models for each goodness of fit estimators, we carry out a Bayesian analysis and assign a probability Ψ 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 Ψ 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 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.

thumbnail 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.

thumbnail Fig. C.2

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

C.2 Validation of the assumption on the achromatic SPF between the J and the Η 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 Η 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.l, 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 Η is smaller than 1% and does not dominate our uncertainty.

thumbnail Fig. C.3

Comparison of the total intensity SPF in the Η (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.

References

  1. Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948 [Google Scholar]
  2. Anche, R. M., Douglas, E., Milani, K., et al. 2023, PASP, 135, 125001 [NASA ADS] [CrossRef] [Google Scholar]
  3. Arnold, J. A., Weinberger, A. J., Videen, G., & Zubko, E. S. 2019, AJ, 157, 157 [NASA ADS] [CrossRef] [Google Scholar]
  4. Arriaga, P., Fitzgerald, M. P., Duchêne, G., et al. 2020, AJ, 160, 79 [Google Scholar]
  5. Augereau, J. C., & Papaloizou, J. C. B. 2004, A&A, 414, 1153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Augereau, J. C., Lagrange, A. M., Mouillet, D., Papaloizou, J. C. B., & Grorod, P. A. 1999, A&A, 348, 557 [NASA ADS] [Google Scholar]
  7. Bardyn, A., Baklouti, D., Cottin, H., et al. 2017, MNRAS, 469, S712 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bentley, M. S., Schmied, R., Mannel, T., et al. 2016, Nature, 537, 73 [Google Scholar]
  9. Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bhowmik, T., Boccaletti, A., Thébault, P., et al. 2019, A&A, 630, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Booth, M., Wyatt, M. C., Morbidelli, A., Moro-Martín, A., & Levison, H. F. 2009, MNRAS, 399, 385 [NASA ADS] [CrossRef] [Google Scholar]
  12. Chen, C., Mazoyer, J., Poteet, C. A., et al. 2020, ApJ, 898, 55 [CrossRef] [Google Scholar]
  13. Choquet, E., Hagan, J. B., Pueyo, L., et al. 2014, in IAU Symp., eds. M. Booth, B. C. Matthews, & J. R. Graham, 299, 30 [NASA ADS] [Google Scholar]
  14. Choquet, É., Perrin, M. D., Chen, C. H., et al. 2016, ApJ, 817, L2 [NASA ADS] [CrossRef] [Google Scholar]
  15. Choquet, É., Milli, J., Wahhaj, Z., et al. 2017, ApJ, 834, L12 [NASA ADS] [CrossRef] [Google Scholar]
  16. Choquet, É., Bryden, G., Perrin, M. D., et al. 2018, ApJ, 854, 53 [Google Scholar]
  17. Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
  18. de Boer, J., Langlois, M., van Holstein, R. G., et al. 2020, A&A, 633, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Delorme, P., Schmidt, T., Bonnefoy, M., et al. 2017, A&A, 608, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Dobrica, E., Engrand, C., Duprat, J., et al. 2009, Meteor. Planet. Sci., 44, 1643 [CrossRef] [Google Scholar]
  21. Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531 [Google Scholar]
  22. Dominik, C., Min, M., & Tazaki, R. 2021, Astrophysics Source Code Library [record ascl:2104.010] [Google Scholar]
  23. Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [NASA ADS] [CrossRef] [Google Scholar]
  24. Duchêne, G., McCabe, C., Pinte, C., et al. 2010, ApJ, 712, 112 [CrossRef] [Google Scholar]
  25. Engler, N., Schmid, H. M., Quanz, S. P., Avenhaus, H., & Bazzon, A. 2018, A&A, 618, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Engler, N., Boccaletti, A., Schmid, H. M., et al. 2019, A&A, 622, A192 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Engler, N., Milli, J., Gratton, R., et al. 2023, A&A, 672, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Esposito, T. M., Duchêne, G., Kalas, P., et al. 2018, AJ, 156, 47 [Google Scholar]
  29. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  30. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, AJ, 154, 7 [Google Scholar]
  32. Hagan, J. B., Choquet, É., Soummer, R., & Vigan, A. 2018, AJ, 155, 179 [NASA ADS] [CrossRef] [Google Scholar]
  33. Henning, T., & Stognienko, R. 1996, A&A, 311, 291 [NASA ADS] [Google Scholar]
  34. Hughes, A. M., Duchêne, G., & Matthews, B. C. 2018, ARA&A, 56, 541 [Google Scholar]
  35. Jackson, A. P., Wyatt, M. C., Bonsor, A., & Veras, D. 2014, MNRAS, 440, 3757 [NASA ADS] [CrossRef] [Google Scholar]
  36. Jennings, J., Booth, R. A., Tazzari, M., Rosotti, G. P., & Clarke, C. J. 2020, MNRAS, 495, 3209 [Google Scholar]
  37. Jenniskens, P. 1993, A&A, 274, 653 [NASA ADS] [Google Scholar]
  38. Kasdin, N. J., Bailey, V. P., Mennesson, B., et al. 2020, SPIE Conf. Ser., 11443, 114431U [Google Scholar]
  39. Kiselev, N., Rosenbush, V., Levasseur-Regourd, A.-C., & Kolokolova, L. 2015, in Polarimetry of Stars and Planetary Systems (Cambridge University Press), 379 [CrossRef] [Google Scholar]
  40. Kolokolova, L., Das, H. S., Dubovik, O., Lapyonok, T., & Yang, P. 2015, Planet. Space Sci., 116, 30 [NASA ADS] [CrossRef] [Google Scholar]
  41. Kral, Q., Thébault, P., Augereau, J. C., Boccaletti, A., & Charnoz, S. 2015, A&A, 573, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Kral, Q., Clarke, C., & Wyatt, M. C. 2018, Circumstellar Discs: What Will Be Next?, eds. H. Deeg, & J. Belmonte (Springer, Cham), 165 [Google Scholar]
  43. Krivov, A. V., & Wyatt, M. C. 2021, MNRAS, 500, 718 [Google Scholar]
  44. Langlois, M., Dohlen, K., Augereau, J.-C., et al. 2010, SPIE Conf. Ser., 7735, 77352U [NASA ADS] [Google Scholar]
  45. Lasue, J., Levasseur-Regourd, A. C., Hadamcik, E., & Alcouffe, G. 2009, Icarus, 199, 129 [NASA ADS] [CrossRef] [Google Scholar]
  46. Lebreton, J., Augereau, J.-C., Thi, W.-F., et al. 2012, A&A, 539, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Lecavelier Des Etangs, A., Vidal-Madjar, A., & Ferlet, R. 1996, A&A, 307, 542 [Google Scholar]
  48. Levasseur-Regourd, A. C., Zolensky, M., & Lasue, J. 2008, Planet. Space Sci., 56, 1719 [NASA ADS] [CrossRef] [Google Scholar]
  49. Li, A., & Greenberg, J. M. 1998, A&A, 331, 291 [NASA ADS] [Google Scholar]
  50. Ma, J., Schmid, H. M., & Tschudi, C. 2023, A&A, 676, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, SPIE Conf. Ser., 9908, 990834 [Google Scholar]
  52. Marino, S., Matrà, L., Stark, C., et al. 2016, MNRAS, 460, 2933 [Google Scholar]
  53. Markkanen, J., Agarwal, J., Väisänen, T., Penttilä, A., & Muinonen, K. 2018, ApJ, 868, L16 [NASA ADS] [CrossRef] [Google Scholar]
  54. Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
  55. Marshall, J. P., Milli, J., Choquet, É., et al. 2018, ApJ, 869, 10 [NASA ADS] [CrossRef] [Google Scholar]
  56. Marshall, J. P., Milli, J., Choquet, E., et al. 2023, MNRAS, 521, 5940 [CrossRef] [Google Scholar]
  57. 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 [Google Scholar]
  58. Mazoyer, J., Arriaga, P., Hom, J., et al. 2020, SPIE Conf. Ser., 11447, 1144759 [NASA ADS] [Google Scholar]
  59. McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, Nat. Astron., 7, 431 [NASA ADS] [CrossRef] [Google Scholar]
  60. Milli, J., Mouillet, D., Lagrange, A.-M., et al. 2012, A&A, 545, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Milli, J., Mawet, D., Pinte, C., et al. 2015, A&A, 577, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Milli, J., Vigan, A., Mouillet, D., et al. 2017, A&A, 599, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Milli, J., Engler, N., Schmid, H. M., et al. 2019, A&A, 626, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Min, M., Hovenier, J. W., & de Koter, A. 2005, A&A, 432, 909 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Min, M., Kama, M., Dominik, C., & Waters, L. B. F. M. 2010, A&A, 509, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Miret-Roig, N., Galli, P. A. B., Brandner, W., et al. 2020, A&A, 642, A179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Montesinos, B., Eiroa, C., Krivov, A. V., et al. 2016, A&A, 593, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Mouillet, D., Lagrange, A.-M., Beuzit, J.-L., & Renaud, N. 1997, A&A, 324, 1083 [NASA ADS] [Google Scholar]
  69. Mulders, G. D., Min, M., Dominik, C., Debes, J. H., & Schneider, G. 2013, A&A, 549, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Muñoz, O., Frattin, E., Jardiel, T., et al. 2021, ApJS, 256, 17 [CrossRef] [Google Scholar]
  71. Olofsson, J., Samland, M., Avenhaus, H., et al. 2016, A&A, 591, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Olofsson, J., Milli, J., Bayo, A., Henning, T., & Engler, N. 2020, A&A, 640, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Olofsson, J., Thébault, P., Kral, Q., et al. 2022, MNRAS, 513, 713 [NASA ADS] [CrossRef] [Google Scholar]
  74. Pawellek, N., Wyatt, M., Matrà, L., Kennedy, G., & Yelverton, B. 2021, MNRAS, 502, 5390 [NASA ADS] [CrossRef] [Google Scholar]
  75. Perrin, M. D., Duchene, G., Millar-Blanchaer, M., et al. 2015, ApJ, 799, 182 [Google Scholar]
  76. Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Pollack, J. B., Hollenbach, D., Beckwith, S., et al. 1994, ApJ, 421, 615 [Google Scholar]
  78. Purcell, E. M., & Pennypacker, C. R. 1973, ApJ, 186, 705 [NASA ADS] [CrossRef] [Google Scholar]
  79. Quirico, E., Moroz, L. V., Schmitt, B., et al. 2016, Icarus, 272, 32 [CrossRef] [Google Scholar]
  80. Ray, A., & Srivastava, D. C. 2008, J. Struct. Geol., 30, 1593 [NASA ADS] [CrossRef] [Google Scholar]
  81. Ren, B., Choquet, É., Perrin, M. D., et al. 2019, ApJ, 882, 64 [Google Scholar]
  82. Ren, B., Choquet, É., Perrin, M. D., et al. 2021, ApJ, 914, 95 [CrossRef] [Google Scholar]
  83. Ren, B. B., Rebollido, I., Choquet, É., et al. 2023, A&A, 672, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Rodigas, T. J., Stark, C. C., Weinberger, A., et al. 2015, ApJ, 798, 96 [Google Scholar]
  85. Rouleau, F., & Martin, P. G. 1991, ApJ, 377, 526 [NASA ADS] [CrossRef] [Google Scholar]
  86. Ruane, G., Ngo, H., Mawet, D., et al. 2019, AJ, 157, 118 [NASA ADS] [CrossRef] [Google Scholar]
  87. Sauvage, J.-F., Fusco, T., Petit, C., et al. 2016, J. Astron. Telesc. Instrum. Syst., 2, 025003 [Google Scholar]
  88. Schmid, H. M., Bazzon, A., Roelfsema, R., et al. 2018, A&A, 619, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Schneider, G., Silverstone, M. D., Hines, D. C., et al. 2006, ApJ, 650, 414 [Google Scholar]
  90. Schneider, G., Grady, C. A., Hines, D. C., et al. 2014, AJ, 148, 59 [NASA ADS] [CrossRef] [Google Scholar]
  91. Schneider, G., Grady, C. A., Stark, C. C., et al. 2016, AJ, 152, 64 [NASA ADS] [CrossRef] [Google Scholar]
  92. Singh, G., Bhowmik, T., Boccaletti, A., et al. 2021, A&A, 653, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Smart, W. M. 1930, MNRAS, 90, 534 [NASA ADS] [Google Scholar]
  94. Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28 [Google Scholar]
  95. Soummer, R., Perrin, M. D., Pueyo, L., et al. 2014, ApJ, 786, L23 [Google Scholar]
  96. Stark, C. C., Schneider, G., Weinberger, A. J., et al. 2014, ApJ, 789, 58 [Google Scholar]
  97. Strubbe, L. E., & Chiang, E. I. 2006, ApJ, 648, 652 [NASA ADS] [CrossRef] [Google Scholar]
  98. Takeuchi, T., & Artymowicz, P. 2001, ApJ, 557, 990 [NASA ADS] [CrossRef] [Google Scholar]
  99. Tazaki, R., Ginski, C., & Dominik, C. 2023, ApJ, 944, L43 [NASA ADS] [CrossRef] [Google Scholar]
  100. Thebault, P., & Kral, Q. 2019, A&A, 626, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Thébault, P., & Wu, Y. 2008, A&A, 481, 713 [Google Scholar]
  102. Thebault, P., Kral, Q., & Augereau, J. C. 2014, A&A, 561, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Thebault, P., Olofsson, J., & Kral, Q. 2023, A&A, 674, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. van Holstein, R. G., Snik, F., Girard, J. H., et al. 2017, SPIE Conf. Ser., 10400, 1040015 [Google Scholar]
  105. van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020a, A&A, 633, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020b, Astrophysics Source Code Library [record ascl:2004.015] [Google Scholar]
  107. Wyatt, M. C. 2006, ApJ, 639, 1153 [NASA ADS] [CrossRef] [Google Scholar]
  108. Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321 [Google Scholar]
  109. Zubko, E., Kimura, H., Shkuratov, Y., etal. 2009, J. Quant. Spec. Radiat. Transf., 110, 1741 [NASA ADS] [CrossRef] [Google Scholar]
  110. Zubko, E., Muinonen, K., Videen, G., & Kiselev, N. N. 2014, MNRAS, 440, 2928 [NASA ADS] [CrossRef] [Google Scholar]
  111. Zubko, E., Videen, G., Arnold, J. A., et al. 2020a, ApJ, 895, 110 [NASA ADS] [CrossRef] [Google Scholar]
  112. Zubko, E., Zheltobryukhov, M., Chornaya, E., et al. 2020b, MNRAS, 497, 1536 [NASA ADS] [CrossRef] [Google Scholar]

1

ESO programme 099.C-0147(B).

2

We use an α and wsmooth parameters of 1.03 and 10−4, respectively. The extracted profile and in particular the peak location were not very sensitive to these choices.

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).

4

Marino et al. (2016) found a width of 23 au in the mm, which corresponds to a relative width of ~25%.

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

6

The relative change between the Qϕ 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.

All Tables

Table 1

Discs with the degree of linear polarisation extracted over a range of scattering angles.

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.

Table 3

Deprojected ellipse parameters.

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).

Table 5

Parameters for the 16 640 models generated.

Table C.1

Grid of parameters for the 13 000 models generated.

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.

All Figures

thumbnail Fig. 1

Raw coronagraphic frames of HD181327 and HD202917. Left: first (top) and last (bottom) frames of the science target HD181327. Right: two examples of frames of the reference star HD202917, observed immediately after HD181327, and with parallactic angles matching that of the first (top) and last (bottom) raw frame of HD181327. 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).

In the text
thumbnail Fig. 2

Final images of the intensity (Stokes I, left image), azimuthal Stokes Qϕ (middle) and Uϕ (right), calibrated in mJy arcsec−2 (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.

In the text
thumbnail Fig. 3

Marginal probability distribution of the parameters of the best ellipse fitting ring: the E–W and N–S offset the ellipse centre x0 and y0, 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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 5

Interpretation of the IRDIS polarised intensity. Left: polarised intensity Qϕ Middle: best scattered light model explaining the data (unconvolved). Right: residuals after subtraction of the best model.

In the text
thumbnail 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.

In the text
thumbnail 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).

In the text
thumbnail 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.

In the text
thumbnail 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).

In the text
thumbnail Fig. 10

SPF and pSPF on the same vertical scale (in arbitrary units) as a function of the scattering angle.

In the text
thumbnail Fig. 11

Polarised fraction as a function of the scattering angle.

In the text
thumbnail 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 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.001i and 60% of more absorbing material 2.43 + 0.59i in a size distribution parameterised with v = −2.3, and provides a reduced of 10.3.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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 etal. 2015).

In the text
thumbnail Fig. A.1

Posterior distribution of the parameters describing the polarised intensity of the disc.

In the text
thumbnail Fig. A.2

Posterior distribution of the parameters describing the total intensity of the disc with NICMOS.

In the text
thumbnail 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.

In the text
thumbnail 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).

In the text
thumbnail 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.

In the text
thumbnail Fig. C.2

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

In the text
thumbnail Fig. C.3

Comparison of the total intensity SPF in the Η (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.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.