# The near-infrared degree of polarization in debris disks. Toward a self-consistent approach to model scattered light observations

J. Olofsson, P. Thébault, A. Bayo, Th. Henning, J. Milli

► To cite this version:

J. Olofsson, P. Thébault, A. Bayo, Th. Henning, J. Milli. The near-infrared degree of polarization in debris disks. Toward a self-consistent approach to model scattered light observations. Astronomy & Astrophysics - A&A, 2024, 688, 10.1051/0004-6361/202450100. insu-04836838

HAL Id: insu-04836838

https://insu.hal.science/insu-04836838v1

Submitted on 14 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The near-infrared degree of polarization in debris disks

## Toward a self-consistent approach to model scattered light observations ⭐, ⭐☆

J. Olofsson<fnref n="1" />,<fnref n="2" />✉, P. Thébault<fnref n="3" />, A. Bayo<fnref n="1" />✉, Th. Henning<fnref n="2" />, and J. Milli<fnref n="4" />

Received 23 March 2024 / Accepted 28 May 2024

## ABSTRACT

**Context.** Debris disks provide a unique opportunity to probe the properties of small µm-sized particles, allowing us to peer into the constituents of their parent bodies, namely the young analogs of comets and the asteroids of our Solar System.

**Aims.** In the past, studies of the total intensity phase function – that is, the brightness of the disk as a function of the scattering angle – have proven powerful in constraining the main characteristics of the dust particles in debris disks. Nonetheless, there can remain some degeneracies in the modeling, which can be alleviated when considering polarized intensity observations.

**Methods.** We obtained new near-infrared scattered-light observations of four young debris disks, and used state-of-the-art algorithms to recover the total intensity and linear polarimetric images of the disks. These images allow us to constrain the degree of linear polarization as a function of the scattering angle.

**Results.** All four debris disks are detected in polarized intensity, and three are also recovered in total intensity. We measured a peak degree of polarization of ≤40% for all three disks. For the disk around HD 129590, we are also able to determine the degree of polarization in the radiation-pressure-driven halo. To reproduce the observed polarization fractions, we find that the particles must consist of highly refractive and absorbing material. For HD 129590, by measuring the polarization fraction beyond the birth ring, we constrain the width of the size distribution to be increasingly small toward greater radii, which is compatible with the effect of radiation pressure. We put these findings to the test and present a self-consistent approach to producing synthetic images, assuming different profiles for the radiation pressure strength, and accounting for the presence of unbound grains. We find the contribution of these grains to be especially critical in order to reproduce the increasing degree of polarization with stellocentric distance.

**Conclusions.** Some of our results – namely a very small blow-out size and very large (n, k) values for the optical constants, which are required to reproduce the observed degree of polarization – might seem difficult to reconcile with our understanding of cosmic dust. Similar results have been obtained for other disks and we discuss the current limitation of available light-scattering models as well as possible avenues to alleviate these limitations.

**Key words.** instrumentation: high angular resolution – instrumentation: polarimeters – Kuiper belt: general – zodiacal dust – circumstellar matter

## 1. Introduction

Near-infrared (NIR) scattered-light and millimeter observations of debris disks probe micrometer (µm)- and millimeter(mm)- sized particles, the small end of the size distribution. The lifetime of these particles is smaller than the age of the central star, as they can be removed by radiation pressure, Poynting-Robertson drag, and collisions (Wyatt 2008). This implies that there must be constant replenishment in a collisional cascade initiated from the larger bodies in the disk (Krivov & Wyatt 2021). Characterizing the properties of such planetesimals is key in order to understand how they grew and the process of planet formation in general. Unfortunately, they cannot be observed directly and we therefore

need to rely on indirect methods to try to infer some of their main characteristics. One possible avenue to achieve this goal is to better constrain the properties of the small dust particles that are released in destructive collisions and to try to make inferences as to the characteristics of the parent bodies. For instance, Olofsson et al. (2022a) showed that the µm-sized particles in the debris disk around HD 32297 have optical properties compatible with fluffy, highly porous aggregates. By comparing with Solar System objects, this result indicates that the parent bodies in the disk around HD 32297 are likely pristine comet-like objects that collide infrequently, thus avoiding compaction of the particles.

A powerful tool to constrain some of the properties of the smallest dust particles in a debris disk is the study of the phase function, which is the brightness of the disk as a function of the scattering angle<fnref n="1" />. By modeling the phase function with

<footnote n="1">¹ European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany
e-mail: johan.olofsson@eso.org</footnote>

<footnote n="2">² Max Planck Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="3">³ LESIA-Observatoire de Paris, UPMC Univ. Paris 06, Univ. Paris-Diderot, France</footnote>

<footnote n="4">⁴ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote>* The reduced images are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https:// cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/688/A42</footnote>

<footnote>** Based on observations made with ESO Telescopes at the Paranal Observatory under programs ID 105.20GP.001 and 109.237K.001.</footnote>

<footnote n="1">¹ The complementary angle between the star, the particle, and the observer, between 0° and 180°.</footnote>
light-scattering models (e.g., Mie theory, Mie 1908 for compact spherical grains), one can infer for instance the typical sizes of the grains. The phase functions can be estimated either in total intensity (Stokes *I*) or in linear polarimetry (a combination of Stokes *Q* and *U*) and the disk appears quite different in the two cases. In total intensity, the front side of the disk, along its minor axis, is usually best revealed thanks to the strong forward scattering peak. On the other hand, in linear polarimetric observations, the major axis of the disk is brighter, as the polarized phase function usually peaks (or shows a plateau) at scattering angles close to 90°. Taken in isolation, either approach can reliably provide first-order estimates of the typical sizes and porosity of the particles, but there is a strong synergy in combining both total intensity and polarimetry. As outlined in the conclusions of Min et al. (2016), for aggregates, the shape of the phase function in total intensity is mostly governed by the size of the aggregates, while the polarized phase function is more closely related to the properties of the individual monomers constituting the aggregates. Furthermore, the ratio between polarized flux and total intensity provides the degree of polarization (ideally also as a function of the scattering). As this is a ratio between two images, the dependency on the dust density distribution, for example, is naturally removed, allowing us to further characterize the properties of the dust particles.

Furthermore, NIR observations of debris disks offer the rather unique opportunity to probe different grain sizes, solely depending on the stellocentric distance. Indeed, several drag forces must be accounted for when modeling the dynamics of dust particles, and these forces are size dependent. The main force is radiation pressure, which is parametrized by the $\beta$ ratio between the radiation pressure and gravitational forces. As summarized in Krivov (2010), radiation pressure increases the eccentricity of the smaller particles as soon as they are released, redistributing them in an extended halo beyond the birth ring (Thebault et al. 2023). The eccentricity of the particles depends on $\beta$ (as $e \sim \beta/(1-\beta)$, which in turn depends on particle size $s$ (beyond a few $\mu$m, $\beta$ becomes proportional to $\sim 1/s$). If the particles are released in the birth ring at a distance $a_0$ to the star, then their semi-major axis $a$ will be such that their pericenters lie at $a_0$, meaning $a(1-e) = a_0$. We can thus derive the distance of their apocenters as $a(1+e) = a_0(1+e)/(1-e) \sim a_0/(1-2\beta)$. This means that as the stellocentric distance increases, the maximum grain size of the size distribution decreases. Accordingly, we would expect to find a gradient in the optical properties as a function of the distance to the star.

Over recent decades, the degree of polarization has been estimated for about a dozen debris disks. Using Subaru CIAO K-band observations, Tamura et al. (2006) estimated a 10% degree of polarization for the disk around $\beta$ Pictoris, and found negligible variation across the major axis of the disk between $\sim 2.5''$ and $6''$. At R-band wavelengths, Gledhill et al. (1991) previously reported a slightly larger polarization degree of close to $\sim 17\%$ at larger separations ($15-30''$). Tamura et al. (2006) interpreted the low polarization fraction as the presence of particles with sizes comparable to the wavelength of observation (a few $\mu$m). Graham et al. (2007) presented Hubble Space Telescope (HST) observations of the disk around the low-mass star AU Mic (V-band). These authors found that the degree of linear polarization increases as a function of the projected distance, from 5% in the inner regions to 40% at 80 au. To reproduce the observations, the authors concluded that the dust particles are likely highly porous ($\sim 90\%$) and that ballistic cluster-cluster aggregation is a promising avenue to explain their findings.

Similar results of the degree of polarization increasing with stellocentric distance were obtained for several other disks (HIP 79977, Thalmann et al. 2013; HD 111520, Draper et al. 2016; HD 32297, Asensio-Torres et al. 2016; HD 15115, Engler et al. 2019). As these disks are highly inclined, if not perfectly edge-on, the polarized fraction is estimated along the projected major axis, and information about the dependency on the scattering angle is lost. These studies all found similar results, with a low ($\le 10\%$) degree of polarization, linearly increasing up to $\sim 30-45\%$ at larger separations (only reaching 15-20% for HD 15115). The only face-on disk for which the degree of polarization has been estimated is the disk around the low-mass star TWA 7 (Ren et al. 2021). Combining HST and SPHERE observations, these latter authors estimated a high polarization fraction in the main ring and the faint third and outermost ring of 85% and 75%, respectively. As the disk is close to face-on, this corresponds to scattering angles of close to 90°.

More recently, the degree of polarization as a function of the scattering angle was measured for a handful of disks (HD 35841, Esposito et al. 2018; HD 191089, Ren et al. 2019; HR 4796, Arriaga et al. 2020; HD 114082, Engler et al. 2023). NIR observations reveal variety in the shape of the polarization fraction, which usually peaks at $\sim 90°$ (except for HR 4796, where the peak is closer to $40°$), as well as in the maximum peak values, which range from $\sim 10\%$ up to $\sim 50\%$. In all cases, the degree of polarization as a function of the scattering angle could only be determined for the birth ring, and could not be extracted for the extended halos of the disks.

By design, polarimetric observations do not suffer from significant artifacts when removing the contribution of the star to reveal the faint disks. In summary, the stellar photons are expected to be largely unpolarized, while any photon that has been scattered off by some dust particle will show linear polarization in a preferential direction. Subtracting images obtained with orthogonal polarization directions is therefore an efficient way to remove the stellar photons while keeping photons that have been scattered by the disk. This means that the main challenge in measuring the degree of polarization lies in recovering the image in total intensity and minimizing well-known self-subtraction artifacts (e.g., Milli et al. 2012). This has been a very active field of research over the past few years, and there are several algorithms designed to tackle this kind of issue, such as DI-sNMF (data imputation using sequential non-negative matrix factorization; Ren et al. 2020), MAYONNAISE (Pairet et al. 2021), REXPACO (Flasseur et al. 2021), and mustard (Juillard et al. 2023), among others. On top of these techniques, which are designed to improve the post-processing of the observations, alternative approaches have focused on building large libraries of reference images that can be used to best reproduce the observations and perform reference star differential imaging (RDI, Xie et al. 2022).

In this paper, we present new observations of four young debris disks obtained using the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE, Beuzit et al. 2019) instrument at the Very Large Telescope. The aim here is to measure the degree of polarization of the debris disks and attempt to constrain some of the properties of the dust particles. In Sect. 2 we present the observations, which are analyzed in Sect. 3. Our results are presented in Sect. 4. In Sect. 5 we present a more complex, self-consistent model, with the aim being to reproduce the observations of HD 129590 in order to provide additional context for our results, before discussing our findings in Sect. 6.
<table><caption>**Table 1.** Stellar and disk properties.</caption><thead><tr><th>Star</th><th>SpT</th><th>d<sub>★</sub><sup>(a)</sup><br>(pc)</th><th>f<sub>disk</sub><sup>(b)</sup><br>(10<sup>-3</sup>)</th></tr></thead><tbody><tr><td>HD 191089</td><td>F5</td><td>50.11 ± 0.05</td><td>1.5</td></tr><tr><td>HD 157587</td><td>F5</td><td>99.87 ± 0.23</td><td>3.2</td></tr><tr><td>HD 115600</td><td>F2</td><td>109.04 ± 0.25</td><td>2.3</td></tr><tr><td>HD 129590</td><td>G3</td><td>136.32 ± 0.44</td><td>7.0</td></tr></tbody></table>

Notes. <sup>(a)</sup>Distances from Gaia Collaboration (2016, 2021). <sup>(b)</sup>Values from Esposito et al. (2020).

## 2. Observations and data reduction

The selection of the four disks was performed with the following criteria: to be around young stars, bright, sufficiently inclined to maximize the chances of a detection in total intensity (but not too inclined either so that the front and near sides of the disk can be distinguished, that is, $60^{\circ} \le i \le 85^{\circ}$), and to have been previously imaged at near-IR wavelengths (to ensure that a detection will be likely). The stars had to be observable during a given ESO period (in that case odd-numbered, P105 and later P109). Since at the time of the P105 deadline the star-hopping mode was new, we settled for a relatively small sample of disks that had not been previously observed with SPHERE in polarimetry. In the end, among all possible candidates, the following science targets were retained: HD 191089 (Ren et al. 2019; Esposito et al. 2020), HD 157587 (Millar-Blanchaer et al. 2016; Esposito et al. 2020), HD 115600 (Currie et al. 2015; Gibbs et al. 2019; Esposito et al. 2020; Olofsson et al. 2022b), and HD 129590 (Matthews et al. 2017; Esposito et al. 2020; Olofsson et al. 2022b, 2023), with the following reference stars: HD 191131, HD 158018, HD 117255, and HD 129280. Table 1 summarizes some of the stellar and disk properties.

The disks were observed using the SPHERE/IRDIS instrument (Dohlen et al. 2008), making use of the star-hopping technique (Wahhaj et al. 2021). This mode allows to hop back and forth between the science target and a reference calibrator (on average 15 min per cycle), which has no known companion nor disk. The observations were performed using the dual-beam polarimetric mode (DPI, de Boer et al. 2020; van Holstein et al. 2020), in pupil-tracking, with the BB_H filter. Because the degree of polarization is smaller than unity, there is however the risk that a disk can be detected in total intensity but remains undetected (or poorly recovered) in polarized intensity. To mitigate this risk, we also observed each star individually, in DPI, without observing the reference star (thus allowing to observe the science target longer within the 1–1.5 h allocated observing block).

Table A.1 summarizes the observations used in this study, along with some of the atmospheric conditions during the observations. For the second column of the Table, “SCI” and “CAL” refer to the science and reference stars of the star-hopping sequence, respectively, while “DPI” refers to the stand-alone DPI observations.

### 2.1. Linear polarimetric images

The reduction of the stand-alone DPI observations is done using the IRDAP<fnref n="2" /> package (version 1.3.4, van Holstein et al. 2020), but for the star-hopping sequence, some preprocessing is required.

A star-hopping sequence usually consists of a concatenation containing several OBs, alternating between the science and calibrator stars. At the moment, IRDAP only handles files that have the same OBJECT keyword in the headers. For a given concatenation, we therefore grouped the fits files of the same target together (SCI or CAL) and reduced them independently using IRDAP.

The IRDAP package provides several outputs. First, it returns two images, $Q_{\phi}$ and $U_{\phi}$, the former revealing the polarized signal from the disk (if there is any), the latter can be used as a proxy for the uncertainties as it contains no astrophysical signal. Second, IRDAP also provides a cube of $N_{f}$ frames (the left and right sides of the IRDIS detector are summed together) and a list of parallactic angles for each of the $N_{f}$ frames. This cube will be used to derive the total intensity image. Even though the detector integration time (DIT) is the same for all pairs of SCI-CAL, we normalized all the frames by their corresponding DITs. The left column of Fig. 1 shows the polarimetric $Q_{\phi}$ images, with a square root scale, where all four disks are detected.

### 2.2. Total intensity images

To retrieve the disk in total intensity, we performed data imputation with sequential non-negative matrix factorization, using the DI-sNMF<fnref n="3" /> package (Ren et al. 2020). For each science target we only used the star-hopping sequences and did not try to use the DPI stand-alone observations. We first used the cube (produced by IRDAP, see above) of the associated calibrator to build the NMF components (a non-negative, nonorthogonal basis), applying a central numerical mask with a radius of 8 pixels. These components are then used to perform RDI on the science frames, applying a mask of 1" in radius (0.55" for HD 115600) to make sure the disk signal is not included. For each science frame, a model of the point spread function is constructed from the components, and subtracted to the original image. Afterwards, each frame is de-rotated by its corresponding parallactic angle, and the cube is median-collapsed to produce the final image in total intensity. For all the science targets in our sample, we used 5 NMF components. The resulting images are shown in the central column of Fig. 1, with a square root scaling. The disks around HD 157587, HD 11560, and HD 129590 are well detected in total intensity, while the disk around HD 191089 cannot be recovered. This disk is the faintest in our sample (Table 1) and from the modeling of the polarimetric data, this is also the disk with the smallest inclination, $i \sim 61^{\circ}$ (see later), which is always more challenging for total intensity observations (Milli et al. 2012).

In Appendix A, we further investigate the strength of self-subtracted effects which can severely impair the analysis of total intensity observations. This is done by using the final total intensity image as a “model” of the disk and following an approach similar to forward modeling. We find that when using DI-sNMF on a star-hopping sequence, self-subtraction is not significantly impacting neither the total intensity images nor the degree of polarization.

### 2.3. Merging of the different datasets

As reported in Table A.1, for each science target, several datasets are available. At the time of the call for proposals the star-hopping mode was a recent addition to the SPHERE instrument, and the performance in total intensity for extended emission remained unknown. We had therefore requested two star-hopping sequences of 1.5 h each for each science target to maximize the

<footnote n="2">² Available at https://github.com/robvanholstein/IRDAP</footnote>

<footnote n="3">³ Available at https://github.com/seawander/nmf_imaging</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Observations of HD 191089, HD 157587, HD 115600, and HD 129590 (from top to bottom, respectively). The left subpanel shows the polarimetric observations, the central subpanel shows the total intensity, while the right subpanel shows the degree of polarization (polarimetry over total intensity). The scaling is linear for the right panels and between 0 and 1. For the left and center panels, the scaling is in square root. On the left panel, the horizontal bar represents 1" and the distance in au is reported below. North is up, east is to the right, as indicated by the compass on the leftmost panels.</figcaption></figure>
on-source integration time, on top of the stand-alone DPI obser-
vations. It should also be noted that the star-hopping sequences,
being performed in DPI, can also be used to produce the Q<sub>φ</sub> and
U<sub>φ</sub> images.

Figure A.1 shows all the reductions from all the differ-
ent epochs that we obtained. The four leftmost columns show
the results in total intensity, while the four rightmost columns
show the Q<sub>φ</sub> images. In total intensity we only show the post-
processing of the star-hopping sequences, but for the Q<sub>φ</sub> images,
we highlight in the lower right corner if the image comes from a
star-hopping or stand-alone DPI sequence.

For the DPI images, we median-combined all the different
epochs in which the disks were detected (given that all observa-
tions were done with the same DITs), and this median-combined
image will be used for the rest of this study. Regarding the
total intensity images, we tried combining the different dat-
acubes together, along with their respective parallactic angles,
but this did not yield significant improvement. In the remainder
of this study, we therefore used data coming from the best of the
star-hopping sequence. The last column of Table A.1 shows the
datasets that were used. If there is a tick mark for a “SCI” object,
but not for the corresponding “CAL” entry, this means that we
only used this dataset for the DPI Q<sub>φ</sub> image.

2.4. Degree of polarization

From the polarimetric and total intensity images, we are then
able to estimate the degree of linear polarization, by simply
dividing the former by the latter. The degree of polarization is
displayed in the right column of Fig. 1 using a linear scaling
between 0 and 1.

Simply dividing the two images one with the other is suf-
ficient since no additional astrophysical effects need to be
accounted for. For instance, Olofsson et al. (2020) highlighted
the impact that a non-negligible vertical scale height can have
on the determination of the phase function. For highly inclined
and vertically thick disks, there is a column density enhancement
along the semi-major axis of the disk, solely due to projection
effects. Since this depends on the dust density distribution, this
effect affects the polarimetric and total intensity images the same
way, and therefore cancels out when computing the degree of
polarization.

As mentioned previously, the disk around HD 191089 is not
detected in total intensity, and therefore for the rest of this study,
we focus on the remaining three targets (we refer the reader
interested in the disk around HD 191089 to Ren et al. 2019 and
Esposito et al. 2020). For HD 157587, we successfully recover
the degree of polarization along the projected semi-minor axis
of the disk, in the north east direction. For HD 115600, on the
other hand, the minor axis of the disk is hidden behind the coro-
nagraph, and the degree of polarization can only be measured
along the major axis of the disk. For both these disks, we do not
recover strong signal beyond the birth ring, and cannot estimate
the degree of polarization in the extended halo. This is only pos-
sible for the fourth object, HD 129590, for which we can recover
the degree of polarization along the minor and major axis, as
well as in the halo beyond the birth ring. On the lower right
panel of Fig. 1, the fact that the degree of polarization along the
minor axis is close to 0 does not mean it is unconstrained. At this
location, the disk is well detected both in polarimetry and total
intensity, meaning that the degree of polarization can reliably be
estimated to a few percents (see Sect. 4.3 for a further analysis).
Interestingly, along the major axis, the degree of polarization is
increasing with the stellocentric distance.

Even though the disk around HD 129590 is very bright both
in total intensity and polarimetry, the other two disks around
HD 157857 and HD 115600 are fainter. Therefore, to increase
the signal to noise in the images, we binned the data on an
hexagonal grid (and not a squared grid to avoid possible alias-
ing). For each star, each image (total intensity and polarimetry),
and each hexagon, we compute the median value of the flux
(as well as the median absolute deviation for the uncertainties),
distance to the star, and scattering angle. The latter two quanti-
ties are computed in the midplane of the disk, assuming a given
inclination and position angle for the disk (see Fig. B.2, the scat-
tering angle being the arccosine of the dot product between the
line of sight and the coordinates at the disk midplane). These
angles are obtained from modeling the polarimetric data, which
is described in Appendix B and the results are presented in
Table B.1. The binned images are shown in Fig. 2, with a similar
structure as Fig. 1. The phase functions and degree of polariza-
tion as a function of the scattering presented in the rest of this
study are computed from these binned images.

Overall, we measure peak degrees of polarization that are rel-
atively low, below ~50%, comparable to values reported in the
literature for other debris disks. As further discussed in Sect. 4,
for HD 157857, the maximum degree of polarization is in the
range 40–50%. For HD 115600, the peak value is most likely
below 40%. Finally, for HD 129590, the degree of polarization
is very low in the birth ring, below 20%, but increases to values
closer to 50% as the stellocentric distance increases.

It should be noted that the three disks for which we deter-
mined the degree of polarization were shown to display some
degree of asymmetry by Crotts et al. (2024). They reported that
though the disk around HD 115600 appears axisymmetric and
does not show color or brightness asymmetries, the analysis of
the vertical profile tentatively suggests a possible warp. When
modeling the SPHERE observations of this disk, we also find
a large opening angle (ψ ~ 0.13, Table B.1), suggesting that the
vertical structure of this disk is unusual. However, this should
not affect the inclination and position angle determination (cru-
cial for the scattering angle calculation). For HD 129590, Crotts
et al. (2024) noted a brightness asymmetry that is only detected
in their K1 band data but not in the H band observations. How-
ever, the disk does not appear to have strong color asymmetries
and our model, which assumes a circular disk (Appendix B), can
account for most of the signal in the SPHERE H band observa-
tions. We therefore cannot further comment on the origin of the
asymmetry detected in the K1 GPI observations. Lastly, Crotts
et al. (2024) found a significant color asymmetry for the disk
around HD 157587, the East side of the disk being relatively bluer
than the West side. A possible explanation being that the distri-
bution of small dust particles is not uniform throughout the disk.
This will be further discussed in Sect. 4.

3. Analysis

The main objective is to constrain the properties of the dust parti-
cles in the three debris disks for which the degree of polarization
can be computed. Because we targeted disks with relatively
large inclinations, combined with the angular resolution of the
SPHERE instrument, we can measure the degree of polarization
for a wide range of scattering angles.

3.1. Modeling approach

For a given dust model (see next subsections), we can compute
the goodness of fit for the degree of polarization phase function,
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Same as Fig. 1, but the images were binned onto a hexagonal grid instead of having the native SPHERE resolution (12.26 mas). On the right column, pixels saturating above unity are represented in orange.</figcaption></figure>

but as discussed in Milli et al. (2024) this is not necessarily the
complete picture. For instance, it may well be that the total inten-
sity and polarized intensity phase functions are both completely
off with respect to the observations, but that their ratio is still a
good match to the data. Therefore, there is nonredundant infor-
mation in the total intensity phase function that can be included
in the modeling approach.

A first approach to finding the best fitting model would be
to simply sum the χ² values calculated for the total intensity
phase function and for the degree of polarization as a function
of the scattering angle. Nonetheless, there is the risk that one

may dominate the other, yielding for instance a good match to the
total intensity data and a bad fit to the degree of polarization. For
this reason, we chose to normalize the χ² values coming from the
two observables and computed the ratio between both χ² values
when the models are set to null. This ratio is then used to weight
down one of the χ² value. In practice, the “null” χ² for the total
intensity was always larger than for the degree of polarization,
and the former was weighted down.

Additionally, the dust density distribution, the flux normal-
ization of the observations, the magnitude of the star, all these
parameters do not affect the degree of polarization since it is a
ratio between two quantities that are affected the same way by the aforementioned parameters. This is however not the case for the total intensity phase function taken alone (on top of the geomet- ric effects discussed in Olofsson et al. 2020). For a given dust model, comparing the total intensity phase function directly to the observations is not straightforward. Therefore, prior to com- puting the total intensity $\chi^2$, we first find the scaling factor that minimizes the residuals between the model and the observations (see Eq. (7) of Olofsson et al. 2016). By doing so, we are effec- tively losing information related to the scattering efficiencies of the particles. Given that this part of the modeling is less “com- plete” compared to the modeling of the degree of polarization, the associated $\chi^2$ is weighted down by a factor $\frac{1}{2}$.

In the following, we first describe some of our attempts to model the observations and how they failed to reproduce the data before presenting the approach that we retained for the rest of this study.

## 3.2. Preliminary attempts

To analyze the variation of the degree of polarization as a function of the scattering angle, we first attempted to use the *AggScatVIR*<fnref n="4" /> library (Tazaki & Dominik 2022; Tazaki et al. 2023). It provides the phase functions in total and polarized intensity, at different wavelengths (in our case 1.63 µm), for dif- ferent particle shapes. There are two main families for the shape of the particles, namely, aggregates or irregular grains, and both can have different sizes. For the irregular grains, the sizes, or rather, volume-equivalent radius, range between 0.2 and 1.6 µm, and there are two compositions to choose from. For the aggre- gates family, several aggregation models are available (from compact to fractal aggregates), with different total number *N* of spherical monomers (effectively increasing the particle volume). The monomers themselves can have different sizes ($s_{min}$ = 100, 200, and 400 nm) and compositions. There are however some limitations, for instance, due to the significant computational cost, optical properties for very large aggregates (e.g., *N* = 1024 or 2048) are not always available, and this is especially true for highly fractal particles (e.g., FA1.1, see Fig. 2 of Tazaki et al. 2023).

Given that the models are precomputed and made avail- able in the library, the comparison between the models and the observations is extremely fast. Unfortunately, there are no combinations of particle shapes, sizes, or compositions that can adequately reproduce the measured degree of polarization. This is especially the case in the birth rings of the three disks, where the maximum degree of polarization is ≤40–50%. One possi- ble explanation is that the birth ring should host a wide particle size distribution all the way up to mm-sized grains, which are not included in the *AggScatVIR* library (the largest particle size being of a few µm).

The second approach to modeling the observations is to use optical constants of known material with various compositions, use the Mie theory or the Distribution of Hollow Sphere (DHS, Min et al. 2005) model, and compute the phase function over a size distribution, varying the minimum and maximum grain sizes ($s_{min}$ and $s_{max}$, respectively), the slope of the size distri- bution, and the porosity of the grains. To compute the total and polarized intensity phase functions, we used the *optool*<fnref n="5" /> package (Dominik et al. 2021). We tried different dust compositions, mixing pyroxene (Dorschner et al. 1995), amorphous carbon

(Zubko et al. 1996), and crystalline water ice (Warren & Brandt 2008) in different proportions. Unfortunately, this approach did not yield any promising results, as it was especially challenging to reproduce low degrees of polarization.

## 3.3. Make up your dust

In the end, we settled for an approach similar to the one described in Arriaga et al. (2020, see also Milli et al. 2024). Instead of relying on measured optical constants, they (the real and imagi- nary parts) become free parameters ($n$ and $\log_{10}(k)$). We used the DHS model, assuming the maximum filling factor $f_{max}$ = 0.8. There are two additional free parameters related to the size dis- tribution, which are the minimum and maximum grain sizes ($s_{min}$ and $s_{max}$, respectively). Preliminary tests suggested that the slope of the size distribution and porosity remain largely uncon- strained, and we therefore left these values fixed to –3.5 and 25%. It should be noted that we here assume that the minimum grain size (one of the free parameter) corresponds to the radia- tion pressure blow-out size $s_{blow-out}$. Otherwise, one would need to account for a break in the size distribution, of a given ampli- tude, for a given size ($s_{min} \le s_{blow-out} \le s_{max}$). This is beyond the scope of the exercise and the implications of this assumption are further explored and discussed in Sects. 5 and 6. To identify the most probable solution, we use the *MultiNest* algorithm (Feroz et al. 2009) and the Python package *PyMultiNest* (Buchner et al. 2014).

# 4. Results

Table 2 shows the best fitting parameters for the three debris disks studied here. Results for each individual star are discussed in the rest of this section.

## 4.1. HD 157587

Figure 3 shows the observations and best fit models for the degree of polarization (top) and total intensity (bottom), in an annulus between projected separations of 0.7" ≤ *r* < 1.1" (cor- responding to 70 < *r* < 110 au, tracing the birth ring of the disk). There is a significant dispersion for the degree of polarization, with large uncertainties, especially for scattering angles larger than 40°. This is most likely caused by the overall faintness of the debris disk, especially in total intensity, and most likely not related to the observing conditions. The star was observed with an average seeing of ~0.55", a coherence time of 5.5–6 ms, and the on-sky rotation was not negligible (~50°). The disk around HD 157857 is the one with the second smallest inclination of the sample (*i* ~ 70°), while the disk with the smallest inclina- tion (HD 191089, *i* ~ 61°) is not recovered in total intensity. This might be indicative that recovering the total intensity images of optically thin debris disks remains a challenge if the inclination is smaller than ~75°, even when using state-of-the-art observing and post-processing techniques. Alternatively, the dispersion in the degree of polarization for scattering angles larger than ~40° could also be the consequence of different size distributions on either sides of the disk. As mentioned previously, Crotts et al. (2024) reported a color asymmetry between the east and west sides of the disk. Since Fig. 3 shows the degree of polarization of both the east and west sides, an over-abundance of small grains on one side could result in a different (most likely larger) degree of polarization. This would translate in a larger dispersion close to 90° scattering angle. However, given the low S/N of the detec- tion in both polarized and total intensity and that the brightness

<footnote n="4">⁴ Available at https://github.com/rtazaki1205/AggScatVIR</footnote>

<footnote n="5">⁵ Available at https://github.com/cdominik/optool</footnote>
<table><caption>**Table 2.** Best-fitting parameters to model the total intensity phase function and degree of polarization as a function of the scattering angle.</caption><thead><tr><th>Star</th><th>Region</th><th>s<sub>min</sub><br>(µm)</th><th>s<sub>max</sub><br>(µm)</th><th>n</th><th>log<sub>10</sub>(k)</th></tr></thead><tbody><tr><td>HD 157587</td><td>0.70" ≤ r &lt; 1.10"</td><td>0.03<sup>+0.06</sup><sub>-0.02</sub></td><td>97<sup>+375</sup><sub>-84</sub></td><td>3.1<sup>+0.5</sup><sub>-0.8</sub></td><td>0.69<sup>+0.07</sup><sub>-0.10</sub></td></tr><tr><td>HD 115600</td><td>0.35" ≤ r &lt; 0.50"</td><td>0.11<sup>+0.06</sup><sub>-0.07</sub></td><td>0.49<sup>+0.14</sup><sub>-0.13</sub></td><td>2.0<sup>+1.1</sup><sub>-0.7</sub></td><td>1.6<sup>+0.2</sup><sub>-0.3</sub></td></tr><tr><td>HD 129590</td><td>0.30" ≤ r &lt; 0.50"</td><td>0.36<sup>+0.01</sup><sub>-0.01</sub></td><td>110<sup>+75</sup><sub>-52</sub></td><td>4.3<sup>+0.1</sup><sub>-0.0</sub></td><td>-2.69<sup>+0.40</sup><sub>-2.28</sub></td></tr><tr><td>HD 129590</td><td>0.50" ≤ r &lt; 0.70"</td><td>0.16<sup>+0.02</sup><sub>-0.03</sub></td><td>1.9<sup>+0.4</sup><sub>-0.3</sub></td><td>3.4<sup>+0.1</sup><sub>-0.1</sub></td><td>-0.04<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>HD 129590</td><td>0.70" ≤ r &lt; 0.90"</td><td>0.48<sup>+0.01</sup><sub>-0.01</sub></td><td>0.96<sup>+0.03</sup><sub>-0.03</sub></td><td>3.2<sup>+0.1</sup><sub>-0.1</sub></td><td>0.03<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>HD 129590</td><td>0.90" ≤ r &lt; 1.10"</td><td>0.25<sup>+0.21</sup><sub>-0.04</sub></td><td>1.2<sup>+0.2</sup><sub>-0.4</sub></td><td>2.4<sup>+0.1</sup><sub>-0.1</sub></td><td>-0.11<sup>+0.07</sup><sub>-0.04</sub></td></tr></tbody></table>

**Notes.** The second column shows the results for different annuli for HD 129590 (see text for details).

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Observations and best-fit model for HD 157587 (orange circles and solid black line, respectively). *Top*: degree of polarization. *Bottom*: total intensity phase function. The inner and outer radii for the stellocentric distances are indicated in the upper right of the bottom panel.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Same as Fig. 3 but for HD 115600.</figcaption></figure>

asymmetry is only seen in the *J* band GPI data (Crotts et al. 2024), we opt not to further investigate the origin of this large dispersion.

Still, the best fit model can successfully reproduce the overall shape of both the degree of polarization and total intensity as a function of the scattering angle, even though the strong forward scattering peak is slightly under-estimated in the bottom panel. We find that the minimum grain size has to be very small, but with significant uncertainties, suggesting $s_{min} \le 0.1$ µm, and that the maximum grain size is most likely larger than ~20 µm. Since the spatial region in which the modeling is performed encompasses the birth ring of the disk ($0.7'' \le r < 1.1''$), it is expected that we need a wide range of grain sizes, all the way up to at least a few tens of µm. Regarding the real and imaginary parts of refraction, we obtain values of $n \sim 3.1$ and $k \sim 4.9$, rather large values (see Sect. 6).

## 4.2. HD 115600

Figure 4 shows the observations and best fitting model for the disk around HD 115600, where the degree of polarization is estimated between $0.35'' < r < 0.5''$ (corresponding to $38 < r < 54.5$ au). Overall, the results for this target are unfortunately not very reliable. The disk is the second faintest of our sample, it is more compact that the one around HD 157587, and its inclination is larger. Combined with a smaller on-sky rotation during the star-hopping sequence, we are not able to accurately recover the projected semi-minor axis of the disk, especially in total intensity. As shown in the bottom panel of Fig. 4, it is very likely that we are missing the forward scattering peak, casting doubts on the robustness of the results. It seems unlikely that we are probing a very narrow size distribution ($s_{min} \sim 0.1$ µm and $s_{max} \sim 0.5$ µm) in the birth ring of the disk, on top of the unrealistic imaginary part for the optical constants ($k \sim 36$). Nonetheless, we can still reliably constrain the degree of polarization to be $\le 40\%$ along the major axis of the disk.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Same as Fig. 3 but for HD 129590, and this time for four concentric annuli (see text and annotations in the upper right corner of the bottom panels).</figcaption></figure>

### 4.3. HD 129590

The disk around HD 129590 is better recovered in total intensity compared to HD 115600. According to Table A.1, the seeing and coherence time were not too dissimilar for the two datasets, and in fact, the observing conditions were on average better for HD 115600 compared to HD 129590. Besides the fact that the disk around HD 129590 is brighter than the one around HD 115600 (Table 1), another main difference is the on-sky rotation during the observations, which is a factor two larger for HD 129590$^6$, significantly helping to recover the disk in total intensity. The halo beyond the birth ring is well detected both in total and polarized intensity, and as mentioned before, this is the only disk in our sample for which we are able to constrain the degree of polarization in the halo.

Since radiation pressure and gas drag (cold gas has been detected around HD 129590, Kral et al. 2020) are size-sorting processes (see Olofsson et al. 2022b for the latter), there should be a size segregation at different stellocentric distances (Thebault et al. 2014). We can therefore attempt to estimate the degree of polarization as a function of the scattering angle for different regions beyond the birth ring. We proceed the same way as for HD 157587 and HD 115600, but instead of focusing only on the birth ring, we compute the degree of polarization in four concentric rings: [0.3", 0.5"] (encompassing the birth ring, $41 < r < 68$ au), [0.5", 0.7"], [0.7", 0.9"], and [0.9", 1.1"], corresponding to [68, 95] au, [95, 123] au, and [123, 150] au, respectively.

Figure 5 shows the results for the four regions. The top panels show the best fit models to the degree of polarization as a function of the scattering angle, and the bottom panels show the total intensity phase functions for the four different regions. According to the results reported in Table 2, all four regions require a minimum grain size that is below 0.5 µm (between 0.16 and

0.48 µm). Regarding the optical constants, the value of *n* and *k* are all fairly comparable to each other outside of the birth ring, even if the different regions are modeled independently from each other (2.4 < *n* < 3.3 and 0.77 < *k* < 1.1). In the birth ring, the best fit model requires a larger *n* value and smaller *k* imaginary part to reproduce the low degree of polarization. Interestingly, we see that the width of the size distribution is overall decreasing was we venture outside the birth ring. In the birth ring we find *s*<sub>max</sub> ~ 110 µm, and the maximum grain sizes in the next three rings are 1.9, 1.0, and 1.2 µm, from the inner to the outer regions, respectively. Figure 4 of Olofsson et al. (2023) shows the radiation pressure ratio *β* as a function of the particle size *s* for HD 129590, and they found that the blow-out size for this solar-type star would be close to 1.2 µm. The comparison is not one-to-one though because they assumed a different dust composition (a mix of pyroxene and carbon). Nonetheless, the differences between the small minimum grain size derived here and the blow-out size reported in Olofsson et al. (2023) will be further discussed in Sect. 6. It should however be noted that with the approach that we followed here, we cannot constrain the density of the dust particles, and therefore cannot compute the variation of *β* as a function of *s*. Still, it seems that we are probing a wide range of sizes in the birth ring, but as we probe farther and farther out in the extended halo, the width of the size distribution decreases while the minimum size remains (relatively) constant, the expected behavior of radiation pressure.

## 5. Toward a self-consistent model

To provide additional context to our findings, we describe here an attempt at a more rigorous modeling approach. Since we are able to measure the degree of polarization in the extended halo for the disk around HD 129590, we use this star as a reference. Thebault & Kral (2019) showed that the contribution of grains that are on unbound orbits (*e* ≥ 1 because of radiation pressure) can be significant. Since our results put very stringent constraints on the

<footnote>$^6$ It should be noted that the declination of −59° of HD 115600 makes it a challenge to significantly increase the parallactic rotation with typical 1–1.5 h observing blocks.</footnote>
minimum and maximum grain sizes, our goal here is to assess
if and how much these unbound grains can affect the degree
of polarization. In the following, we first give a brief descrip-
tion of how synthetic images are computed, before describing
the starting hypotheses for four fiducial cases.

**5.1. Model description**

The approach follows the methodology of the DyCoSS code (see Thébault 2012), which is an iterative one. In summary, we first draw $n_{\text{part}}$ particles, following a size distribution $dn(s) \propto s^{-3.5}ds$, between 0.01 µm and 1 mm. As will be further discussed in Sect. 5.2, each particle is assigned a $\beta$ value. The central star has a mass $M_{\star} = 1.3 M_{\odot}$ (as reported in Kral et al. 2020 for HD 129590, the star chosen as an example in this section), and each particle feels a central mass of $M_{\star}(1-\beta)$. The $n_{\text{part}}$ are released all at once, with their semi-major axis drawn from a normal distribution, centered at 48 au with a standard deviation of 1.5 au (at a distance of 136.32 pc this is meant to mimic the birth ring of the disk around HD 129590). The opening angle of the disk is set to be constant with $h/r = 0.03$, and the mean anomaly is uniformly drawn between $-\pi$ and $\pi$. We then run a first N-body simulation, with a fine timestep (to ensure that we actually register the unbound grains), and at each timestep we save the position of all particles. We let the simulation evolve for $n_{\text{iter}}$ timesteps and eventually stop the simulation. From the positions of the particles, and their sizes, we can compute an optical depth radial profile, which maximum$^7$ is normalized to $7 \times 10^{-3}$ (the fractional luminosity of the disk, Esposito et al. 2020). This first simulation cannot be used as is, since the result will depend on the total duration of the N-body integration. Therefore, we then iterate and run another simulation, with the same timestep, and saving the positions of the particles at each timestep. The main difference being that collisions are now accounted for (making use of the optical depth profile from the previous simulation), and the simulation will last until 99.9% of the initial $n_{\text{part}}$ have been destroyed (more details can be found in Olofsson et al. 2022b). Since the unbound grains will never come back inside the birth ring, we set another condition for them to be destroyed. If they reach distances larger than 5 000 au, they are removed from the simulation. Once 99.9% of the grains have been destroyed, we re-evaluate the optical depth profile, normalize it, and we can iterate once more, until two consecutive simulations have converged (in that case we check that the optical depth profiles are similar). The implicit assumption for this approach to be valid is that the disk is in a steady state, meaning that it will not become brighter or fainter over time.

Once we have iterated a sufficient number of times, we can
produce images in scattered light. With the (x, y, z) positions
of the particles, saved at all timesteps, we can rotate them to
account for the inclination and position angle of the disk (i = 82°
and φ = -60.6°, see Table B.1), and calculate the value of the
scattering angle. We used optool to compute the polarized and
total intensity phase functions, at the wavelength of the obser-
vations, and assumed n = 4.25 and k = 2 × 10⁻³ (as derived for
the birth ring of HD 129590, Table 2, with a porosity of 25% and
using fmax = 0.8). The contribution of each particle is multiplied
by its cross-section s², the value of the phase function (either in
polarized or total intensity), and divided by the squared distance

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Radiation pressure $\beta$ ratio as a function of grain size $s$ for four different simulations</figcaption></figure>

to the central star to account for illumination effects. The two
images in polarized and total intensity are then convolved with
the telescope point-spread function from the observations, and
are then used to compute the degree of polarization image.

**5.2. Fiducial cases**

It should first be noted that running one simulation can take
up to several hours, producing very large files (tens of GB,
depending on the integration timestep). This approach is there-
fore (presently) not very suitable for directly fitting observations;
hence the choice of very specific cases in this section.

To better illustrate the impact of $\beta$ as a function of the size $s$, we investigate four fiducial cases using the different $\beta(s)$ curves shown in Fig. 6 (all the other parameters remaining the same otherwise). The motivation is to use some of the results we presented in Sect. 4 as inputs for the models, compute synthetic images, and qualitatively compare them to the original observations. For Run 1, there is a clear cut-off at the size $s_{\text{blow-out}}$ ($\le 0.5$ µm, see Sect. 6.1 for further discussion) and all the grains smaller than this size are set on unbound orbits. This blow-out size was chosen to be quite close to the minimum grain size derived for the birth ring of HD 129590 (0.35 µm) For Run 2, $\beta(s)$ crosses the threshold $\beta = 0.5$ twice and therefore the smallest grains are bound to the star, but there is a small interval of sizes for which the particles will be unbound. For Run 3, none of the grains are unbound, but $\beta$ reaches a maximum value of 0.47, and finally, Run 4 is quite similar to Run 1, except that the blow-out size is larger, with $s_{\text{blow-out}} \sim 1.85$ µm.

**5.3. Diagnostics and run comparison**

It is first interesting to check a few diagnostics plots to bet-
ter understand the differences between the three runs. The left
panels of Fig. 7 shows the final optical depth profiles as a func-
tion of the stlocentric distances, while the right panels show
the cross-section as a function of the grain size s, estimated
in four concentric annuli (not normalized with respect to each
other). The regions are at the same stlocentric distances as the
regions used for the analysis of HD 129590, assuming a distance
of 136.32 pc.

Starting with Run 1, the optical depth beyond the main
ring (centered at 48 au) follows a power-law in -1.5, which is

<footnote n="7">⁷ A more exact approach would be to compute the resulting spectral energy distribution and normalize the optical depth over the range of distances that contribute most to it, but the approach used here is sufficient for first-order estimates.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Diagnostics for Runs 1 to 4, shown from top to bottom, respectively. **Left:** optical depth profile as a function of stellocentric distance. The dashed line shows a profile in $r^{-1.5}$. For the bottom panel, the dashed-dotted line follows a profile in $r^{-1}$. **Right:** cross-section as a function of grain size $s$. The different colors show the cross-section measured in several concentric rings.</figcaption></figure>

expected for the steady-state evolution of a debris disk (Thébault & Wu 2008). Regarding the contribution of different grain sizes to the geometrical cross-section for Run 1, the blow-out size is clearly identifiable by the sharp drop at around 0.4 µm, visible for all four regions. For sizes that are smaller than this critical size, the cross-section increases again for smaller $s$, as the size distribution of the released particles in $\propto s^{-3.5}$ is very top-heavy. It is interesting to note that the farther the region considered is, the narrower the distribution is. This is the expected behavior for radiation pressure; particles with a given $\beta$ value, released from a distance $a_0$ can only reach separations up to $\sim a_0/(1 - 2\beta)$.

For Run 2, for which there is only a narrow range of sizes with $\beta > 0.5$, the optical depth profile differs from the one of Run 1. There is a significant “bump” at $\sim 200$ au. This bump is caused by the small-end part of the size distribution. Grains with $s \lesssim 0.03$ µm all have similar $\beta$ values close to 0.4 (the plateau on the left side of Fig. 6). The bump shows the location of the apocenters of these, very numerous and long-lived (Thébault & Wu 2008) particles that all have similar $\beta$ value ($a_0/(1 - 2\beta) \sim 48/(1 - 2 \times 0.4) = 240$ au). Despite this local increase in optical depth, the profile otherwise follows a power-law compatible with $r^{-1.5}$. Looking at the size-dependent cross-section plot, the

range of sizes for which the grains are unbound is easily identified by the strong dip between 0.1 and 0.3 µm. The decrease of the maximum size as a function of the stellocentric distance is quite comparable to the first Run.

For Run 3, there are no unbound grains, and for the bound grains, the $\beta$ value reaches a maximum of 0.47. This explains the sharp drop in optical depth at $\sim 800$ au (indeed $a_0/(1 - 2\beta) = 48/(1 - 2 \times 0.47) = 800$ au). The “bump” at $\sim 120$ au is caused by the plateau of $\beta(s)$ for very small particle sizes. The bump is closer to the star, because the $\beta$ values flatten at $\beta \sim 0.3$ in this case (compared to 0.4 for Run 2). For the cross-section, since there are no unbound particles, there is no discontinuity on the small end of the size distribution. We still see the expected behavior of radiation pressure, as we probe farther and farther out of the birth ring, the maximum grain size decreases.

Run 4 being in nature very similar to Run 1, the differences between the two runs are quite minute. The most notable difference, besides the size for which the cross-section is dropping, concerns the optical depth profile at large separation from the star. We can see that the slope deviates from the $r^{-1.5}$ profile and flattens. This is because the unbound grains are much more numerous in this example ($s_{min}$ remaining the same but $s_{blow-out}$
being larger), and their contribution rather follows a slope in $r^{-1}$, meaning that their relative contribution to the total flux (compared to bound grains) increases with the distance to the star (e.g., Fig. 2 of Thebault et al. 2023). This can best be seen on the lower left panel of Fig. 7, showing that the contribution of the unbound grains takes over the one from bound grains at a separation of 600–700 au.

## 5.4. Synthetic images

Figure 8 shows the observations for HD 129590 at the top, and the images for Runs 1 to 3. From left to right, we show the polarized intensity, total intensity, and degree of polarization. We can likely rule out the presence of bound very small particles in the disk around HD 129590. Indeed for Runs 2 and 3, the major axis of the disk is too bright in polarimetry, and the back side of the disk is too bright in total intensity. The degree of polarization also seems to be too large compared to the observations. This is because very small dust particles have a polarized phase function peaking at 90° scattering angles, an almost isotropic total intensity phase function, and in general a large degree of polarization.

Conversely, the results of Run 1 are a validation of our modeling approach presented in Sect. 3. By using the results from the modeling of the birth ring of HD 129590 (optical constants *n* and *k*) and using a $\beta(s)$ function that yields a subµm blow-out size below which all grains are unbound, we obtain a good match to the observations not only in the birth ring but also in the extended halo. Indeed, in the polarimetric images, the back side of the disk is equally faint as in the observations, while the front side appears with similar brightness. In total intensity we are able to reproduce the strong forward scattering along the minor axis, as well as the apparent bulge associated with it. The most notable difference between the observations and the model being that the halo might be slightly fainter in the model compared to the observations, and the arc reported in Olofsson et al. (2023) is not as visible in total intensity. This could be due to our parametrization of the birth ring, using a normal distribution. Using an asymmetric distribution that extends farther out might help with the transition between the birth ring and the halo. Even though we are not fitting the observations, the degree of polarization seems to agree quite well with the observations. Its value increases as a function of stlocentric distances, though not as much as in the observations, only from ~20% in the birth ring, up to ~30% in the outer regions. A possible explanation could be that the blow-out size for Run 1 is close to 0.5 µm (Fig. 7), while we found a minimum grain size of 0.35 µm when fitting the observations (Table 2).

## 5.5. Unbound grains, minimum grain size, and blow-out size

Until now we only investigated the effect the shape of the $\beta(s)$ curve can have on the final images in polarized and total intensity as well as for the degree of polarization. If one wished to ignore the contribution of unbound grains, there are simpler approaches to compute images (e.g., Lee & Chiang 2016; Olofsson et al. 2022a). It is therefore worth investigating how much unbound grains impact the final scattered light images and degree of polarization (see also Thebault & Kral 2019 for a more in-depth analysis). For Run 1 we recomputed scattered light images, but this time ignoring any particle with $\beta \ge 0.5$. The resulting images are shown in the top panels of Fig. 9, along side the previous images for Run 1 (including the contribution of unbound grains) in the middle panels.

At first glance, the differences seem to be quite minute, but upon closer inspection we can see a negative branch in the polarimetric image, and that the front side of the disk appears bulkier compared to the original Run 1. The biggest difference lies in the degree of polarization. It appears more structured, because of the negative polarization branch, but most importantly, the maximum degree of polarization does not vary with increasing stlocentric distances and remains overall quite low (~15%, see later).

We can test the hypothesis mentioned in Sect. 3, where we equated the minimum size to the blow-out size. We just discussed the importance of unbound grains and it poses the question of whether their contribution might “contaminate” our estimation of $s_{min}$ when modeling the observations as we did previously. Run 4 was computed to begin addressing this point, as the blow-out size is larger (1.85 µm versus 0.5 µm), but the contribution of unbound grains is included. The bottom panels of Fig. 9 show that the total intensity image is quite comparable to the one of Run 1 but the degree of polarization shows slightly larger values (~30% in the birth ring). A positive note is that in this case as well the degree of polarization is increasing with the stlocentric distance, in line with the observations. But the polarimetric image appears quite different than the observed one: there is a strong negative polarization branch for scattering angles close to ~50°. Since the image shows the absolute value of the pixels, this causes this apparent “null branch” along the minor axis of the disk, which is not seen in the observations.

This is further quantified in Fig. 10, showing the degree of polarization as a function of the scattering angle, in a similar way as it was measured for the observations, in four different concentric regions. The black dotted line shows the results for Run 1 without the contribution of unbound grains, the red solid line corresponds to Run 1 with the unbound grains, and the dashed blue line is for Run 4, including the unbound grains as well. For the birth ring (leftmost panel), the best agreement is found for Run 1 without the unbound particles. This is because this is the closest match to the fitting results obtained in Sect. 4. But as soon as we start venturing outside of the birth ring, the maximum degree of polarization remains constant and fails to reproduce the observations. When the contribution of unbound grains is included, we see that the degree of polarization is slightly too large in the birth ring (though not in total disagreement). However, as we go farther and farther, the degree of polarization continuously increases with separation. This is even more so the case for Run 4, but in all four regions, the degree of polarization is too large compared to the observations. Overall, this underlines the importance of including unbound grains to best reproduce the observations.

The main findings of this section can be summarized as follows: first, our results greatly highlight the importance of not only obtaining observations in polarized and total intensity but also determining the degree of polarization in debris disks, since this latter can help in disentangling otherwise relatively similar images. Second, this also shows that despite their smaller contribution to the cross section (right side of Fig. 7), the contribution of unbound grains does help reproduce the observations and should therefore not be disregarded, which strengthens the results of Thebault & Kral (2019). As discussed in Thebault et al. (2023) the relative contribution of unbound grains increases with the stlocentric distance, and we hypothesize that this is the main driver for the degree of polarization to increase with stlocentric distance. As we probe farther and farther out, the average grain size is effectively decreasing. Indeed, small unbound grains are contributing more and more, while the
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison of the observations (top panels) with the first three models that include the contribution of unbound grains. The left and middle panels show images in polarized and total intensity, respectively, while the right panels show the degree of polarization.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Images for Run 1 without (top) or with (middle) the contribution of unbound grains. The bottom images show the results for Run 4 for which the blow-out size is larger (the contribution of unbound grains is included).</figcaption></figure>

width of the size distribution becomes narrower due to radia-
tion pressure. Since smaller particles are expected to have larger
degrees of polarization this provides a natural explanation for
the increase in the degree of polarization with stlocentric dis-
tance. Nonetheless, considering only the increase of the degree
of polarization is not sufficient to fully explain the observations
if its absolute value is not correct (e.g., ~15% versus ~30% in
the birth ring). This exercise shows that indeed the contribution
of unbound grains does contaminate slightly the fitting approach
described in Sect. 4 (Run 1 with or without unbound grains, left-
most panel of Fig. 10), but that (i) this contamination is relatively
marginal in the birth ring and (ii) that unbound grains must be
included to reproduce the increasing degree of polarization in

the extended halo. Overall, the minimum grain size that we infer
when modeling the observations has to be close to the blow-out
size.

5.6. *Future perspectives*

In its current implementation, the approach described in this
Sect. cannot easily be used to perform complex modeling
given the time it takes to obtain the final images. That being
said, given its overall simplicity, there are reasons to be opti-
mistic, especially when considering additional observations that
could help decrease the number of free parameters. There
are already several ALMA datasets at high angular resolution
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Degree of polarization as a function of the scattering angle for four concentric regions. Observations for HD 129590 are shown in orange circles. The different lines correspond to the degree of polarization measured from the modeled images for Run 1 (with and without the contribution of unbound grains) and Run 4.</figcaption></figure>

(e.g., Marino et al. 2016 for HD 181327; Kennedy et al. 2018 for HR 4796; Daley et al. 2019 and Vizgan et al. 2022 for AU Mic), and more are coming with the ALMA ARKS large program. From these high angular resolution, one can place stringent constraints on some geometric parameters (inclination, position angle, and quite possibly the opening angle of the disk) as well as on the surface density profile of the birth ring. The remaining free parameters should be related to the dust optical properties, either a composition (or mix of compositions), or in the “worst case scenario”, values for *n*, *k*, and a curve to relate *s* to *β*. This could open the possibility of performing at least a coarse exploration of the parameter space and provide novel constraints on the dynamics of the dust particles.

# 6. Discussion

## 6.1. Radiation pressure and blow-out size

Even though the results of Run 1 seem to agree with the observations of HD 129590, we must acknowledge possible pitfalls or limitations of the approach. Some of the modeling results obtained in Sect. 4 might seem challenging or unrealistic, especially when it comes to the inferred blow-out size (the optical constants will be discussed in Sect. 6.2). Leaving aside the results for HD 115600 (since the projected minor axis of the disk is not well detected), we have one F- and one G-type star for which we derive $s_{min} \le 0.1$ µm and $s_{min} = 0.36$ µm, respectively. These minimum sizes appear to be very small and, if representative of the blow-out sizes, would severely challenge our understanding of the strength of radiation pressure.

We here argue that this apparent challenge is not a real issue and that it is most likely due to limitations of available light scattering models that can properly describe the morphology of the dust particles. We want to emphasize that from a dynamical point of view what really matters is the distribution of $\beta$ values. This is the quantity that governs the final dust density distribution as a function of the stellocentric distance. Associating a given size $s$ to a given $\beta$ value is, in the end, nothing more than a “necessary evil” to account for the phase function (even the cross-section could in the end be approximated by $1/\beta^2$ for sizes large enough).

For the remainder of this discussion, we focus on the case of HD 129590, since we are able to measure the degree of polarization in the halo, and were able to propose a self-consistent model that explains the observations reasonably well. The lesson learned from the exercise presented in Sect. 5 is that if we have a good description of the dust properties in the birth ring, then we are able to reproduce the observations quite well (even

beyond the birth ring). The only caveat being that we need a $\beta(s)$ curve in agreement with the inferred dust properties; in other words, a blow-out size close to the minimum size derived from fitting the phase function and degree or polarization (as discussed previously in Sect. 5).

The apparent challenge of finding very small minimum grain sizes therefore all boils down to the choice of the light scattering model used to fit the observed profiles as a function of the scattering angles. As discussed in Sect. 3.2, we tried several approaches, either trying mixtures of known materials, or pre-computed optical properties of complex geometries. We could only obtain a good match to the observations using the DHS model, combined with rather unusual optical constants (see Sect. 6.2). We were not able to use the library of phase functions provided in the AggScatVIR library mostly because they failed to reproduce the low degree of polarization observed in the birth rings of the disks. We hypothesized that this is because the largest particle size in the library only goes as far as a few µm while we should expect to probe a much wider size distribution in the birth ring. It would be very informative to test how the optical properties integrated over a wider size distribution, including aggregates at least for the small-end of the distribution, would look like and what the inferred minimum grain size would be in such case. Unfortunately, this cannot be tested at the moment.

Still, aggregate grains remain an interesting avenue to further explore. For instance, Arnold et al. (2022) obtained promising results when modeling observations of AU Mic using agglomerated debris particles. Despite the challenging edge-on configuration of the debris disk, they show that using aggregates does help reproducing the degree of polarization along the major axis of the disk. Furthermore, Tazaki & Dominik (2022) showed for instance that some aggregates can yield maximum degree of polarization as low as ~40% at a wavelength of 1.6 µm. Porous aggregates, even with an equivalent radius of 2 µm can have relatively low maximum degree of polarization if the size of the individual monomer is sufficiently large (0.4 µm). Interestingly, this is quite close to the minimum grain size that we infer for the birth ring of HD 129590, and as the authors mention, the polarization properties of the aggregate as a whole are more correlated to the properties of the monomers themselves, because of the additional internal scattering events that can happen inside the aggregate itself. This was also discussed in Min et al. (2016), who argued that polarization depends on the size of the monomers, while total intensity observations are more sensitive to the size of the aggregates themselves.
We therefore suggest that our observations could be mostly sensitive to the size of the individual monomers. Because we are using the DHS model, we cannot at the same time have the optical properties of a much larger aggregate particle, leading us to actually under-estimate the minimum grain size in the distribution. To some extent, this is further supported by the fits to the total intensity phase functions, where the forward scattering peak is in some cases under-estimated. Larger aggregates should in principle lead to stronger forward scattering peaks, while preserving the polarized flux, in better agreement with the observations. If correct, this means that the minimum grain size we obtain is in fact a lower limit for the sizes of the particles.

## 6.2. Optical constants

Our modeling results of the degree of polarization and total intensity phase function from Sect. 4 suggest relatively large *n* and very large *k* values for the particles optical constants. Other studies have already reported similar issues in the past, for instance Duchêne et al. (2020) for the disk around HD 32297 (*n* ~ 3.8, log₁₀*k* ~ -1.4), Arriaga et al. (2020) for HR 4796 (*n* = 3.4 and *k* = 3.7 when fitting jointly polarized and total intensity observations), or Milli et al. (2024) for HD 181327 (*n* = 3.4, *k* = 1), all of these studies having used DHS (and also compared with the Mie theory as well). All three studies summarized their findings by comparing the derived optical constants with those from known materials (e.g., astronomical silicates, iron, carbon, among others), highlighting that the inferred real and imaginary parts are quite unusual, and we refer the interested reader to these works for further information.

Instead, we focus here on the possible limitations of the light scattering theory that we used to model our observations, to further pursue the discussion initiated in Sect. 6.1. Muñoz et al. (2021) derived the total intensity phase function and degree of polarization for several samples of forsterite particles, for various narrow size distributions. Their sample “XS” has a size distribution between 0.1 and 1.0 µm, while the “XL” sample has a size distribution between ~20 and 100 µm. They also present the results for three other samples that have intermediate distributions between these two (“S”, “M”, and “L”, respectively). They derive maximum degree of polarization smaller than 20% for all five samples, but also see a strong negative branch for the samples “XS” and “S” that is not seen for the other three samples. The total intensity phase functions show strong forward scattering peaks, for all five samples, as well as some backward scattering. Overall, their results for samples “M”, “L”, or even “XL” might compare well (visually at least) with our results obtained for HD 129590 (for which we have the most reliable measurements and the widest range of scattering angles), these samples having size distributions in the range ~1–10 or 20–100 µm.

Of greater interest for this discussion, Muñoz et al. (2021) also investigated the dust properties that would be inferred by modeling their laboratory measurements with the Mie theory, in other words, performing the same exercise we did in Sect. 4, but with full knowledge of the ground truth. They show that to be able to reproduce the total intensity phase functions of their low-absorbing, irregular forsterite grains (*k* ~ 10⁻⁵), they need to artificially increase the absorptivity of the spherical grains and require *k* values of 10⁻²–10⁻¹. They performed a similar exercise on the degree of polarization as a function of the scattering angle. Only changing the imaginary part *k* was not sufficient to match the observations, and they therefore used a size distribution of particles. Their results indicate that the inferred

grain sizes are severely under-estimated when assuming spherical grains. They also find that the imaginary part of the optical constants are over-estimated, as well as the real part *n*. Overall, even though they used the Mie theory, implying a spherical geometry for the particles, while we use the DHS model more akin to irregular particles, this echoes resoundingly our findings of (very) large (*n*,*k*) values, and small minimum sizes. Furthermore, Muñoz et al. (2021) did not attempt to simultaneously fit the degree and total intensity phase functions as we attempted here.

## 6.3. Future perspectives

The bottleneck clearly seems to be related to how the optical properties of dust particles are computed (the phase functions, degree of polarization but also the variation of *β* as a function of *s*). Despite the slightly pessimistic tone of the discussion in Sects. 6.1 and 6.2, we show in Sect. 5 that if we can find a reasonable match to the optical properties in the birth ring, regardless how “unrealistic” they might seem, then the rest of the pieces fall in place together. With the adequate *β*(*s*) function, and accounting for the contribution of the unbound particles, we are able to reproduce the observations on a wide range of spatial scales. Even if the absolute values remain highly uncertain, the relative values for *n* and *k* for instance can still be used to better understand the kind of cosmic dust debris disks are made of. It is indeed quite interesting that we find similar values as for the disk around HD 32297, HR 4796, and HD 181327, most of these studies being performed using different models, assumptions, and approaches. Even though we are clearly in the small statistics regime, regardless of the exact composition of the dust particles, it seems to be relatively similar around different stars, of different spectral types. Constraining the degree of polarization (and not only the polarized or total intensity phase functions) for a larger sample of young debris disks is therefore crucial to better understand the building blocks of planetesimals.

Additionally, one could try to step away from depending on light scattering models when trying to reproduce observations. Recently, both Lawson et al. (2021) and Hom et al. (2024) took steps in this direction, as both studies fixed the total intensity phase function, based on measurements of solar system bodies or other debris disks. However, the phase function they used still remains independent of the grain size in both studies and it would be the same in the birth ring and halo of the disk. Nonetheless, we can still think of other approaches to depend less on the most commonly used light scattering models. The first one would be to pursue the work on optical properties of aggregates to increase the maximum size of the particles, possibly setting some limits on the fractal dimension or “fluffiness” of the largest grains. Another approach would be to directly make use of laboratory measurements such as the ones presented in Muñoz et al. (2021). Granted, there is “only” one composition analyzed, but there are several samples – of a wide range of sizes – that are almost fully characterized, the main ingredient missing being the radiation pressure efficiency in order to compute *β*(*s*). Otherwise, one could interpolate the phase function and degree of polarization over the different sizes, and such a model would provide stringent constraint on the strength of radiation pressure, *β*(*s*) being one of the few free parameters of the model (assuming the composition is indeed representative of dust in debris disks).

Another approach, which needs maturing, would be to rely mostly on the observations. With high angular resolution ALMA observations, the geometric parameters of the birth ring can be derived quite accurately. Using this prior knowledge it might be possible to work our way back to the optical properties, in a data
driven way. First the degree of polarization and polarimetric (or total intensity) images are highly complementary. Indeed, the former does not depend on the dust density distribution, meaning that we have direct access to some of the optical properties, which are not entirely disconnected from what is measured in the polarimetric (or total intensity) images. Furthermore, using cross-section plots such as the one presented in the right panels of Fig. 7, it might be possible to retrieve the dust properties outside-in. Indeed, the outermost regions probe a relatively narrow range of grain sizes. As we probe closer in, the maximum grain size of the distribution will increase. In Olofsson et al. (2020), we presented another data-driven approach to modeling polarimetric observations in which the polarized phase function is an output of the modeling, and it might be possible to expand on a similar approach, informed by the dynamics imposed by radiation pressure.

## 7. Summary

In this paper, we present new SPHERE observations of four debris disks, and attempt to measure the degree of polarization. We show that this remains a challenge, despite significant progress in instrumentation and post-processing techniques; HD 191089 is not detected in total intensity, and the results for HD 115600 cannot be exploited to fully constrain the dependence with the scattering angle (as the minor axis is not well recovered). Nonetheless, we were able to measure the degree of polarization for the birth ring of the disk around HD 157857, and more significantly, we measured the degree of polarization in the birth ring as well as in the extended halo of the disk surrounding HD 129590.

For the three disks for which we are able to determine the degree of polarization in the birth ring, we find small a maximum degree of polarization, of namely $≤40\%$ (with some dispersion for HD 157857), which remains challenging to model assuming compact spherical grains. Reproducing the observations using a precomputed library of optical properties of aggregates was also found to be challenging, which we interpreted as due to the contribution of larger particles, as expected in the birth rings of debris disks. We were only able to find a satisfactory match to the observations by leaving the optical constants ($n, k$) as free parameters. Nonetheless, for the disk around HD 129590, we are also able to constrain the degree of polarization in the halo, and fitting the data independently, we do find that the results are relatively consistent with each other (except for the maximum grain size), which gives us confidence in our approach. For this disk, and for the regions outside of the birth ring in particular, we find that the maximum grain size is much smaller compared to the inner regions, which is naturally explained by the stellar radiation pressure; only small particles are set on eccentric orbits and can venture outside, where they are released.

We tested our findings by using them as inputs for a self-consistent model, including the contribution of unbound particles, and present our results in Sect. 5. The main conclusion of this exercise is that we can reproduce the observations relatively well only if unbound grains are included in the final images and if the blow-out size is close to the minimum grain size inferred from the modeling. The driving argument to include the contribution of unbound particles is to reproduce the increase in the degree of polarization with stellocentric distances, which we argue is the consequence of a decreasing average effective size with increasing separation.

Our modeling results presented in Sect. 4 might appear unusual in several respects; in particular the inferred minimum

grain size, as well as the optical constants. As discussed at length in Sect. 6, the root cause of these (apparent) inconsistencies is most likely related to the light scattering model that we use (which is the only one that can provide a good fit to the data). In spite of the challenges underlined in Sect. 6, we discuss possible avenues to circumvent some of these issues and further our understanding of cosmic dust in young debris disks. The new approach presented in Sect. 5 to produce synthetic NIR scattered-light images requires very few free parameters, as the dust density distribution in the halo is solely governed by the effect of radiation pressure. We show that if we have a suitable solution for the dust properties in the birth ring, we are then able to (at least qualitatively) reproduce the extended halo beyond the birth ring. The power of this approach will be greatly enhanced by the high-angular-resolution ALMA observations in the near future.

**Acknowledgements.** We thank the anonymous referee for a constructive, timely, and helpful report. The comments helped strengthen the results of this paper, especially (but not only) with respect to self-subtraction effects, as well as improving the overall clarity of the paper. We are thankful to Ryo Tazaki for pointing out interesting references that helped shape the discussion of this paper. J.M. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon Europe research and innovation program (grant agreement No. 101053020, project Dust2Planets, PI F. Ménard). SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF-Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RI3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). We also acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU in France. This work has also been supported by a grant from the French Labex OSUG@2020 (Investissements d'avenir – ANR10 LABX56). The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). It has also been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). M.R.M., H.M.S., and S.D. are pleased to acknowledge this financial support of the SNSF. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Astropy<fnref n="8" />, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020), and Numba (Lam et al. 2015).

## References

Arnold, J. A., Weinberger, A. J., Videen, G., & Zubko, E. S. 2022, ApJ, 930, 123
Arriaga, P., Fitzgerald, M. P., Duchêne, G., et al. 2020, AJ, 160, 79
Asensio-Torres, R., Janson, M., Hashimoto, J., et al. 2016, A&A, 593, A73
Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33
Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123
Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155
Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
Crotts, K. A., Matthews, B. C., Duchêne, G., et al. 2024, ApJ, 961, 245
Currie, T., Lisse, C. M., Kuchner, M., et al. 2015, ApJ, 807, L7
Daley, C., Hughes, A. M., Carter, E. S., et al. 2019, ApJ, 875, 87
de Boer, J., Langlois, M., van Holstein, R. G., et al. 2020, A&A, 633, A63
Dohlen, K., Saisse, M., Origne, A., et al. 2008, SPIE Conf. Ser., 7018, 701859

<footnote n="8">⁸ http://www.astropy.org</footnote>
Dominik, C., Min, M., & Tazaki, R. 2021, OpTool: Command-line driven tool for creating complex dust opacities, Astrophysics Source Code Library [record ascl:2104.010]

Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503

Draper, Z. H., Duchêne, G., Millar-Blanchaer, M. A., et al. 2016, ApJ, 826, 147

Duchêne, G., Rice, M., Hom, J., et al. 2020, AJ, 159, 251

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

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

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

Esposito, T. M., Kalas, P., Fitzgerald, M. P., et al. 2020, AJ, 160, 24

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Flasseur, O., Thé, S., Denis, L., Thiébaut, É., & Langlois, M. 2021, A&A, 651, A62

Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Gibbs, A., Wagner, K., Apai, D., et al. 2019, AJ, 157, 39

Gledhill, T. M., Scarrott, S. M., & Wolstencroft, R. D. 1991, MNRAS, 252, 50

Graham, J. R., Kalas, P. G., & Matthews, B. C. 2007, ApJ, 654, 595

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hom, J., Patience, J., Chen, C. H., et al. 2024, MNRAS, 528, 6959

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

Juillard, S., Christiaens, V., & Absil, O. 2023, A&A, 679, A52

Kennedy, G. M., Marino, S., Matrà, L., et al. 2018, MNRAS, 475, 4924

Kral, Q., Matrà, L., Kennedy, G. M., Marino, S., & Wyatt, M. C. 2020, MNRAS, 497, 2811

Krivov, A. 2010, Research in Astronomy and Astrophysics, 10, 383

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

Lam, S. K., Pitrou, A., & Seibert, S. 2015, in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM 15 (New York, NY, USA: Association for Computing Machinery)

Lawson, K., Currie, T., Wisniewski, J. P., et al. 2021, AJ, 162, 293

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

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

Matthews, E., Hinkley, S., Vigan, A., et al. 2017, ApJ, 843, L12

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

Millar-Blanchaer, M. A., Wang, J. J., Kalas, P., et al. 2016, AJ, 152, 128

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

Milli, J., Choquet, E., Tazaki, R., et al. 2024, A&A, 683, A22

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

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

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

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

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

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

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

Olofsson, J., Thébault, P., Bayo, A., et al. 2023, A&A, 674, A84

Pairet, B., Cantalloube, F., & Jacques, L. 2021, MNRAS, 503, 3724

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

Ren, B., Pueyo, L., Chen, C., et al. 2020, ApJ, 892, 74

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

Tamura, M., Fukagawa, M., Kimura, H., et al. 2006, ApJ, 641, 1172

Tazaki, R., & Dominik, C. 2022, A&A, 663, A57

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

Thalmann, C., Janson, M., Buenzli, E., et al. 2013, ApJ, 763, L29

Thébault, P. 2012, A&A, 537, A65

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

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

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

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

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

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Vizgan, D., Hughes, A. M., Carter, E. S., et al. 2022, ApJ, 935, 131

Wahhaj, Z., Milli, J., Romero, C., et al. 2021, A&A, 648, A26

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

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

Xie, C., Choquet, E., Vigan, A., et al. 2022, A&A, 666, A32

Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321
# Appendix A: Observations and data processing

## A.1. Observing log

Table A.1 summarizes all the observations obtained in Pro-grammes 105.20GP.001 and 109.237K.001.

## A.2. Gallery of final products

Figure A.1 shows a gallery of all the reduced observa-tions obtained within Programme IDs 105.20GP.001 and 109.237K.001. The left half shows observations in total intensity and the right half shows the polarimetric $Q_{\phi}$ images.

## A.3. Assessment of self-subtraction effects

Self-subtraction remains a significant challenge for many post-processing algorithms and can prevent us from properly measuring the total intensity phase function of the disks, as well as their morphology (e.g., Milli et al. 2012). Since the degree of polarization depends on the total intensity images, we here want to estimate if they are severely impacted by self-subtraction effects. Instead of relying on a model-dependent approach, our goal is to check whether our final total intensity images are an accurate representation of the disks. We therefore duplicated the original SPHERE data cube and for each frame, we subtracted the final total intensity image to it (rotated by the corresponding parallactic angle). We then run the DI-sNMF pipeline on this newly created cube, using the same reference star to build the components, and compute a residual map. This is equivalent to performing forward modeling of the disk, an approach that has been routinely used in the past decade. This exercise is then repeated twice more, with a disk image that is 5% and 10% fainter to estimate how stringent the constraints are. The results are presented in Fig. A.2. While there is still some residual signal close to the coronagraph, the disks are overall very well removed. For all three stars we can see some signal along the trace of the disks on the rightmost panels (10% flux decrease). For HD 129590 and HD 115600 this is also the case with a 5% decrease in flux, while for HD 157857 it is more difficult to assert whether we start seeing some signal coming from the disk. In all cases, when subtracting the unaltered $I_{\text{disk}}$ image to the data-cube, there is no clear trace of any scattered light. This means that the final total intensity images (left panels) are indeed reliable representation of the surface brightness of the disks and that self-subtraction effects are well mitigated with this approach (at the < 5 – 10% level).
<table><caption>**Table A.1.** Log of the IRDIS BB_H pupil-tracking DPI observations.</caption><thead><tr><th>Star</th><th>Type</th><th>Date</th><th>DIT [s]</th><th>N<sub>f</sub></th><th>&lt;Airmass&gt;</th><th>&lt;Seeing&gt; ["]</th><th>&lt;τ&gt; [ms]</th><th>ΔPA [°]</th><th>OB grade</th><th>Used?</th></tr></thead><tbody><tr><td>HD 191089</td><td>SCI</td><td>2021-07-20</td><td>32</td><td>64</td><td>1.10</td><td>0.85</td><td>3.7</td><td>10.45</td><td>A &amp; B</td><td>✓</td></tr><tr><td>↔ HD 191131</td><td>CAL</td><td></td><td>32</td><td>24</td><td>1.10</td><td>0.86</td><td>3.7</td><td>10.64</td><td></td><td></td></tr><tr><td>HD 191089</td><td>DPI</td><td>2021-07-21</td><td>32</td><td>96</td><td>1.45</td><td>0.96</td><td>3.8</td><td>5.11</td><td>A</td><td>✓</td></tr><tr><td>HD 191089</td><td>SCI</td><td>2021-07-21</td><td>32</td><td>69</td><td>1.07</td><td>0.74</td><td>6.5</td><td>10.54</td><td>Problem</td><td>✗</td></tr><tr><td>↔ HD 191131</td><td>CAL</td><td></td><td>32</td><td>24</td><td>1.06</td><td>0.74</td><td>5.0</td><td>14.04</td><td></td><td></td></tr><tr><td>HD 191089</td><td>SCI</td><td>2021-09-04</td><td>32</td><td>64</td><td>1.09</td><td>0.53</td><td>4.5</td><td>11.71</td><td>B</td><td>✓</td></tr><tr><td>↔ HD 191131</td><td>CAL</td><td></td><td>32</td><td>24</td><td>1.09</td><td>0.58</td><td>5.2</td><td>11.78</td><td></td><td></td></tr><tr><td>HD 191089</td><td>SCI</td><td>2022-05-04</td><td>32</td><td>64</td><td>1.05</td><td>0.77</td><td>5.3</td><td>18.76</td><td>A</td><td>✓</td></tr><tr><td>↔ HD 191131</td><td>CAL</td><td></td><td>32</td><td>24</td><td>1.04</td><td>0.81</td><td>5.4</td><td>14.60</td><td></td><td>✓</td></tr><tr><td>HD 191089</td><td>DPI</td><td>2022-05-17</td><td>32</td><td>16</td><td>1.35</td><td>1.03</td><td>2.8</td><td>0.88</td><td>C</td><td>✗</td></tr><tr><td>HD 157587</td><td>DPI</td><td>2021-07-16</td><td>64</td><td>52</td><td>1.18</td><td>0.67</td><td>3.0</td><td>2.53</td><td>C</td><td>✓</td></tr><tr><td>HD 157587</td><td>DPI</td><td>2021-09-29</td><td>64</td><td>8</td><td>1.16</td><td>0.99</td><td>2.8</td><td>0.40</td><td>X</td><td>✗</td></tr><tr><td>HD 157587</td><td>SCI</td><td>2022-05-01</td><td>64</td><td>24</td><td>1.01</td><td>0.62</td><td>2.9</td><td>354.54</td><td>C</td><td>✗</td></tr><tr><td>↔ HD 158018</td><td>CAL</td><td></td><td>64</td><td>12</td><td>1.01</td><td>0.57</td><td>2.9</td><td>327.77</td><td></td><td>✗</td></tr><tr><td>HD 157587</td><td>SCI</td><td>2022-05-10</td><td>64</td><td>8</td><td>1.01</td><td>1.02</td><td>4.6</td><td>17.67</td><td>Incomplete</td><td>✗</td></tr><tr><td>↔ HD 158018</td><td>CAL</td><td></td><td>64</td><td>4</td><td>1.01</td><td>1.01</td><td>4.0</td><td>5.12</td><td></td><td>✗</td></tr><tr><td>HD 157587</td><td>DPI</td><td>2022-05-11</td><td>64</td><td>52</td><td>1.17</td><td>1.22</td><td>2.6</td><td>2.73</td><td>C</td><td>✓</td></tr><tr><td>HD 157587</td><td>DPI</td><td>2022-05-14</td><td>64</td><td>52</td><td>1.10</td><td>0.77</td><td>2.9</td><td>7.20</td><td>A</td><td>✓</td></tr><tr><td>HD 157587</td><td>SCI</td><td>2022-08-01</td><td>64</td><td>32</td><td>1.03</td><td>0.56</td><td>5.5</td><td>53.23</td><td>A</td><td>✓</td></tr><tr><td>↔ HD 158018</td><td>CAL</td><td></td><td>64</td><td>12</td><td>1.02</td><td>0.53</td><td>6.0</td><td>48.90</td><td></td><td>✓</td></tr><tr><td>HD 115600</td><td>SCI</td><td>2021-07-11</td><td>64</td><td>11</td><td>1.25</td><td>0.52</td><td>4.2</td><td>13.91</td><td>Problem</td><td>✗</td></tr><tr><td>HD 115600</td><td>DPI</td><td>2021-07-16</td><td>64</td><td>52</td><td>1.58</td><td>0.61</td><td>4.2</td><td>18.08</td><td>A</td><td>✓</td></tr><tr><td>HD 115600</td><td>SCI</td><td>2022-02-10</td><td>64</td><td>8</td><td>1.26</td><td>0.74</td><td>6.9</td><td>4.13</td><td>Problem</td><td>✗</td></tr><tr><td>HD 115600</td><td>SCI</td><td>2022-02-14</td><td>64</td><td>32</td><td>1.25</td><td>0.45</td><td>6.8</td><td>27.54</td><td>A</td><td>✓</td></tr><tr><td>↔ HD 117255</td><td>CAL</td><td></td><td>64</td><td>12</td><td>1.25</td><td>0.49</td><td>8.2</td><td>26.17</td><td></td><td>✓</td></tr><tr><td>HD 115600</td><td>SCI</td><td>2022-04-11</td><td>64</td><td>32</td><td>1.22</td><td>1.08</td><td>1.8</td><td>29.03</td><td>C</td><td>✓</td></tr><tr><td>↔ HD 117255</td><td>CAL</td><td></td><td>64</td><td>12</td><td>1.22</td><td>1.07</td><td>2.4</td><td>27.25</td><td></td><td></td></tr><tr><td>HD 115600</td><td>DPI</td><td>2022-04-17</td><td>64</td><td>64</td><td>1.66</td><td>0.60</td><td>10.2</td><td>21.11</td><td>A</td><td>✓</td></tr><tr><td>HD 115600</td><td>SCI</td><td>2022-05-01</td><td>64</td><td>32</td><td>1.24</td><td>0.73</td><td>4.0</td><td>30.34</td><td>A</td><td>✓</td></tr><tr><td>↔ HD 117255</td><td>CAL</td><td></td><td>64</td><td>12</td><td>1.24</td><td>0.75</td><td>4.4</td><td>25.93</td><td></td><td></td></tr><tr><td>HD 115600</td><td>SCI</td><td>2022-05-13</td><td>64</td><td>8</td><td>1.28</td><td>1.03</td><td>3.9</td><td>3.99</td><td>Incomplete</td><td>✗</td></tr><tr><td>HD 129590</td><td>DPI</td><td>2021-07-16</td><td>64</td><td>52</td><td>1.45</td><td>0.66</td><td>4.4</td><td>12.38</td><td>B</td><td>✓</td></tr><tr><td>HD 129590</td><td>SCI</td><td>2022-04-01</td><td>64</td><td>32</td><td>1.06</td><td>0.66</td><td>4.8</td><td>64.05</td><td>A</td><td>✓</td></tr><tr><td>↔ HD 129280</td><td>CAL</td><td></td><td>64</td><td>11</td><td>1.06</td><td>0.68</td><td>5.6</td><td>53.76</td><td></td><td>✓</td></tr><tr><td>HD 129590</td><td>DPI</td><td>2022-04-04</td><td>64</td><td>52</td><td>1.08</td><td>0.92</td><td>2.6</td><td>36.79</td><td>B</td><td>✓</td></tr><tr><td>HD 129590</td><td>SCI</td><td>2022-04-28</td><td>64</td><td>32</td><td>1.07</td><td>0.78</td><td>3.3</td><td>42.48</td><td>A</td><td>✓</td></tr><tr><td>↔ HD 129280</td><td>CAL</td><td></td><td>64</td><td>12</td><td>1.07</td><td>0.76</td><td>3.7</td><td>33.97</td><td></td><td></td></tr><tr><td>HD 129590</td><td>SCI</td><td>2022-05-14</td><td>64</td><td>34</td><td>1.04</td><td>0.62</td><td>4.5</td><td>64.01</td><td>A &amp; ?</td><td>✓</td></tr><tr><td>↔ HD 129280</td><td>CAL</td><td></td><td>64</td><td>12</td><td>1.04</td><td>0.65</td><td>4.5</td><td>52.38</td><td></td><td></td></tr></tbody></table>

**Notes.** We report the type of observations (SCI, CAL, or DPI), the observing date, the detector integration time (DIT), the number of total frames ($N_f$), the average airmass, seeing, coherence time ($\tau$), the range of parallactic angle during the sequence, the grades of the Observing Blocks, and whether the observations were used in this study. For clarity, star-hopping sequences are highlighted in light gray and the calibrator is marked with an arrow (↔).
<figure><img src="image_12.png" /><figcaption>**Fig. A.1.** Gallery of the final data products for all the epochs listed in Table A.1. The four columns on the left are for total intensity observations, the four on the right are for polarized intensity.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. A.2.** Assessing the strength of self-subtraction effects. From left to right: Total intensity image computed using DI-sNMF, which serve as input for the "forward modeling"-like exercise ($I_{\text{disk}}$). The next three panels show the results when first subtracting $I_{\text{disk}} \times 0.95$, and $I_{\text{disk}} \times 0.9$ to the original datacube and performing DI-sNMF. The color scale is the same for all panels of one row. From top to bottom we show the results for HD 157857, HD 115600, and HD 129590, respectively.</figcaption></figure>
## Appendix B: Geometric modeling of the polarimetric data

To determine the position angle and inclination of the disks studied here we modeled the $Q_\phi$ observations similarly to what was described in Olofsson et al. (2022b). As a matter of fact, for HD 129590 and HD 115600, we used the results presented in the aforementioned paper. As discussed in Olofsson et al. (2023), polarimetric observations are best suited to derive the morphological parameters of the disk compared to total intensity observations. Figure B.1 shows the observations, best fit results, and residuals, and Table B.1 shows the best fit parameters. The interested reader is referred to Olofsson et al. (2022b) for further detail on the modeling approach, but in short, the volumetric dust density distribution $N_{\text{dens}}$ of the disk follows

$$N_{\text{dens}}(r, z) \propto \left[ \left( \frac{r}{a_0} \right)^{-2\alpha_{\text{in}}} + \left( \frac{r}{a_0} \right)^{-2\alpha_{\text{out}}} \right]^{-1/2} \times \exp \left[ - \left( \frac{|z|}{\tan(\psi)r} \right)^{\gamma} \right], \quad (B.1)$$

where $r$ is the stellocentric distance, $z$ the height above the midplane, $a_0$ a reference radius, $\alpha_{\text{in}}$ and $\alpha_{\text{out}}$ two indices for the dust radial distribution. The vertical structure of the disk is parametrized by an opening angle $\psi$ and an exponent $\gamma$ to control the fall-off along the $z$ direction.

<table><caption>**Table B.1.** Results of the geometric modeling of the observations, ordered by increasing distance from earth.</caption><thead><tr><th>Star</th><th>d<sub>★</sub><br>[pc]</th><th>a<sub>0</sub><br>["]</th><th>i<br>[°]</th><th>φ<br>[°]</th><th>α<sub>in</sub></th><th>α<sub>out</sub></th><th>ψ<br>[10<sup>-3</sup> rad]</th><th>γ</th><th>χ<sup>2</sup><sub>r</sub></th></tr></thead><tbody><tr><td>HD 191089</td><td>50.11 ± 0.05</td><td>0.88 ± 0.02</td><td>61.3 ± 0.3</td><td>-108.0 ± 0.3</td><td>8.7 ± 0.7</td><td>-3.4 ± 0.2</td><td>46 ± 19</td><td>5.21 ± 2.74</td><td>0.23</td></tr><tr><td>HD 157587</td><td>99.87 ± 0.23</td><td>0.79 ± 0.03</td><td>70.1 ± 0.3</td><td>-50.7 ± 0.3</td><td>9.7 ± 0.3</td><td>-2.9 ± 0.1</td><td>15 ± 11</td><td>5.79 ± 2.93</td><td>0.25</td></tr><tr><td>HD 115600</td><td>109.04 ± 0.25</td><td>0.45 ± 0.02</td><td>76.6 ± 0.3</td><td>-156.0 ± 0.2</td><td>2.2 ± 0.2</td><td>-6.0 ± 0.2</td><td>132 ± 4</td><td>8.32 ± 1.52</td><td>1.04</td></tr><tr><td>HD 129590</td><td>136.32 ± 0.44</td><td>0.35 ± 0.02</td><td>82.0 ± 0.2</td><td>-60.6 ± 0.3</td><td>32.4 ± 2.5</td><td>-1.8 ± 0.1</td><td>1 ± 1</td><td>9.70 ± 0.23</td><td>1.27</td></tr></tbody></table>

**Notes.** The columns show the stellar names, the distance, the reference radius $a_0$, the inclination $i$, position angle $\phi$, inner and outer slopes of the density distribution ($\alpha_{\text{in}}$ and $\alpha_{\text{out}}$, respectively), opening angle $\psi$, the exponential fall-off $\gamma$, and the reduced $\chi^2$.
<figure><img src="image_14.png" /><figcaption>**Fig. B.1.** Modeling results for the four debris disks discussed in this study. The $Q_{\phi}$ image (the ellipse represents where the goodness of fit is estimated), best fit model, and residuals are shown from top to bottom. From left to right, the stars are HD 191089, HD 157587, HD 115600, and HD 129590 (same as in Olofsson et al. 2022b for the latter two stars). The scaling is linear and the same for each column (colorbar shown in the middle panel). The horizontal bar in the bottom panel represents 1".</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. B.2.** Stellocentric distance in arcsec and scattering angle in degrees (top and bottom, respectively) in the midplane, binned to the hexagonal grid, page 24 of 24</figcaption></figure>