© ESO 2020

1 Introduction

Dust grains are the building blocks of planets, but there are relatively few ways to accurately characterise their properties. Studies of solar system bodies provide the strongest constraints on the constituents of comets or asteroids (e.g. Frattin et al. 2019; Bertini et al. 2019), but they do not inform us directly on what is taking place during the planet formation stage. Observations of discs around young (≤100 Myr old) stars are sensitive to grain sizes that are smaller than typically 100 μm or a few millimetres. Debris discs, discs of second generation dust, are ideal targets to study dust grains. As those discs are optically thin at all wavelengths, we have direct access to the absorption and scattering properties of the grains, without having to account for non-trivial optical depth effects or multiple scattering events. There are two possible ways to characterise dust properties in debris discs, first, via their thermal emission by measuring the spectral slope at (sub-) millimetre wavelengths (Draine 2006; MacGregor et al. 2016, constraining the slope of the grain size distribution or the maximum grain size) or modelling their emission features in the mid-infrared (MIR, Olofsson et al. 2012, informing us about the dust composition). The second avenue is to study how stellar light is scattered off of the dust grains (in total intensity or polarised light and at optical or near-infrared wavelengths), either by measuring the colour of the disc between different bands (Debes et al. 2008; Rodigas et al. 2015) or studying the phase function (e.g. Olofsson et al. 2016; Milli et al. 2017, 2019; Ren et al. 2019). Both approaches can bring constraints on the typical grain sizes as well as their porosity. The phase function informs us how efficiently the light is scattered as a function of the scattering angle between the star, the dust grain, and the observer. This approach requires the disc to be spatially resolved and, therefore, became more popular in the past years with the availability of high angular resolution instruments such as the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE, Beuzit et al. 2019) or the Gemini Planet Imager (GPI, Perrin et al. 2015), but pioneering works were led with Hubble Space Telescope observations as well (e.g. Graham et al. 2007; Stark et al. 2014).

There are two different methods to estimate the phase function of debris discs, either from total intensity or polarised light observations, and in this study, we mostly focus on polarimetric observations. Out of the two approaches, the first one consists of fitting a model to the observations and the phase function is a free parameter of the modelling, either using scattering theory such as Mie or more complex ones (and in that case, the free parameter(s) mostly govern the grain size distribution or the porosity), or using parametrised approximations such as the Henyey-Greenstein one (Henyey & Greenstein 1941). The second approach being to measure the surface brightness of the disc as a function of the scattering angle, without the use of a model. The first approach requires a disc model for the dust density distribution as well as a scattering theory (or an approximation) that is able to reproduce the true phase function, while the second approach does not account for changes in column density at different azimuthal angles (as is the case along the semi-major axis of inclined discs). We note that in this method, the final phase function is not exactly equal to the surface brightness as several correction factors haveto be applied (e.g. illumination effects, point spread function dilution, and aperture shape, Milli et al. 2019).

The phase function is best measured for discs with a significant inclination (but not perfectly edge-on as most of the azimuthal information is lost), as a wide range of scattering angles can be probed. Inclinations around 75° are ideal, but asa consequence, there is an increase in column density along the major axis, if the disc is not infinitely flat. This increase in column density is illustrated in Fig. 1 showing disc models computed with different opening angles (increasing from left to right, see Sect. 2.2 for how the images are computed). The models shown in Fig. 1 all have an isotropic phase function and therefore directly trace the dust column density. As the opening angle increases, the major axis of the disc becomes brighter compared to the minor axis. Therefore when retrieving the phase function by measuring the surface brightness as a function of the scattering angle, we have to account for the column density variations along the disc. We here present a novel approach to retrieve the phase function in debris discs, in a non-parametric way, but that remains model-dependent on the density distribution throughout the disc.

Fig. 1 Convolved synthetic images of a debris disc with increasing opening angle. From left to right, ψ = arctan(h∕r) = 0.01, 0.03, 0.05, 0.07 (see text forfurther details).

2 Determining the phase function

In this section, we describe our new approach to determine the phase function in a non-parametric way, and briefly present the observationsused to test it beforehand.

2.1 Observations

Because of the brightness of the disc around HR 4796 A, we use the SPHERE/Zurich Imaging Polarimeter (ZIMPOL) polarimetric observations presented in Milli et al. (2019) and Olofsson et al. (2019) to illustrate our method. HR 4796 A is a young (8 ± 2 Myr old, Stauffer et al. 1995) nearby (71.9 ± 0.7 pc, Gaia Collaboration 2018) A-type star surrounded by one of the brightest debris disc (L_disc∕L_⋆ ~ 5 × 10⁻³, Moór et al. 2006). We used the Q_ϕ and U_ϕ images presented in Milli et al. (2019) as the signal to noise ratio is larger close to the semi-minor axis, at the expense of some small artefacts along the semi-major axis. The observations were obtained without a coronagraph and the reader is referred to Milli et al. (2019) for more information on the observing sequence and data processing. For all the calculations described below, we mask all the points that are within 0.′′ 35 from the estimated position of the star and all the pixels that are outside of an elliptical mask with a semi-major axis of 1.′′ 25 (see left panel of Fig. 2) are excluded when computing the goodness of fit.

2.2 Synthetic disc images

We use an updated version of the code presented in Olofsson et al. (2016, 2018) which can quickly produce images of (eccentric) debris discs which are not infinitely flat¹. To summarise briefly, the density distribution is computed as $$n(r,z)\propto {\left[{\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{\text{in}}}+{\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{\text{out}}}\right]}^{-1/2}\times {\text{e}}^{-z^2/2h^2},$$(1)

where n is the volumetric density, r₀ is a reference radius, α_in and α_out are the slopes of the density distribution and h is the vertical height of the disc (parametrised with the opening angle ψ such as tan ψ = h∕r). For eccentric discs, parametrised with two free parameters, the eccentricity e and the argument of pericenter ω, the reference radius r₀ depends on the azimuthal angle γ such as $$r_0=\frac{a(1-e^2)}{1+e\text{cos}(\omega +\gamma )},$$(2)

where a is the semi-major axis, γ = arctan2(y, x) is the azimuthal angle in the midplane², and x and y are the pixels coordinates (with the origin at the centre of the image) after projecting for the inclination i and rotating according to the position angle ϕ. To produce images, the code first defines a bounding box, sufficiently large for the volume density to be negligible at the borders. For each pixel of the image, the entry and exit points of the bounding box are calculated, and that “column” is divided into m = 100 equal parts of the same volume V. For each cell, the volume density is computed following Eq. (1), the scattering angle θ is computed from the dot product between the unit vector along the line of sight and the 3D coordinates at the centre of the cell, and the flux is the product of the volume density, V, and the scattering S₁₁ (θ) (or polarised S₁₂ (θ)) phase function. For each pixel, the final flux is the sum over the m cells. The user can provide a 2D array for S₁₂ (or S₁₁), a 1D array for θ (between 0 and π) and the code interpolates at the proper scattering angle when computing an image. The array for the phase function is 2D so that there can be one phase function for the north side and a different one for the south side. Both sides are identified based on the sign of the azimuthal angle (which ranges between − π and π). We note that the code is not flux-calibrated, therefore, the total dust mass or the polarisation degree are not mandatory input parameters.

To compare the synthetic images to the observations, we follow the approach explained in the appendix of Engler et al. (2018). To summarise briefly, the modelled Q_ϕ image is decomposed in Q and U images, according to the polar coordinates on the detector. Then both images are convolved with a 2D normal distribution with a full width at half maximum (FWHM) of 34 mas (σ = 2 pixels, comparable to the observing conditions as reported in Milli et al. 2019). The convolved Q and U images are then combined to obtain the final Q_ϕ image (see Engler et al. 2018 for further details).

Fig. 2 Observations, best-fit model with an isotropic phase function, final best-fit model with the revised phase function, and residuals to the SPHERE/ZIMPOL observations of HR 4796 A, with the same linear stretch (from left to right). The goodness of fit is estimated between the inner circle and the elliptical mask shown in the leftmost panel. In the central right panel, the location of the pericenter is marked with a circle.

2.3 Description of the approach

The motivation of our approach is to decouple the geometric parameters of the disc (e.g. radius, inclination, eccentricity) from the dust properties (the phase function in this case). This is achieved by comparing a first model computed using an isotropic phase function, which solely traces the dust density distribution, to the observations. By computing the ratio of surface brightness between this first model and the observations we can infer a revised phase function that should account for most of the differences. This revised phase function can then be injected in a new model to be compared to the observations.

We follow a three-step approach; for a given set of parameters we compute a first model with an isotropic phase function for both north and south sides (with the same number of pixels, and the same pixel scale as the observations, 7.2 mas per pixel) to trace the density distribution as a function of the scattering angle. Because the code is not flux-calibrated, we first find the best scaling factor S_scale that minimises the difference between the model and the observations, as $$S_{\text{Scale}}=\frac{{\displaystyle \sum \frac{I_{\text{obs}}\times I_{\text{model}}}{{\sigma }^2}}}{{\displaystyle \sum {\left(\frac{I_{\text{model}}}{\sigma }\right)}^2}},$$(3)

where I_model and I_obs are the images of the model and the Q_ϕ image, while σ are the uncertainties estimated from the U_ϕ image (see next section, and see Sect. 4.5 for a discussion on total intensity observations). The purpose of this first scaling is to try to account for most of the (unknown) multiplicative factors that govern the total flux of the disc (e.g. total dust mass, albedo).

Then, for each side of the disc (north and south), we estimate the surface brightness of both the model and the observations as a function of the scattering angle in the midplane of the disc. Given that the resulting distribution can be quite noisy for the observations we apply a numerical mask, selecting only the pixels where the model is brighter than 0.45 times its peak brightness. Given that we use an isotropic phase function for this first model, the whole ring is recovered, for a wide range of scattering angles.

The second step of the approach is to bin the surface brightness as a function of the scattering angle, for both the model and the observations. The binning is performed over 50 linearly spaced bins, and for each bin we compute the median value. For the observations, we also compute the median absolute deviation σ_i in each bin. Afterwards, we average the resulting distributions using a running mean over 5 neighbouring bins, to smooth them. For the observed profile, the corresponding uncertainties are estimated as $${\sigma }_{\text{binned}}={\left({\sum }_{i=0}^N\frac{1}{{\sigma }_i^2}\right)}^{-1/2},$$(4)

where N = 5 is the number of neighbouring bins.

This process is illustrated in Fig. 3 (the left panel shows an asymmetry as the disc model is slightly eccentric). Given the mask that excludes the central 0.′′ 35, we are able to probe scattering angles in the range [22°, 163°] on both the north and south sides. The polarised phase function is estimated from the ratio between the functions representingthe observations and the model (propagating the uncertainties at the same time), in other words, the red curve of theright panel divided by the red curve in the left panel of Fig. 3. The scaling factor S_scale of the isotropic model accounts for most of the differences with the observations but this model is not necessarily a good match, and as a consequence the y-axis of Fig. 3 may be different. When computing the ratio between the two averaged surface brightness, there may still be some contribution from unconstrained “flux-calibration” factors and those factors are incorporated in the resulting phase function, which is therefore a scaled (up or down) version of the true phase function³.

We then compute the final model with the newly evaluated phase function. Once the second model is computed we need to find the new best scaling factor S_Scale that minimises the χ² following Eq. (3). We note that computing the second model is only necessary to compare the model to the observationsand estimate a goodness of fit to find the best parameters of the disc model.

From Fig. 3, we can see that by estimating the phase function from the observed surface brightness, for a non-flat disc, we over-estimate it quite significantly at scattering angles near 90° (or under-estimate it a smaller angles), especially for highly inclined discs (e.g. Olofsson et al. 2016; Milli et al. 2019).

Fig. 3 Surface brightness of the north side of the disc, as a function of the scattering angle for the model with an isotropic phase function and for the observations (left and right, respectively). Each black circle represents a pixel. The red curves correspond to the median in each bin of θ.

3 Modelling and results

To determine most accurately the shape of the polarised phase function, we modelled the observations, trying to constrain the most relevant parameters of the disc. We put a strong emphasis on free parameters that can have an effect on the local increase in column density along the major axis. Therefore, the free parameters are the semi-major axis a, the inclination i, the outer slope of the density distribution α_out (α_in being fixed to + 25), the eccentricity e, the argument of periapsis ω, and the opening angle ψ. The position angle of the disc has already been well constrained for this dataset and we therefore use a value of − 152.1° following the results presented in Olofsson et al. (2019). The uncertainties are estimated from the U_ϕ image, computing the standard deviation in concentric annuli of 2 pixel width. Neither the Q_ϕ nor the U_ϕ images are convolved. Overall, since our model is convolved with a point spread function representative of the observations (following the approach outlined in Engler et al. 2018), that illumination effects are naturally accounted for in the model, and that we do not use an aperture to measure the surface brightness profiles, we do not need to apply correction factors such as the ones mentioned in the introduction and described in Milli et al. (2019).

To identify the most probable solution, we use the MultiNest algorithm (Feroz et al. 2009) interfaced to Python via the PyMultiNest package (Buchner et al. 2014). The probability distributions are plotted using the corner package (Foreman-Mackey 2016) and are presented in Fig. A.1. The best-fit values and their uncertainties are reported in Table 1, while the best-fit model and the residuals are presented in the centre right and right panels of Fig. 2 (the model with the isotropic phase function is shown in the centre left panel). We note that the uncertainties reported in Table 1 are most likely under-estimated and should be taken with some caution. This is most likely due to the uncertainties derived from the U_ϕ image used to compute the goodness of fit. By measuring the standard deviation in concentric annulii, a strategy commonly used in direct imaging studies, we may be under-estimating the true uncertainties, yielding larger χ² values. As a consequence, the Monte-Carlo algorithm may explore a narrower range of values, leading to narrow probability distributions. Finally, the polarised phase functions for the north and south sides of the disc are shown in Fig. 4, in black and dashed red, respectively.

We find that the semi-major axis is well constrained at a =1.′′066 ± 0.001, the inclination is i = 77.60° ± 0.06, the outer slope of the density distribution is α_out = −11.8 ± 0.2, the eccentricity e = 0.026 ± 0.002, the argument of periapsis is ω = −147.8° ± 5.2 (and its location is shown in the middle panel of Fig. 2), and finally, we find that the opening angle is ψ = 0.035 ± 0.001. From Fig. A.1, the degeneracy between e and ω is clearly noticeable, which explains why we find a smaller value for e compared tothe results of Milli et al. (2017, 2019), Olofsson et al. (2019) who found the argument of periapsis to be closer to the projected semi-minor axis of the disc.

Fig. 4 Polarised phase function for the north and south sides (black and dashed red, respectively) as a function of the scattering angle.

4 Discussion

4.1 Residuals and caveats

The residuals image overall shows that most of the signal from the disc has been removed, especially in the south side of the disc. There is still some signal left in the northern side, and this implies that the smoothed surface brightness profile that we estimated from the observations may not accurately capture the true surface brightness distribution of the northern side. As pointed out in Olofsson et al. (2019), the brightness asymmetry along the two sides cannotbe solely explained by pericenter glow (Wyatt et al. 1999) and that there may be an over-density of small dust grains at the north side of the disc. Our modelling approach is blind to local increase in dust density at different azimuthal angles, and there is no easy way around this issue. Nonetheless, the fact that the phase functions are quite similar between the north and south sides is quite reassuring, as we would not expect to have very different dust grains in different places of the disc (i.e., not surviving half an orbit). If indeed, there is an over-density of small dust grains along the north side as suggested in Olofsson et al. (2019), then the slight bump at 90° may not be real, and the true phase function may be more similar to the one of the south side.

4.2 Other attempt at determining the phase function

The approach presented in this study relies on some assumptions, for instance, how to estimate the surface brightness profiles of the model and of the observations. We originally attempted at least another approach that would have circumvented some of those issues, and we briefly discuss it here.

With the parameters of the disc fixed (a, i, etc.), we sampled the phase function over a small number of angles (10), and tried to fit the actual values of the phase function, without any prior nor parametrisation. The idea being that the code should find the shape of the phase function that minimises the χ². Unfortunately, the fitting never really converged. We postulate that the main issue with this approach is that we are trying to minimise second order effects. The χ² is mostly dominated by the geometric shape of the ring, and small changes in the shape of the phase function yields very small changes in the final χ² values. One possible work-around could be to work with relative χ² values, but fine-tuning the evaluation of the goodness of fit may not be that trivial, and overall, we deemed this possible solution out of the scope of this paper.

Another alternative possibility, as mentioned in the introduction, would be to assume a parametric form for the phase function with a handful of free parameters (e.g. weighted sums of several Henyey-Greenstein functions, or a polynomial form). However, the challenge of such an approach would be to estimate when to increase the complexity of the form and when to stop. With our method, for a given set of disc parameters (e.g. a, e, i), the phase function that we retrieve is the best phase function that would minimise the goodness of fit (but it does not necessarily mean that it is a good solution).

4.3 Impact of the inclination

To quantify when the column density variations make a significant difference, we compute a grid of models using an isotropic phase function. We then compute the synthetic surface brightness as a function of the scattering angle similarly to the left panel of Fig. 3. For each model we estimate the ratio between the maximum and minimum values of the profile. In the grid, we explore the two parameters that have the most important impact on the density increase; the inclination i and opening angle ψ. The inclination ranges between 30 and 85°, and the opening angle ranges between 0.01 and 0.06 radians. The semi-major axis, position angle, α_in, and α_out are set to the same values as before. To simplify the problem, we set e = 0 (therefore ω no longer matter). Figure 5 shows the ratio of density enhancement between the major and minor axis of the disc for the grid. For discs that have an opening angle of 0.04, the effect starts to become significant for inclinations larger than ~60°.

4.4 The vertical height of the disc around HR 4796 A

Milli et al. (2019) injected the phase function they inferred into a disc model and obtained residuals that are comparable to the ones of Fig. 2, but assumed a vertically thin disc (vertical height of 1 au at a reference radius of 77.4 au). Following the description laid out in Augereau et al. (1999), this would translate into an opening angle of ~0.009 for the code used in this work. Therefore, the true shape of the phase function depends on whether the disc is vertically flat or not.

Kennedy et al. (2018) modelled ALMA observations but could not firmly conclude on the vertical height of the disc. They mentioned that the disc could be vertically resolved, with a typical height of ~10 au at a radius of 80 au, but that a flat disc was also consistent with their observations. Given that a dynamical cold disc would be vertically thin (Thébault & Wu 2008, but see also Thébault 2009), and would explain its narrowness, the authors remained cautious and the actual vertical thickness of the disc remains a matter of debate.

Nonetheless, Thébault (2009) argued that debris discs should have a minimum aspect ratio H∕r of 0.04 ± 0.02 (where H is the local half width at half maximum) and that discs are most likely stratified for different grain sizes; the smallest dust grains having a larger aspect ratio compared to larger grains. Therefore, debris discs may appear vertically thicker in scattered light observations than at millimetre wavelengths. They also argue that the disc vertical height cannot directly be related to its dynamical excitation. Converting our best-fit value for the opening angle to the aspect ratio as defined in Thébault (2009), we obtain a value of 0.041, consistent with their results.

For an inclined disc, the width measured along the major and minor axis should depend on the vertical thickness; if the disc is vertically flat the minor axis should appear narrower than the major axis due to the inclination, as illustrated in Fig. 1. To better quantify this, we computed several models, with the same spatial resolution as the observations, varying the following two free parameters, the opening angle ψ and the outer slope of the dust density distribution α_out (which governs the width of the disc). For those models the semi-major axis and the inclination are the ones of the best fit model, and we used an isotropic phase function (in Appendix A.2 we repeat the same exercise for a different phase function to test the impact of this choice). Each image is convolved as explained in the previous section. We measure the FWHM along the major and minor axis of the disc model and compute their ratio. Figure 6 shows how the ratio varies as a function of the two free parameters. While the slope of the dust density distribution has a small effect on the ratio of FWHM, the opening angle has the strongest impact. Overall, this suggests that the angular resolution of the observations is sufficient to constrain the height of the disc, and that this can in principle be done by measuring the width of the disc as a function of the azimuthal angle.

However, this approach cannot be easily applied to our observations. The innermost regions are affected by strong noise, making the determination of the FWHM of the minor axis difficult. Furthermore, the background level is not the same along the minor and major axis of the disc, which may bias the peak values of both profiles, and hence the values of the FWHM. That being said, the test presented in Fig. 6 strongly suggests that the width of the disc at different azimuthal angles informs us about its vertical structure. This supports our findings that the distribution of small dust grains in the disc around HR 4796 A is most likely vertically extended.

Fig. 5 Map showing the surface brightness ratio between the major and minor axis of disc models with an isotropic phase function, as a function of the inclination and opening angle of the disc.

Fig. 6 Map showing the ratio of the FWHM along the semi-major and semi-minor axis of discs models, as a function of the outer slope of dust density distribution and the opening angle of the disc.

4.5 Angular differential observations

We presented a new approach to estimate the phase function measured in polarimetric observations, but it would be extremely valuable to also be able to measure the phase function in total intensity from angular differential imaging (ADI) observations. The main challenge when modelling ADI observations is that self-subtraction effects cannot easily be dealt with (Milli et al. 2012). After median-collapsing the de-rotated cube (after performing principal component analysis, or any other algorithms), we cannot measure the surface brightness of the disc free of biases. Therefore we cannot properly estimate the surface brightness of the disc to correct the phase function for column density effects.

One way to avoid biases induced by self-subtraction is to perform “forward modelling” by subtracting a disc model in the cube (and we then need another free parameter to scale up or down the model before subtracting it). Having the scattered light phase function as a free parameter (as mentioned in a previous sub section) runs into the same shortcomings. One possible, but time costly, approach would be to run a disc model with an isotropic total intensity phase function, find the best scaling factor for the given set of disc parameters, measure the residuals in the final image and modify the shape of the phase function accordingly before re-evaluating the scaling factor.

Nonetheless, we here attempt to roughly quantify the changes to the phase function measured on total intensity observationsof HR 4796 A. Milli et al. (2017) presented SPHERE/IRDIS observations of the disc, and measured the phase function from the surface brightness distribution, correcting for self-subtraction effects based on a disc model. This model has parameters that are compatible with our best fit solution, with a slightly different value for e (0.06) and therefore forω as well (− 105.7° in our reference frame) given the degeneracy between the two parameters. They then fitted two weighted Henyey-Greenstein functions f_HG to the measured phase function, as w × f_HG(g₀) + (1 − w) × f_HG(g₁), where f_HG takes the following form $$f_{\text{HG}}(g)=\frac{1}{4\pi }\frac{1-g^2}{{(1+g^2-2g\text{cos}\theta )}^{3/2}}.$$(5)

Milli et al. (2017) found that g₀ = 0.99, g₁ = −0.14, and w = 0.83. If the disc is not flat, the phase function they inferred does not take into account column density variations. We therefore used their analytical form, and applied an additional correction factor based on the dust density variation as a function of the scattering angle for a non-flat disc. The revised phase function being the original phase function divided by the surface brightness of the best fit model with an isotropic phase function. Both phase functions are shown in Fig. 7, normalised to unity at 90°. The most notable difference is that backward scattering becomes much more significant when accounting for the column density variations due to the inclination. It is worth noting, however, that this result remains model dependent for both the self-subtraction and the density variation corrections (and the discussion of Sect. 4.4).

Interestingly, as a side note, Ren et al. (2019) measured the surface brightness of the disc and halo around HD 191089 and measured strong backward scattering for the halo but not for the main ring. If the halo is vertically very thin, extending mostly from the densest regions, that is, the midplane, then we should not expect significant column density changes due to the inclination of the system (~59°) and therefore the measurement of the phase function from the halo would not suffer from the same biases described in this paper.

Fig. 7 Fit to the total intensity phase function derived in Milli et al. (2017) using two Henyey-Greenstein functions (black) “corrected” for projection effects (dashed red). Both functions are normalised to unity at 90°.

4.6 Dust properties in the disc around HR 4796 A

With the revised phase functions, we can attempt to revise the dust properties inferred in Milli et al. (2019). We computed a grid of polarised phase functions, using the OpacityTool (Woitke et al. 2016, Toon & Ackerman 1981). The code can compute absorption and scattering properties, as well as six elements of the scattering matrix, including S₁₁ and S₁₂, assuming a distribution of hollow spheres (DHS, Min et al. 2005). We assume that the grain size distribution is a differential power-law of the form dn(s) ∝ s^−3.5ds, where s is the grain size, and we set s_max = 1 mm. We compute the polarised phase function, integrated over the size distribution, varying s_min between 0.01 μm and 110 μm, and the porosity between 0 and 99%. The optical constant are taken from Dorschner et al. (1995, amorphous silicates) with a 85% mass fraction and Zubko et al. (1996, amorphous carbon) with a 15% mass fraction. The maximum filling factor is set to 0.8. To estimate whether we can reproduce the inferred phase functions, we computed a grid of models, and the χ² maps for both the north and south sides are shown in Fig. 8 (left and right, respectively). For both panels, the insets show the inferred phase functions and their best-fit model. For both sides, we find that the best model is obtained for s_min = 0.01 μm. For the north side, we find that theporosity should be 16%, and 0% for the south side, with significant uncertainties as illustrated in Fig. 8. In Appendix A.3 we show in more detail the effects of porosity and minimum grain size on the shape of the phase function.

While the best-fit model reproduces rather well the phase function for the south side, it fails to capture the shape of the north side. We have to keep in mind that, as discussed before, the northern phase function may be biased by a possible over-density of small dust grains, but overall, we find that the phase functions suggest the presence of very small dust grains, with rather low porosity values.

With the exercise described above, the slope of the grain size distribution is set to − 3.5, to minimise the number of free parameters, but which may be a strong hypothesis. The grain size distribution can show some wavy structures (e.g. Thébault & Augereau 2007) or there can be an over-abundance of small dust grains in bright debris discs (as is the case for HR 4796 A, Thebault & Kral 2019). Therefore, we repeated the same exercise as before, but integrating the size distribution between s and s + δ_s, keeping the slope as − 3.5. The motivation being to identify the characteristic size that can best explain the phase functions. We use PyMultiNest to find the best fit model and for the north side, we obtain s_min = 0.03 ± 0.02μm, δ_s = 0.29 ± 0.15μm, and a porosity of 16 ± 1%. For the south side, we find s_min = 0.02 ± 0.01 μm, δ_s = 0.36 ± 0.29μm, and a porosity compatible with 0%. The shape of the best fit model is quite similar to the previous best fit models. The results of this approach are quite similar to the previous ones, the phase functions are best reproduced by very small dust grains.

Using the OpacityTool, we can also compute the asymmetry parameter g_sca and compute the unitless β ratio between the stellar radiation pressure and the gravitational force for different grain sizes, as $$\beta (s)=\frac{3L_{\star }}{16\pi G{c}^2{M}_{\star }}\frac{Q_{\text{pr}}(s)}{\rho s},$$(6)

where L_⋆ and M_⋆ are the stellar luminosity and mass (25.75L_⊙ and 1.31M_⊙, respectively, Olofsson et al. 2019), G the gravitational constant, c the speed of light, ρ the dust density, and Q_pr the radiation pressure efficiency (equal to Q_ext(λ, s) − g_sca(s) × Q_sca(λ, s)) averaged over the stellar spectrum. For porosity values of 0 and 16% we find that the blow-out sizes (for which β ≤ 0.5, assuming the parent bodies are on circular orbit) are ~13.5 and 17μm (for larger porosity values, the blow-out size increases even more). All the grains that are smaller no longer are bound to the star and would be removed from the system rapidly. We therefore reach the same conclusions as the ones presented in Milli et al. (2017, 2019), that the Mie or DHS theory cannot adequately explain the full picture. Indeed, Augereau et al. (1999) found that to reproduce the spectral energy distribution of the disc, the minimum grain size should be close to 10μm, which is rather compatible with the aforementioned blow-out size but would fail to reproduce the measured phase function. Relatively large aggregates composed ofsub-μm sized monomers may be a viable alternative to explain the observations. As explained in Min et al. (2016), the polarisation properties of aggregates are intimately related to the size of the individual monomers and not to the overall size of the aggregate itself.

Fig. 8 Maps of the χ2 when trying to reproduce the shape of the phase functions using DHS for the north and south sides (left and right, respectively). The two free parameters are the porosity and the minimum grain sizes. The insets show the observations and the best-fitting model.

5 Conclusions

In this paper, we presented an alternative approach to estimate the phase function from polarised observations of debris discs with an emphasis on discs that have a non-negligible vertical scale height. While our method remains model-dependent it does not require a parametrised form for the phase function (e.g. Henyey-Greenstein). The total flux depends both on the local density and the true phase function, but when the discs are highly inclined, and not infinitely flat, there are variations in column density along the major axis, due to projection effects, and those variations have to be taken into account.

We presented an approach to account for those column density variation effects, and find that the inferred phase function is quite different from previous estimates. We tested our model to SPHERE/ZIMPOL observations of HR 4796 A and derived phase functions for both the north and south sides, both being relatively similar. We reach similar conclusions as the ones outlined in Milli et al. (2019), that is, we cannot fully reconcile all key aspects with a single scattering theory (e.g. phase function and blow-out size). Our modelling results suggest that the disc is not vertically flat, with an opening angle of ψ ~ 0.035. The vertical scale height can successfully be constrained by the model based on how the width of the disc varies as a function of the azimuthal angle.

We also note that our modelling approach remains blind to any local increase of the dust density, and that it cannot readily be applied to ADI observations. We remark however that similar biases are probably occurring when deriving the total intensity phase function, which may lead to an under-estimation of backward scattering.

Acknowledgements

We are grateful to the referee, Kees Dullemond, for providing helpful comments, especially with respect to the determination of the vertical scale height, as well as several pointers to help clarify the paper. This research has made use of the SIMBAD database (operated at CDS, Strasbourg, France). This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013). J.O. and A.B. acknowledge support from the ICM (Iniciativa Científica Milenio) via the Nucleo Milenio de Formación planetaria grant. J.O. acknowledges support from the Universidad de Valparaíso and from Fondecyt (grant 1180395). A.B. acknowledges support from Fondecyt (grant 1190748).

Appendix A Miscellaneous

A.1 Corner plot

Figure A.1 shows the corner plot for the modelling, with the density plots and the projected probability distributions for each parameter.

Fig. A.1 Corner plot of the posterior density distributions for the modelling.

A.2 The impact of the phase function on the apparent width of the disc

In Sect. 4.4 we computed the ratio of FWHM along the major and minor axis of disc models to assess whether the width of thedisc at different azimuthal angles can inform us about its vertical scale height. In Fig. 6, we used an isotropic phase function, but the phase function can change the intensity along the major and minor axis, and as a consequence the measureof the FWHM. Therefore, we here repeat the same exercise but with a different phase function. We computed models in total intensity (and not polarised light) using the Henyey-Greenstein approximation and with a g value of 0.9. This choice is motivated by the fact that the phase function strongly peaks at small scattering angles, significantly enhancing the intensity along the minor axis. Therefore, this phase function and the isotropic one are very complementary ones, allowing us to further assess the robustness of our findings. The fact that we are computing total intensity images and not polarimetric images is not relevant for the interpretation of the results. We proceeded in the same way as described in Sect. 4.4, computing the models, convolving them and measuring the two FWHM. The results are presented in Fig. A.2 and the results are compatible with Fig. 6 with comparable values for the ratio. We therefore conclude that while the choice of the phase function as a small impact on the appearance of the disc, it is of second order compared to the effect of the density distribution.

Fig. A.2 Same as Fig. 6 using a different phase function for the models.

A.3 Porosity and grain size

To complement Fig. 8, Fig. A.3 shows the effect of the porosity (left panel) and minimum grain size (right panel) on the phase function, compared to the phase function of the south side of the disc, to highlight how the shape varies as a function of those two parameters.

Fig. A.3 Different phases functions compared to the phase function estimated for the south side of the disc, for different values of the porosity (left) and minimum grain size (right). The color bars are in units of percent and μm (left and right, respectively).

References

All Tables

All Figures

	Fig. 1 Convolved synthetic images of a debris disc with increasing opening angle. From left to right, ψ = arctan(h∕r) = 0.01, 0.03, 0.05, 0.07 (see text forfurther details).
In the text

Fig. 2 Observations, best-fit model with an isotropic phase function, final best-fit model with the revised phase function, and residuals to the SPHERE/ZIMPOL observations of HR 4796 A, with the same linear stretch (from left to right). The goodness of fit is estimated between the inner circle and the elliptical mask shown in the leftmost panel. In the central right panel, the location of the pericenter is marked with a circle.

In the text

	Fig. 3 Surface brightness of the north side of the disc, as a function of the scattering angle for the model with an isotropic phase function and for the observations (left and right, respectively). Each black circle represents a pixel. The red curves correspond to the median in each bin of θ.
In the text

	Fig. 4 Polarised phase function for the north and south sides (black and dashed red, respectively) as a function of the scattering angle.
In the text

	Fig. 5 Map showing the surface brightness ratio between the major and minor axis of disc models with an isotropic phase function, as a function of the inclination and opening angle of the disc.
In the text

	Fig. 6 Map showing the ratio of the FWHM along the semi-major and semi-minor axis of discs models, as a function of the outer slope of dust density distribution and the opening angle of the disc.
In the text

	Fig. 7 Fit to the total intensity phase function derived in Milli et al. (2017) using two Henyey-Greenstein functions (black) “corrected” for projection effects (dashed red). Both functions are normalised to unity at 90°.
In the text

	Fig. 8 Maps of the χ2 when trying to reproduce the shape of the phase functions using DHS for the north and south sides (left and right, respectively). The two free parameters are the porosity and the minimum grain sizes. The insets show the observations and the best-fitting model.
In the text

	Fig. A.1 Corner plot of the posterior density distributions for the modelling.
In the text

	Fig. A.2 Same as Fig. 6 using a different phase function for the models.
In the text

	Fig. A.3 Different phases functions compared to the phase function estimated for the south side of the disc, for different values of the porosity (left) and minimum grain size (right). The color bars are in units of percent and μm (left and right, respectively).
In the text