# Azimuthal asymmetries in the debris disk around HD61005 A massive collision of planetesimals?

J. Olofsson, M. Samland, H. Avenhaus, C. Caceres, Th. Henning, A. Moor, J.

Milli, H. Canovas, S. P. Quanz, M. R. Schreiber, et al.

► To cite this version:

J. Olofsson, M. Samland, H. Avenhaus, C. Caceres, Th. Henning, et al., Azimuthal asymmetries in the debris disk around HD61005 A massive collision of planetesimals?. Astronomy & Astrophysics - A&A, 2016, 591, 10.1051/0004-6361/201628196, hal-01440164

HAL Id: hal-01440164

https://hal.science/hal-01440164v1

Submitted on 9 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Azimuthal asymmetries in the debris disk around HD 61005

## A massive collision of planetesimals? ★★★

J. Olofsson¹,²,³, M. Samland³, H. Avenhaus⁴,², C. Caceres¹,², Th. Henning³, A. Moór⁵, J. Milli⁶, H. Canovas¹,², S. P. Quanz⁷, M. R. Schreiber¹,², J.-C. Augereau⁸,⁹, A. Bayo¹,², A. Bazzon⁷, J.-L. Beuzit⁸,⁹, A. Boccaletti¹⁰, E. Buenzli⁷, S. Casassus⁴,², G. Chauvin⁸,⁹, C. Dominik¹¹, S. Desidera¹², M. Feldt³, R. Gratton¹², M. Janson³,¹³, A.-M. Lagrange⁸,⁹, M. Langlois¹⁴,¹⁵, J. Lannier⁸,⁹, A.-L. Maire¹²,³, D. Mesa¹², C. Pinte¹⁶,⁸, D. Rouan¹⁰, G. Salter¹⁵, C. Thalmann⁷, and A. Vigan¹⁵,⁶

*(Affiliations can be found after the references)*

Received 26 January 2016 / Accepted 16 April 2016

## ABSTRACT

**Context.** Debris disks offer valuable insights into the latest stages of circumstellar disk evolution, and can possibly help us to trace the outcomes of planetary formation processes. In the age range 10 to 100 Myr, most of the gas is expected to have been removed from the system, giant planets (if any) must have already been formed, and the formation of terrestrial planets may be on-going. Pluto-sized planetesimals, and their debris released in a collisional cascade, are under their mutual gravitational influence, which may result into non-axisymmetric structures in the debris disk.

**Aims.** High angular resolution observations are required to investigate these effects and constrain the dynamical evolution of debris disks. Furthermore, multi-wavelength observations can provide information about the dust dynamics by probing different grain sizes.

**Methods.** Here we present new VLT/SPHERE and ALMA observations of the debris disk around the 40 Myr-old solar-type star HD 61005. We resolve the disk at unprecedented resolution both in the near-infrared (in scattered and polarized light) and at millimeter wavelengths. We perform a detailed modeling of these observations, including the spectral energy distribution.

**Results.** Thanks to the new observations, we propose a solution for both the radial and azimuthal distribution of the dust grains in the debris disk. We find that the disk has a moderate eccentricity (e ~ 0.1) and that the dust density is two times larger at the pericenter compared to the apocenter.

**Conclusions.** With no giant planets detected in our observations, we investigate alternative explanations besides planet-disk interactions to interpret the inferred disk morphology. We postulate that the morphology of the disk could be the consequence of a massive collision between ~1000 km-sized bodies at ~61 au. If this interpretation holds, it would put stringent constraints on the formation of massive planetesimals at large distances from the star.

**Key words.** circumstellar matter – zodiacal dust – techniques: high angular resolution

## 1. Introduction

Debris disks are the leftovers of star and planetary formation processes (see Wyatt 2008; Krivov 2010; Matthews et al. 2014, for recent reviews). Departure from photospheric emission at infrared (IR) wavelengths was first discovered around Vega (Aumann et al. 1984) using the Infrared Astronomical Satellite (IRAS). This excess emission was originally thought to be the remnant of the cloud out of which Vega formed. Several decades later, we now know that the dusty “debris” responsible for the IR emission are arranged in a disk comparable to the Edgeworth-Kuiper belt in the solar system. The dust grains, with sizes between a few µm to a few millimeters, located at tens of au

from the central star are heated by the stellar radiation and re-emit at mid-IR and mm wavelengths. Since the original discovery, and mostly thanks to space-based missions such as the Infrared Space Observatory, *Spitzer*, and *Herschel*, several hundred of main sequence stars are known to harbor debris disks (e.g., Eiroa et al. 2013; Chen et al. 2014).

Recent decades have seen incredible progress in the field of disk observations (e.g., Augereau et al. 1999; Kalas et al. 2005; Buenzli et al. 2010; Lebreton et al. 2012; Millar-Blanchaer et al. 2015). Observations with ever improving spatial resolution have revealed asymmetric disks (e.g., Lagrange et al. 2016; Kalas et al. 2015) as well as complex, moving, small-scale structures (Boccaletti et al. 2015). Nonetheless, spatially resolved observations of debris disks remain relatively rare and even though there are theoretical works focusing on the dynamical evolution of debris disks (e.g., Dominik & Decin 2003; Kenyon & Bromley 2006), they still need to be confronted with the observations. The current paradigm is that the primordial gas-rich proto-planetary disks are thought to have a half-life time of about 2–3 Myr (Hernández et al. 2007). As a disk evolves Pluto-sized planetesimals can form (Johansen et al. 2015) along with giant planets which may accrete their mass from the gas reservoir (core accretion or gravitational instability scenarios).

<footnote>* Based on observations made with ESO Telescopes at the Paranal Observatory under programs ID 095.C-0298 and 095.C-0273. Based on *Herschel* observations, OBSIDs: 1342270977, 1342270978, 1342270979, 1342270989, and 1342255147. *Herschel* is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.</footnote>

<footnote>** The reduced images as FITS files, and data of Fig. 1 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/591/A108</footnote>
<table><caption>**Table 1.** Log for the VLT/SPHERE and ALMA observations.</caption><thead><tr><th colspan="8">VLT/SPHERE</th></tr><tr><th>Observing date<br>[YYYY-MM-DD]</th><th>Prog. ID</th><th>Instrument Mode</th><th>Filter</th><th>Seeing<br>["]</th><th>Airmass</th><th>Coherence time<br>[ms]</th></tr></thead><tbody><tr><td>2015-02-03</td><td>95.C-0298</td><td>IRDIFS</td><td>H<sub>2</sub>H<sub>3</sub>/YJ</td><td>0.67</td><td>1.01</td><td>22.0</td></tr><tr><td>2015-03-30</td><td>95.C-0298</td><td>IRDIFS_EXT</td><td>K<sub>1</sub>K<sub>2</sub>/YH</td><td>1.24</td><td>1.04</td><td>1.7</td></tr><tr><td>2015-05-01</td><td>95.C-0273</td><td>IRDIS DPI</td><td>B_H</td><td>1.16</td><td>1.29</td><td>1.7</td></tr></tbody><thead><tr><th colspan="8">ALMA</th></tr><tr><th>Observing date<br>[YYYY-MM-DD]</th><th>Prog. ID</th><th>Mode</th><th>Resolution<br>[kHz]</th><th>Frequency range<br>[GHz]</th><th>PWV<br>[mm]</th><th>Integration<br>[s]</th></tr></thead><tbody><tr><td rowspan="2">2014-03-20</td><td rowspan="2">2012.1.00437.S</td><td rowspan="2">Continuum<br>Gas</td><td>31 250.00</td><td>211.91–228.97</td><td rowspan="2">0.79</td><td rowspan="2">120.96</td></tr><tr><td>488.28</td><td>230.05–230.99</td></tr></tbody></table>

After a few Myr, the gaseous content is quickly removed from the disk by efficient processes such as photo-evaporation (e.g., Alexander et al. 2006; Owen et al. 2011). Only already formed planets, planetesimals and dust grains will thus remain while the disk enters its debris disk phase. After a short phase of runaway growth, terrestrial planets may form in the inner regions of the disk, with Pluto-sized bodies in the outer regions, via chaotic growth of these oligarchs, on a timescale of 10–100 Myr (Kenyon & Bromley 2006, 2008, 2010). The time evolution of the entire system then becomes more regular and less chaotic. The km-sized bodies, arranged in one or more planetesimal belt(s), evolve under their mutual gravitational influence. Through collisions, they continuously release small particles in a collisional cascade. Small dust grains, in turn, are removed from the system either by radiation pressure or Poynting-Robertson drag. Therefore, one can consider that after ~100 Myr, planetary formation has stopped. The system is left with one (or more) planetesimal belt(s) and, quite possibly, with planets of various masses. By observing systems in the range 10–100 Myr, one can therefore study the time evolution of debris disks. Spatially resolved images can provide constraints on the radial and azimuthal distribution of the dust, giving us insight about the dynamics at stake in these systems.

Here, we present Very Large Telescope (VLT) Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) and Atacama Large Millimeter/submillimeter Array (ALMA) observations of the debris disk around the solar type star HD 61005 (G8V), located at a distance of 35.4 ± 1.1 pc (van Leeuwen 2007). The age of the system is believed to be within 40±10 Myr old, based on membership of the Argus association (Desidera et al. 2011; De Silva et al. 2013; Elliott et al. 2014). The uncertainties for the age are mostly related to the dispersion in ages reported in the literature for both the Argus association and the IC 321 super cluster. In the last ten years, it has been spatially resolved on multiple occasions with several instruments; Hines et al. (2007, Hubble Space Telescope HST/NICMOS); Maness et al. (2009, HST/ACS); Buenzli et al. (2010, VLT/NaCo); Ricarte et al. (2013, Submillimeter Array SMA); and Schneider et al. (2014, HST/STIS). The disk earned its nickname of The Moth because of the swept-back wings first revealed in the HST observations of Hines et al. (2007). The wings may originate from the interaction between the interstellar medium (ISM) and the disk itself; as the star moves through the local ISM, (small) dust grains are set on eccentric orbits, drifting away from the central star.

While Hines et al. (2007) and Maness et al. (2009) mostly studied the wings, the study presented in Buenzli et al. (2010)

resolved the debris disk as a ring, thanks to the better angular resolution provided by the NaCo instrument. They found the disk to be almost edge-on (i = 84.3 ± 1°), to have a semi-major axis of 61.25 ± 0.85 au with an eccentricity of e = 0.045 ± 0.015, which translates into an offset of 2.75 ± 0.85 au of the star with respect to the disk. A brightness asymmetry between the two ansae is observed, which cannot be fully explained by the off-centering. No planets were detected by Buenzli et al. (2010), for observations that should, in principle, have detected planets with masses starting from a few Jupiter masses. The integrated luminosity of the disk is large at IR wavelengths ($L_{disk}/L_{\star} \sim 3 \times 10^{-3}$) and is best modeled by two spatially separated dust belts. The main dust belt is the one located at ~60 au presented in Buenzli et al. (2010), but a warm component is usually required to reproduce the spectral energy distribution (SED). Ricarte et al. (2013) argue that this additional belt is mandatory to match the SED at about 20 µm, however its location remains unconstrained. In this paper, we focus on constraining the properties of the debris disk which is resolved at unprecedented angular resolution at both near-IR and mm wavelengths. Studying the properties and origin of the swept-back wings is beyond the scope of this paper, since they are marginally detected with these observations.

## 2. Observations, data processing, and stellar parameters

Table 1 summarizes the VLT/SPHERE and ALMA observations presented in this study.

### 2.1. VLT/SPHERE IRDIS observations and data reduction

The star HD 61005 was observed with the VLT/SPHERE (Beuzit et al. 2008), within the guaranteed time consortium. The observations were obtained in different instrumental set-ups in February, March, and May 2015, using the dual-band imager (IRDIS, Dohlen et al. 2008; Vigan et al. 2010), the integral field spectrograph (IFS, Claudi et al. 2008), and the dual-polarization imager (IRDIS DPI, Langlois et al. 2014).

#### 2.1.1. IRDIS dual-band observations

The February observations were conducted in the H₂H₃ dual band (centered on 1.59 and 1.67 µm) for IRDIS and the YJ-band (0.95–1.35 µm, at spectral resolution R ~ 54) for IFS. The March observations used the K₁K₂ filters (centered on 2.11 and 2.25 µm) for IRDIS and the YH-band (0.95–1.65 µm, R ~ 33) for IFS. All of these observations were performed using an
apodized Lyot coronagraph, consisting of a focal mask with a diameter of 185 milli-arcsec (N\_ALC\_YJH\_S) and a corresponding pupil mask. Coronagraphic observations were performed in pupil stabilized mode to use angular differential imaging (ADI) post-processing (Marois et al. 2006) to attenuate residual speckle noise. The observation strategy can be summarized as follows: 1) photometric calibration: imaging of star offset from coronagraph mask to obtain PSF for relative photometric calibration; 2) centering: imaging with star behind mask with four artificially induced satellite spots for centering; 3) science: coronagraphic sequence; 4) centering: same as point two; 5) photometric calibration: same as point one; 6) sky background observation using same DIT as coronagraphic sequence. Finally, true north and plate scale are determined using astrometric calibrators as part of the SPHERE GTO survey for each run (Maire et al. 2016).

Basic reduction of the IRDIS data (background subtraction, flat fielding, centering) was performed using the SPHERE Data Reduction Handling (DRH) pipeline (Pavlov et al. 2008, version 15.0). The output consists of cubes for each filter, re-centered onto a common origin using the satellite spot reference. The cubes are then corrected for the true north position determined from the astrometric calibrations and for distortion. We collapsed the two filters together to increase the signal-to-noise ratio (S/N), and from now on we will refer to the $H_2H_3$ and $K_1K_2$ datasets as H and K observations, respectively. The datacubes in both bands were then processed using a principal component analysis (PCA, using the implementation of the scikit-learn Python package, Pedregosa et al. 2011) approach from all frames within the data-cube.

## 2.1.2. IRDIS dual polarization observations

On May 1, 2015, the target was observed using IRDIS in DPI mode in H-band. The same coronagraph was used for these observations. IRDIS DPI splits the light into two perpendicular polarization directions imaged at the same time on the same detector. Full cycles of half-wave plate (HWP) positions (0, 22.5, 45, and 67.5°) were taken to construct the Stokes Q and U vectors. The strategy of the observations was to take as long exposures as possible without saturating the detector just outside the coronagraph to achieve the best possible inner working angle as well as the best S/N for the outer part of the disk. Two integration times (DIT = 64 s with a total integration time of 768 s and DIT = 16 s with a total integration time of 3008 s) were used. IRDIS suffers from a de-polarization effect at certain detector position angles, which depend on the parallactic angle at the time of observation. Because the parallactic angle changed rapidly during the observations, we updated the detector angle at regular intervals.

The DPI data were reduced using a custom pipeline that is different to the DRH pipeline, which closely follows the processes described in Avenhaus et al. (2014), using the double-difference method (see also Canovas et al. 2011) to construct Stokes Q and U vectors from the data. The pipeline has been adapted to suit the IRDIS instrument. The data were centered using the centering frames taken just before and after the science observations. A de-rotation was applied to bring all files to the same orientation (see section above). Furthermore, the frames were corrected to take account of the fact that the IRDIS pixel scale differs slightly (~0.6%) in the two principal detector directions. The files were then corrected for true north as determined for IRDIS.

Because scattered light in an optically thin (debris) disk is expected to be polarized perpendicular to the line between the star

and the image point in question, we then construct local Stokes Q and U vectors, denoted $Q_\phi$ and $U_\phi$ (see also Benisty et al. 2015). In case of single scattering of the stellar light, $Q_\phi$ is expected to contain the disk signal, while $U_\phi$ is expected to contain no signal (with possible exceptions when the optical depth is large, Canovas et al. 2015) but noise on the same level as the $Q_\phi$ image and can serve as a noise estimator. $Q_\phi$ and $U_\phi$ can be calculated as:

$$ 
\begin{aligned}
Q_{\phi} &= +Q \cos(2\Phi) + U \sin(2\Phi) \\
U_{\phi} &= -Q \sin(2\Phi) + U \cos(2\Phi),
\end{aligned}
\qquad (1) $$

where $\phi$ refers to the azimuth in polar coordinates, and $\Phi$ is the position angle of the location of interest (x, y), with respect to the stellar location $(x_0, y_0)$ as:

$$ \Phi = \arctan \frac{x - x_0}{y - y_0} + \theta, \qquad (2) $$

where $\theta$ corrects for instrumental effects such as a small misalignment of the half-wave plate.

During the data reduction process, one HWP cycle equivalent to 256 s of data (DIT = 16 s) was taken out because the telescope had lost tracking for a short amount of time, rendering this data unusable. The result are two pairs of $Q_\phi$ and $U_\phi$ images, one for the DIT = 16 s and DIT = 64 s observations each. These were then combined with a weighted average to produce the final $Q_\phi$ and $U_\phi$ images.

## 2.1.3. IFS observations

The IFS data proved difficult to be properly reduced at the time of this analysis, mainly because of centering problems. Presenting and analyzing these observations will be postponed for a future study.

## 2.1.4. Processed images

Figure 1 shows the final reduction for our dataset (with a central mask of radius 0.15"). The top row displays the IRDIS ADI data in H and K-bands (left and middle panels, respectively), and the IRDIS DPI $Q_\phi$ image (right), all with a linear stretch. For each image, the bottom row shows a S/N map, with a linear stretch between [-3$\sigma$, 3$\sigma$]. The noise map is calculated from the reduced images and represents the standard deviation in concentric annuli, centered on the star, with a constant width (2 pixels). We did not mask out the disk when computing the noise maps, hence the uncertainties might be slightly over-estimated. For the DPI observations, the noise map is computed from the $U_\phi$ image which does not seem to contain any signal from the disk (Fig. C.1 shows the $U_\phi$ image with the same linear stretch as the $Q_\phi$ image in Fig. 1).

We note that the disk is not homogeneously detected at high S/N. The east side is detected at larger S/N in all datasets and it appears brighter than the west side. Such asymmetry, already reported in Buenzli et al. (2010), are further investigated in this paper. For the sake of simplicity, in the rest of the paper, we refer to the ADI and DPI datasets as “scattered” and “polarized” observations, respectively. Strictly speaking, this is not correct as polarized photons must have been scattered by dust grains.

## 2.2. ALMA observations

HD 61005 was observed with ALMA in Band 6 (PI: David Rodriguez, program 2012.1.00437.S), in the frequency range
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Reduced SPHERE observations of HD 61005 used in the analysis. North is up, East is left. *From left to right:* IRDIS ADI in *H* and *K*-bands (PCA with 6 components), and IRDIS DPI *Q*<sub>Φ</sub> in *H*-band. *Top row* shows the data in linear stretch, with a central mask of 0.15", and the *bottom row* shows estimated signal-to-noise maps (see text for detail), with a stretch between [-3σ, 3σ].</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** *From left to right:* observations, best-fit model and residuals of the ALMA data. For all panels the color map is a linear stretch between -0.27 and 1.23 mJy/beam. The standard deviation estimated in an empty region of the observations is of 0.09 mJy/beam. For both the observations and the model the contours are set at [3, 5, 7.5, 10]σ, and [-σ, σ] for the residuals (no residuals beyond the 2σ level). The beam size is shown in the lower left corner of each panel.</figcaption></figure>

211.97–230.99 GHz. The target was observed several times, with precipitable water vapor ranging from 5.19 to 0.79 mm. We only kept the observations performed on the 20 of March 2014 for which the water vapor was minimum. Out of the four spectral windows, three were used to derive the continuum emission (211.97–228.97 GHz, with a 31 250 kHz resolution), while the last one was used to search for CO gas emission (230.05–230.99 GHz, with a 488.28 kHz resolution), which was not detected. Data processing was performed within CASA using the standard scripts provided by the observatory. We only

kept the spectral windows used for the continuum observations and averaged the complex visibilities along the 128 different spectral channels (while flagging points with negative weights). Figure C.2 shows the (*u*, *v*) plane coverage for the continuum observations, with minimum and maximum baselines of 11.8 and 334.9 m, respectively. The reconstructed image (with so-called briggs weighting and pixel size of 0.13") is shown in the left panel of Fig. 2 (sensitivity of 0.09 mJy/beam). The beam size is 1.36" × 0.73" (48 au × 26 au) with a position angle of -86.5°. We do not attempt to measure the total flux of the disk directly
<table><caption>**Table 2.** Broadband photometric measurements of HD 61005, and the equivalent widths of the far-IR filters (see text for details).</caption><thead><tr><th>λ<br>[μm]</th><th>F<sub>ν</sub><br>[mJy]</th><th>σ<br>[mJy]</th><th>EW<br>[μm]</th><th>Instrument</th></tr></thead><tbody><tr><td>0.428</td><td>895.17</td><td>14.02</td><td></td><td>TYCHO B</td></tr><tr><td>0.534</td><td>1810.23</td><td>18.34</td><td></td><td>TYCHO V</td></tr><tr><td>1.235</td><td>2753.74</td><td>65.94</td><td></td><td>2MASS J</td></tr><tr><td>1.662</td><td>2440.48</td><td>103.40</td><td></td><td>2MASS H</td></tr><tr><td>2.159</td><td>1738.75</td><td>38.43</td><td></td><td>2MASS Ks</td></tr><tr><td>3.353</td><td>819.28</td><td>31.69</td><td></td><td>WISE W1</td></tr><tr><td>4.603</td><td>453.05</td><td>8.76</td><td></td><td>WISE W2</td></tr><tr><td>11.56</td><td>78.40</td><td>1.08</td><td></td><td>WISE W3</td></tr><tr><td>22.09</td><td>44.28</td><td>1.55</td><td></td><td>WISE W4</td></tr><tr><td>68.92</td><td>717.00</td><td>5.33</td><td>21.41</td><td>PACS Blue</td></tr><tr><td>97.90</td><td>703.58</td><td>6.84</td><td>31.29</td><td>PACS Green</td></tr><tr><td>153.94</td><td>472.65</td><td>14.58</td><td>69.76</td><td>PACS Red</td></tr><tr><td>251.50</td><td>235.6</td><td>13.5</td><td>67.61</td><td>SPIRE PSW</td></tr><tr><td>352.83</td><td>118.9</td><td>7.6</td><td>95.75</td><td>SPIRE PMW</td></tr><tr><td>511.60</td><td>49.8</td><td>5.2</td><td>185.67</td><td>SPIRE PLW</td></tr><tr><td>1300.0<sup>a</sup></td><td>4.6</td><td>0.7</td><td>105.40</td><td>ALMA Band 6</td></tr></tbody></table>

**Notes.** (a) Results from the modeling of the ALMA data (Sect. 3).

from these observations, but we do estimate it when modeling the complex visibilities (Sect. 3).

## 2.3. Spectral energy distribution

The star HD 61005 was observed by Herschel (Pilbratt et al. 2010) with the Photodetector Array Camera and Spectrometer instrument (PACS, Poglitsch et al. 2010), within the program OT2_tcurrie_1. The observation numbers (OBSID) are the two pairs 1342270977, 1342270978 and 1342270979, 1342270980 for the 70 µm and 100 µm observations, respectively. The 160 µm map used the four OBSID combined. The data were processed using the HIPE software (build 12.0.2083, Ott 2010), the very same way as described in Olofsson et al. (2013). HD 61005 was also observed with the Spectral and Photometric Imaging Receiver instrument (SPIRE, Griffin et al. 2010) in small scan map mode (OBSID: 1342255147 within the program OT2_kstape01_1). We used the Timeline Fitter task in HIPE to derive SPIRE photometry for our target. Calibration errors (~5.5%, Bendo et al. 2013) are included in the uncertainties. We also gathered photometric observations using VOSA<fnref n="1" /> (Bayo et al. 2008) and the dataset used to build the SED can be found in Table 2. The meaning of the third column is explained in Sect. 5.

Finally, we downloaded the *Spitzer/IRS* spectrum from the Cornell Atlas of *Spitzer/IRS* Sources database<fnref n="2" /> (Lebouteiller et al. 2011).

## 2.4. Stellar parameters

The stellar photospheric model is taken from the ATLAS9 Kurucz library (Castelli et al. 1997) with an effective temperature of $T_{\star} = 5500$ K (Casagrande et al. 2011). With the dilution factor used to scale the photospheric model to the optical and near-IR

photometric measurements, at a distance of 35.4 pc, we find a radius $R_{\star} = 0.84 R_{\odot}$. We derived a luminosity of $L_{\star} = 0.58 L_{\odot}$. To derive the stellar mass, which will become important when discussing the dust properties and the effect of radiation pressure on dust grains, we use isochrones from Siess et al. (2000), for an age of 40 Myr and effective temperature of 5500 K. We find that the stellar mass must be of about 1.1 $M_{\odot}$ (the corresponding luminosity matching our estimated $L_{\star}$). We find a slightly smaller mass (1 $M_{\odot}$) when using the isochrones from Baraffe et al. (2015), but the differences may arise from different model prescriptions (e.g. overshooting). In the following, we adopt a mass of 1.1 $M_{\odot}$. The SED with the broadband photometric measurements, the *Spitzer/IRS* spectrum as well as the photospheric model are shown in Fig. 7 of Sect. 5.

## 2.5. Preamble on the modeling strategy

In this study, we aim to model observations from different facilities, at different wavelengths, using different techniques (interferometry and direct imaging). Therefore, before detailing the modeling strategy for each individual dataset, we provide a quick preamble on the methodology.

We first model the ALMA data (Sect. 3) assuming a circular disk, fitting the reference radius ($r_0$, where the dust density peaks), the outer slope for the dust density distribution ($\alpha_{out}$), the position angle ($\phi$), the inclination ($i$), and the total flux at 1.3 mm ($f_{1300}$).

Prior to the modeling of the SPHERE observations, we attempt to constrain some of the dust properties, to limit the number of free parameters. We use the best fit results for the inclination and position angle from the modeling of the ALMA data to derive the polarized intensity as a function of the azimuthal angle from the $Q_{\phi}$ image. We constrain the minimum and maximum grain sizes ($s_{min}$ and $s_{max}$, respectively) as well as the porosity fraction of the dust grains (Sect. 4.1). This enables us to reduce the pool of free parameters when modeling the SPHERE DPI observations. For the SPHERE ADI observations, we use the Henyey-Greenstein approximation for the phase function, which disconnects the modeling process from the aforementioned dust properties. Thanks to the great complementarity between the ADI and DPI observations, we model them simultaneously to best constrain the azimuthal and radial dust density distribution (Sect. 4).

Finally, in Sect. 5, we model the SED of HD 61005, using the results inferred from the modeling of the SPHERE data on the location of the disk to derive stringent constraints on the dust properties (minimum grain size, dust composition, and total dust mass).

# 3. The parent planetesimal belt: constraints from the ALMA observations

Roughly speaking, different wavelengths trace different grain sizes. Therefore, we chose to model the ALMA observations independently of the SPHERE ones, the latter probing the small dust grains in the debris disk while the millimeter observations most likely trace a population of larger grains that more closely follow the parent planetesimals' belt.

## 3.1. Modeling strategy

The modeling of the ALMA observations is performed in the Fourier space, attempting to reproduce both the real and

<footnote n="1">¹ http://svo2.cab.inta-csic.es/theory/vosa/</footnote>

<footnote n="2">² The Cornell Atlas of *Spitzer/IRS* Sources is a product of the Infrared Science Center at Cornell University, supported by NASA and JPL. http://cassis.sirtf.com/atlas/query.shtml</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** One- and two-dimensional (diagonal and lower triangle, respectively) projections of the posterior probability distributions for the results of the modeling of the ALMA observations.</figcaption></figure>

imaginary parts of the complex visibilities (averaged along the spectral dimension). In Appendix A, we explain how we generate synthetic images and the different notations are summarized in Table C.1. From a synthetic image at the wavelength of 1.3 mm, we first scale the total flux of the image to the free parameter $f_{1300}$ (in mJy) before computing the Fourier transform of the image. We then interpolate the Fourier transform at the spatial frequencies of the observations. The goodness of fit is the sum of the weights (estimated in CASA<fnref n="3" />) times the squared difference between the observed and modeled complex visibilities (the weights are proportional to $1/\sigma^2$). We consider the following free parameters: the inclination $i$, the position angle $\phi$, the total flux of the disk at 1.3 mm $f_{1300}$, the reference radius<fnref n="4" /> $r_0$, and the outer power-law slope for the dust distribution $\alpha_{\text{out}}$, which is parametrized as

$$n \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}, \quad (3)$$

where $r$ is the distance from the star and $n$ the number density. The inner power-law slope is set to $\alpha_{\text{in}} = 5$, preliminary tests indicating this parameter is poorly constrained by the observations. To find the most probable solution, we use an affine invariant ensemble sampler Monte-Carlo Markov Chain, implemented in the `emcee` package, using 200 walkers, a burn-in phase of 500 iterations and a total length of the chains of 2000 iterations after the burn-in phase. At the end of the run, we find that the mean acceptance fraction (the mean fraction of steps accepted for each walker within the chain) is of 0.48 (a good sign of convergence and stability, Gelman & Rubin 1992). The maximum auto-correlation time for all the parameters is of 60 steps, indicating that the chains should have stabilized by the end of the simulations.

<table><caption>**Table 3.** Best fit results for the modeling of the ALMA observations.</caption><thead><tr><th>Parameter</th><th>Uniform prior</th><th>σ<sub>kde</sub></th><th>Best-fit value</th></tr></thead><tbody><tr><td>r<sub>0</sub> [au]</td><td>[40, 80]</td><td>0.2</td><td>66.4<sup>+6.1</sup><sub>-8.7</sub></td></tr><tr><td>α<sub>out</sub></td><td>[-15, -1.5]</td><td>0.1</td><td>-6.6<sup>+1.4</sup><sub>-6.1</sub></td></tr><tr><td>φ [°]</td><td>[55, 85]</td><td>0.1</td><td>70.7<sup>+1.9</sup><sub>-2.3</sub></td></tr><tr><td>i [°]</td><td>[79, 89]</td><td>0.1</td><td>84.5<sup>+2.9</sup><sub>-2.5</sub></td></tr><tr><td>f<sub>1300</sub> [mJy]</td><td>[1, 20]</td><td>0.1</td><td>4.6<sup>+0.7</sup><sub>-0.6</sub></td></tr></tbody></table>

## 3.2. Results

The projected posterior probability distributions are displayed in Fig. 3, for the different free parameters (using the `triangle` Python package, Foreman-Mackey et al. 2014). To derive the best-fit values as well as the uncertainties, we smooth the distributions with a kernel density estimator (the width of the Gaussian kernel $\sigma_{kde}$ are reported in Table 3), and the best-fit value is the peak position of the distribution. The confidence intervals ($a_1, a_2$) for the parameter $a$ are estimated as follows:

$$\int_{a_{\text{min}}}^{a_1} p(a) \mathrm{d}a = \int_{a_2}^{a_{\text{max}}} p(a) \mathrm{d}a = \frac{1 - \gamma}{2}, \quad (4)$$

where $\gamma = 0.68$ and $p(a)$ is the smoothed posterior probability distribution (integral normalized to 1) for parameter $a$ (e.g., Pinte et al. 2008). We note that not all the distributions reach zero on each side of their maximum (especially for $\alpha_{\text{out}}$ and $i$) and, therefore, these uncertainties should be treated carefully. Table 3 summarizes our results for the modeling of the Band 6 observations, and Fig. 2 shows the observations, the best-fit model, and the residuals (from left to right). The synthetic image of the best fit model is processed through CASA (using the `ft` method with the same antenna configuration as the observations) and the image is reconstructed with the `clean` algorithm with the same parameters as for the observations. We note that there are some residuals on the east side that may suggest that the disk is brighter on one side, even at mm wavelength. However, these residuals are below $3\sigma$, therefore we cannot conclude they are significant. The apparent brightness asymmetry could be due to the asymmetric $(u, v)$ coverage of the observations.

Overall, we find that most of the parameters are well constrained, except for the outer power-law slope of the dust density distribution, for which we can safely exclude slopes shallower than $\alpha_{\text{out}} = -4$. This is explained by the beam size of the observations which is larger than the debris disk for steep values of $\alpha_{\text{out}}$. Otherwise, we find the reference radius of the disk to be $r_0 \sim 66$ au, the position angle $\phi \sim 70.7^\circ$, the inclination $i \sim 84.5^\circ$, and the flux at 1.3 mm $f_{1300} \sim 4.6$ mJy. These results agree well with the parameters reported in Buenzli et al. (2010, $i = 84.3 \pm 1^\circ$, $\phi = 70.3 \pm 1^\circ$, $r_0 = 61.25 \pm 0.85$ au) and Ricarte et al. (2013, $r_0 = 67 \pm 2$ au, $\phi = 71.5 \pm 5^\circ$). The relatively large beam size of the ALMA observations can explain the slight discrepancy for $r_0$ between the modeling of the ALMA data and the value inferred by Buenzli et al. (2010). This value will be revisited when modeling the SPHERE observations (Sect. 4). Ricarte et al. (2013) obtained a total flux of $7.2 \pm 0.3$ mJy at 1.3 mm with their SMA observations (Steele et al. 2016 obtained $8.0 \pm 0.8$ mJy analyzing the same observations), while we find the total flux to be well constrained at $4.6 \pm 0.7$ mJy (within the SMA and ALMA respective $3\sigma$ uncertainties). The shortest baselines being of $B_{\text{min}} \sim 11$ and 16 m (for the ALMA and SMA observations, respectively), the largest scales the observations

<footnote n="3">³ The absolute values of the weights derived within CASA may be inaccurate but their relative values are not.</footnote>

<footnote n="4">⁴ Here we assume the disk is circular.</footnote>
are sensitive to are of the order of 14.2" and 10.1", respectively (0.6λ/Bmin), much bigger than the disk. It is therefore unlikely that flux from the disk is filtered out by the interferometers. The differences between the SMA and ALMA data may arise from missing frequencies, differences in beam sizes, or calibration uncertainties (the uncertainties reported by Steele et al. 2016 being more conservative than the ones reported by Ricarte et al. 2013).

Finally, we note that the famous wings responsible for the disk's nickname are not detected in the ALMA data. Despite a good angular resolution, a disk model (without wings) can successfully reproduce the observations, and we see no trace of the wings in the residuals (with an rms of ~0.09 mJy/beam). Our results therefore agree with the ones of Ricarte et al. (2013); only small dust grains are likely present in the wings.

## 4. Constraining the dust radial distribution from the SPHERE observations

We model the SPHERE images by producing synthetic images at the central wavelength of the H-band observations, $\lambda_c = 1.63~\mu\text{m}$. Given the low S/N of the K-band ADI observations, preliminary attempts to model these data showed that the dust distribution cannot be better constrained than with the H-band observations. We therefore focus the modeling effort on the H-band ADI and DPI data.

Figure 1 highlights the complementarity of both the scattered and polarized light images; the ADI and DPI data have very different S/N at the ansae and along the semi-minor axis of the disk. Combining both datasets, therefore, offers the opportunity to study the dust distribution in great detail.

The modeling process has a high dimensionality with many possible free parameters, and regions of intermediate to low S/N. Therefore, to obtain novel yet reliable constraints on the dust distribution we choose to perform a prior analysis on the observations to reduce the number of free parameters. Having a proper description of the polarized phase function prior to the modeling of the DPI observations greatly helps reducing the dimensionality of the modeling (e.g., the minimum and maximum grain sizes as well as the porosity of the dust grains). In this section, we first describe how the polarized phase function is derived and modeled, then we present the modeling strategy and summarize the results we obtain.

### 4.1. Phase function of the polarized light

To compute the polarized phase function from the DPI dataset, we define an elliptical mask with the following parameters: inner and outer radii ($r_{in}$ and $r_{out}$), the inclination $i$, and the position angle $\phi$. For each pixel within the elliptical mask, we compute the scattering angle as the dot product of the line of sight and the location of the pixel with respect to the star. We divide the elliptical mask in two, for the east and west sides and, for each side, we compute the minimum and maximum scattering angles. We then divide each side of the mask into 30 smaller regions corresponding to different bins of the phase function (see Fig. C.3 for an illustration). Each pixel is multiplied by its squared distance to the star, to account for illumination effect. The measured phase function is found by averaging the flux in the observations, in each individual region. Since the observed uncertainties $\sigma_i$ are not the same for each pixel within a given intersection, the “average”

<figure><img src="image_5.png" /><figcaption>Fig. 4. Azimuthal dependency of the polarized intensity in the DPI H-band observations, in black and red for the East and West sides, respectively. The shaded curves are representative of the uncertainties estimated from the noise map and the shaded gray area denotes low S/N regions estimated from the S/N map in Fig. 1. The thin solid lines display examples for different grain sizes for the same dust composition as the best fit.</figcaption></figure>

uncertainty in the intersection is computed as follows:

$$ \sigma = \left( \sum_{i=1}^{N} \frac{1}{\sigma_i^2} \right)^{-1/2}, \qquad (5) $$

where $N$ is the number of pixels in the considered region. The phase function is then normalized to its maximum value (east and west sides are divided by the same value).

For the elliptical mask, we use the best-fit results of the modeling of the ALMA observations (see Table 3) for $i$ and $\phi$ and choose $r_{in}$ and $r_{out}$ small and large enough (50 and 72 au, respectively) so that they encompass the trace of the debris disk. Figure 4 shows the phase function for the DPI H-band observations, the east side being in black, and the west side in red. The uncertainties are shown in a shaded color, and the gray area indicates regions of low S/N, which were estimated from the S/N map derived from the $U_{\phi}$ data. As can be seen in Fig. 1 (bottom right panel), the disk is only detected for a narrow range of azimuthal angles on the west side. This strongly suggests that we underestimated the uncertainties calculated using Eq. (5).

The difference in brightness between the two sides of the disk is striking in Fig. 4 and the east side appears almost twice as bright as the west side for most of the scattering angles. Even though the S/N for the west side is relatively low it seems that the polarized phase function peaks at an angle compatible with the east side.

To alleviate the number of free parameters, we aim to constrain the grain size distribution for the modeling of the DPI observations directly from the phase function displayed in Fig. 4. We emphasise that we do not have access to the polarization degree $(\sqrt{Q^2+U^2})/I$ as we do not have an unbiased measure of the total intensity $I$ prior to the modeling. Indeed, the ADI process introduces self-subtraction effects (e.g., Milli et al. 2012), which can eventually be quantified with a model that describes the observations well (which we do not yet have). Since the west side suffers from low S/N, we perform the modeling on the east side's phase function. To reproduce the phase function, we assume that the signal in the $Q_{\phi}$ image is proportional to the size dependent
$S_{12}(s)$ element of the Müller matrix. The matrix enables us to compute the Stokes vectors $I$ and $Q$ for the scattered and polarized light, respectively. Assuming single scattering event (reasonable assumption in low density environment such as debris disks), the scattered light will be the product of the first diagonal element of the matrix $S_{11}$ times the stellar intensity $I_0$. The polarized intensity will be proportional to the second element of the first column of the matrix $S_{12}$ times $I_0$. We compute $S_{12}$ for different grain sizes between $s_{min}$ and $s_{max}$, using the Mie theory, and we average $S_{12}$ over a grain size distribution with a slope $p < 0$ ($dn(s) \propto s^p ds$) as follows,

$$S_{12}^{\text{avg}} = \frac{\int_{s_{\text{min}}}^{s_{\text{max}}} S_{12}(s) \times s^p ds}{\int_{s_{\text{min}}}^{s_{\text{max}}} s^p ds}. \qquad (6)$$

We then compute the best scaling factor $f_{S_{12}}$ (to be multiplied to $S_{12}^{\text{avg}}$) that will minimize the difference between the profiles $S_{12}^{\text{avg}}$ and $S_{12}^{\text{obs}}$ (with uncertainties $\sigma$)

$$f_{S_{12}} = \frac{\sum \left( \frac{S_{12}^{\text{obs}} \times S_{12}^{\text{avg}}}{\sigma^2} \right)}{\sum \left( \frac{S_{12}^{\text{avg}}}{\sigma} \right)^2}. \qquad (7)$$

We perform a simple grid search over the following parameters: $s_{min}$, $s_{max}$, and the optical properties of the dust grains. We fix $p = -3.5$, as expected for a collisional cascade in a debris disk (Dohnanyi 1969). For the dust composition, we consider a base medium of amorphous silicate with olivine stoichiometry (MgFeSiO$_4$, Dorschner et al. 1995) to which we can add some porosity using the Bruggeman mixing rule. The minimum grain size can vary between 0.01 and 10 µm. To constrain the maximum grain size, we vary the quantity $\Delta s$ ($=s_{max} - s_{min}$) between 0.01 and 100 µm (in log space), and finally the porosity can change by steps of 10%. We find that the polarized phase function at 1.63 µm is best reproduced by small spherical dust grains, with typical sizes in the range $0.3 \le s \le 0.35$ µm and a porosity fraction of ~80%. The best-fit solution is shown with a thick cyan line in Fig. 4, and reproduces well both the overall shape and the peak position of the observed phase function. Also shown in Fig. 4 are two examples for different grain size distributions around 0.25–0.3 µm and 0.35–0.4 µm, to illustrate how sensitive the phase function is with respect to the grain sizes. We note that, within such a narrow range of sizes, the value of the slope $p$ of the grain size distribution remains unconstrained. In this section, we assumed that the polarized flux is directly proportional to $S_{12}$, while it is also related to the dust density in the disk. In the rest of this paper, we try to determine the azimuthal distribution of the dust in the disk. Therefore, this prior analysis must be regarded as a first order approximation. The motivation of this prior analysis was to reduce the dimensionality of the modeling, but it also comes at a slight cost in the interpretation of the grain size distribution. The main conclusion of this section is that it does not seem that we observed large dust grains (which would have a stronger $S_{12}$ signal at small scattering angles). In Sect. 6.1 we discuss this result further, but it could also be the consequence of low S/N along the semi-minor axis of the disk.

## 4.2. Modeling strategy

Since both the ADI and DPI datasets were taken at the same wavelength, for a given set of parameters we compute

two images: an unpolarized light image using the Henyey-Greenstein (HG) analytical prescription of phase function $S_{11}$ (Henyey & Greenstein 1941) and a polarized light image using the Mie theory. Using the HG approximation, which is parametrized by the anisotropic scattering factor $g$ ($-1 \le g \le 1$), gives us more control when trying to reproduce the observed phase function (hence less free parameters to be considered during the modeling). To reproduce the DPI observations, we use the grain properties derived previously (Sect. 4.1). The absorption and scattering efficiencies, as well as the Müller matrix element $S_{12}$, are computed with the Mie theory, which is valid for compact spherical grains. This approach may not appear self-consistent, but it is a way to disentangle the modeling of unpolarized and polarized light that may not be well accounted for by spherical grains (e.g., Milli et al. 2015).

The pool of free parameters includes the reference radius $r_0$, the inclination $i$, the position angle $\phi$, the opening angle $\psi$, and the outer power-law slope $\alpha_{out}$ for the dust density distribution. Preliminary tests indicated that we can hardly constrain the inner power-law slope, so we fixed $\alpha_{in} = 5$. This enables us to focus on other parameters that describe the geometry of the debris disk. For instance, Buenzli et al. (2010) conclude that the eccentricity of the disk was not enough to explain the brightness asymmetry between the east and the west sides, and we aim to address this interpretation of previous observations. Therefore, we include the eccentricity $e$, the rotation angle $\Phi_e$, the density damping $\eta$, and its azimuthal shape (via the width $w$ of the Gaussian profile) and its reference angle $\Phi_\eta$ (see Appendix A). The azimuthal profile has the shape of a Gaussian profile with a $\sigma = w$, a peak of 1 for the azimuthal angle $\Phi_\eta$ and a minimum value of $\eta \ge 0$ (hence an amplitude of $1 - \eta$).

Because we now consider azimuthal variations for the dust density distribution, it is parametrized slightly differently. The semi-major axis of the disk is defined as $r_0/(1 - e^2)$. Although the formal denomination of $r_0$ is the so-called semi latus rectum, we will simply refer to it as the reference radius in the rest of the analysis. The radius at which the dust density peaks now depends on the azimuthal angle $\theta$. This radius $r_\theta$ is defined as $r_0/(1 + e \times \cos\theta)$, and the azimuth-dependent dust density distribution $n(\theta)$ is parametrized similarly to Eq. (3), replacing $r_0$ with $r_\theta$ for each azimuthal angle.

The dust mass $M_{\text{dust}}$ is not varied, but the synthetic images are scaled during the fitting process. For the polarized images, the modeled image is the absolute value of the Stokes $Q_\phi$ parameter (we are assuming that there are no multiple scattering events in the debris disk). We scale the modeled image by a factor $f_{\text{DPI}}$, which is found analytically to minimize the residuals (similarly to Eq. (7)). For the ADI dataset, this approach is not possible because the post-processing of the data cube can introduce self-subtraction effects. We therefore opted for a forward modeling strategy, similar to the one described in Thalmann et al. (2014). For a given set of parameters, we compute one synthetic image at the wavelength 1.63 µm. We then produce a cube of 64 images, each one rotated to match the parallactic angles of each frame of the observations (total rotation of 93.3°). Each frame of the cube of synthetic images is multiplied by $f_{\text{ADI}}$ and is then subtracted from the corresponding observed frame. The PCA process is performed, keeping only the six main components. Another approach would be to perform the PCA on the observations, save the coefficients, and apply them to the modeled cube. But the observations contain signal from the disk that may be accounted for in some of the principal components (even though they should be similarly subtracted in the modeled cube). Therefore, to ensure as little subtraction as possible, we chose to perform the PCA on
<table><caption>**Table 4.** Best fit results for the ADI and DPI H-band observations.</caption><thead><tr><th>Parameter</th><th>Uniform prior</th><th>σ<sub>kde</sub></th><th>Best-fit value</th></tr></thead><tbody><tr><td>r<sub>0</sub> [au]</td><td>[40, 80]</td><td>0.1</td><td>60.4<sup>+0.8</sup><sub>-0.5</sub></td></tr><tr><td>i [°]</td><td>[75, 88]</td><td>0.1</td><td>84.1<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>α<sub>out</sub></td><td>[-10, -1.75]</td><td>0.01</td><td>-2.70<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>e</td><td>[0, 0.6]</td><td>0.0025</td><td>0.093<sup>+0.018</sup><sub>-0.014</sub></td></tr><tr><td>Φ<sub>e</sub> [°]</td><td>[0, 360]</td><td>1.</td><td>127.3<sup>+12.7</sup><sub>-4.3</sub></td></tr><tr><td>φ [°]</td><td>[65, 75]</td><td>0.1</td><td>70.6<sup>+0.2</sup><sub>-0.3</sub></td></tr><tr><td>η</td><td>[0, 1]</td><td>0.0025</td><td>0.47<sup>+0.02</sup><sub>-0.04</sub></td></tr><tr><td>Φ<sub>η</sub> [°]</td><td>[0, 360]</td><td>1.0</td><td>138.1<sup>+7.2</sup><sub>-5.6</sub></td></tr><tr><td>w [°]</td><td>[5, 180]</td><td>1.0</td><td>51.7<sup>+8.8</sup><sub>-3.2</sub></td></tr><tr><td>ψ</td><td>[0.02, 0.12]</td><td>0.005</td><td>0.058<sup>+0.001</sup><sub>-0.002</sub></td></tr><tr><td>g</td><td>[0, 0.95]</td><td>0.001</td><td>0.54<sup>+0.01</sup><sub>-0.02</sub></td></tr><tr><td>f<sub>ADI</sub></td><td>[3, 8.5]</td><td>0.01</td><td>6.20<sup>+0.09</sup><sub>-0.07</sub></td></tr></tbody></table>

the model-subtracted cube. The end goal being to minimize the flux in the final image.

The goodness of the fit is the sum of the squared ratio between the final images (residuals for the scattered and polarized data) and the noise map (bottom row of Fig. 1). For the modeling of both the ADI and DPI observations, we therefore have a total of 12 free parameters: $r_0$, $i$, $\alpha_{out}$, $e$, $\Phi_e$, $\phi$, $\eta$, $\Phi_\eta$, $w$, $\psi$, $g$, and $f_{ADI}$. Here, we also use the emcee Python package, with 100 so-called walkers in the chain, and first burn in 500 runs for each of these walkers. We then run the chain for 2000 iterations in total (see Millar-Blanchaer et al. 2015 for a similar modeling approach). To speed up the process, we cropped and re-sampled the SPHERE images. The size of one pixel in the new image (or cube) is resampled to be 2.25 times bigger than in the original images, the size of the new image being 200 × 200 pixels, and we use a central mask of radius 0.15". We chose not to convolve the synthetic images by a PSF because the cropped and down-sampled images have a pixel size of 0.028", while the approximation of the instrumental PSF with a Gaussian profile would have a width of 0.024" (approximating the PSF as an Airy disk for an 8 m telescope). At the end of the run, we find that the mean acceptance fraction (the mean fraction of steps accepted for each walker within the chain) is of 0.35, with a maximum auto-correlation time of 79 steps.

### 4.3. Results

The projected posterior probability distributions are displayed in Fig. B.1, for the different free parameters. To derive the best-fit values as well as the uncertainties, we proceed similarly as in Sect. 3.2. The results are presented in Table 4 and Fig. 5 shows the observations, the residuals and the best-fit models (from left to right) for the ADI and the DPI data (top and bottom rows, respectively).

All parameters seem to be well constrained. The most probable solution has a semi-major axis of 60.9 au ($r_0/[1-e^2]$) for an eccentricity of 0.093. With a rotation angle of $\Phi_e \sim 127^\circ$ the pericenter is located slightly toward the observer, on the east side<fnref n="5" />. An azimuthal density variation seems to be necessary to

reproduce the observations with a damping factor $\eta$ of ~0.47 with a reference angle of ~138° (hence almost co-located with the pericenter of the eccentric disk). The azimuthal variation of the dust density distribution is a Gaussian profile with a width of ~50°. The position angle and inclination are consistent with the results from Buenzli et al. (2010). The aspect ratio of $\psi \sim h/r \sim 0.06$ agrees well with numerical simulations of the vertical structures of debris disks (e.g., Thebault 2009). Finally, the phase function is anisotropic for low scattering angles with $g \sim 0.54$. The dust density distribution for the best-fit model, viewed from above the disk, is displayed in Fig. B.2.

Nonetheless, as shown in the residuals of the ADI observations, our best-fit model does not manage to remove all the signal along the semi-minor axis of the disk. It successfully suppresses most of the signal at larger scattering angles, but the best-fit model seems to fail at properly describing the scattering at smaller angles. We tried to implement a phase function with two weighted HG functions, but could not significantly improve the residuals. We included all the “basic” parameters related to the geometry of the disk ($i$, $\phi$, $r_0$, $\psi$) and yet failed to perfectly match the observations. Possible explanations can be related to the phase function (the HG phase function remains an approximation), the radial segregation of the grain size distribution, or the azimuthal dust density distribution (this will be further discussed in Sect. 6.1). We crudely assumed a Gaussian profile for the azimuthal distribution, but the actual distribution could be skewed in one or another direction which could explain the brighter region along the semi-minor axis. Another (highly speculative but interesting) explanation could be that we may be seeing an inner disk. Because of the high inclination of the disk, an inner dust belt may appear to the observer as if it was merging with the main belt. The main challenge with this explanation is that we see no indication of this type of a belt in the DPI observations, but the noise is larger in the innermost regions of this dataset. One possible way to address this point would be to detect gas, which could trace velocities compatible with a radius smaller than 61 au. Yet no CO was detected in the ALMA observations (Sect. 6.2).

## 5. Constraining the dust mineralogy from the SED

The SED of the debris disk around HD 61005 is constructed from the fluxes reported in Table 2 and the Spitzer/IRS spectrum. Modeling an SED from unresolved observations is a degenerate problem. We are basically trying to find the adequate temperature of the grains, which can be changed either by their radial distances, their sizes, or their nature. The modeling of the SPHERE images provide strong constraints on the geometric parameters of the disk, and we can therefore focus on better constraining the dust properties.

### 5.1. Modeling the thermal emission from the disk

To model the SED of HD 61005, for a given set of parameters, the goodness of fit is computed as

$$ \chi^2 = \sum_i \omega_i \times \left[ \frac{F_{\text{obs}}(\lambda_i) - F_{\star}(\lambda_i) - F_{\text{model}}(\lambda_i)}{\sigma_i} \right]^2, \quad (8) $$

where $F_{\text{obs}}$ is the observed flux with its associated uncertainty $\sigma_i$, $F_{\star}$ the stellar contribution, and $F_{\text{model}}$ the modeled thermal emission of the debris disk. The $\omega_i$ values are weights to the observed data points at wavelength $i$. This weighting is designed

<footnote n="5">⁵ For $\Phi_e = 180^\circ$, the pericenter would be along the semi-major axis on the east side. Smaller angles would move the pericenter towards the observer, while larger angles would move the pericenter towards the back side of the disk.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Left to right: observed, residuals, and best-fit models for the ADI and DPI datasets (top and bottom, respectively).</figcaption></figure>

to account for the fact that we simultaneously model broad-band photometric observations (e.g., Herschel/PACS) and spectroscopic observations (Spitzer/IRS): two adjacent points in the IRS spectrum do not have the same significance as PACS 70 and 100 µm observations. We therefore opt for a similar strategy to the one described in Ballering et al. (2013). We compute the average spectral resolution of the IRS data and compare it to the equivalent widths of the broadband filters. One IRS point will have a weight $\omega = 1$ and broadband photometric point will have a weight equal to the number of IRS points that would fit in the corresponding equivalent width of the filter. In Table 2, we report the equivalent widths for the various instruments used to model the thermal emission in the far-IR. These values were taken from the Spanish Virtual Observatory filter profile service<fnref n="6" />. For the ALMA Band 6 observations, we assumed the equivalent width to be equal to the spectral range (105 µm).

To account for different mineralogical components, we use amorphous silicate grains of olivine stoichiometry (MgFeSiO₄, Dorschner et al. 1995, $\rho = 3.5 \text{ g cm}^{-3}$) and amorphous water ices (Li & Greenberg 1998, $\rho = 1.2 \text{ g cm}^{-3}$). To mimic porosity, we also add the fraction of vacuum to the pool of free parameters. We mix the optical constants of the different dust components, using the Bruggeman mixing theory, and the Mie theory to compute the absorption coefficients $Q_{abs}$. Since computing $Q_{abs}$ for large grain sizes can be the bottle-neck when computing thousands of models, we first created a library of opacity files ($s_{min} = 0.01 \text{ µm}$ and $s_{max} = 5 \text{ mm}$, with steps of 10% for the

water ices, and porosity fractions). During the modeling process, the appropriate $Q_{abs}$ are drawn and interpolated for the proper grain size distribution and composition.

The free parameters include the grain size distribution ($s_{min}$ and $p$, we fix $s_{max} = 5 \text{ mm}$), the relative fractions for the fractions of amorphous water ices and porosity, and the inner and outer slopes for the dust density ($\alpha_{in}$ and $\alpha_{out}$, respectively). We fix the reference radius $r_0$ to the best-fit results of the modeling of the SPHERE observations. For each set of parameters the dust mass $M_{dust}$ is found by scaling the thermal emission of the model to the observed SED (similarly to Eq. (7)). The values for the dust mass are saved for posterior estimation of their corresponding uncertainties. To speed up the modeling of the SED, we neglect the non-azimuthal dust density distribution derived from the analysis of the SPHERE data and assume the disk to be a circular ring. We use the `emcee` package to search for the most probable model, using 100 walkers, burning the first 500 runs for each walker, and then running the chain for 1000 more iterations. The acceptance fraction at the end of the run is of 0.28 (and a maximum auto-correlation time across all parameters of 72 steps).

## 5.2. Results

Table 5 summarizes the best fit results and the derived uncertainties and projected posterior probability distributions for each parameter are displayed in Fig. 6. The SED of HD 61005, along with the photometric and spectroscopic observations, the photospheric model, and the best-fit model are shown in Fig. 7. The shape of the mid- to far-IR excess is well reproduced by

<footnote n="6">⁶ http://svo2.cab.inta-csic.es/theory/fps3/</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Posterior probability distributions for the modeling of the SED. The solid lines indicate the most probable value, and the vertical dashed lines indicate the derived uncertainties.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Spectral energy distribution of HD 61005. Photometric measurements are shown in red, along with the *Spitzer/IRS* spectrum (uncertainties, shown in gray, are $3\sigma$, most of them being smaller than the symbol). The dashed gray line is the photospheric model, the black and cyan lines are the best fit models (including the stellar contribution, but with and without the inner component, see text for details).</figcaption></figure>

<table><caption>**Table 5.** Best-fit results for the modeling of the SED.</caption><thead><tr><th>Parameter</th><th>Uniform prior</th><th>σ<sub>kde</sub></th><th>Best-fit value</th></tr></thead><tbody><tr><td>α<sub>out</sub></td><td>[-10, -0.5]</td><td>0.04</td><td>-1.92<sup>+0.06</sup><sub>-0.04</sub></td></tr><tr><td>p</td><td>[-5.0, -2.0]</td><td>0.02</td><td>-3.84<sup>+0.02</sup><sub>-0.05</sub></td></tr><tr><td>s<sub>min</sub> [μm]</td><td>[0.01, 10]</td><td>0.1</td><td>1.9<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>f<sub>ice</sub></td><td>[0, 0.95]</td><td>0.1</td><td>0.3<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>f<sub>porosity</sub></td><td>[0, 0.95]</td><td>0.1</td><td>0.3<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>M<sub>dust</sub> [10<sup>-1</sup> M<sub>☉</sub>]</td><td>—</td><td>0.1</td><td>3.7<sup>+0.4</sup><sub>-0.1</sub></td></tr></tbody></table>

our model, but the turn-off point, where the excess emission starts, is not well matched (near $\lambda \sim 20 \mu$m). This hints towards the presence of an additional component, an inner disk that has also been postulated by many authors in the literature (e.g., Morales et al. 2011; Ballering et al. 2013; Ricarte et al. 2013; Chen et al. 2014). Even when using detailed dust properties (instead of Planck functions at single temperatures) to model the excess, the best-fit model still underestimates the flux in the mid-IR beyond $3\sigma$. To constrain the properties of an additional belt, we fit a Planck function to the residuals (the only free parameter being the temperature, the scaling to the residuals being done similarly to Eq. (7)). We find the temperature of the Planck function to best reproduce the residuals to be ~220 K. In Fig. 7, the

Planck function is shown as a dotted-dashed black line and the final best-fit model (stellar contribution plus the inner and outer belts) as the solid black line. To assess the relevance of adding the inner disk we follow a similar approach as in Moór et al. (2015b, and references therein) and use the Akaike Information Criterion (AIC). We find that the addition of the inner component significantly improves the final model (the AIC including the inner belt is about 80% of the AIC without it), and therefore is deemed necessary to reproduce the whole SED of the disk around HD 61005.

Our modeling results suggest that the outer slope of the dust density distribution is relatively shallow, which would indicate that most of the grains responsible for the emission are located in an extended disk. The minimum grain size is found to be ~1.9 $\mu$m, and the grain size distribution has a steep slope of $p \sim -3.84$. The composition of the grains would be a mixture of amorphous silicate grains with small fraction of water ices and porosity. The total dust mass (between $s_{min}$ and $s_{max}$) is of the order of 0.37 $M_{\oplus}$.

# 6. Discussion

In the following, we try to put together the results of the modeling of the ALMA, SPHERE datasets as well as of the SED. We discuss what can be inferred regarding the dust properties, the gas content, the morphology of the disk to better characterize what is happening in this system. Because we modeled several datasets (ALMA, SPHERE, SED), Table 6 summarizes which parameters are constrained with which datasets, to clarify the discussion.

## 6.1. Dust properties

### 6.1.1. Self-subtraction corrected phase function

With a satisfying model for the scattered light, we can further investigate the phase function in the SPHERE ADI H-band observations. Indeed, the model enables us to estimate the self-subtraction introduced by the PCA analysis. To do so, we first perform the PCA on the raw data cube and save the resulting coefficients. We then take the best-fit model, produce a data cube at the proper parallactic angles, project it into the new orthogonal basis found for the observations, and multiply this projection with the same coefficients as for the raw observations (e.g., Soummé et al. 2012). This provides us with an estimate of the signal from the disk that is subtracted when doing the PCA. The attenuation is estimated by the ratio $r_{corr}$ between this estimate and the original model. When deriving the phase function for the ADI dataset (similarly to Sect. 4.1), each pixel is multiplied by the quantity $1/(1 - r_{corr})$. Figure 8 shows the phase function
<table><caption>**Table 6.** Summary of the modeling strategy adopted in this study.</caption><thead><tr><th>Parameter</th><th>ALMA data [1]</th><th>Polarized intensity [2]</th><th>SPHERE DPI [3]</th><th>SPHERE ADI [4]</th><th>SED [5]</th></tr></thead><tbody><tr><td>r<sub>0</sub> [au]</td><td>Fitted</td><td>—</td><td>Fitted</td><td>Fitted</td><td>Fixed to [3,4]</td></tr><tr><td>α<sub>in</sub></td><td>5</td><td>—</td><td>5</td><td>5</td><td>5</td></tr><tr><td>α<sub>out</sub></td><td>Fitted</td><td>—</td><td>Fitted</td><td>Fitted</td><td>Fitted</td></tr><tr><td>φ [°]</td><td>Fitted</td><td>Fixed to [1]</td><td>Fitted</td><td>Fitted</td><td>—</td></tr><tr><td>i [°]</td><td>Fitted</td><td>Fixed to [1]</td><td>Fitted</td><td>Fitted</td><td>—</td></tr><tr><td>f<sub>1300</sub> [mJy]</td><td>Fitted</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>s<sub>min</sub></td><td>0.1 µm</td><td>Fitted</td><td>Fixed to [2]</td><td>—</td><td>Fitted</td></tr><tr><td>s<sub>max</sub></td><td>5 mm</td><td>Fitted</td><td>Fixed to [2]</td><td>—</td><td>5 mm</td></tr><tr><td>f<sub>porosity</sub></td><td>0</td><td>Fitted</td><td>Fixed to [2]</td><td>—</td><td>Fitted</td></tr><tr><td>e</td><td>0</td><td>—</td><td>Fitted</td><td>Fitted</td><td>—</td></tr><tr><td>Φ<sub>e</sub> [°]</td><td>—</td><td>—</td><td>Fitted</td><td>Fitted</td><td>—</td></tr><tr><td>η</td><td>1</td><td>—</td><td>Fitted</td><td>Fitted</td><td>—</td></tr><tr><td>Φ<sub>η</sub> [°]</td><td>—</td><td>—</td><td>Fitted</td><td>Fitted</td><td>—</td></tr><tr><td>w [°]</td><td>—</td><td>—</td><td>Fitted</td><td>Fitted</td><td>—</td></tr><tr><td>ψ</td><td>0.05</td><td>—</td><td>Fitted</td><td>Fitted</td><td>—</td></tr><tr><td>g</td><td>—</td><td>—</td><td>—</td><td>Fitted</td><td>—</td></tr><tr><td>f<sub>ADI</sub></td><td>—</td><td>—</td><td>—</td><td>Fitted</td><td>—</td></tr><tr><td>p</td><td>-3.5</td><td>-3.5</td><td>-3.5</td><td>-3.5</td><td>Fitted</td></tr><tr><td>f<sub>ice</sub></td><td>0</td><td>0</td><td>—</td><td>—</td><td>Fitted</td></tr><tr><td>M<sub>dust</sub> [M<sub>☉</sub>]</td><td>—</td><td>—</td><td>—</td><td>—</td><td>Fitted</td></tr></tbody></table>

**Notes.** Even though the SPHERE DPI and ADI data are fitted simultaneously, we separated the entries to better see dependencies. A dash denotes that the parameter is irrelevant during the fitting of a given dataset.

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Phase function derived from the ADI H-band observations, corrected for self-subtraction effects (see text for details).</figcaption></figure>

derived for the east and west sides (black and red, respectively). The bump at ~120° on the west side is an artifact introduced by $r_{corr}$ (low signal in the disk model leads to an erroneous correction factor). In Fig. 8, we also show the phase function derived during the modeling (for $g = 0.55$), scaled to the phase function of the east for scattering angles larger than 40°. The inferred HG function appears to be a relatively good representation of the corrected phase function for scattering angles larger than ~50°, but severely fails to reproduce it at smaller angles. This would be an explanation for the strong residuals along the semi-minor axis in Fig. 5. Overall, this suggests quite extreme forward scattering in the disk around HD 61005 – a result that remains model-dependent since the PCA attenuation is estimated from our results (see also Milli et al. 2016, for a similar discussion

for the disk around HR 4796 A). Nonetheless, the strong peak of scattering for small angles (down to ~10°) strongly points towards the presence of large dust grains in the debris disk. However, this seems in contradiction with the prior analysis of the DPI dataset where we had to actually reject the presence of large dust grains ($s \sim 0.3 \mu m$). For the SED, we need a whole range of sizes (from $\mu m$- to mm-sized grains), but the emission is actually dominated by the smaller grains (discussed further in this section). Numerous studies of debris disks around different stars (e.g., Rodigas et al. 2015; Lebreton et al. 2012; Milli et al. 2015, 2016) already reported and discussed similar issues in reconciling the modeling results of different kind of observations. Throughout this paper, for their computational merits, we have either used the Mie theory or the HG prescription while both of them may not be accurate descriptions of the nature of the dust grains.

### 6.1.2. Reconciling different observations

To try to reconcile the different dust properties inferred from our modeling results, the radial extent of the disk is interesting to discuss. The ALMA observations point towards a narrow belt ($\alpha_{out} \le -5$), while the SED, ADI, and DPI observations suggest that the disk is extended ($\alpha_{out} \sim -2$ or $-3$). To assess the typical grain sizes we probe with the SED, we computed the integrated infrared luminosity for each grain size bins between $s_{min}$ and $s_{max}$, and built the cumulative distribution function. We find that more than 98% of the disk infrared luminosity is accounted for by dust grains smaller than $5 \mu m$. The rest of the (very steep) grain size distribution contributes marginally to the total emission. This means that for the DPI observations and the SED, we are observing an extended disk containing small dust grains. We computed the $\beta$ ratio between the radiation pressure and gravitational forces as in Burns et al. (1979), assuming $L_* = 0.58 L_⊙$ and $M_* = 1.1 M_⊙$. For the dust composition found from the SED modeling, we find that grains larger than
$s_{\text{blow}} \sim 0.8 \, \mu\text{m}$, for which $\beta \le 0.5$, should remain on bound orbits (assuming they were released on circular orbits). There is a small discrepancy between $s_{\text{blow}}$ and the minimum grain size we found ($s_{\text{min}} \sim 2 \, \mu\text{m}$) but as discussed in Pawellek et al. (2014), this could be caused by the assumption of spherical grains when modeling the SED. Assuming our inferred dust composition is not too far off compared to the real composition of the dust grains, it seems that radiation pressure should be efficient for (sub-) $\mu\text{m}$-sized grains. Consequently, the parent belt (inferred from the ALMA data) would be fairly narrow, but small grains placed on eccentric orbits (or even unbound) owing to radiation pressure make the disk look more extended at near-IR wavelengths.

The strong forward scattering peak derived from the corrected phase function from the ADI observations suggests the presence of large grains in an extended disk, which seems in contradiction with the aforementioned scenario. Nonetheless, according to Min et al. (2016), large dust aggregates behave like large grains (the size of the whole aggregate) in scattered light and like small grains (the size of the individual monomers) for polarized light. Such aggregates also display mild backward scattering that we do not detect with our observations but may be the consequence of low S/N along the back side of the disk.

Therefore, we hypothesize that the parent planetesimals are arranged in a narrow belt (traced by ALMA) and that the small-size end of the grain size distribution consists of dust aggregates that are pushed away by radiation pressure. This makes the disk appear extended in the SPHERE observations and when modeling the SED. For future studies of this system, it would be interesting to implement a radial segregation of the grain size distribution (varying $dn(s)$ as a function of $r$, see for instance Stark et al. 2014 for a discussion of the phase function depending on the radial distance).

## 6.2. Gas mass upper limits

No CO emission was detected in the ALMA data, and we follow a similar procedure to that of Matrà et al. (2015) and Moór et al. (2015a) to derive upper limits for the CO gas mass. We extract the spectrum from the datacube reduced within CASA, centered around 230.5 GHz ($^{12}$CO 2–1). The upper limit for the CO mass was determined as

$$M_{\text{CO}} = \frac{4\pi m d_{\star}^2}{h\nu_{2-1}A_{2-1}} \frac{S_{2-1}}{x_{2-1}}, \quad (9)$$

where $m$ is the mass of the CO molecule, $d_{\star}$ the distance of the star, $\nu_{2-1}$ the frequency of the transition, $A_{2-1}$ the Einstein coefficient$^7$, $x_{2-1}$ the fractional population of the upper level, and $S_{2-1}$ the observed integrated line flux. The fractional level was calculated assuming local thermal equilibrium, using the Boltzmann equation, with a temperature of 20 K (as discussed in Moór et al. 2015a, $M_{\text{CO}}$ is not strongly dependent on the temperature). $S_{2-1}$ was computed from the spectrum as $S_{\text{rms}}\Delta v \sqrt{N}$, where $\Delta v$ is the channel velocity width, $N$ the number of channels over a width of 15 km s$^{-1}$ (centered at the local standard of rest velocity, $v_{\text{LSR}}$), and $S_{\text{rms}}$ the standard deviation of the spectrum within that velocity range. From the heliocentric radial velocity of 22.5 km s$^{-1}$ (Desidera et al. 2011), we obtain a $v_{\text{LSR}}$ of 3.68 km s$^{-1}$. Assuming optically thin gas, this leads to an upper limit of $6.9 \times 10^{-7} M_{\oplus}$. As discussed in Matrà et al. (2015), the LTE hypothesis may not hold in low-density environment and

our upper limit could underestimate the CO gas mass. In the ISM the CO/H$_2$ abundance is $10^{-4}$ and with our inferred dust mass of $\sim 0.37 M_{\oplus}$ (for our grain size distribution), this would lead to a gas-to-dust ratio that is much smaller than unity. However, it is unlikely that this ratio will remain the same in a circumstellar disk because of effects such as photo-dissociation. Overall, with the available observations leading to a non-detection, and given the several assumptions made (LTE, CO/H$_2$ ratio) we can hardly conclude on the gas-to-dust ratio on the debris disk. This question should be addressed by future, deeper CO observations (or other species than the CO molecule).

## 6.3. Stirring of the planetesimals

In the last few years, theoretical works have aimed at characterizing the time evolution of debris disks (e.g., Kenyon & Bromley 2006, 2008). One of the purpose of these studies has been to better understand the growth of planetesimals to sizes of about a thousand km, in the framework of terrestrial planetary formation. It is believed that once these large bodies have been formed, they stir the population of smaller planetesimals. The stirring increases the relative velocities of these planetesimals and increases their chances of collisions, therefore initiating a collisional cascade. According to these models, collisions between planetesimals become destructive only once a Pluto-sized object is formed. Since the timescale for the growth of planetesimals scales with the distance $r$ to the star, self-stirring is thought to be an inside-out process. Here, we follow a similar approach to that described in Moór et al. (2015b), who studied several debris disks spatially resolved with the Herschel observatory. They examined if the disks in their sample are consistent with a self-stirring scenario. To achieve this, they compared the timescale for forming a 1000 km-sized object in a self-stirred debris disk with the age of the stars. To quantify this timescale $t_{1000}$ (in Myr), we used Eq. (41) from Kenyon & Bromley (2008), which we recall below:

$$t_{1000} = 145 x_{\text{m}}^{-1.15} (r/80)^3 (2M_{\odot}/M_{\star})^{3/2} \text{ [Myr]}, \quad (10)$$

where $r$ is the distance of the dust belt, and $x_{\text{m}}$ is a scaling factor to parametrize the disk’s initial surface density (the higher, the more massive). Mustill & Wyatt (2009) argue that $x_{\text{m}}$ values larger than 10 are highly unlikely since the initial disk would have been gravitationally unstable. Assuming that the timescale for the growth of the planetesimals is the age of the system (i.e., $t_{1000} = 40$ Myr), we find that $x_{\text{m}}$ should be of the order of about 3.3 (for $r = 61$ au). This value would suggest that either the primordial disk was $\sim 3$ times more massive than the minimum-mass solar nebula, or that other sources of stirring (e.g., induced by a planet) need to be taken into consideration for this system. Debris disks with known planets, such as the ones studied in Moór et al. (2015b, e.g., HR 8799, $\beta$ Pictoris, HD 95086), have much larger $x_{\text{m}}$ values ($\ge 10$), while most of the debris disks usually have values for $x_{\text{m}}$ smaller than 3. It is therefore reasonable to assume that the disk around HD 61005 is bright because self-stirring by planetesimals has recently started. It does not imply that alternative stirring mechanisms should not be considered, but based on the distance of the dust belt and the age of the system, it is not mandatory to invoke the presence of a massive planet to explain what is currently observed.

## 6.4. An eccentric and asymmetric disk

With the new SPHERE observations, we confirm the findings of Buenzli et al. (2010) that the disk displays a strong brightness

<footnote>$^7$ Taken from the SPLATALOGUE catalog at http://www.cv.nrao.edu/php/splat/</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** *From left to right:* observations, models including some of the results of the SPHERE modeling, and residuals. *Top row* is for $\eta = 0.47$ and the *bottom row* for $\eta = 0.6$. For all panels the color map is a linear stretch between -0.27 and 1.23 mJy/beam. For both the observations and the model the contours are set at [3, 5, 7.5, 10]$\sigma$, and [-2, -1, 1, 2]$\sigma$ for the residuals. The beam size is shown in the lower left corner of each panel.</figcaption></figure>

asymmetry. Both the ADI and DPI observations show that the eastern side is brighter than the western side at an unprecedented angular resolution. The asymmetry is detected both in scattered light along the semi-minor axis (ADI dataset) and in scattered and polarized light along the semi-major axis (NACO observations of Buenzli et al. 2010 and DPI observations, respectively). When modeling the observations, we had to include a density damping factor along the azimuthal direction to reproduce the SPHERE observations. From preliminary tests, the DPI dataset (in which the brightness asymmetry is most visible) cannot be properly reproduced without these azimuthal variations. Therefore, we argue it can hardly be a scattering artifact and that it is, indeed, most likely related to a dust density enhancement in the disk.

On the other hand, the ALMA observations do not display such a strong asymmetry at first glance. To further investigate if the best fit model to the SPHERE observations is compatible with the ALMA data, we generated models including eccentricity and azimuthal variations. We took the best fit model to the ALMA data (Table 3), and the following additional parameters from Table 4: $e = 0.093$, $\Phi_e = 127^\circ$, $w = 52^\circ$, $\Phi_\eta = 138^\circ$, and two values for $\eta$ (0.47 and 0.6). Figure 9 shows the comparison between these models and the observations. For $\eta = 0.47$ (derived from the SPHERE observations), the brightness of the western side is slightly underestimated but the residuals barely reach the $2\sigma$ level. For $\eta = 0.6$, the model becomes more similar to the observations and the residuals do not reach the $2\sigma$ level. Consequently, we cannot exclude that the disk shows departure

from centro-symmetry in the ALMA data, but the sensitivity of the observations does not allow us to put meaningful constraints on the azimuthal variations for the population of large grains. We can exclude damping factors of $\eta \sim 0.5$ to the $2\sigma$ level, but a value of 0.6 cannot be ruled out based on the available observations.

### 6.4.1. Possible planet-disk interactions?

Massive planets are often invoked to explain the eccentricity of a debris disk since the planet may shepherd it. The case of Fomalhaut is a good example of these types of studies. The disk shows an eccentricity comparable to the one of HD 61005 ($e \sim 0.11 \pm 0.1$, Kalas et al. 2005) and a candidate companion was detected by Kalas et al. (2008). Studies such as the one of Beust et al. (2014) investigate in detail the possible interactions between the planet and the planetesimal in the disk. In the last two decades, many numerical works have been led to investigate how planets can shape debris disks (e.g., Wyatt et al. 1999; Wyatt 2006; Nesvold et al. 2013; Pearce & Wyatt 2014; Faramaz et al. 2014), or how the planetesimals evolve within their mutual gravitational interactions (e.g., Kral et al. 2015; Jackson et al. 2014).

Via secular interactions, an eccentric planet can shape the debris disk. Because of differential precession, the planet-disk interactions can result in a (short-lived) spiral arm (Mouillet et al. 1997; Wyatt 2005). Then, when taking into account the effect of collisions in the disk, the spiral quickly dissipates
(Nesvold et al. 2013), and the disk becomes eccentric (and apse-
aligned with the planet's orbit). Because of the eccentricity, the
disk is expected to be brighter at the pericenter (the so-called
pericenter glow effect). Nonetheless, this may not be what is
happening in the disk around HD 61005. First, we did not de-
tect any point sources in neither the H- nor the K-band ADI
datasets, even with their large parallactic rotations. The exact de-
tection limits on the masses of possible planets will be published
in the survey paper of the SPHERE consortium, but prelimi-
nary tests using the PynPoint package (Amara & Quanz 2012)
shows that we should have been able to detect planets with a con-
trast of about ~10 mag in the near-IR. At an age of about 40 Myr,
we can therefore rule out any planets with a mass larger than
~7 Jupiter mass at a distance of 0.2". Second, as already dis-
cussed in Buenzli et al. (2010, and confirmed in this study), the
eccentricity of the disk is not sufficient for the pericenter glow to
solely explain the observed brightness asymmetry. Our modeling
results suggest that there is a ~50% reduction in density between
the pericenter and the apocenter, with a peak density closer to
the observer (Φ<sub>e</sub> ~ Φ<sub>η</sub> ~ 130°). Given that both Φ<sub>e</sub> and Φ<sub>η</sub> are
free parameters and that their most probable values are so close
to each other, it strongly points towards an increase of the dust
density distribution at the pericenter of the eccentric disk. In-
terestingly, the simulations of Pearce & Wyatt (2014) show the
opposite effect: the density is lower at the pericenter when in-
cluding an eccentric planet in the disk (the particles spend more
time at the apocenter than at the pericenter).

Wyatt (2006) describes the interactions between a massive
planet with the dust belt, in different configurations. Whether
the planetesimals are in 2:1 or 3:2 resonances with the planet,
we can expect the disk to display asymmetries not only in scat-
tered light but also at millimeter wavelengths. The large grains
traced by mm observations are expected to closely follow the
planetesimals that are in resonance with the planet. As assessed
previously, with the current sensitivity of the mm observations,
we cannot confirm nor rule out that large dust grains display de-
parture from centro-symmetry.

Therefore, if a planet were to be detected (at other wave-
lengths or with other facilities), it would be interesting to con-
front our observations and modeling results with the orbital pa-
rameters of the planet. But in the following, we hypothesize
that no massive planet was formed around HD 61005 and we
try to explain the inferred morphology of the disk under this
assumption.

6.4.2. The impact of gas in a debris disk

Lyra & Kuchner (2013) show that narrow, eccentric debris disks can be the consequence of the presence of gas. While debris disks are generally thought of as depleted of gas, several recent studies have demonstrated it is not necessarily the case, even at ages as old as a few tens of Myr (e.g., Hughes et al. 2008; Dent et al. 2014; Kospal et al. 2013; Moór et al. 2015a). Lyra & Kuchner (2013) suggest that the dust grains can heat up the surrounding gas via photoelectric heating. This will locally increase the gas pressure and, since the grains follow the pressure gradient, the pressure maximum will concentrate them even further. In the end, the disk may appear as narrow and eccentric, even though no planets are shaping it (although the original disk must be narrow). In Sect. 6.2, we concluded that the gas-to-dust ratio remains highly uncertain because of several unconstrained parameters (grain size distribution, CO/H₂ abundance ratio, LTE assumption). Overall, this scenario may be a viable solution to

explain our findings, but it cannot be verified at the moment (the
azimuthal density variation would have to be addressed as well).

6.4.3. The aftermath of a collision between Pluto-sized bodies?

An alternative possibility would be that we are witnessing an
impact between two massive bodies in the debris disk. This so-
lution is mostly interesting because we find a peak of density at
the pericenter of a slightly eccentric disk. In Sect. 6.3, we con-
cluded that self-stirring by massive planetesimals may be tak-
ing place in the 40 Myr old debris disk. At this age, and a ra-
dius of 60 au, a few massive oligarchs may have been formed,
which are disrupting the orbits of other planetesimals. The colli-
sional cross-sections of these bodies would have increased along
with the probabilities of intercepting each others' orbits. Numer-
ical studies of these types of impact (e.g., Jackson et al. 2014;
and Kral et al. 2015) suggest that our results may be compati-
ble with the scenario of an impact between two massive bod-
ies. This type of collision would release significant amount of
debris, and all of it would have to pass through the same loca-
tion where the collision initially occurred. This would enhance
the density and therefore the chances of subsequent collisions
at the location of this pinch point. In which case, the scattering
cross-section will locally increase along the eastern side com-
pared to the western side. Similar scenarios have been postulated
for the highly asymmetric disk around HD 15115 (Mazoyer et al.
2014), β Pictoris (Dent et al. 2014), and HD 181327 (Stark et al.
2014). According to the work of Jackson et al. (2014), the debris
disk around HD 61005 could maintain this asymmetry for about
~0.45 Myr at 60 au (1000 orbital periods). However, in these nu-
merical simulations, the eccentricity and overall complexity of
the resulting disk highly depends on the initial conditions of the
impact. Furthermore, Leinhardt & Stewart (2012) find that the
grain size distribution of the fragments released by the collision
of gravity-dominated bodies is best described by a power-law
slope in -3.85. Our modeling results of the SED are in excellent
agreement with this value, supporting the scenario that the disk
mainly consists of fragments from a catastrophic collision.

Nonetheless, a challenging aspect to this postulate is the
dust mass inferred from the SED modeling. Simulations of the
stochastic phase of planetesimal growth by Stewart & Leinhardt
(2012) suggest that collision of massive bodies should release
about ~15% of the total mass in debris. Assuming “debris” cor-
responds to grains with sizes between 2 µm and 5 mm (used to
model the SED), this would correspond to a total mass of about
2.5 M<sub>⊕</sub> for the catastrophic impact, which would make this type
of event an outliner in most theoretical works about the forma-
tion of planetesimals. For a 40 Myr old solar-mass star, the mod-
els of Kenyon & Bromley (2008, 2010) predict a dust mass pro-
duction of the order of ~0.05 M<sub>⊕</sub> (for 1 µm ≤ s ≤ 1 mm) for a
disk that is initially 3 times more massive than the MMSN (to
be compared to 0.25 M<sub>⊕</sub> interpolating our results to the same
grain size distribution). If the primordial disk around HD 61005
was indeed several times more massive than the MMSN (see
Sect. 6.3), it would evolve more rapidly and massive bodies
would have had the time to form. Overall, if this interpretation
of the multi-wavelength observations of the Moth holds, it would
suggest that massive bodies can still be formed at a large distance
from the star at an age of ~40 Myr. This would provide valuable
constraints on the formation and evolution of planetesimals in
the first 100 Myr.

The presence of spatially separated belts (inner and outer
belts) is often discussed in the context of planet-disk interactions
(e.g., Su et al. 2013). Massive planets may scatter away any close-by debris, hence opening a dust-free region in the disk. Planets may also be able to scatter debris from the outer regions and send them inwards, which would result in a bi-modal distribution of the dust (e.g., Bonsor & Wyatt 2012). In the case of HD 61005, we do not have strong constraints on the location of the inner belt. But, given that the evolution of a debris disk is an inside-out process (Kenyon & Bromley 2008), it is unlikely the inner belt is primordial. Given the age of the star, the inner regions should already have been depleted of any detectable amount of small dust grains. Still, assuming that there are no massive planets around HD 61005, is it possible to reconcile the presence of an analog of the asteroid belt with the above interpretation? If it is indeed the case that a few oligarchs have formed in the main belt at 60 au, the disk would be in an active and chaotic state. Some of these massive bodies may have collided in the main belt, but some may have scattered debris inwards. However, simulations, such as the ones presented in Jackson et al. (2014) and Kral et al. (2015), do not display cleared gaps. Consequently, under our assumption that there are no giant planets, it seems that our interpretation of the dynamical evolution of the outer regions is not directly related with what is taking place in the inner regions. It may well be that telluric planets have formed in the inner regions and further investigation remains necessary to better constrain the location of the innermost belt.

## 7. Conclusion

In this study, we presented new VLT/SPHERE and ALMA observations of the debris disk around the 40 Myr-old solar-type star HD 61005. We modeled both the radial and azimuthal distribution of the dust grains. We find that the disk is slightly eccentric (e ~ 0.1), with a dust density two times smaller at the apocenter compared to the pericenter. The sensitivity of the ALMA observations cannot confirm nor rule out an asymmetric disk for large mm-sized dust grains. We confirm that the eccentricity of the disk is not sufficient to explain the known brightness asymmetry in the near-IR. Thanks to the multi-wavelength analysis, we highlight a possible radial segregation of the grain size distribution since the disk appears more extended at near-IR wavelengths compared to far-IR and mm wavelengths. The parent planetesimals are confined within a narrow belt (with a semi-major axis of about 61 au), while small dust grains (or dust aggregates) would be pushed away by radiation pressure. No CO gas was detected with the ALMA observations and we did not detect any giant planets around the star. For this reason, we attempt to explain the morphology of the disk without invoking planet-disk interactions. We propose that we could be witnessing the aftermath of collision(s) between massive bodies in the main belt. The debris disk may be bright because self-stirring recently started to become significant in the belt, and it triggered a massive collision. The location of the collision became a pinch point where all debris must go through at each orbit, hence locally increasing the dust density and the scattering cross-section, as well as the chances of subsequent collisions. Nonetheless, this scenario would require a massive primordial disk to reconcile numerical simulations of debris disk evolution with the inferred dust mass. Finally, the presence of an inner belt would be consistent with the modeling of the SED, and we conclude its presence and dynamical evolution is unrelated to the evolution of the outer belt.

*Acknowledgements.* We would like to thank the anonymous referee for the valuable feedback we received. The comments improved the discussion of the paper,

mostly with respect to the scenario of a massive collision. They also helped with the presentation of the modeling strategy for all the different datasets presented in this paper. This paper makes use of the following ALMA data: ADS/JAO.ALMA 2012.1.00437.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. J. Olofsson, H. Avenhaus, C. Caceres, H. Canovas, M. R. Schreiber, A. Bayo, and S. Casassus acknowledge support from the Millennium Nucleus RC130007 (Chilean Ministry of Economy). J. Olofsson acknowledges support from ALMA/CONICYT Project 31130027, H. Avenhaus the support from FONDECYT 2015 Postdoctoral Grant 3150643, C. Caceres the support from CONICYT FONDECYT grant 3140592, H. Canovas the support from ALMA/CONICYT (grants 31100025 and 31130027), and A. Bayo the financial support from the Proyecto Fondacyt Iniciación 11140572. This work was supported by the Momentum grant of the MTA CSFK Lendület Disk Research Group. A. Moór acknowledges support from the Bolyai Research Fellowship of the Hungarian Academy of Sciences. J. Olofsson would like to thank Daniel Foreman-Mackey for his help with the parallelization of the code with emcee and Mark Booth for valuable discussions about ALMA data. This research made use of Astropy, a community-developed core Python package for Astronomy (Astronomy Collaboration 2013). This research has made use of the SIMBAD database (operated at CDS, Strasbourg, France). The authors have used the TOPCAT software (Taylor 2005) in this work.

## References

Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 216

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

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

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

Aumann, H. H., Beichman, C. A., Gillett, F. C., et al. 1984, ApJ, 278, L23

Avenhaus, H., Quanz, S. P., Schmid, H. M., et al. 2014, ApJ, 781, 87

Ballering, N. P., Rieke, G. H., Su, K. Y. L., & Montiel, E. 2013, ApJ, 775, 55

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

Bayo, A., Rodrigo, C., Barrado Y Navascués, D., et al. 2008, A&A, 492, 277

Bendo, G. J., Griffin, M. J., Bock, J. J., et al. 2013, MNRAS, 433, 3062

Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6

Beust, H., Augereau, J.-C., Bonsor, A., et al. 2014, A&A, 561, A43

Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, in SPIE Conf. Ser., 7014, 18

Boccaletti, A., Thalmann, C., Lagrange, A.-M., et al. 2015, Nature, 526, 230

Bonsor, A., & Wyatt, M. C. 2012, MNRAS, 420, 2990

Buenzli, E., Thalmann, C., Vigan, A., et al. 2010, A&A, 524, L1

Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icarus, 40, 1

Canovas, H., Rodenhuis, M., Jeffers, S. V., Min, M., & Keller, C. U. 2011, A&A, 531, A102

Canovas, H., Ménard, F., de Boer, J., et al. 2015, A&A, 582, L7

Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138

Castelli, F., Gratton, R. G., & Kurucz, R. L. 1997, A&A, 318, 841

Chen, C. H., Mittal, T., Kuchner, M., et al. 2014, ApJS, 211, 25

Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, in SPIE Conf. Ser., 7014, 3

De Silva, G. M., D'Orazi, V., Melo, C., et al. 2013, MNRAS, 431, 1005

Dent, W. R. F., Wyatt, M. C., Roberge, A., et al. 2014, Science, 343, 1490

Desidera, S., Covino, E., Messina, S., et al. 2011, A&A, 529, A54

Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in SPIE Conf. Ser., 7014, 3

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

Dominik, C., & Decin, G. 2003, ApJ, 598, 626

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

Dullemond, C. P., & Dominik, C. 2008, A&A, 487, 205

Eiroa, C., Marshall, J. P., Mora, A., et al. 2013, A&A, 555, A11

Elliott, P., Bayo, A., Melo, C. H. F., et al. 2014, A&A, 568, A26

Faramaz, V., Beust, H., Thébault, P., et al. 2014, A&A, 563, A72

Foreman-Mackey, D., Price-Whelan, A., Ryan, G., et al. 2014, Python Package, triangle.py v0.1.1

Gelman, A., & Rubin, D. B. 1992, Statist. Sci., 7, 457

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3

Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70

Hernández, J., Hartmann, L., Megeath, T., et al. 2007, ApJ, 662, 1067

Hines, D. C., Schneider, G., Hollenbach, D., et al. 2007, ApJ, 671, L165

Hughes, A. M., Wilner, D. J., Kamp, I., & Hogerheijde, M. R. 2008, ApJ, 681, 626

Jackson, A. P., Wyatt, M. C., Bonsor, A., & Veras, D. 2014, MNRAS, 440, 3757
Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, *Science Advances*, **1**, 1500109

Kalas, P., Graham, J. R., & Clampin, M. 2005, *Nature*, **435**, 1067

Kalas, P., Graham, J. R., Chiang, E., et al. 2008, *Science*, **322**, 1345

Kalas, P. G., Rajan, A., Wang, J. J., et al. 2015, ApJ, **814**, 32

Kenyon, S. J., & Bromley, B. C. 2006, AJ, **131**, 1837

Kenyon, S. J., & Bromley, B. C. 2008, ApJS, **179**, 451

Kenyon, S. J., & Bromley, B. C. 2010, ApJS, **188**, 242

Kospal, A., Moor, A., Juhasz, A., et al. 2013, ApJ, **776**, 77

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

Krivov, A. V. 2010, RA&A, **10**, 383

Lagrange, A.-M., Langlois, M., Gratton, R., et al. 2016, A&A, **586**, L8

Langlois, M., Dohlen, K., Vigan, A., et al. 2014, in SPIE Conf. Ser., **9147**, **91471R**

Lebouteiller, V., Barry, D. J., Spoon, H. W. W., et al. 2011, ApJS, **196**, 8

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

Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, **745**, 79

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

Lyra, W., & Kuchner, M. 2013, *Nature*, **499**, 184

Maire, A.-L., Bonnefoy, M., Ginski, C., et al. 2016, A&A, **587**, A56

Maness, H. L., Kalas, P., Peek, K. M. G., et al. 2009, ApJ, **707**, 1098

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

Matrà, L., Panić, O., Wyatt, M. C., & Dent, W. R. F. 2015, MNRAS, **447**, 3936

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

Mazoyer, J., Boccaletti, A., Augereau, J.-C., et al. 2014, A&A, **569**, A29

Millar-Blanchaer, M. A., Graham, J. R., Pueyo, L., et al. 2015, ApJ, **811**, 18

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

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

Milli, J., Vigan, A., Mouillet, D., et al. 2016, A&A, submitted

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

Moór, A., Henning, T., Juhász, A., et al. 2015a, ApJ, **814**, 42

Moór, A., Kóspál, Á., Ábrahám, P., et al. 2015b, MNRAS, **447**, 577

Morales, F. Y., Rieke, G. H., Werner, M. W., et al. 2011, ApJ, **730**, L29

Mouillet, D., Larwood, J. D., Papaloizou, J. C. B., & Lagrange, A. M. 1997, MNRAS, **292**, 896

Mustill, A. J., & Wyatt, M. C. 2009, MNRAS, **399**, 1403

Nesvold, E. R., Kuchner, M. J., Rein, H., & Pan, M. 2013, ApJ, **777**, 144

Olofsson, J., Juhász, A., Henning, T., et al. 2012, A&A, **542**, A90

Olofsson, J., Henning, T., Nielbock, M., et al. 2013, A&A, **551**, A134

Ott, S. 2010, in Astronomical Data Analysis Software and Systems XIX, eds. Y. Mizumoto, K.-I. Morita, & M. Ohishi, ASP Conf. Ser., **434**, 139

Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, **412**, 13

Pavlov, A., Möller-Nilsson, O., Feldt, M., et al. 2008, in SPIE Conf. Ser., **7019**, 39

Pawellek, N., Krivov, A. V., Marshall, J. P., et al. 2014, ApJ, **792**, 65

Pearce, T. D., & Wyatt, M. C. 2014, MNRAS, **443**, 2541

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Mach. Learn. Res., **12**, 2825

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, **518**, L1

Pinte, C., Padgett, D. L., Ménard, F., et al. 2008, A&A, **489**, 633

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, **518**, L2

Ricarte, A., Moldvai, N., Hughes, A. M., et al. 2013, ApJ, **774**, 80

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

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

Siess, L., Dufour, E., & Forestini, M. 2000, A&A, **358**, 593

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

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

Steele, A., Hughes, A. M., Carpenter, J., et al. 2016, ApJ, **816**, 27

Stewart, S. T., & Leinhardt, Z. M. 2012, ApJ, **751**, 32

Su, K. Y. L., Rieke, G. H., Malhotra, R., et al. 2013, ApJ, **763**, 118

Taylor, M. B. 2005, in Astronomical Data Analysis Software and Systems XIV, eds. P. Shopbell, M. Britton, & R. Ebert, ASP Conf. Ser., **347**, 29

Thalmann, C., Mulders, G. D., Hodapp, K., et al. 2014, A&A, **566**, A51

Thebault, P. 2009, A&A, **505**, 1269

van Leeuwen, F. 2007, A&A, **474**, 653

Vigan, A., Moutou, C., Langlois, M., et al. 2010, MNRAS, **407**, 71

Wolf, S., & Hillenbrand, L. A. 2005, Comput. Phys. Commun., **171**, 208

Wyatt, M. C. 2005, A&A, **440**, 937

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

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

Wyatt, M. C., Dermott, S. F., Telesco, C. M., et al. 1999, ApJ, **527**, 918

e-mail: johan.olofsson@uv.cl

<footnote n="1">¹ Instituto de Física y Astronomía, Facultad de Ciencias, Universidad de Valparaíso, Av. Gran Bretaña 1111, Playa Ancha, Valparaíso, Chile</footnote>

<footnote n="2">² ICM nucleus on protoplanetary disks, Protoplanetary discs in ALMA Early Science, Universidad de Valparaíso, Valparaíso, Chile</footnote>

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

<footnote n="4">⁴ Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile</footnote>

<footnote n="5">⁵ Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences, PO Box 67, 1525 Budapest, Hungary</footnote>

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

<footnote n="7">⁷ Institute for Astronomy, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland</footnote>

<footnote n="8">⁸ Univ. Grenoble Alpes, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG, UMR 5274), 38000 Grenoble, France</footnote>

<footnote n="9">⁹ CNRS, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG, UMR 5274), 38000 Grenoble, France</footnote>

<footnote n="10">¹⁰ LESIA, Observatoire de Paris, CNRS, Université Pierre et Marie Curie 6 and Université Denis Diderot Paris 7, 5 place Jules Janssen, 92195 Meudon, France</footnote>

<footnote n="11">¹¹ Astronomical Institute Anton Pannekoek, University of Amsterdam, PO Box 94249, 1090 GE Amsterdam, The Netherlands</footnote>

<footnote n="12">¹² INAF Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote n="13">¹³ Department of Astronomy, Stockholm University, AlbaNova University Center, 10691 Stockholm, Sweden</footnote>

<footnote n="14">¹⁴ CNRS, Centre de Recherche Astrophysique de Lyon, Observatoire de Lyon, École Normale Supérieure de Lyon, Université Lyon 1, 9 avenue Charles André, Saint-Genis Laval, 69230, France</footnote>

<footnote n="15">¹⁵ Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388 Marseille, France</footnote>

<footnote n="16">¹⁶ UMI-FCA, CNRS/INSU France, and Departamento de Astronomía, Universidad de Chile, Casilla 36-D Santiago, Chile</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Schematic view of how a disk, seen edge-on, is described in the DDiT code. Different shells ($n_r = 10$) are represented in different fading colors, while dust belts are represented without transparency (see text for details).</figcaption></figure>

# Appendix A: Description of the code

## A.1. Computing the spectral energy distribution

The DDiT code is mostly based on the code already described in Olofsson et al. (2012, 2013, which was based on the “Debris disk radiative transfer simulation tool” described in Wolf & Hillenbrand 2005). However, as opposed to the work presented in Olofsson et al. (2012, 2013), only one dust species per model is considered in this study. The dust properties are defined by the optical constants, the density ($\rho$), the minimum and maximum grain sizes ($s_{min}$ and $s_{max}$, respectively), and the slope $p$ of the grain size distribution ($dn(s) \propto s^p ds$, $p < 0$). The absorption efficiencies, $Q_{abs}(\lambda, s)$, are computed using the Mie theory for $n_g = 100$ grain sizes logarithmically spaced between $s_{min}$ and $s_{max}$.

Concerning the disk model itself, this is defined by the following parameters: a reference radius $r_0$, two power-law indexes $\alpha_{in} > 0$ and $\alpha_{out} < 0$, the opening angle of the disk ($\psi = \arctan(h/r)$, where $h$ is the height from the midplane), and the total dust mass $M_{disk}$. We first define the inner and outer radii of the disk by finding $r_{in} = r_0 \times C^{(1/\alpha_{in})}$ and $r_{out} = r_0 \times C^{(1/\alpha_{out})}$, where $C$ is a threshold value for the minimum density that we choose to be $C = 10^{-2}$. The code is 1D and the disk is divided into $n_r$ spherical shells centered at $r_j$, logarithmically spaced between $r_{in}$ and $r_{out}$ (see Fig. A.1).

The dust number density distribution for the shell at distance $r_j$ and the grain size $s_k$, $N_{dens}(r_j, s_k)$, is defined as

$$N_{\text{dens}}(r_j, s_k) = \left[ \left( \frac{r_j}{r_0} \right)^{-2\alpha_{\text{in}}} + \left( \frac{r_j}{r_0} \right)^{-2\alpha_{\text{out}}} \right]^{-1/2} \times N_s(s_k), \quad (A.1)$$

where $N_s(s_k)$ is a function depending solely on the grain size (similarly to Eq. (2) of Dullemond & Dominik 2008)

$$N_s(s_k) = \left( \frac{s_k}{s_{\min}} \right)^p \times s_k \times \Delta\log(s)_k, \quad (A.2)$$

where $\Delta\log(s)_k$ is the width of each bin in logarithmic space of $s$. Since the $n_g$ grain sizes are logarithmically spaced, each $\Delta\log(s)_k$ is the same, except the first and last ones which are half of that value. The mass density distribution is then given by $M_{\text{dens}}(r_j, s_k) = N_{\text{dens}}(r_j, s_k) \times 4/3\pi s_k^3 \rho$. Both $N_{\text{dens}}(r_j, s_k)$ and $M_{\text{dens}}(r_j, s_k)$ are then normalized so that the integral of $M_{\text{dens}}(r_j, s_k) \times V_j$ (where $V_j$ is the volume of the shell j) over $s$ and $r$ equals the total input dust mass $M_{\text{disk}}$.

The radially integrated optical depth at the wavelength of the peak of the stellar emission ($\tau_{\lambda,\text{peak}}$) is computed at the outer edge $r_{\text{out}}$, assuming a spherical symmetry of the disk. This quantity is then divided by $f_{\psi} = \sin(\psi)$ to obtain the final optical depth of the model, for a given $\psi$. This $1/f_{\psi}$ factor accounts for the fact that the code is 1D, i.e., that the disk is not a sphere but has a maximum vertical height $h$ from the midplane (hence can be seen as dust belts). In Fig. A.1, the shells are shown in fading colors, while the dust belts are represented without transparency. The $f_{\psi}$ factor is the ratio of the volume of a sphere of radius $R$ to which we subtract two spherical cones (up and down) of height $R \times (1 - \sin(\psi))$, divided by the volume of a sphere of the same radius. One spherical cone has a volume of $2\pi/3R^3(1 - \sin(\psi))$, which simplifies the volume ratio to $\sin(\psi)$. By doing so, we can keep $M_{\text{disk}}$ constant and leave the output of the code unchanged while simply updating the optical depth. Once the model has been computed, a check is performed to see if $\tau_{\lambda,\text{peak}}$ is smaller than unity for the given dust mass and $\psi$.

When computing the SED, we consider that the contribution of scattered light is negligible at mid to far-IR wavelengths compared to the thermal emission. Furthermore, since debris disks are optically thin at all wavelengths, geometrical parameters such as the inclination have no impact on the final SED. Finally, as the observed SED is constructed from spatially unresolved observations, disk or aperture sizes do not have to be taken into consideration.

For each shell, the temperature of a dust grain is computed as a function of its size $s$ and distance $r$ to the star. Assuming radiative equilibrium, one can express $r$ as

$$r(s_k, T_{\text{dust}}) = \frac{R_{\star}}{2} \times \sqrt{\frac{\int_{\lambda} F_{\star}(\lambda) Q_{\text{abs}}(\lambda, s_k) d\lambda}{\int_{\lambda} \pi B_{\lambda}(\lambda, T_{\text{dust}}) Q_{\text{abs}}(\lambda, s_k) d\lambda}}, \quad (A.3)$$

where $R_{\star}$ is the stellar radius, $F_{\star}$ the stellar surface flux, and $B_{\lambda}$ the Planck function at the dust temperature $T_{\text{dust}}$. By inverting numerically the above equation, we obtain the radial dependency of the temperature. Afterwards, the flux $F_{\text{th}}$ (the thermal radiation of a single dust grain) received by an observer at distance $d_{\star}$ is computed as

$$F_{\text{th}}(\lambda, r_j, s_k) = \left( \frac{4\pi s_k^2}{4\pi d_{\star}^2} \right) \pi B_{\lambda}(\lambda, T_{\text{dust},j}) Q_{\text{abs}}(\lambda, s_k). \quad (A.4)$$

For a given shell, the total thermal emission is obtained by multiplying $F_{\text{th}}(\lambda, r_j, s_k)$ by the normalised $N_{\text{dens}}(r_j, s_k)$. The final SED is the sum of the thermal emission from all $n_r$ shells.

## A.2. Computing synthetic images

At optical or near-IR wavelengths, stellar light scattered by dust grains can no longer be considered negligible and has to be accounted for. We follow a parallel approach as for the thermal emission and the stellar light scattered by a single dust grain is computed as follows

$$F_{\text{scat}}(\lambda, r_j, s_k, \theta) = \frac{F_{\star}(\lambda)}{4\pi r_j^2} \pi s_k^2 \left( \frac{1}{4\pi d_{\star}^2} \right) Q_{\text{sca}}(\lambda, s_k) S_{11}(\lambda, s_k, \theta), \quad (A.5)$$
<figure><img src="image_12.png" /><figcaption>Fig. A.2. Sketch of a disk (inclined by i) as seen from the side to illustrate the entry and exit points ($D_{En}$ and $D_{Ex}$) when computing images. The cut shows the y and z in the disk frame axis and is made at $x = 0$ (the x axis pointing toward the reader).</figcaption></figure>

where $Q_{sca}$ is the scattering efficiencies calculated with the Mie theory and $S_{11}(\theta)$ is the scattering phase function.

In the model parametrization, we define two quantities for the image itself, the number of pixels $n_p$ (images are squares) and the size $l_p$ of one pixel in units of au. Then, for each of the $n_p^2$ pixels in the image, we first apply adequate projection and rotation to account for the geometrical parameters i and $\phi$ (inclination and position angle, respectively). For $\phi = 0^{\circ}$, the semi-major axis is along the south-north direction, and the front side of the disk is chosen to be on the east side. For increasing $\phi$ the disk rotates counter-clockwise, towards the east.

Since the disk has a non-zero opening angle $\psi$, it is not infinitely flat in the vertical direction. As in other codes (e.g. GRATER, Augereau et al. 1999), we assume a vertical Gaussian profile for the dust distribution, with a $\sigma$ width equal to $\psi$. This avoids having blocky images, especially for highly inclined disks. Figure A.2 shows a crude sketch of a cut in the (y, z) reference plane of the disk, for $x = 0$. The disk has an inclination i and for each pixel, we trace rays along the line of sight and search for intersections with the disk. The disk can be described by four main geometrical shapes: inner and outer shells, and upper and lower planes. Each of them can be either “entry” or “exit” points for the rays (blue and yellow lines in Fig. A.2, respectively). One should note that only slightly different intersection conditions are to be considered for cases where $\pi/2 - i < \psi$ (e.g., edge-on disk) and that the code can produce images for such disks. We can therefore obtain the coordinates of both the entry and exit points ($D_{En}$ and $D_{Ex}$, respectively). Since a pixel has a size $l_p$, the column between $D_{En}$ and $D_{Ex}$ (modulus of $|l|$) has a volume of $l_p^2 \times |l|$. Since this column can probe a significant range of distances $r$ to the star, we divide it into $n_{los}$ boxes, along the line of sight. For each of these boxes, we evaluate the distance to the star $r$, the volume of the box $V_{box}$, and the scattering angle (dot product between the line of sight and the position in the disk with respect to the star). We then interpolate $F_{th}$, $F_{scat}$, and $N_{dens}$ at the distance $r$ of the center of the box, and interpolate the phase function at the scattering angle. We then sum the product $F_{th/scat} \times N_{dens} \times V_{box}$ for each box, to obtain both the scattered light and thermal emission for this given pixel.

### A.2.1. Polarized images

The code can also produce polarized light images on top of the total intensity images. Neglecting multiple scattering events, a reasonable assumption in optically thin disks, we compute $F_{pol}(\lambda, r_j, s_k, \theta)$ similarly to Eq. (A.5) simply replacing $S_{11}$

<figure><img src="image_13.png" /><figcaption>Fig. A.3. Sketch of an eccentric disk. The color-coding represents the density distribution. The rotation angle is indicated by $\Phi_e$, the reference angle for the azimuthal variation by $\Phi_\eta$, and $\eta = 0.3$ in this example.</figcaption></figure>

by $S_{12}$. The second element of the Müller matrix $S_{12}$ is calculated using the Mie theory (while $S_{11}$ can be approximated by the Henyey-Greenstein function). We can therefore produce synthetic images along the Stokes vector $Q$, which can be used to model polarimetric observations.

### A.2.2. Eccentric disks

The code can produce images of eccentric disks. The inner and outer shells as described in Fig. A.2 are modified to inner and outer edges with an eccentricity $e$ parameter and a rotation angle $\Phi_e$ (see Fig. A.3 for a sketch of a face-on disk). For each pixel in the image, at position $(x, y)$ (in the reference frame rotated by $\Phi_e$), we check if it is in between the inner and outer edges, which are defined as,

$$ \text{Edge}_{\text{inner, outer}} = \frac{r_{\text{in, out}}}{1 + e \times \cos[a + \Phi_e]}, \qquad (A.6) $$

where $a = \text{atan2}(y, x)$. For $\Phi_e = \phi = 0^{\circ}$, the pericenter is located on the south side.

### A.2.3. Azimuthal density variations

It is also possible to introduce azimuthal dust density variations in the synthetic images. The density distribution is parametrized with three parameters: a damping factor $\eta$ ($0 \le \eta \le 1$), a reference angle $\Phi_\eta$, and the azimuthal variation follows a Gaussian profile of $\sigma = w$ that peaks at 1 for $\text{atan2}(y, x) = \Phi_\eta$ and has an amplitude of $1 - \eta$. The attenuation factor is of the form

$$ N_{\text{att}}(x, y) = \eta + (1 - \eta) \times \exp\left(-\frac{(\text{atan}(y, x) - \Phi_{\eta})^2}{2w^2}\right). \qquad (A.7) $$

For $\Phi_\eta = \phi = 0^{\circ}$, the density maximum is located on the south side. For a given azimuthal angle, $N_{dens}$ is multiplied by $N_{att}$ before computing the image.
Appendix B: Results

<figure><img src="image_14.png" /><figcaption>Fig. B.1. Projected posterior probability distributions for the combined modeling of the ADI and DPI H-band observations.</figcaption></figure>
<figure><img src="image_15.png" /></figure>

**Fig. B.2.** Top view of the dust density distribution (normalized to its maximum) for the best-fit model to the SPHERE observations. The direction to the observer is indicated by the arrow and the spatial scale in the lower right corner.

## Appendix C: Miscellaneous

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

**Fig. C.1.** *U*<sub>φ</sub> component of the *H*-band DPI observations.

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

**Fig. C.2.** Coverage of the (*u*, *v*) plane for the ALMA observations used in this paper.

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

**Fig. C.3.** Elliptical mask binned over different scattering angles to depict how the phase function is calculated. The second bin of the phase function would be the mean flux in the region highlighted in white. The color coding shows the scattering angle for a given position in the elliptical mask. The central circle shows the numerical mask to hide the central region.
<table><caption>**Table C.1.** Notations used in this study.</caption><thead><tr><th>Notations</th><th>Units</th><th>Explanation</th></tr></thead><tbody><tr><td>L<sub>★</sub></td><td>[L<sub>☉</sub>]</td><td>Stellar luminosity</td></tr><tr><td>M<sub>★</sub></td><td>[M<sub>☉</sub>]</td><td>Stellar mass</td></tr><tr><td>R<sub>★</sub></td><td>[R<sub>☉</sub>]</td><td>Stellar radius</td></tr><tr><td>T<sub>★</sub></td><td>[K]</td><td>Stellar effective temperature</td></tr><tr><td>d<sub>★</sub></td><td>[pc]</td><td>Stellar distance</td></tr><tr><td>i</td><td>[°]</td><td>Inclination of the disk</td></tr><tr><td>φ</td><td>[°]</td><td>Position angle of the disk (disk aligned North-South for φ = 0°, front side is towards the East)</td></tr><tr><td>r<sub>0</sub></td><td>[au]</td><td>Semi latus rectum, the semi-major axis of the disk being r<sub>0</sub>/(1 - e<sup>2</sup>)</td></tr><tr><td>e</td><td></td><td>Eccentricity of the disk</td></tr><tr><td>η</td><td></td><td>Dust density damping factor</td></tr><tr><td>Φ<sub>e</sub></td><td>[°]</td><td>Reference angle to define the pericenter (along the semi-major axis, at the southern side for Φ<sub>e</sub> = φ = 0°)</td></tr><tr><td>Φ<sub>η</sub></td><td>[°]</td><td>Reference angle for the peak density (along the semi-major axis, at the southern side for Φ<sub>η</sub> = φ = 0°)</td></tr><tr><td>w</td><td>[°]</td><td>Width of the Gaussian profile for the azimuthal dust density variation</td></tr><tr><td>g</td><td></td><td>Asymmetry parameter for the HG phase function</td></tr><tr><td>α<sub>in</sub></td><td></td><td>Slope for the inner dust density distribution</td></tr><tr><td>α<sub>out</sub></td><td></td><td>Slope for the outer dust density distribution</td></tr><tr><td>ψ</td><td>[rad]</td><td>Opening angle of the disk</td></tr><tr><td>β</td><td></td><td>Radiation pressure to gravitational forces ratio</td></tr><tr><td>s<sub>blow</sub></td><td>[μm]</td><td>Dust grain size for which β = 0.5</td></tr><tr><td>s<sub>min</sub></td><td>[μm]</td><td>Minimum grain size used in the model</td></tr><tr><td>s<sub>max</sub></td><td>[mm]</td><td>Maximum grain size</td></tr><tr><td>p</td><td>—</td><td>Grain size distribution exponent</td></tr><tr><td>f<sub>1300</sub></td><td>[mJy]</td><td>Flux at 1.3 mm</td></tr><tr><td>f<sub>ice</sub></td><td></td><td>Fraction of water ice mixed to the dust opacity</td></tr><tr><td>f<sub>porosity</sub></td><td></td><td>Fraction of vacuum mixed to the dust opacity</td></tr><tr><td>M<sub>dust</sub></td><td>[M<sub>⊕</sub>]</td><td>Dust mass</td></tr><tr><td>f<sub>ADI</sub></td><td>[arbitrary]</td><td>Scaling factor for the ADI images</td></tr></tbody></table>