© J. Olofsson et al. 2022

1. Introduction

The second generation dust that we observe in debris disks is continuously replenished from a collisional cascade of large planetesimals (see Krivov 2010; Hughes et al. 2018), or, in some rarer cases, by transient events such as the violent breakup of larger planetary embryos (Jackson et al. 2014; Kral et al. 2015). The physical properties of the particles released from such collisions still remain quite elusive, even though it is becoming clear that the assumption of compact spherical grains fails to reproduce contemporary observations (e.g., Milli et al. 2017, 2019). Debris disks are interesting targets to try and better characterize the optical dust properties for two reasons. First, they are optically thin at all wavelengths and we do not have to account for multiple scattering events or estimating the temperature in an optically thick medium, simplifying (to some extent) the modeling of the observations. Second, the dynamics of the particles strongly depend on their sizes. For most stellar spectral types, radiation pressure, which can be parametrized by the unitless ratio β between the radiation pressure and gravitational forces (β ∝ 1/s for grains larger than a few μm), naturally results in a different spatial extent for different grain sizes¹. The smaller particles are set on highly eccentric orbits, while larger ones remain on orbits very similar to the ones of the parent bodies. Since observations at different wavelengths probe different grain sizes, we can take advantage of this natural spatial segregation (Thebault et al. 2014).

Millimeter (mm) observations inform us about the spatial distribution of the large dust particles, and therefore best trace the location of the birth ring of planetesimals where the collisions are taking place. On the other hand, near-infrared (IR) scattered and polarized light observations are sensitive to the other extreme of the size distribution and trace μm-sized dust grains. Consequently, there is a great synergy between the two observational techniques; by constraining the location of the birth ring from mm observations, we can then have a better description of where the small dust grains should be launched from on highly eccentric orbits. This leads to a more accurate description of the spatial distribution of“observable” dust particles (with sizes s ≲ 1 mm), which can be used to better constrain their optical properties, such as the absorption efficiencies and near-IR phase function.

The debris disk around the A0 star HD 32297 is an ideal candidate to investigate the spatial distribution of second generation dust using multi wavelengths observations. It has been spatially resolved in the near-IR (e.g., Kalas 2005; Rodigas et al. 2014; Schneider et al. 2014; Bhowmik et al. 2019; Esposito et al. 2020). The disk is seen almost perfectly edge-on, and displays extended swept-back wings (best seen in the Hubble Space Telescope observations of Schneider et al. 2014). The disk has also been observed at mm wavelengths, and spatially resolved along the major-axis (MacGregor et al. 2018; Cataldi et al. 2020). Interestingly, MacGregor et al. (2018) reported the presence of a halo in the ALMA Band 6 (1.3 mm) dataset and concluded that this additional component cannot reasonably arise from μm-sized dust grains. Nonetheless, the detection of such a halo, possibly composed of mm-sized particles, is surprising as these grains should in principle not venture very far away from the birth ring, justifying further investigation using multi-wavelength observations. The almost edge-on configuration of the disk does come at a cost since information is lost due to projection effects. However, highly inclined disks allow for a wider range of scattering angles to be probed to constrain the phase function. Furthermore, due to the same projection effects, extended halos are easier to detect at all wavelengths as we probe larger column densities at all distances. HD 32297, having high angular resolution observations at both near-IR and sub-mm wavelengths and displaying a rather unique halo in ALMA observations, is therefore an interesting target to further study.

In the literature, observations (either near-IR or mm) are usually reproduced using geometric models that do not necessarily include any physics in them, and the interpretation of the physical processes at play is done a posteriori (Augereau et al. 1999; Marino et al. 2016; Engler et al. 2017; Kennedy et al. 2018; Daley et al. 2019; Matrà et al. 2019; Milli et al. 2019; Ren et al. 2019; Olofsson et al. 2020, among others). In some cases, more complex models are used to model the observations, taking for instance into account the effect of radiation pressure or stellar winds on the small particles (e.g., Esposito et al. 2016; Olofsson et al. 2019; Adam et al. 2021), or the intrinsic width of the parent planetesimal belt (Kennedy 2020). We here propose a novel approach at modeling near-IR and sub-mm observations, belonging to this second“family” of models. This approach can provide more stringent constraints on the properties of the dust particles. The philosophy is to identify the spatial scales, which are intrinsically related to typical grain sizes, that are most representative of the observations. By computing images for different β intervals, we can then identify the intervals (i.e., spatial extent) best suited to match the observations. This information can then be used to derive the properties of the dust grains (e.g., absorption efficiencies, phase function, porosity). In this paper, we first describe the multi-wavelength observations of the disk around HD 32297, how we model them, and discuss the results, before concluding.

2. Observations

2.1. SPHERE observations

To probe the population of small dust particles, we used the“Spectro-Polarimetric High Contrast Exoplanet REsearch” (SPHERE, Beuzit et al. 2019) dual-beam polarimetric imaging (DPI, Dohlen et al. 2008) J-band observations of HD 32297, which were first presented in Bhowmik et al. (2019). We reduced the data using the IRDAP² package presented in van Holstein et al. (2020, version 1.3.1). The outputs of the pipeline are the Q_ϕ and U_ϕ images, the former containing the polarized signal from the disk (left panel of Fig. 1), while the latter is free of astro-physical signal (assuming single scattering events) and is used as a proxy for the uncertainties.

2.2. ALMA observations

To constrain the location of the birth ring, we used the Band 8 (615 μm) ALMA observations published in Cataldi et al. (2020, program ID 2017.1.00201.S), with naturally weighted beam size of 0.63″ × 0.54″ (marginally smaller than the observations presented in MacGregor et al. 2018 with a beam of 0.76″ × 0.51″). The data was reduced using the“Common Astronomy Software Applications” package (CASA, version 5.6, McMullin et al. 2007). Once the reduction was performed using the script provided by the observatory, we re-evaluated the weights (∝1/σ²) using the statwt task, and used the split task to average the complex visibilities to a single channel per spectral window, and a time bin of 30 s, before exporting them using the uvplot package (Tazzari 2017). The bottom left panel of Fig. 1 shows the observations (using tclean and natural weighting) and the beam size in the lower left corner. The noise level was estimated to σ = 0.3 mJybeam⁻¹, from the clean image, in regions where there is no disk signal.

3. Geometrical fit to the observations

As mentioned before, the outward extent of the disk is governed by the radiation pressure, parametrized by the β ratio. By“launching” thousands of particles with different β values, similarly to Lee & Chiang (2016), we can therefore probe different spatial scales. Furthermore, by working with β values, we can make minimal assumptions on the exact sizes of the dust particles or their composition. Our goal is to identify which spatial regions contribute the most to best match the observations, taking into account that particles with different β values probe different spatial scales.

3.1. The dynamics of dust particles

For a given model³, we“launch” n_dust particles with different β values (see next paragraph), representing dust grains of various sizes, assuming that their parent bodies have a semi-major axis a, eccentricity e, and argument of periapsis ω. The semi-major axis a is drawn from a normal distribution⁴ centered at a₀ with a standard deviation δ_a. The launch point is uniformly distributed in mean anomaly between –π and π, and is then converted into a true anomaly v by solving the Kepler equation for e. The initial velocity of the particle is assumed to be the Keplerian velocity of the parent body. Following, Wyatt et al. (1999); Wyatt (2006), and Lee & Chiang (2016), for each particle (β given and v), the new orbital elements are computed as

$$\begin{array}{l}{a}_{\text{n}}=\frac{a\left(1-\beta \right)}{1-2\beta \frac{1+e\cos \left(v\right)}{1-e^2}},\hfill \\ e_{\text{n}}=\frac{1}{1-\beta }\times \sqrt{e^2+2\beta e\cos \left(v\right)+{\beta }^2,}\hfill \\ {\omega }_{\text{n}}=\omega +\arctan \left[\frac{\beta \sin \left(v\right)}{\beta \cos \left(v\right)+e}\right].\hfill \end{array}$$(1)

For each particle, we draw a longitude of ascending nodes uniformly between −π and π and an inclination following a normal distribution centered at 0 radians, with a standard deviation ψ. We then draw a value for the mean anomaly (∈] −π, π]) and compute the corresponding true anomaly (using the eccentricity e_n). The (x; y; z) positions of the particle, as well as its distance to the star $r=\sqrt{x^2+y^2+z^2}$, are then calculated using all the orbital parameters. To limit the number of free parameters, we set the eccentricity of the parent bodies to e = 0 (in that case, the value of ω no longer matters, removing another free parameter). The disk around HD32297 is not known for showing any signs of eccentricity and the edge-on configuration is the most challenging to properly constrain this parameter. Implicitly, we are making the assumption that all the parent bodies have small eccentricities (to ensure enough collisions are taking place), but that on average, the parent belt is circular.

We consider n_dust particles, with β values between β_min and β_max. To account for the fact that particles on highly eccentric orbits will spend most of their time in low density regions, and therefore survive longer, we also compute the following “correction factor” α = (1 − β)^3/2/[1 − e² − 2β × (1 + ecos(v))]^3/2 for each particle (Strubbe & Chiang 2006; Thébault & Wu 2008). Up until now, the free parameters are a₀, δ_a, and ψ.

Fig. 1 From left to right: observations, best fit model, and residuals, for the SPHERE and ALMA observations (top and bottom, respectively). The scaling is linear and the same for all horizontal panels. For the ALMA observations, the beam is shown in the leftmost panel, and the contours on the right panel are at [–1,1,2] σ (σ = 0.3 mJy beam−1).

3.2. Synthetic images

To create images to be compared with the observations, there are three additional parameters, the pixel size, as well as the inclination of the disk i and its position angle ϕ. Each of the (x; y; z) position is projected and rotated to account for i and ϕ, respectively. We then find the image pixel closest to the new values, accounting for the pixel size.

For scattered light images, we compute the scattering angle θ (the angle between the star, the particle, and the observer) using the dot product between a unit vector along the line of sight and the rotated 3D coordinates of the particle. The contribution of a particle to a pixel of the image will then be αS ₁₂(θ)=(rβ)², where S ₁₂ is the polarized phase function (see later), and θ the scattering angle. The 1/r² accounts for the illumination factor, and the 1/β² represents the cross section of the particle.

For mm images, from the distance r (in au), we estimate the size-independent temperature of the grains T_dust as $T_{\text{dust}}=278.3L_{★}^{1/4}/\sqrt{r}$ (Wyatt 2008), compute the Planck function B_v(λ, T_dust) at the wavelength of interest λ, and the contribution of a single particle will be αB_v(λ, T_dust)/β², therefore assuming that the particles behave as perfect blackbodies (the 1/β² accounts for the surface area of the particles).

3.3. Words of caution

The main benefit of the approach described above is that it relies as little as possible on the optical properties of the dust particles, therefore decreasing the number of free parameters, such as the dust composition, the porosity, or the light scattering theory used. While this speeds up the computational time, this also comes at a cost. For instance, both the phase function S₁₂ and the temperature T_dust should in fact depend on the size of the particles, and the contribution to the thermal emission (polarized scattered light) should also depend on the absorption (scattering) efficiencies of the particles.

3.4. Modeling strategy

To find the best fit solution we used the emcee affine-invariant ensemble sampler (Foreman-Mackey et al. 2013), with the following four free parameters: a₀, δ_a, and two standard deviations ψsphere, and ψ_ALMA, for the near-IR and mm observations, respectively. The distribution of β values will be evaluated independently for each realization of the emcee sampler (see later), which is parametrized to use 20“walkers” and a chain length of 10 000 for each of them. The need for two different scale heights for the SPHERE and ALMA data is motivated by the work presented in Olofsson et al. (2022) where they showed that the scale height can vary as a function of the wavelength (disks being flatter at near-IR wavelengths compared to mm wavelengths) if the disk harbors gas, which is the case for HD 32297 (Greaves et al. 2016; MacGregor et al. 2018; Cataldi et al. 2020). Vertical stratification is also expected when the disk is devoid of gas (Thebault 2009), but in this case the disk should be seen flatter at mm wavelengths compared to near-IR observations. The inclination and position angle of the disk are not free parameters in the modeling, since the ALMA observations do not have a comparable angular resolution compared to the SPHERE ones. Furthermore, Cataldi et al. (2020) noted a possible discrepancy when comparing the inclination derived from near-IR scattered light observations and the C I emission line. The authors found a smaller inclination (i ~ 77.9°) in the latter case, but the origin of this discrepancy is not clear. We therefore opted to fix both parameters to i = 86.9° and ϕ = – 132.3° (from the modeling of the same SPHERE observations, presented in Olofsson et al. 2022), thus reducing the amount of free parameters to four.

The other (fixed) parameters required for the fitting are the following; L_✶ = 7.61 L_⊙ (Olofsson et al. 2022), d_✶ = 129.71 pc (Gaia Collaboration 2016, 2021), n_dust = 10000000, n_β = 10, and λ = 1.22 and 615 μm (SPHERE and ALMA, respectively). β_min and β_max are set to 0.01 and 0.49, respectively, to capture a wide range of spatial scales. As mentioned in Strubbe & Chiang (2006), the contribution of unbound particles to the surface brightness profiles should be marginal at best, and we therefore opted for a cut-off at 0.49. The pixel scales for the synthetic images are 12.26 and 5 mas for the SPHERE and ALMA observations, respectively.

During the fitting process, the sampler will draw values for each of the four free parameters, and synthetic images will be computed for both SPHERE and ALMA. Since our goal is to identify the spatial scales that are most representative of the observations, rather than creating one single image by summing over the whole range of β, we instead create n_β images, sampling n_β intervals of β values (linearly spaced between β_min and β_max). For each set of free parameters (a₀, δ_a, and both ψ values), we then find the linear combination of these n_β images that best reproduces the observations (using the lmfit package, Newville et al. 2021, which also provides the confidence interval for each weight). In practice, this means that we are relaxing the assumption of a fixed grain size distribution from the cross-section model, and are re-evaluating the contributions of each β interval.

For the ALMA observations, the lmfit modeling is performed in the u – υ space. Since the Fourier transform is a linear operation, we can compute the visibilities for each n_β images and find the best combination to fit the observations by scaling them. Following Cataldi et al. (2020), when computing the images, we also included an offset to the model, of δ_RA = –0.06″ and δ_Dec = 0.04″.

For the SPHERE observations, we still need to account for the polarized phase function. Instead of using a parametrized one (e.g., Henyein-Greenstein), we follow the approach outlined in Olofsson et al. (2020), where the phase function is an output of the modeling process (and therefore does not depend on the optical properties of the particles). For a given set of free parameters, we run a first model with an isotropic without any polarized phase function and find the best linear combination of the n_β frames using lmfit. Since the phase function is not included, the resulting image will therefore trace the dust density distribution. For both the observations and the model, we then compute the scattering angle for each pixel, and estimate the brightness profile as a function of θ. The phase function that will best reproduce the observations is estimated as the brightness profile of the observations divided by the one of the model (as detailed in Olofsson et al. 2020 only one iteration is required to derive the most suitable phase function). We then compute a new model, this time using the phase function that was just estimated for all the n_β images. We again find the linear combination of those images minimizing the differences compared to the SPHERE DPI observations.

Once a model has been computed for both the SPHERE and ALMA observations, we compute the joint goodness of fit. Because the ALMA and SPHERE observations do not have the same degrees of freedom⁵, one dataset may dominate the total X² over the other. Therefore, prior to starting the modeling, for both datasets we first compute the X² for the null hypothesis, and compute a scaling factor so that the largest of the two values is equal to the other one. In practice, we multiply the X² from the ALMA observations by a factor 0.22. Finally, the likelyhood for this realization of free parameters is returned to the emcee sampler, and another set of free parameters can be drawn. It would have been ideal to avoid the two-fold approach (emcee combined with lmfit) and include the re-weighting of the size distribution as free parameters in the emcee fitting. Unfortunately, computing one model is quite expensive, and ensuring the convergence of the fitting process over a 14 dimension space would be too costly.

4. Results and discussion

4.1. Disk geometry

Table 1 lists the free parameters, the uniform linear priors, as well as the best-fit results along with their uncertainties. The density distributions are displayed in Fig. A.l, using the corner package (Foreman-Mackey 2016) and show that while both a₀ and δ_a are well constrained, the probability density distributions for the two scale heights ψ_sphere and ψ_alma are most likely upper limits, meaning that the disk is not vertically resolved in either dataset (as already mentioned in Olofsson et al. 2022 regarding the SPHERE observations). For these two parameters, for the best fit model we took the value at the peak of the distributions (ψ_sphere = 0.005 and ψ_alma = 0.007 rad). For the confidence interval, we estimated the bin of the histogram for which the integral of the distribution up to that bin reaches 84% of the total integral (≤0.014 and ≤0.071, respectively). We find that the birth ring of the disk is best described by a normal profile centered at ~l″(129.7au) with a standard deviation of 0.22″(28.5 au), in agreement with the modeling results of MacGregor et al. (2018), Bhowmik et al. (2019), and Cataldi et al. (2020). From a₀ and the full width at half maximum of the radial profile, we find a fractional width (FWHM/a₀) of 0.52. Even though the almost edge-on inclination is not ideal to constrain the width of the disk, this is in good agreement with other disks observed with ALMA (median width of 0.74, Marino 2021, 2022). When summing the total flux in the best-fit model for the ALMA observations, we derive a total flux of 23.4 ± 0.3 mJy, compatible with the value reported in Cataldi etal. (2020, 22.0± 2 mJy).

The upper panels of Fig. 1 show the observations, best-fit model, and residuals for the SPHERE DPI observations. Most of the signal from the disk is removed, though some residuals remain, both positive and negative along the major axis of the disk (similarly to Duchêne et al. 2020 and Olofsson et al. 2022). The lower panels of the same Figure show the results for the ALMA observations. To convolve the best-fit model, we replaced the observed visibilities by the ones of the model, and used tclean to compute the image. For the residuals, we proceeded the same way, but this time subtracting the complex visibilities of the model to the observed ones⁶. The contours on the residuals are [–1,1,2] σ (and there is nothing above 3 σ). There is still some signal, at the 1 –1.5 σ level above and below the major axis of the disk, which most likely explains why Cataldi et al. (2020) derived a smaller inclination when modeling the observations.

4.2. The nature of the mm halo

This sub-section focuses on the halo that was reported in ALMA Band 6 observations by MacGregor et al. (2018). We first demonstrate that the halo is also detected in the Band 8 observations, before discussing whether is it composed of mm- or μm-sized dust grains. Afterwards, we discuss what this implies for the optical dust properties and their emissivity at long wavelengths, and discuss its possible origin and detectability.

Fig. 2 Normalized flux contributions as a function of β for the ALMA observations. The red line shows the contributions derived from the cross-section model. The open symbols show the contributions of the β intervals for the best fit model when relaxing the assumption on the grain size distribution. The values that are below 2.5 × 10−2 are considered as upper limits, and are represented by downward triangles. The horizontal bars show the width of the β intervals, and the vertical bars are the lσ confidence interval. The other solid lines show the contributions for models that account for the optical properties of the particles (see text for details).

4.2.1. Recovering the halo in Band 8

Before further describing the results of our modeling, we compute a reference model using the best-fit parameters listed in Table 1. But for this model, instead of evaluating the weights of each β intervals, the size distribution is fixed to the canonical differential power-law dn(s) ∝ s^{−3 5}ds (assuming β ∝ s⁻¹ it is equivalent to dn(β) ∝ β^3/2dβ, Lee & Chiang 2016). This model, referred to as the cross-section model since it should account for most of the effects (size, distribution, cross section, extended lifetime), will serve as a guide to interpret our findings. Figure 2 shows the relative contributions to the thermal emission in the sub-mm for each of the n_β images, normalized to their maximum. We first focus on the red line which shows the flux contribution of each β interval obtained for this cross-section model, which does not take into account the optical properties of the dust grains and has a fixed grain size distribution (the other solid lines will be discussed later on). It should first be noted that this profile remains relatively flat for β ≳ 0.2 despite the apocen-ter of those grains being farther and farther away from the star and their temperature thus becoming lower. For the integrated flux, this effect is counter-balanced by (i) the top-heavy size distribution and (ii) the correction factor α, accounting for these grains’ extended lifetime. The combination of those two factors results in a flux multiplied by a factor [β(1 – β)/(1 – 2β)]^3/2, which increases significantly for large β values.

The results when relaxing the constraints on the size distribution (linear combination obtained using lmfit, instead of imposing the size distribution) are shown as open black symbols in Fig. 2. These symbols show the flux contributions required to reproduce the ALMA observations as a function of β (contributions below 2.5 × 10⁻² are represented by downward triangles and are considered upper limits). The ALMA observations can be best reproduced using only two intervals of β values, β₀ = 0.034 and β₈ = 0.42, these two intervals contributing to ~90 ± 4 and ~10 ± 1% of the total flux, respectively. While it is not surprising that the total flux is dominated by the main belt of the disk (large grains having low β values), the observations are best explained when also considering particles with large β values. Since those particles have large eccentricities, they contribute to a spatially extended emission, whose contribution is not negligible (~ 10%). This requirement is purely driven by the radial extent of the disk. This means that the halo reported in MacGregor et al. (2018) is therefore also detected in the B and 8 observations.

4.2.2. The size of the dust particles

Up until now, we only made very few assumptions that directly relate β to the grain size s. The only occasions where we needed to use the grain size were when computing the images as we need the cross-section and surface area of the grains (for scattered light and mm observations, respectively). If we want to further discuss the properties of the dust grains, we need to relate the β values with typical sizes s. For instance, up until now, we cannot discard the possibility that the high-beta particles contributing to ~10% of the total flux are in fact large mm-sized grains that acquired high eccentricities by another mechanism than radiation pressure. If mm-sized grains are responsible for the halo (particles with β = 0.42), our modeling results suggest that they need to have an eccentricity of β₈/(1 – β₈) ~ 0.72 and a semi-major axis (1 – β₈)/(1 – 2β₈) ~ 3.5 times larger than the semi-major axis of the parent body they were released from. Since we find the parent belt to be at 1″, corresponding to 129.7 au, this translates into an apocenter distance of ~800 au for the mm-sized particles. This would therefore require a physical mechanism able to excite the eccentricities of large mm-grains to such high values.

As mentioned previously, the presence of CO and C I gas has been reported for the disk around HD 32297, and MacGregor et al. (2018) suggested that gas drag could be a possible mechanism to alter the dynamics of the large grains. Indeed, Olofsson et al. (2022) showed that gas drag can have an impact on the radial distribution of particles in debris disks. However, they also showed that large grains are largely unaffected by the gas drag force and do not migrate outward. As a matter of fact, large particles are expected to migrate inward (see also Krivov et al. 2009). Alternative mechanisms mentioned in MacGregor et al. (2018) could be interactions with the interstellar medium or planet-disk interactions. However, as discussed by the authors, the former should predominantly act on small particles, while the latter has mostly been studied for scattered light observations (e.g., Thebault et al. 2012; Lee & Chiang 2016) and not for large grains. It is therefore quite challenging to identify a mechanism that could excite the eccentricity of mm-sized grains up to ~0.7. Consequently, following the principle of Occam’s razor, we now assume that large β values do correspond to small dust particles, as radiation pressure can naturally and efficiently increase the eccentricity of those grains.

4.2.3. Optical properties of the dust particles

The assumption that the halo that is observed at mm wavelength is composed of small particles is however not free of challenges, as noted by MacGregor et al. (2018), as small dust grains are usually poor emitters at long wavelengths. To further investigate this, we used the optool⁷ package (Dominik et al. 2021, Toon & Ackerman 1981) to compute the absorption and scattering efficiencies (Q_abs and Q_sca, respectively), and the asymmetry parameter ɡ_sca for grains of different sizes. We used a mixture of amorphous silicates (optical constants from Dorschner et al. 1995) and amorphous carbon (Zubko et al. 1996). The mixing is done with volume fractions of 60% amorphous silicates, 15% amorphous carbon, and 25% porosity (the standard“DIANA” setup, Woitke et al. 2016, commonly used for circumstellar disks). When computing the optical properties of dust particles, we used the Mie theory (Mie 1908) for compact spheres, and the“Distribution of Hollow Spheres” (DHS, Min et al. 2005) to mimic irregularly shaped grains, with different values for the maximum filling factor f_max (to parametrize the departure from spherical grains). We used optool to compute the properties for a single grain size s, and compute the corresponding β value as

$$\beta \left(s\right)=\frac{3L_{\star }}{16\pi G{c}^2{M}_{\star }}\times \frac{Q_{\text{pr}}\left(s\right)}{\rho s},$$(2)

where L_⋆ and M_⋆ are the stellar luminosity and mass (7.61 L_⊙ and 1.78 M_⊙, respectively), G the gravitational constant, and ρ the dust density (evaluated by optool depending on the input parameters). The radiation pressure efficiency Q_pr is equal to Q_abs + (1 – ɡ_sca)Q_sca averaged over the stellar spectrum. Since the relationship between β and s will depend on the shape of the grains, we sampled a range of sizes large enough to ensure that it will cover the range of β between 0 and 0.5 after converting s to β using Eq. (2). Afterwards, we interpolated the Q_abs values for the values of β used in the modeling. The colored lines (other than the red one) in Fig. 2 show the expected flux distribution from the cross-section model (which assumes a size distribution in s^−3.5) multiplied by the absorption efficiencies, the only ingredient that was missing from the previous analysis. The first thing to be noted is that irregularly shaped grains are overall more efficient than compact spheres at contributing to the thermal emission as discussed in Min et al. (2016) and Tazaki & Tanaka (2018). Even though small grains are not efficient emitters at mm wavelengths, considering grains with a maximum filling factor f_max = 0.8 helps increase their contribution to the mm flux, compared to spherical grains.

The second interesting result is that the expected flux distribution of the cross-section model (red solid line) clearly over-predicts the flux required to reproduce the halo (interval corresponding to β = 0.42) by a factor ~4, but that when accounting for the values of Q_abs, those models under-predict the required flux by a factor ~2–3 (other colored solid lines). Ideally, this could be used to better constrain the values of Q_abs of particles of different size⁸. Indeed, the cross-section model accounts for most of the effects, the size distribution, the correction factor, the surface area, and (at first order) the temperature of the grains. Therefore, dividing the flux contributions needed to model the observations by the contributions expected from the cross-section model should in principle yield a robust estimate of the Q_abs values as a function of β. Unfortunately, this is not straightforward, since we do not have a continuous distribution of fluxes as the output of our modeling strategy. The fitter“prefers” to use two discrete bins of β rather than a continuous distribution as a function of β. It may be the case that the contribution of particles with β values in the range [0.3,0.5] is all accounted for in the single β₈ interval, thus artificially increasing its contribution to the thermal emission. A possible explanation being that the observations do not have infinite spatial resolution and are not free of noise, and therefore, the two intervals with β₀ and β₈ are sufficient to explain most of the signal. Moving forward, possible improvements to the modeling approach would be to penalize significant variations of flux between two adjacent bins of β, assuming that they should be correlated with each other (e.g., Jennings et al. 2020), but this out of the scope of the present paper.

4.2.4. The origin and detectability of the halo

Regarding the origin of the halo, Olofsson et al. (2022) showed that gas drag can alter the radial distribution of the dust particles (on top of the effect of radiation pressure). Since the disk around HD 32297 harbors some amount of gas, it begs the question of how this additional drag force impacts our results. But because our modeling approach is focused on the different spatial scales, the exact mechanism responsible for the segregation of the particles does not really matter, whether it is radiation pressure or a combination of gas drag and radiation pressure, as long as the dependence on the size is similar. In both cases, it is the smallest particles that are pushed the farthest away from the star, while the larger grains remain in the birth ring (or even migrate inward). As a matter of fact, Olofsson et al. (2022) showed that in the presence of gas, the surface brightness of the disk (dominated by the small grains) can only become shallower than in the gas-free case, which can only help creating and maintaining the halo. In that case, the correction factor α would be under-estimated as the small grains can survive longer in the disk compared to the gas-free case. Regardless, even if the segregation of the particles is driven by gas drag and not solely radiation pressure, it supports our conclusion that the halo is caused by small dust particles.

It should be noted though that the detectability of the halo highly depends on the spatial resolution (the beam sizes of the Band 6 and 8 observations being comparable). While large grains should be confined to the birth ring (a rather compact region), the small grains are set on highly eccentric orbits, and therefore their contribution to the total flux is more spatially diluted. Even if their contribution to the total flux is significantly boosted by the correction factor [(1 – β)/(1 – 2β)]^3/2, at higher angular resolution, the contribution of the small particles will nonetheless be spread over several elements of resolution and may not be easily detectable above the noise level (see also, Lynch & Lovell 2021). Finally, it is not surprising that the two disks for which a mm halo has been confirmed (HD 32297 and HD 61005, MacGregor et al. 2018) are edge-on disks. The edge-on configuration is indeed more favorable to detect such halos as the line of sight goes through more material, even though Marino et al. (2016) also reported the possible detection of a halo (or outer belt) for the face-on disk around HD 181327.

4.3. The polarized phase function

To model the near-IR SPHERE observations, several h2@ intervals are required, ranging between ~0.1 up to ~0.45, as shown

in Fig. 3. Similarly to Fig. 2, the cross-section model is shown with a red line, whereas the flux contributions required to fit the SPHERE observations when relaxing the constraint on the size distribution are shown with open black symbols. The other colored lines are the contributions from the cross-section model, multiplied this time by the absorption efficiencies Q_sca (also computed using DHS with optool). To successfully reproduce the observations, four intervals of β have to contribute above 10% of the total flux, β₁ = 0.08, β₃ = 0.18, β₇ = 0.37, and β₉ = 0.47. In relative terms (since all the profiles are normalized to their maximum), the impact of accounting for Q_sca is less pronounced than the impact the Q_abS values had for thermal emission. Nonetheless, departure from spherical compact grains help to increase the contribution of the smallest particles to the total flux in scattered light. That being said, Fig. 3 suggests that the models (including scattering efficiencies or not) under-predict the required flux contribution of small dust grains. As mentioned before, this could be related to the use of only a handful of β intervals in the modeling, or it could be related to the presence of gas in the system. The cross-section model does not include the additional effect of outward migration due to the gas drag (which is more efficient on small particles, Olofsson et al. 2022), possibly resulting in an under-estimation of the contributions for large β values. Lastly, the normalization of the profiles could slightly mislead the interpretation of the Figure. All the profiles for the models are indeed normalized to unity for the interval β (their maximum), while the contributions for the best fit model are normalized to the last bin β₉. If instead, all the profiles were normalized to, for instance, the bin β₁, this will shift down all the open symbols by a factor ~1.8, bringing the last point closer to the profiles of the models. As mentioned for the modeling of the ALMA observations, introducing a penalty for loosely correlated flux contributions between adjacent bins could be an interesting improvement to the approach.

The polarized phase function derived from the best fit solution is shown as a black solid line in Fig. 4 and the shaded area correspond to the uncertainties estimated during the fitting process. Even though the disk is almost perfectly edge-on (i = 86.9°), we cannot probe the full range of scattering angles, since the projected semi-minor axis of the disk lies behind the coronagraph, masking both small and large scattering angles. For this reason, and as noted in Olofsson et al. (2020), we cannot normalize the phase function over 4 steradians, and therefore the absolute values of the phase function cannot be properly calibrated. We restrict the range where the shape of the phase function can be trusted to the interval 25° ≤ θ ≤ 160°.

We then fitted the polarized phase function using optool, with the same composition as before, using the DHS model with f_max =0.8. The two free parameters are the size of the grain s (between 0.01 and 5 μm) as well as the porosity fraction $\mathcal{P}$ (between 0 and 0.75, in volume). Since the observed phase function is not normalized over 4π steradians, for each model, we find the scaling factor that best minimizes the differences between the model and the observations (see Eq. (7) of Olofsson et al. 2016). Using emcee we found that the best solution is for s = 0.29 ± 0.01 μm and $\mathcal{P}=0.4\pm 0.2$ (mostly unconstrained), and this solution is shown as a blue dashed line in Fig. 4. While the shape at scattering angles larger than ~90° matches the observed profile relatively well, the best fit solution cannot reproduce the plateau at smaller scattering angles, as the peak of the DHS model is too narrow.

The DHS model is ideal to compute the optical properties of irregular grains with low porosity. But given that we cannot find a satisfying fit to the observations, we tried another model, the Modified Mean Field Theory (MMF, Tazaki & Tanaka 2018), which mimics high porosity grains or aggregates of monomers. We used the same composition as in Sect. 4.2 and for MMF, there are three main parameters; the size of the monomer s₀, the size s of the aggregate, and the fractal dimension D_f. For s₀ we used the default value of 0.1 μm (e.g., Tazaki et al. 2016) and varied s between 0.5 and 5 μm (s cannot be smaller than s₀). The second free parameter is D_f that we vary between 1.1 and 2.1. Unfortunately, there are some limitations when using MMF, especially for small wavelengths and large D_f values, for which the scattering matrix (as well as the asymmetry parameter) cannot be determined reliably⁹ as multiple scattering events cannot be neglected. We therefore could not explore the range 2.1 < D_f ≤ 3. As noted in Tazaki (2021), the range of D_f between 1.7–2.1 would correspond to fluffy aggregates formed by ballistic cluster cluster aggregation (BCCA), and D_f can be in the range 1.1–1.4 for non ballistic CCA clusters (analogous to a linear chain of monomers). Even though we could not compute the phase function for those cases, larger values of D_f (~3) would correspond to a ballistic particle cluster aggregation process (BPCA), leading to nearly homogeneous aggregates with high porosity (Tazaki 2021 and references therein). Using emcee, we found that the best solution is obtained for D_f = 1.73 ± 0.01 and s ≥ 2.5 μm, and the corresponding polarized phase function is shown as a solid red line in Fig. 4. The range of possible values for D_f is narrow, but our best fit value for s is a lower limit, as grains larger than 2.5 μm also provide an equally good fit to the derived phase function. The value D_f ~ 1.7 suggests a“fluffy” structure for the aggregates, which should have been formed by collisions between clusters of comparable sizes (Tazaki et al. 2016).

To better understand the origin of those fluffly aggregates we can turn to studies of comets in the Solar System. As outlined in the review by Levasseur-Regourd et al. (2018), the current paradigm is that cometary dust is mostly composed of aggregates of various“compactness”. Most recently, the Rosetta mission provided valuable insights on the constituents of the 67P/Churyumov-Gerasimenko (67P hereafter) comet. Analyzing data from the Grain Impact Analyzer and Dust Accumulator (Colangeli et al. 2007) instrument, Fulle et al. (2015) identified two different populations of particles; compact particles and fluffy aggregates, the latter having filling factors (the fraction of the particle volume occupied by monomers) as low as 10⁻³. A further analysis of the data by Fulle et al. (2016) predicted that the fractal dimension of the aggregates should be close to or smaller than ~1.87, and that those particles should amount to ~15% of the non volatile volume. In parallel, Mannel et al. (2016) analyzed observations taken by the Micro-Imaging Dust Analysis System (Riedler et al. 2007) on board the Rosetta spacecraft, and even though most of the particles analyzed are quite compact, they detected one fluffy aggregate with an estimated fractal dimension of D_f = 1.7.

Our current understanding of grain growth is that fluffy aggregates must have formed during the proto-planetary phase. The small relative velocities of the sub-|am-sized monomers leads to the formation of aggregates, whose fractal dimension is smaller than 2 (see Blum & Wurm 2008 for a review). Grain growth continues in a“hit-and-stick” regime until we reach the bouncing barrier, which should happen for sizes between 1 cm to 1 mm depending on the location in the disk (Zsom et al. 2010; Lorek et al. 2018), leading to the compaction of the particles. Consequently, fractal particles must be pristine, and must have survived the compaction phase during the formation of planetes-imals. Such particles have most likely been incorporated within the voids between larger pebbles when forming larger bodies. Based on the Rosetta mission observations, Fulle & Blum (2017) suggested that 67P cannot have experienced any catastrophic collision in its lifetime, which would have otherwise resulted in compaction and destruction of the fluffy aggregates. The age of HD 32297 is not well determined, but the star should be older than 15Myr (Rodigas et al. 2014) but younger than ~30Myr (Kalas 2005), a possible analog to the young Solar System. It is therefore plausible that the second dust generation we observe around HD 32297 is the result of collisions of the progenitors of cometary bodies similar to 67P. They would contain a fraction of pristine fractal aggregates within their nuclei, and may have evolved over several Myr without suffering any major collisions (thus avoiding compaction of the aggregates), until they are eventually destroyed to release the particles that we observe. The smallest of those particles are then set on high eccentricity orbit where they can survive for long periods of time. Nonetheless, this explanation is not free of hurdles; Thébault & Augereau (2007) found that for a disk as luminous as the one around HD 32297, the collisional timescale of the particles should be much shorter than the age of the star. It is therefore a timescale that is hard to reconcile with the possible pristine nature of the particles. Measuring the total intensity phase function and deriving the degree of polarization (the ratio between polarized and total intensity) would greatly help confirming that the particles are indeed akin to fluffy aggregates. Unfortunately, we cannot readily use the angular differential imaging observations presented in Bhowmik et al. (2019) as the post-processing of the observations strongly biases the determination of the total intensity phase function, especially near the projected minor axis of the disk. Alternative post-processing techniques, such as reference star differential imaging, making use of a calibration star, might alleviate some of those issues.

It should also be noted though that the typical sizes of the particles detected in the vicinity of 67P are smaller than ~ 1 mm, larger than the size of the aggregates required to reproduce the polarized phase function of HD 32297 (even though our best fit value is a lower limit). Measuring the degree of polarization of the debris disk would provide additional constraints on the properties of the dust particles. As a final remark, as mentioned previously, the modeling approach assumes that the phase function is the same for all the particles, regardless of their sizes, a strong but necessary assumption. Additional observations, both in total intensity and polarization, at optical and near-IR wavelengths (other than J band), might help in that regard.

Fig. 3 Same as Fig. 2, but for the SPHERE observations.

Fig. 4 Polarized phase function derived from the modeling of the SPHERE observations (dark line and shaded area), as a function of the scattering angle θ. The best fit solutions using DHS and MMF are shown in dashed blue and solid red, respectively.

5. Summary

In this paper we presented a new approach to model multi-wavelength spatially resolved observations of debris disks. The emphasis is put on the spatial extent of particles that experience different strength of the radiation pressure force. By parametrizing the model using different intervals of β values, we are able to alleviate the modeling of the relationship between the grain size and ias much as possible (but not entirely). The characterization of the optical properties can then be done a posteriori.

By modeling simultaneously near-IR SPHERE and mm ALMA observations of the disk around HD 32297, we take advantage of the synergy between the different wavelengths, better constraining the location of the parent planetesimals to model both datasets. We find that the disk is best described by a reference radius of 1″ and a standard deviation for the width of 0.22″ and the disk is not resolved vertically in either dataset. The ALMA observations are best reproduced by a combination of two intervals of i values (small and large i), while the SPHERE observations require several intervals of β values to be fitted adequately.

We confirm the presence of extended emission in the ALMA Band 8 observations, corresponding to the halo reported in MacGregor et al. (2018). We show that this halo does not necessarily originate from the thermal emission of large mm-sized grains (which would require a mechanism to excite their eccentricities above 0.7), but rather that small grains just above the blow-out size can explain this extended emission (that amounts to about ~10% of the total flux). This is due to the larger number of such particles, in part because of the top-heavy size distribution, but most importantly, because of their eccentric orbits significantly increasing their lifetime. Both effects contribute to compensating for their poorer emissivity at long wavelengths. Because of the large eccentricity involved, the detectability of the halo should also depend on the spatial resolution of the observations, making it more difficult to detect at larger resolution. We also discuss the fact that irregularly shaped grains can increase the total flux by an order of magnitude compared to compact spherical grains, for the same total dust mass. However, even for irregular grain shapes, models that include the dust opacity values Q_abs are still a factor 2–3 below the measured ~10)% contribution of the halo to the total flux in the mm, which thus remains an open question to be further investigated.

To model the polarized phase function determined from the SPHERE observations, the Modified Mean Field Theory model provides the most satisfactory fit to the observations. We find that it is best modeled by aggregates with a size larger than 2.5 μm, with a fractal dimension of 1.7, suggesting that the particles we detect at near-IR wavelengths are comparable to fluffy aggregates. Such particles have been detected in comets of the Solar System (e.g., 67P) and are considered to be pristine remnants of its infancy. We hypothesize that the particles in the disk around HD 32297 were released from the collisions of similar cometary bodies, that have not suffered many collisions before being eventually destroyed.

Acknowledgements

We thank the referee for their helpful comments, that helped clarify several parts of the paper, especially regarding the assumptions that are implicitly made in the modeling approach. J.O. and A.B. acknowledge support by ANID, – Millennium Science Initiative Program – NCN19_171. J.O. acknowledges support from the Universidad de Valparaíso, and from Fondecyt (grant 1180395). G.M.K. is supported by the Royal Society as a Royal Society University Research Fellow. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00201.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF-Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). We also acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stel-laire (PNPS) of CNRS-INSU in France. This work has also been supported by a grant from the French Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). It has also been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). MRM, HMS, and SD are pleased to acknowledge this financial support of the SNSF. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/ consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020), and Numba (Lam et al. 2015).

Appendix A Corner plots

Figure A.1 shows the 2D density distributions as well as the projected probability density distributions, for the free parameters when modeling the SPHERE and ALMA observations.

Fig. A.1 Corner plot for the modeling results of the SPHERE and ALMA observations

References

All Tables

All Figures

	Fig. 1 From left to right: observations, best fit model, and residuals, for the SPHERE and ALMA observations (top and bottom, respectively). The scaling is linear and the same for all horizontal panels. For the ALMA observations, the beam is shown in the leftmost panel, and the contours on the right panel are at [–1,1,2] σ (σ = 0.3 mJy beam−1).
In the text

Fig. 2 Normalized flux contributions as a function of β for the ALMA observations. The red line shows the contributions derived from the cross-section model. The open symbols show the contributions of the β intervals for the best fit model when relaxing the assumption on the grain size distribution. The values that are below 2.5 × 10−2 are considered as upper limits, and are represented by downward triangles. The horizontal bars show the width of the β intervals, and the vertical bars are the lσ confidence interval. The other solid lines show the contributions for models that account for the optical properties of the particles (see text for details).

In the text

	Fig. 3 Same as Fig. 2, but for the SPHERE observations.
In the text

	Fig. 4 Polarized phase function derived from the modeling of the SPHERE observations (dark line and shaded area), as a function of the scattering angle θ. The best fit solutions using DHS and MMF are shown in dashed blue and solid red, respectively.
In the text

	Fig. A.1 Corner plot for the modeling results of the SPHERE and ALMA observations
In the text