The challenge of measuring the phase function of debris disks. Application to HR\,4796
AAstronomy & Astrophysics manuscript no. phasefunction c (cid:13)
ESO 2020June 16, 2020
The challenge of measuring the phase function of debris disks
Application to HR 4796
J. Olofsson , , J. Milli , A. Bayo , , Th. Henning , and N. Engler Instituto de Física y Astronomía, Facultad de Ciencias, Universidad de Valparaíso, Av. Gran Bretaña 1111, Playa Ancha, Valparaíso,Chilee-mail: [email protected] Núcleo Milenio Formación Planetaria - NPF, Universidad de Valparaíso, Av. Gran Bretaña 1111, Valparaíso, Chile Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France Max Planck Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany ETH Zurich, Institute for Particle Physics and Astrophysics, Wolfgang-Pauli-Strasse 27, CH-8093 Zurich, SwitzerlandJune 16, 2020
ABSTRACT
Context.
Debris disks are valuable systems to study dust properties. Because they are optically thin at all wavelengths, we have directaccess to the absorption and scattering properties of the dust grains. One very promising technique to study them is to measure theirphase function, i.e., the scattering e ffi ciency as a function of the scattering angle. Disks that are highly inclined are promising targetsas a wider range of scattering angles can be probed. Aims.
The phase function (polarized or total intensity) is usually either inferred by comparing the observations to synthetic diskmodels assuming a parametrized phase function, or estimating it from the surface brightness of the disk. We argue here that the latterapproach can be biased due to projection e ff ects leading to an increase in column density along the major axis of a non flat disk. Methods.
We present a novel approach to account for those column density e ff ects. The method remains model dependent, as onestill requires a disk model to estimate the density variations as a function of the scattering angle. This method allows us however toestimate the shape of the phase function without having to invoke any parametrized form. Results.
We apply our method to SPHERE / ZIMPOL observations of HR 4796 and highlight the di ff erences with previous measure-ments using the surface brightness only, the main di ff erences being at scattering angles smaller than ∼ ◦ . Our modelling resultssuggests that the disk is not vertically flat at optical wavelengths, result supported by comparing the width along the major and minoraxis of synthetic images. We discuss some of the caveats of the approach, mostly that our method remains blind to real local increaseof the dust density, and that it cannot yet be readily applied to angular di ff erential imaging observations. Conclusions.
We show that the vertical thickness of inclined ( ≥ ◦ ) debris disks can a ff ect the determination of their phase functions.Similarly to previous studies on HR 4796, we still cannot reconcile the full picture using a given scattering theory to explain the shapeof the phase function, the blow-out size due to radiation pressure and the shape of the spectral energy distribution, a long lastingproblem for debris disks. Nonetheless, we argue that similar e ff ects as the ones highlighted in this study can also bias the determinationof the phase function in total intensity. Key words.
Stars: individual (HR 4796 A) – circumstellar matter – Techniques: high angular resolution – Scattering
1. Introduction
Dust grains are the building blocks of planets, but there are rela-tively few ways to accurately characterize their properties. Stud-ies of solar system bodies provide the strongest constraints onthe constituents of comets or asteroids (e.g., Frattin et al. 2019,Bertini et al. 2019), but do not inform us directly on what istaking place during the planet formation stage. Observations ofdisks around young ( ≤
100 Myr old) stars are sensitive to grainsizes that are smaller than typically 100 µ m or a few mm. De-bris disks, disks of second generation dust, are ideal targets tostudy dust grains. As those disks are optically thin at all wave-lengths, we have direct access to the absorption and scatteringproperties of the grains, without having to account for non triv-ial optical depth e ff ects or multiple scattering events. There aretwo possible ways to characterize dust properties in debris disks,first, via their thermal emission by measuring the spectral slopeat (sub-) mm wavelengths (Draine 2006, MacGregor et al. 2016,constraining the slope of the grain size distribution or the maxi- mum grain size), or modelling their emission features in the mid-IR (Olofsson et al. 2012, informing about the dust composition). The second avenue is to study how stellar light is scatteredo ff of the dust grains (in total intensity or polarized light, at opti-cal or near-infrared wavelengths), either by measuring the colourof the disk between di ff erent bands (Debes et al. 2008, Rodigaset al. 2015), or studying the phase function (e.g., Olofsson et al.2016; Milli et al. 2017, 2019, Ren et al. 2019). Both approachescan bring constraints on the typical grain sizes as well as theirporosity. The phase function informs us how e ffi ciently the lightis scattered as a function of the scattering angle (between thestar, the dust grain and the observer). This approach requires thedisk to be spatially resolved and therefore became more popularin the past years with the availability of high angular resolutioninstruments such as VLT / SPHERE (Beuzit et al. 2019) or GPI(Perrin et al. 2015), but pioneering works were led with
HubbleSpace Telescope observations as well (e.g., Graham et al. 2007,Stark et al. 2014).
Article number, page 1 of 13 a r X i v : . [ a s t r o - ph . S R ] J un & A proofs: manuscript no. phasefunction ]1.00.50.00.51.0 [ ] ] 1.00.50.00.51.0 [ ] 1.00.50.00.51.0 [ ] Fig. 1.
Convolved synthetic images of a debris disk with increasing opening angle. From left to right, ψ = arctan( h / r ) = .
01, 0 .
03, 0 .
05, 0 .
07 (seetext for further details).
There are two di ff erent methods to estimate the phase func-tion of debris disks, either from total intensity or polarized lightobservations, and in this study, we will focus mostly on polari-metric observations. Out of the two approaches, the first one con-sists of fitting a model to the observations and the phase functionis a free parameter of the modelling, either using scattering the-ory such as Mie or more complex ones (and in that case, thefree parameter(s) mostly govern the grain size distribution orthe porosity), or using parametrized approximations such as theHenyey-Greenstein one (Henyey & Greenstein 1941). The sec-ond approach being to measure the surface brightness of the diskas a function of the scattering angle, without the use of a model.The first approach requires a disk model for the dust density dis-tribution as well as a scattering theory (or an approximation) thatis able to reproduce the true phase function, while the secondapproach does not account for changes in column density at dif-ferent azimuthal angles (as is the case along the semi-major axisof inclined disks). One should note that in this method, the finalphase function is not exactly equal to the surface brightness asseveral correction factors have to be applied (e.g., illuminatione ff ects, point spread function dilution, and aperture shape, Milliet al. 2019).The phase function is best measured for disks with a signif-icant inclination (but not perfectly edge-on as most of the az-imuthal information is lost), as a wide range of scattering anglescan be probed. Inclinations around 75 ◦ are ideal, but as a conse-quence, there is an increase in column density along the majoraxis, if the disk is not infinitely flat. This increase in column den-sity is illustrated in Figure 1 showing disk models computed withdi ff erent opening angles (increasing from left to right, see Sec-tion 2.2 for how the images are computed). The models shown inFig. 1 all have an isotropic phase function and therefore directlytrace the dust column density. One can note that as the open-ing angle increases, the major axis of the disk becomes brightercompared to the minor axis. Therefore when retrieving the phasefunction by measuring the surface brightness as a function ofthe scattering angle, one has to account for the column densityvariations along the disk. We here present a novel approach toretrieve the phase function in debris disks, in a non-parametricway, but that remains model-dependent on the density distribu-tion throughout the disk.
2. Determining the phase function
In this Section, we describe our new approach to determine thephase function in a non-parametric way, and briefly present theobservations used to test it beforehand.
Because of the brightness of the disk around HR 4796, we usethe SPHERE / ZIMPOL polarimetric observations presented inMilli et al. (2019) and Olofsson et al. (2019) to illustrate ourmethod. HR 4796 is a young (8 ± ff er et al. 1995)nearby (71 . ± . L disk / L (cid:63) ∼ × − , Moór et al. 2006). We used the Q φ and U φ images pre-sented in Milli et al. (2019) as the signal to noise ratio is largerclose to the semi-minor axis, at the expense of some small arte-facts along the semi-major axis. The observations were obtainedwithout a coronagraph and the reader is referred to Milli et al.(2019) for more information on the observing sequence and dataprocessing. For all the calculations described below, we mask allthe points that are within 0 (cid:48)(cid:48) .
35 from the estimated position of thestar and all the pixels that are outside of an elliptical mask witha semi-major axis of 1 (cid:48)(cid:48) .
25 (see left panel of Fig. 2) are excludedwhen computing the goodness of fit.
We use an updated version of the code presented in Olofssonet al. (2016, 2018) which can quickly produce images of (ec-centric) debris disks which are not infinitely flat. To summarizebriefly, the density distribution is computed as n ( r , z ) ∝ (cid:32) rr (cid:33) − α in + (cid:32) rr (cid:33) − α out − / × e − z / h , (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 verticalheight of the disk (parametrized with the opening angle ψ suchas tan ψ = h / r ). For eccentric disks, parametrized with two freeparameters, the eccentricity e and the argument of pericenter ω ,the reference radius r depends on the azimuthal angle γ such as r = a (1 − e )1 + e cos( ω + γ ) , (2) Article number, page 2 of 13. Olofsson et al.: The challenge of measuring the phase function of debris disks where a is the semi-major axis, γ = arctan2( y , x ) is the azimuthalangle in the midplane , and x and y are the pixels coordinates(with the origin at the centre of the image) after projecting forthe inclination i and rotating according to the position angle φ .To produce images, the code first defines a bounding box, su ffi -ciently large for the volume density to be negligible at the bor-ders. For each pixel of the image, the entry and exit points of thebounding box are calculated, and that “column” is divided into m =
100 equal parts of the same volume V . For each cell, thevolume density is computed following Eq. 1, the scattering an-gle θ is computed from the dot product between the unit vectoralong the line of sight and the 3D coordinates at the centre ofthe cell, and the flux will be the product of the volume density, V , and the scattering S ( θ ) (or polarized S ( θ )) phase function.For each pixel, the final flux is the sum over the m cells. The usercan provide a 2D array for S (or S ), a 1D array for θ (between0 and π ) and the code will interpolate at the proper scattering an-gle when computing an image. The array for the phase functionis 2D so that there can be one phase function for the north sideand a di ff erent one for the south side. Both sides are identifiedbased on the sign of the azimuthal angle (which ranges between − π and π ). One should note that the code is not flux-calibrated,therefore, the total dust mass or the polarization degree are notmandatory input parameters.To compare the synthetic images to the observations, we fol-low the approach explained in the Appendix of Engler et al.(2018). To summarize briefly, the modelled Q φ image is decom-posed in Q and U images, according to the polar coordinateson the detector. Then both images are convolved with a 2D nor-mal distribution with a full width at half maximum (FWHM) of34 mas ( σ = Q and U imagesare then combined to obtain the final Q φ image (see Engler et al.2018 for further details). The motivation of our approach is to decouple the geometric pa-rameters of the disk (e.g., radius, inclination, eccentricity) fromthe dust properties (the phase function in this case). This isachieved by comparing a first model computed using an isotropicphase function, which solely traces the dust density distribution,to the observations. By computing the ratio of surface bright-ness between this first model and the observations one can infera revised phase function that should account for most of the dif-ferences. This revised phase function can then be injected in anew model to be compared to the observations.We follow a three-step approach; for a given set of parame-ters we compute a first model with an isotropic phase functionfor both north and south sides (with the same number of pixels,and the same pixel scale as the observations, 7 . S scale that minimizes the di ff erence betweenthe model and the observations, as S Scale = (cid:80) I obs × I model σ (cid:80) (cid:16) I model σ (cid:17) , (3)where I model and I obs are the images of the model and the Q φ im-age, while σ are the uncertainties estimated from the U φ image We use arctan2 to place the resulting angle in the proper quadrantwhen both y and x are negative for instance. (see next Section, and see Section 4.5 for a discussion on totalintensity observations). The purpose of this first scaling is to tryto account for most of the (unknown) multiplicative factors thatgovern the total flux of the disk (e.g., total dust mass, albedo).Then, for each side of the disk (north and south) , we es-timate the surface brightness of both the model and the obser-vations as a function of the scattering angle in the midplane ofthe disk. Given that the resulting distribution can be quite noisyfor the observations we apply a numerical mask, selecting onlythe pixels where the model is brighter than 0 .
45 times its peakbrightness. Given that we use an isotropic phase function for thisfirst model, the whole ring is recovered, for a wide range of scat-tering angles.The second step of the approach is to bin the surface bright-ness as a function of the scattering angle, for both the modeland the observations. The binning is performed over 50 linearlyspaced bins, and for each bin we compute the median value. Forthe observations, we also compute the median absolute deviation σ i in each bin. Afterwards, we average the resulting distributionsusing a running mean over 5 neighbouring bins, to smooth them.For the observed profile, the corresponding uncertainties are es-timated as σ binned = N (cid:88) i = σ i − / , (4)where N = (cid:48)(cid:48) .
35, we are able to probescattering angles in the range [22 ◦ , ◦ ] on both the north andsouth sides. The polarized phase function is estimated from theratio between the functions representing the observations and themodel (propagating the uncertainties at the same time), in otherwords, the red curve of the right panel divided by the red curve inthe left panel of Figure 3. The scaling factor S scale of the isotropicmodel accounts for most of the di ff erences with the observationsbut this model is not necessarily a good match, and as a conse-quence the y -axis of Figure 3 may be di ff erent. When computingthe ratio between the two averaged surface brightness there maystill be some contribution from unconstrained “flux-calibration”factors and those factors will be incorporated in the resultingphase function, which is therefore a scaled (up or down) versionof the true phase function .We then compute the final model with the newly evaluatedphase function. Once the second model is computed we need tofind the new best scaling factor S Scale that minimizes the χ fol-lowing Eq. 3. One should note that computing the second modelis only necessary to compare the model to the observations andestimate a goodness of fit to find the best parameters of the diskmodel.From Figure 3, one can see that by estimating the phase func-tion from the observed surface brightness, for a non flat disk, weover-estimate it quite significantly at scattering angles near 90 ◦ (or under-estimate it a smaller angles), especially for highly in-clined disks (e.g., Olofsson et al. 2016, Milli et al. 2019).
3. Modelling and results
To determine most accurately the shape of the polarized phasefunction, we modelled the observations, trying to constrain the One could normalize the phase function over 4 π steradians but sincewe do not probe the full range of scattering angles, this would remainan approximation Article number, page 3 of 13 & A proofs: manuscript no. phasefunction ]1.00.50.00.51.0 [ ] ] 1.00.50.00.51.0 [ ] 1.00.50.00.51.0 [ ] Fig. 2.
Observations, best-fit model with an istropic phase function, final best-fit model with the revised phase function, and residuals to theSPHERE / ZIMPOL observations of HR 4796 A, with the same linear stretch (from left to right). The goodness of fit is estimated between the innercircle and the elliptical mask shown in the leftmost panel. In the central right panel, the location of the pericenter is marked with a circle.
20 40 60 80 100 120 140 160 [ ]2.53.03.54.04.55.05.56.0 S u r f a c e b r i g h t n e ss [ a r b i t r a r y un i t s ]
20 40 60 80 100 120 140 160 [ ]02468
Fig. 3.
Surface brightness of the north side of the disk, as a function of the scattering angle for the model with an isotropic phase function and forthe observations (left and right, respectively). Each black circle represents a pixel. The red curves correspond to the median in each bin of θ . Table 1.
Free parameters and best-fit results.
Parameters Prior Best-fit a [ (cid:48)(cid:48) ] [0 .
90, 1 .
15] 1 . ± . i [ ◦ ] [70, 80] 77 . ± . α out [ − − − . ± . e [0 .
0, 0 .
1] 0 . ± . ω [ ◦ ] [ − − − . ± . ψ [ ◦ ] [0 . .
06] 0 . ± . ff ect on the local increase incolumn density along the major axis. Therefore, the free parame-ters are the semi-major axis a , the inclination i , the outer slope ofthe density distribution α out ( α in being fixed to + e , the argument of periapsis ω , and the opening angle ψ .The position angle of the disk has already been well constrainedfor this dataset and we therefore use a value of − . ◦ followingthe results presented in Olofsson et al. (2019). The uncertaintiesare estimated from the U φ image, computing the standard devia-tion in concentric annuli of 2 pixel width. Neither the Q φ nor the U φ images are convolved. Overall, since our model is convolvedwith a point spread function representative of the observations(following the approach outlined in Engler et al. 2018), that il-lumination e ff ects are naturally accounted for in the model, andthat we do not use an aperture to measure the surface brightnessprofiles, we do not need to apply correction factors such as theones 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). Theprobability distributions are plotted using the corner package(Foreman-Mackey 2016) and are presented in Figure A.1. Thebest-fit values and their uncertainties are reported in Table 1,while the best-fit model and the residuals are presented in thecentre right and right panels of Figure 2 (the model with theisotropic phase function is shown in the centre left panel). Oneshould note that the uncertainties reported in Table 1 are mostlikely 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 thestandard deviation in concentric annulii, a strategy commonly Article number, page 4 of 13. Olofsson et al.: The challenge of measuring the phase function of debris disks S SouthNorth
Fig. 4.
Polarized phase function for the north and south sides (black anddashed red, respectively) as a function of the scattering angle. used in direct imaging studies, we may be under-estimating thetrue uncertainties, yielding larger χ values. As a consequence,the Monte-Carlo algorithm may explore a narrower range of val-ues, leading to narrow probability distributions. Finally, the po-larized phase functions for the north and south sides of the diskare shown in Figure 4, in black and dashed red, respectively.We find that the semi-major axis is well constrained at a = (cid:48)(cid:48) . ± . i = . ◦ ± .
06, the outer slopeof the density distribution is α out = − . ± .
2, the eccentricity e = . ± . ω = − . ◦ ± . ψ = . ± . e and ω is clearly notice-able, which explains why we find a smaller value for e comparedto the results of Milli et al. (2017, 2019), Olofsson et al. (2019)who found the argument of periapsis to be closer to the projectedsemi-minor axis of the disk.
4. Discussion
The residuals image overall shows that most of the signal fromthe disk has been removed, especially in the south side of thedisk. There is still some signal left in the northern side, and thisimplies that the smoothed surface brightness profile that we esti-mated from the observations may not accurately capture the truesurface brightness distribution of the northern side. As pointedout in Olofsson et al. (2019), the brightness asymmetry along thetwo sides cannot be solely explained by pericenter glow (Wyattet al. 1999) and that there may be an over-density of small dustgrains at the north side of the disk. Our modelling approach isblind to local increase in dust density at di ff erent azimuthal an-gles, and there is no easy way around this issue. Nonetheless, thefact that the phase functions are quite similar between the northand south sides is quite reassuring, as one would not expect tohave very di ff erent dust grains in di ff erent places of the disk (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 Olofssonet al. (2019), then the slight bump at 90 ◦ may not be real, and thetrue phase function may be more similar to the one of the southside. The approach presented in this study relies on some assump-tions, for instance, how to estimate the surface brightness pro-files of the model and of the observations. We originally at-tempted at least another approach that would have circumventedsome of those issues, and we briefly discuss it here.With the parameters of the disk fixed ( a , i , etc), we sampledthe phase function over a small number of angles (10), and triedto fit the actual values of the phase function, without any priornor parametrization. The idea being that the code should findthe shape of the phase function that minimizes the χ . Unfortu-nately, the fitting never really converged. We postulate that themain issue with this approach is that we are trying to minimizesecond order e ff ects. The χ is mostly dominated by the geomet-ric shape of the ring, and small changes in the shape of the phasefunction yields very small changes in the final χ values. Onepossible work-around could be to work with relative χ values,but fine-tuning the evaluation of the goodness of fit may not bethat trivial, and overall, we deemed this possible solution out ofthe scope of this paper.Another alternative possibility, as mentioned in the intro-duction, would be to assume a parametric form for the phasefunction with a handful of free parameters (e.g., weighted sumsof several Henyey-Greenstein functions, or a polynomial form).However, the challenge of such an approach would be to estimatewhen to increase the complexity of the form and when to stop.With our method, for a given set of disk parameters (e.g., a , e , i ),the phase function that we retrieve is the best phase function thatwould minimize the goodness of fit (but it does not necessarilymean that it is a good solution). To quantify when the column density variations make a signifi-cant di ff erence, we compute a grid of models using an isotropicphase function. We then compute the synthetic surface bright-ness as a function of the scattering angle similarly to the leftpanel of Figure 3. For each model we estimate the ratio betweenthe maximum and minimum values of the profile. In the grid,we explore the two parameters that have the most important im-pact on the density increase; the inclination i and opening angle ψ . The inclination ranges between 30 and 85 ◦ , and the open-ing angle ranges between 0 .
01 and 0 .
06 radians. The semi-majoraxis, position angle, α in , and α out are set to the same values asbefore. To simplify the problem, we set e = ω nolonger matter). Figure 5 shows the ratio of density enhancementbetween the major and minor axis of the disk for the grid. Fordisks that have an opening angle of 0 .
04, the e ff ect starts to be-come significant for inclinations larger than ∼ ◦ . Milli et al. (2019) injected the phase function they inferred intoa disk model and obtained residuals that are comparable to theones of Figure 2, but assumed a vertically thin disk (verticalheight of 1 au at a reference radius of 77 . Article number, page 5 of 13 & A proofs: manuscript no. phasefunction
30 40 50 60 70 80 i [ ]0.010.020.030.040.050.06 1.11.21.31.41.51.61.7 Fig. 5.
Map showing the surface brightness ratio between the major andminor axis of disk models with an isotropic phase function, as a functionof the inclination and opening angle of the disk.
20 18 16 14 12 10 8 6 out
Fig. 6.
Map showing the ratio of the FWHM along the semi-major andsemi-minor axis of disks models, as a function of the outer slope of dustdensity distribution and the opening angle of the disk. late into an opening angle of ∼ .
009 for the code used in thiswork. Therefore, the true shape of the phase function depends onwhether the disk is vertically flat or not.Kennedy et al. (2018) modelled ALMA observations butcould not firmly conclude on the vertical height of the disk. Theymentioned that the disk could be vertically resolved, with a typi-cal height of ∼
10 au at a radius of 80 au, but that a flat disk wasalso consistent with their observations. Given that a dynamical cold disk would be vertically thin (Thébault & Wu 2008, but seealso Thébault 2009), and would explain its narrowness, the au-thors remained cautious and the actual vertical thickness of thedisk remains a matter of debate.Nonetheless, Thébault (2009) argued that debris disks shouldhave a minimum aspect ratio H / r of 0 . ± .
02 (where H is thelocal half width at half maximum) and that disks are most likelystratified for di ff erent grain sizes; the smallest dust grains havinga larger aspect ratio compared to larger grains. Therefore, debrisdisks may appear vertically thicker in scattered light observa-tions than at millimetre wavelengths. They also argue that thedisk vertical height cannot directly be related to its dynamicalexcitation. Converting our best-fit value for the opening angle tothe aspect ratio as defined in Thébault (2009), we obtain a valueof 0 . ψ and the outer slope of thedust density distribution α out (which governs the width of thedisk). For those models the semi-major axis and the inclinationare the ones of the best fit model, and we used an isotropic phasefunction (in App. A.2 we repeat the same exercise for a di ff er-ent phase function to test the impact of this choice). Each imageis convolved as explained in the previous Section. We measurethe FWHM along the major and minor axis of the disk modeland compute their ratio. Figure 6 shows how the ratio varies as afunction of the two free parameters. While the slope of the dustdensity distribution has a small e ff ect on the ratio of FWHM,the opening angle has the strongest impact. Overall, this sug-gests that the angular resolution of the observations is su ffi cientto constrain the height of the disk, and that this can in principlebe done by measuring the width of the disk as a function of theazimuthal angle.However, this approach cannot be easily applied to our ob-servations. The innermost regions are a ff ected by strong noise,making the determination of the FWHM of the minor axis dif-ficult. Furthermore, the background level is not the same alongthe minor and major axis of the disk, which may bias the peakvalues of both profiles, and hence the values of the FWHM. Thatbeing said, the test presented in Figure 6 strongly suggests thatthe width of the disk at di ff erent azimuthal angles informs usabout its vertical structure. This supports our findings that thedistribution of small dust grains in the disk around HR 4796 ismost likely vertically extended. We presented a new approach to estimate the phase functionmeasured in polarimetric observations, but it would be extremelyvaluable to also be able to measure the phase function in totalintensity from angular di ff erential imaging (ADI) observations.The main challenge when modelling ADI observations is thatself-subtraction e ff ects cannot easily be dealt with (Milli et al.2012). After median-collapsing the de-rotated cube (after per-forming principal component analysis, or any other algorithms),one cannot measure the surface brightness of the disk free of bi-ases. Therefore we cannot properly estimate the surface bright-ness of the disk to correct the phase function for column densitye ff ects. Article number, page 6 of 13. Olofsson et al.: The challenge of measuring the phase function of debris disks
20 40 60 80 100 120 140 [ ]246810 S Original phase functionCorrected phase function
Fig. 7.
Fit to the total intensity phase function derived in Milli et al.(2017) using two Henyey-Greenstein functions (black) “corrected” forprojection e ff ects (dashed red). Both functions are normalized to unityat 90 ◦ . One way to avoid biases induced by self-subtraction is toperform “forward modelling” by subtracting a disk model inthe cube (and we then need another free parameter to scale upor down the model before subtracting it). Having the scatteredlight phase function as a free parameter (as mentioned in a pre-vious sub section) runs into the same shortcomings. One possi-ble, but time costly, approach would be to run a disk model withan isotropic total intensity phase function, find the best scalingfactor for the given set of disk parameters, measure the residualsin the final image and modify the shape of the phase functionaccordingly before re-evaluating the scaling factor.Nonetheless, we here attempt to roughly quantify thechanges to the phase function measured on total inten-sity observations of HR 4796. Milli et al. (2017) presentedSPHERE / IRDIS observations of the disk, and measured thephase function from the surface brightness distribution, correct-ing for self-subtraction e ff ects based on a disk model. This modelhas parameters that are compatible with our best fit solution,with a slightly di ff erent value for e (0 .
06) and therefore for ω as well ( − . ◦ in our reference frame) given the degener-acy between the two parameters. They then fitted two weightedHenyey-Greenstein functions f HG to the measured phase func-tion, as w × f HG ( g ) + (1 − w ) × f HG ( g ), where f HG takes thefollowing form f HG ( g ) = π − g (1 + g − g cos θ ) / . (5)Milli et al. (2017) found that g = . g = − .
14, and w = .
83. If the disk is not flat, the phase function they in-ferred does not take into account column density variations. Wetherefore used their analytical form, and applied an additionalcorrection factor based on the dust density variation as a func-tion of the scattering angle for a non-flat disk. The revised phasefunction being the original phase function divided by the surfacebrightness of the best fit model with an isotropic phase function.Both phase functions are shown in Figure 7, normalized to unity at 90 ◦ . The most notable di ff erence is that backward scatteringbecomes much more significant when accounting for the col-umn density variations due to the inclination. One should keep inmind however that this result remains model dependent for boththe self-subtraction and the density variation corrections (and thediscussion of Section 4.4).Interestingly, as a side note, Ren et al. (2019) measured thesurface brightness of the disk and halo around HD 191089 andmeasured strong backward scattering for the halo but not for themain ring. If the halo is vertically very thin, extending mostlyfrom the densest regions, i.e., the midplane, then one should notexpect significant column density changes due to the inclina-tion of the system ( ∼ ◦ ) and therefore the measurement ofthe phase function from the halo would not su ff er from the samebiases described in this paper. With the revised phase functions, we can attempt to revise thedust properties inferred in Milli et al. (2019). We computeda grid of polarized phase functions, using the
OpacityTool (Woitke et al. 2016, Toon & Ackerman 1981). The code cancompute absorption and scattering properties, as well as six ele-ments of the scattering matrix, including S and S , assuminga distribution of hollow spheres (DHS, Min et al. 2005). We as-sume that the grain size distribution is a di ff erential power-lawof the form d n ( s ) ∝ s − . d s , where s is the grain size, and weset s max = s min between 0 . µ m and 110 µ m, and the porosity between 0 and 99%. The op-tical constant are taken from Dorschner et al. (1995, amorphoussilicates) with a 85% mass fraction and Zubko et al. (1996, amor-phous carbon) with a 15% mass fraction. The maximum fillingfactor is set to 0 .
8. To estimate whether we can reproduce theinferred phase functions, we computed a grid of models, and the χ maps for both the north and south sides are shown in Figure 8(left and right, respectively). For both panels, the insets show theinferred phase functions and their best-fit model. For both sides,we find that the best model is obtained for s min = . µ m. Forthe north side, we find that the porosity should be 16%, and 0%for the south side, with significant uncertainties as illustrated inFig. 8. In Appendix A.3 we show in more detail the e ff ects ofporosity and minimum grain size on the shape of the phase func-tion.While the best-fit model reproduces rather well the phasefunction for the south side, it fails to capture the shape of thenorth side. One has 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 phasefunctions suggest the presence of very small dust grains, withrather low porosity values.With the exercise described above, the slope of the grain sizedistribution is set to − .
5, to minimize the number of free pa-rameters, but which may be a strong hypothesis. The grain sizedistribution can show some wavy structures (e.g., Thébault &Augereau 2007) or there can be an over-abundance of small dustgrains in bright debris disks (as is the case for HR 4796, The-bault & Kral 2019). Therefore, we repeated the same exercise asbefore, but integrating the size distribution between s and s + δ s ,keeping the slope as − .
5. The motivation being to identify thecharacteristic size that can best explain the phase functions. Weuse
PyMultiNest to find the best fit model and for the northside, we obtain s min = . ± . µ m, δ s = . ± . µ m, and aporosity of 16 ± s min = . ± . Article number, page 7 of 13 & A proofs: manuscript no. phasefunction s [ m ] North South
Fig. 8.
Maps of the χ 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 µ m, δ s = . ± . µ m, and a porosity compatible with 0%. Theshape of the best fit model is quite similar to the previous bestfit models. The results of this approach are quite similar to theprevious ones, the phase functions are best reproduced by verysmall dust grains.Using the OpacityTool , we can also compute the asymme-try parameter g sca and compute the unitless β ratio between thestellar radiation pressure and the gravitational force for di ff erentgrain sizes, as β ( s ) = L (cid:63) π Gc M (cid:63) Q pr ( s ) ρ s , (6)where L (cid:63) and M (cid:63) are the stellar luminosity and mass (25 . L (cid:12) and 1 . M (cid:12) , respectively, Olofsson et al. 2019), G the gravita-tional constant, c the speed of light, ρ the dust density, and Q pr the radiation pressure e ffi ciency (equal to Q ext ( λ, s ) − g sca ( s ) × Q sca ( λ, s )) averaged over the stellar spectrum. For porosity val-ues of 0 and 16% we find that the blow-out sizes (for which β ≤ .
5, assuming the parent bodies are on circular orbit) are ∼ . µ m (for larger porosity values, the blow-out sizeincreases even more). All the grains that are smaller no longerare bound to the star and would be removed from the systemrapidly. We therefore reach the same conclusions as the onespresented in Milli et al. (2017, 2019), that the Mie or DHS the-ory cannot adequately explain the full picture. Indeed, Augereauet al. (1999) found that to reproduce the spectral energy dis-tribution of the disk, the minimum grain size should be closeto 10 µ m, which is rather compatible with the aforementionedblow-out size but would fail to reproduce the measured phasefunction. Relatively large aggregates composed of sub- µ m sizedmonomers may be a viable alternative to explain the observa-tions. As explained in Min et al. (2016), the polarization prop-erties of aggregates are intimately related to the size of the in-dividual monomers and not to the overall size of the aggregateitself.
5. Conclusions
In this paper, we presented an alternative approach to estimatethe phase function from polarized observations of debris disks with an emphasis on disks that have a non negligible verticalscale height. While our method remains model-dependent it doesnot require a parametrized form for the phase function (e.g.,Henyey-Greenstein). The total flux depends both on the localdensity and the true phase function, but when the disks are highlyinclined, and not infinitely flat, there are variations in columndensity along the major axis, due to projection e ff ects, and thosevariations have to be taken into account.We presented an approach to account for those column den-sity variation e ff ects, and find that the inferred phase function isquite di ff erent from previous estimates. We tested our model toSPHERE / ZIMPOL observations of HR 4796 and derived phasefunctions for both the north and south sides, both being rela-tively similar. We reach similar conclusions as the ones outlinedin Milli et al. (2019), i.e., we cannot fully reconcile all key as-pects with a single scattering theory (e.g., phase function andblow-out size). Our modelling results suggest that the disk is notvertically flat, with an opening angle of ψ ∼ . Acknowledgements.
We are grateful to the referee, Kees Dullemond, for provid-ing helpful comments, especially with respect to the determination of the verticalscale height, as well as several pointers to help clarify the paper. This researchhas made use of the SIMBAD database (operated at CDS, Strasbourg, France).This research made use of Astropy, a community-developed core Python packagefor Astronomy (Astropy Collaboration et al. 2013). J. O. and A. B. acknowledgesupport from the ICM (Iniciativa Científica Milenio) via the Nucleo Milenio deFormación planetaria grant. J. O. acknowledges support from the Universidadde Valparaíso and from Fondecyt (grant 1180395). A. B. acknowledges supportfrom Fondecyt (grant 1190748).
References
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558,
Article number, page 8 of 13. Olofsson et al.: The challenge of measuring the phase function of debris disks
A33Augereau, J. C., Lagrange, A. M., Mouillet, D., Papaloizou, J. C. B., & Grorod,P. A. 1999, A&A, 348, 557Bertini, I., La Forgia, F., Fulle, M., et al. 2019, MNRAS, 482, 2924Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125Debes, J. H., Weinberger, A. J., & Song, I. 2008, ApJ, 684, L41Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995,A&A, 300, 503Draine, B. T. 2006, ApJ, 636, 1114Engler, N., Schmid, H. M., Quanz, S. P., Avenhaus, H., & Bazzon, A. 2018,A&A, 618, A151Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24Frattin, E., Muñoz, O., Moreno, F., et al. 2019, MNRAS, 484, 2198Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1Graham, J. R., Kalas, P. G., & Matthews, B. C. 2007, ApJ, 654, 595Henyey, L. G. & Greenstein, J. L. 1941, ApJ, 93, 70Kennedy, G. M., Marino, S., Matrà, L., et al. 2018, MNRAS, 475, 4924MacGregor, M. A., Wilner, D. J., Chandler, C., et al. 2016, ApJ, 823, 79Milli, J., Engler, N., Schmid, H. M., et al. 2019, A&A, 626, A54Milli, J., Mouillet, D., Lagrange, A.-M., et al. 2012, A&A, 545, A111Milli, J., Vigan, A., Mouillet, D., et al. 2017, A&A, 599, A108Min, M., Hovenier, J. W., & de Koter, A. 2005, A&A, 432, 909Min, M., Rab, C., Woitke, P., Dominik, C., & Ménard, F. 2016, A&A, 585, A13Moór, A., Ábrahám, P., Derekas, A., et al. 2006, ApJ, 644, 525Olofsson, J., Juhász, A., Henning, T., et al. 2012, A&A, 542, A90Olofsson, J., Milli, J., Thébault, P., et al. 2019, A&A, 630, A142Olofsson, J., Samland, M., Avenhaus, H., et al. 2016, A&A, 591, A108Olofsson, J., van Holstein, R. G., Boccaletti, A., et al. 2018, A&A, 617, A109Perrin, M. D., Duchene, G., Millar-Blanchaer, M., et al. 2015, ApJ, 799, 182Ren, B., Choquet, É., Perrin, M. D., et al. 2019, ApJ, 882, 64Rodigas, T. J., Stark, C. C., Weinberger, A., et al. 2015, ApJ, 798, 96Stark, C. C., Schneider, G., Weinberger, A. J., et al. 2014, ApJ, 789, 58Stau ff er, J. R., Hartmann, L. W., & Barrado y Navascues, D. 1995, ApJ, 454, 910Thébault, P. 2009, A&A, 505, 1269Thébault, P. & Augereau, J.-C. 2007, A&A, 472, 169Thebault, P. & Kral, Q. 2019, A&A, 626, A24Thébault, P. & Wu, Y. 2008, A&A, 481, 713Toon, O. B. & Ackerman, T. P. 1981, Appl. Opt., 20, 3657Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103Wyatt, M. C., Dermott, S. F., Telesco, C. M., et al. 1999, ApJ, 527, 918Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282,1321 Article number, page 9 of 13 & A proofs: manuscript no. phasefunction
Appendix A: Miscellaneous
Appendix A.1: Corner plot
Figure A.1 shows the corner plot for the modelling, with the den-sity plots and the projected probability distributions for each pa-rameter.
Appendix A.2: The impact of the phase function on theapparent width of the disk
In Section 4.4 we computed the ratio of FWHM along the majorand minor axis of disk models to assess whether the width of thedisk at di ff erent azimuthal angles can inform us about its verticalscale height. In Figure 6, we used an isotropic phase function, butthe phase function can change the intensity along the major andminor axis, and as a consequence the measure of the FWHM.Therefore, we here repeat the same exercise but with a di ff erentphase function. We computed models in total intensity (and notpolarized light) using the Henyey-Greenstein approximation andwith a g value of 0 .
9. This choice is motivated by the fact that thephase function strongly peaks at small scattering angles, signif-icantly enhancing the intensity along the minor axis. Therefore,this phase function and the isotropic one are very complementaryones, allowing us to further assess the robustness of our findings.The fact that we are computing total intensity images and not po-larimetric images is not relevant for the interpretation of the re-sults. We proceeded in the same way as described in Section 4.4,computing the models, convolving them and measuring the twoFWHM. The results are presented in Figure A.2 and the resultsare compatible with Figure 6 with comparable values for the ra-tio. We therefore conclude that while the choice of the phasefunction as a small impact on the appearance of the disk, it is ofsecond order compared to the e ff ect of the density distribution. Appendix A.3: Porosity and grain size
To complement Figure 8, Figure A.3 shows the e ff ect of theporosity (left panel) and minimum grain size (right panel) onthe phase function, compared to the phase function of the southside of the disk, to highlight how the shape varies as a functionof those two parameters. Article number, page 10 of 13. Olofsson et al.: The challenge of measuring the phase function of debris disks . . . i [] .
31 2 .
01 1 .
71 1 . o u t . . . . . e [] . . . . a [ ] . . . . . [] . . . i [ ] . . . . out . . . . . e [ ] . . . . . [ ] Fig. A.1.
Corner plot of the posterior density distributions for the modelling. Article number, page 11 of 13 & A proofs: manuscript no. phasefunction
20 18 16 14 12 10 8 6 out
Fig. A.2.
Same as Figure 6 using a di ff erent phase function for the mod-els.Article number, page 12 of 13. Olofsson et al.: The challenge of measuring the phase function of debris disks
20 40 60 80 100 120 140 160 [ ]0.00.51.01.52.02.53.03.5 S S Fig. A.3. Di ff erent phases functions compared to the phase function estimated for the south side of the disk, for di ff erent values of the porosity(left) and minimum grain size (right). The color bars are in units of percent and µµ