© T. Bhowmik et al. 2019

1 Introduction

Debris disks are a class of circumstellar disks in which the dust content is thought to be continuously replenished by collisions of planetesimals, usually distributed in a single belt or in multiple belts. The dust present in these disks can be observed by scattered light imaging and/or thermal emission. Morphological and photo-spectroscopic analysis provide constraints on the spatial distribution and physical properties of grains present in a planetary system. Even though debris disks are believed to be depleted in gas, a number of them are significantly gas rich, very likely because of second-generation production of gas released when planetesimals collide and produce the observed dust (Kral et al. 2017). In addition, debris disks frequently show identifiable structures like blobs (AU Mic, Boccaletti et al. 2018), warps (β Pic, Mouillet et al. 1997), multiple belts (Bonnefoy et al. 2017; Boccaletti et al. 2019, HIP 67497 and NZ Lup), asymmetries (Mazoyer et al. 2014, HD 15115), and other irregularities in the disk. These structures might be sculpted by transient breakups of massive bodies (Jackson et al. 2014; Kral et al. 2015), stellar flybys (Lestrade et al. 2011), stellar companions (Thébault 2012), or interactions with the inter-stellar medium (ISM; Hines et al. 2007; Kalas 2005), or be owing to the presence of gas (Lyra & Kuchner 2013) within the disk. One of the most frequently invoked explanations for the presence of most of these structures is, however, the perturbing effect of a planet (Thebault et al. 2012). In this scenario, structures can be the telltale signature of undetected planets, such as in the case of β Pictoris b, which was first inferred indirectly from a warp in the disk and later imaged directly confirming the prediction (Mouillet et al. 1997; Lagrange et al. 2009).

HD 32297 is an A-type star (Silverstone 2000) (V = 8.14 ± 0.01, J = 7.687 ± 0.024, H = 7.624 ± 0.051, K = 7.594 ± 0.018) located at 133 pc (Gaia DR2, Gaia Collaboration 2018). The age of the star is estimated to be ≳15 Myr (Rodigas et al. 2014) and <30 Myr (Kalas 2005). A fractional luminosity of L_IR ∕L_⋆ ≥ 2.7 × 10⁻³ found with the Infrared Astronomical Satellite (IRAS, Silverstone 2000) provides evidence for the presence of cold dust arranged in a belt. This high fractional luminosity makes it one of the brightest debris disks known to date (e.g., Thebault & Kral 2019). The disk was first resolved in scattered light by the Hubble Space Telescope (HST) NICMOS as an edge-on system and detected up to radial distances of 400 au (Schneider et al. 2005). A significant surface brightness asymmetry (southwest ansae brighter than northeast ansae) was reported and attributed tentatively to the presence of an unseen planet. Kalas (2005) detected the disk at larger separations (400 to 1680 au) from the ground with Keck in the R band and measured a blue color when comparing to the HST near-infrared (NIR) data.

In apparent contradiction to the scattered light observations of Kalas (2005) and Schneider et al. (2005), mid-infrared (MIR) observations show that the northeast side appears brighter than the southwest side at 0.6″ (Fitzgerald et al. 2007a) and beyond 0.75″ (Moerchen et al. 2007), although the impact of the angular resolution and signal to noise ratio (S/N) in these images could be questionable. Maness et al. (2008) and Mawet et al. (2009) observed a similar brightness asymmetry to Schneider et al. (2005) and Kalas (2005) in millimeter emission and K band respectively, but still with low angular resolution.

At large distances, as observed with Keck and confirmed with HST/NICMOS, the disk appears bowed, which has been interpreted as the possible interaction of small dust grains residing at large distances with the ISM (Debes et al. 2009). Rodigas et al. (2014) confirmed the bow shape of the disk at the L′ band. Recently, Lee & Chiang (2016) showed that similar structures could also be the result of an interaction between a planet in an eccentric orbit and planetesimals perturbed by radiation pressure. The presence of small dust grains in a large halo has been confirmed with very deep HST/STIS observations (Schneider et al. 2014).

With the advance of high-contrast imaging, Boccaletti et al. (2012), Currie et al. (2012) and Esposito et al. (2014) were able to resolve the disk in the NIR as an inclined belt (~88°) located at ~130 au (after correcting for the new Gaia distance). According to Boccaletti et al. (2012), there is no significant brightness asymmetry between the two sides of the disk in H and K bands as observed with NACO, while Currie et al. (2012) suggest the southwest side to be brighter than the northeast side at r = 35–80 au.

High-contrast polarimetry combined with total intensity in the NIR was first achieved by Asensio-Torres et al. (2016) with Subaru/HiCIAO observations, and with the aim to break degeneracies on the geometrical parameters and grain properties. They reported a gap in the total intensity H band data, which was not visible in their polarimetric data.

Finally, the most recent observations were carried out with ALMA by MacGregor et al. (2018, hereafter MG18) who concluded that at millimeter wavelengths the disk is composed of a planetesimal belt with an inner edge at 78 au and a outer edge at 122 au, and an extended halo up to 440 au. The presence of millimeter grains in the halo complicates the understanding of the large bow, which is expected to be populated with loosely bound grains placed on very eccentric orbits by stellar radiation pressure (Strubbe & Chiang 2006; Thébault & Wu 2008) that should be sensitive to interactions with the ISM. We note however that there are some other alternatives and the ALMA halo may be explained as the presence of a scattered disk (Geiler et al. 2019).

Furthermore, CO gas emission (Greaves et al. 2016, MG18) was detected with ALMA corresponding to a total mass of ~ 7 × 10⁻² M_⊕ derived from an optically thin CO isotopolog (Moór et al. 2019). This large quantity of gas could interact and drag the dust in this system (Takeuchi & Artymowicz 2001). The smallest bound grains would be pushed backwards and the unbound grains would be slowed down and may accumulate in greater quantity than in gas-depleted systems.

Modeling of the spectral energy distribution (SED) suggests that the main belt is populated with sub-micron grains (Fitzgerald et al. 2007a; Currie et al. 2012), which is in agreement with the color index measured between visible and NIR by Kalas (2005). The presence of such sub-micron grains is unexpected, because they should be smaller than the blowout size due to radiation pressure, and should therefore be ejected on very short timescales (Kral et al. 2013). Moreover, taking into account the size of the belt inferred by Boccaletti et al. (2012) and using photometry from Herschel, Donaldson et al. (2013) were able to draw constraints on the grain composition. These latter authors concluded that the system is made of an inner warm belt at ≥1.1 au, and an outer ring at 110 au populated with ≥2.2 μm, high-porosity grains consistent with cometary-like composition. However, this result was questioned by Rodigas et al. (2014) who favored pure icy grains instead.

In this paper we present new high-contrast imaging observations of HD 32297 with Spectro-Polarimetric High-contrast Exoplanet Research (SPHERE, Beuzit et al. 2019) at VLT in Chile, in the NIR from Y to H band complemented with polarimetric observations in the J band. With the new observations at NIR wavelengths, we resolve the disk at separations as close as 0.15″ and perform a spectral analysis. We briefly discuss the observational technique and the data reduction for total intensity and polarimetric data in Sect. 2. The disk features seen in the images are described in Sect. 3. The procedure of the modeling and parametric study is presented in Sect. 4. In Sect. 5, we describe how the photometry is retrieved from the images in the various spectral channels, for both total intensity and polarimetry. The resulting spectrum is compared to various grain models in Sect. 6. Finally, we discuss the implication of the colors in terms of grain sizes in Sect. 7.

2 Observation and data reduction

2.1 SPHERE

SPHERE is an instrument installed at the VLT for high-contrast direct imaging of giant planets around young nearby stars (Beuzit et al. 2019). Its combination of extreme adaptive optics (AO, Fusco et al. 2014) and advanced coronagraphy has allowed observation of circumstellar disks from the ground. SPHERE consists of three instruments: the Infra-Red Dual-beam Imager and Spectrograph (IRDIS, Dohlen et al. 2008), the Integral Field Spectrograph (IFS, Claudi et al. 2008), and the Zurich IMaging POLarimeter (ZIMPOL, Thalmann et al. 2008). IRDIS is a dual band imager using two narrowband or broadband filters in the Y, J, H, or K bands (Vigan et al. 2010). The IFS is a spectro-imager which delivers 39 simultaneous images across the Y J (IRDIFS mode) or Y JH (IRDIFS-ext mode) bands. These two instruments can be used for parallel observations with the IRDIFS mode.

2.2 Observations

HD 32297 was observed with SPHERE on December 19, 2016, in the IRDIFS mode using pupil stabilization in order to take advantage of angular differential imaging (ADI, Marois et al. 2006)for calibration of stellar residuals during post processing. An apodized Lyot coronagraph N_ALC_YJH_S (Carbillet et al. 2011) with a diameter of 185 mas allows attenuation of the starlight in the AO corrected radius (0.84″ at 1.65 μm).

For IRDIS the observations were performed in the broadband H filter (1.625 μm central wavelength, 0.29 μm filter width). IRDIS has a field of view (FOV) of 11″ × 11″ and pixel size of 12.25 ± 0.02 mas. For IFS we used the YJ mode (0.95– 1.35 μm), which provides a spectral resolution R ≃ 54. The IFS has a FOV of 1.73″ × 1.73″ and a pixel size of 7.46 ± 0.02 mas (Maire et al. 2016). The atmospheric conditions were good with seeing below 0.8″ and correlation time of about 8–9 ms. The observing sequence is as follows: point-spread function (PSF; with the star outside the coronagraphic mask with a neutral density ND2), star center (coronagraphic image with four crosswise replicas of the star created with the deformable mirror to monitor, Langlois et al. 2013), long science coronagraphic exposures, second PSF and sky background. Further details on the observation log are provided in Table 1.

Furthermore, HD 32297 was observed with IRDIS dual polarimetric imaging (DPI) mode a few days apart in field stabilized mode to measure the polarized flux of the disk. Several polarization cycles were taken in the J band (1.25 μm), each consisting of exposures at four orientations of the half wave plate (0°, 22.5°, 45°, 67.5°) to measure the full Stokes parameters. The observing sequence is similar to that of the IRDIFS mode for IRDIS; the two filters of IRDIS are replaced with polarizers to split the polarization into two orthogonal directions (Langlois et al. 2014).

2.3 IRDIFS data reduction

The preliminary data reduction, including flat-field corrections, sky and dark subtractions, star-centering using waffle pattern, bad-pixel removal, distortion correction (Maire et al. 2016), and wavelength calibration, is done at the SPHERE Data Centre¹ (Delorme et al. 2017) using the data reduction and handling (DRH) pipeline (Pavlov et al. 2008; Mesa et al. 2015) for the total intensity data obtained in IRDIFS mode. This step provides a data cube which is further processed with ADI techniques based on several algorithms, such as Karhunen-Loève Image Projection (KLIP, Soummer et al. 2012), classical angular differential imaging (cADI, Marois et al. 2006), and template locally optimized combination of images (TLOCI, Marois 2015), using two pipelines, namely the SpeCal pipeline (Galicher et al. 2018) and another one developed by Boccaletti et al. (2015).

Karhunen-Loève Image Projection is an algorithm based on the principle component analysis (PCA). The algorithm involves reformatting the science data into a covariance matrix to remove redundancy. Eigenvalues and eigenvectors are calculated for the covariance matrix which is used to form the KL (Karhunen Loève) basis, which itself is truncated to a given number of modes (Soummer et al. 2012). The reference image is built by projecting the science data onto the truncated KL basis and then subtracting them out frame by frame. The result is derotated according to the parallactic angle variation as in any other ADI techniques allowing reduction of the stellar halo. The final image is normalized to the maximum of the measured PSF.

In this paper, we work with the KLIP reduced data truncated to ten modes using the pipeline developed by Boccaletti et al. (2015) as it provides the optimal S/N.

2.4 Dual polarimetric imaging data reduction

The stokes vectors Q and U are calculated from the polarization cycle observed through each half-wave plate position. To mitigate the instrumental polarization we use the double subtraction method (Tinbergen 1996). The final Q and U vectors can be written as follows: $$\begin{array}{l}Q=\frac{Q^{+}-Q^{-}}{2},\hfill \\ U=\frac{U^{+}-U^{-}}{2},\hfill \end{array}$$(1)

where Q⁺, U⁺, Q⁻, U⁻ are obtained from observations through the half wave plate positions at 0°, 22.5°, 45°, 67.5°.

We retrieve the Q_ϕ and U_ϕ azimuthal Stokes vectors which are expressed in polar coordinates as explained in Schmid et al. (2006). The Q_ϕ and the U_ϕ vectors can be written as: $$\begin{array}{l}{Q}_{\varphi }=-Q\cos 2\varphi +U\sin 2\varphi ,\hfill \\ U_{\varphi }=-Q\sin 2\varphi -U\cos 2\varphi ,\hfill \end{array}$$(2)

where ϕ is the azimuthal angle with respect to the center of the star. We recover the disk signal in the Q_ϕ vector and the noise in the U_ϕ vector under the assumption that the disk is optically thin and undergoes only single scattering. Further correction to instrumental offset is not done as the U_ϕ image is dominated by noise and has very low true disk signal.

Fig. 1 Top: S/N map of total intensity with IRDIS in BB_H filter. Bottom: S/N map of polarimetric image observed with IRDIS in BB_J filter. The dashed line indicates approximate position of the ansa in both images and the arrows indicate front and back sides of the disk. Both images are rotated to 90° – PA and cropped at 6″ × 1.3″. The color bar shows the intensity in both the images.

3 Morphology of the disk

3.1 General description

The full extension of the disk in total intensity (H band) and polarimetry (J band) is shown in the S/N map in Fig. 1, while a smaller FOV is displayed in Fig. 2 for IRDIS (H band), IRDIS-DPI (J band), and IFS (Y J band, collapsed image). The S/N map is the ratio between the reduced image and the azimuthal standard deviation of the same image.

The disk of HD 32297 is relatively bright compared to the stellar residuals. We detect the disk extending to the same distance as observed by HST-NICMOS (Schneider et al. 2005, up to 3.3″ from the star). The disk is almost edge-on and the shape appears globally symmetrical in both the northeast (NE) and southwest (SW) sides as seen in both the total intensity and polarimetric images. Some local asymmetries between the NE and SW sides of the disk are discussed in detail in the photometric analysis in Sect. 5.

The IRDIS images feature a concavity at large separations towards the northwest, which becomes significant at a stellocentric distance of about 2″. This pattern is observed in both the IRDIS BB_H image and the DPI BB_J image (as well observed in Fig. 1). This concavity or bow shape was also observed with HST/NICMOS as mentioned in Schneider et al. (2014), and attributed to the interaction with the ISM of small particles on very eccentric orbits.

On closer inspection at a smaller scale, inside 1″ the disk takes a half-elliptical shape indicative of an inclined ring (likely the planetesimal belt), as suggested first by Boccaletti et al. (2012) and Currie et al. (2012). The high quality of SPHERE images now allows us to resolve this ellipse, with the innermost part of the disk observed at a stellocentric radial distance as close as 0.15″ (compared to 0.5–0.6″ in Boccaletti et al. 2012). We measure the position of the ring ansa at 0.8– 0.9″. The asymmetry with respect to the major axis is reminiscent of forward scattering suggesting that the bright part is the front side of the disk. Interestingly, the S/N map (Fig. 1) built from the IRDIS BB_H image suggests that the back side of the disk is also detected which is also visible in the IRDIS BB_H intensity image (Fig. 2). This is discussed in Sect. 4.

Even though the photometry of the disk is impacted by the so-called ADI self-subtraction (Milli et al. 2012), the intensity along the disk varies monotonically, and so does not feature any sign of a gap contrary to the observation reported in Asensio-Torres et al. (2016). Aditionally, the stellar residual halo in the SPHERE images is much lower than in HiCIAO observations, confidently ruling out such a gap. This is confirmed with our photometric analysis discussed further in Sect. 5.

The IFS image features a highly symmetrical ring on each side of the minor axis, inside the achievable FOV corresponding to stellocentric distances of <0.9″.

In polarimetry, the intensity along the disk varies differently than in total intensity because the scattering phase function is vastly different. As a result, we observe a peak of intensity located at ~ 0.8″ and associated to the location of the disk ansa. Due to the absence of self-subtraction in DPI, strong signal is observed at small angles <0.15″ (Fig. 2). At this stage, it is not yet clear whether this signal is produced by stellar residual or by the disk itself, although it is moderately visible in the S/N map (Fig. 1). There are two possible explanations for this signal to be a true disk signal, one being a strong forward scattering peak at small angular separation, and another being a detection of a potential inner belt.

Fig. 2 Inner part of the disk with IRDIS in BB_H (KLIP, top), IRDIS in DPI (Qϕ, middle), and IFS (YJ combined, bottom). The top two images are cropped at 3″ × 0.8″ and the IFS image is cropped at 2″ × 0.5″. All images are rotated to 90° – PA. All images are scaled linearly in the range [ − 1 × 10−5, 1 × 10−5 ].

3.2 Position angle of the disk

The position angle (PA) is determined by a method developed for edge-on disks as presented in Lagrange et al. (2012). The science data are first rotated to an initial guess of the PA, and a Gaussian profile is fitted perpendicular to the disk mid-plane in a range of angular separations from which the spine of the disk and the local slope, globally or separately for the NE and SW sides, are derived. This process is repeated until the slope reaches a minimum corresponding to a horizontal disk and providing a measurement of the actual disk PA. The average PA retrieved for both total intensity and polarimetric data is 47.60° ± 0.2° and 47.50° ± 0.15°, respectively. The errors consist of measurement error for the applied method and the additional TN uncertainty of 0.1°.

The local minima in the spine of the disk can determine the position of the ansae (Mazoyer et al. 2014). For higher accuracy, a simple ellipse fit to the spine of the disk can give the position of the ansae corresponding to the semi-major axes of the ellipse. Additionally, the inclination can be derived from the semi-major and semi-minor axes of the fitted ellipse. The ellipse fitting the total intensity data is centered at (−0.06″,−0.008″) and (−0.02″,−0.015″) for the polarimetric data. The plots are presented in Fig. 3. We do this ellipse fit to the spine measured on all the images reduced by the methods, TLOCI (Marois 2015), KLIP 3,5,10 modes (Soummer et al. 2012), classical ADI (Marois et al. 2006), No ADI (Galicher et al. 2018), and DPI Q_ϕ (Schmid et al. 2006), and the dispersion between these measurements is used as error. The position of the ansae measured from the spine is 126.4 ± 12.8 au and the inclination is 88.4° ± 0.6°.

4 Geometrical modeling

4.1 Modeling total intensity with a single phase function

Angular differential imaging techniques induce biases to the disk photometry due to self-subtraction. To overcome this problem and recover the unbiased photometry of the disk we proceed with forward modeling, which involves the generation of synthetic images from a model, given some parameters, undergoing a similar post processing to that used on the original data. We used the GRaTer model to generate synthetic images of debris disks (Augereau et al. 1999). The model assumes a ring of planetesimals releasing dust from a collisional cascade and located at a distance r from the star. The position R₀ is where the dust density peaks corresponding to the position of the ansae and we then impose that the dust density distribution scales radially as r^α_in inwards and r^α_out outwards. The vertical distribution is fixed to be a Gaussian function, the height of which is controlled by the parameter H₀ obtained at a radius R₀.

From this geometrical prescription, GRaTer calculates, for a given model, the resulting scattered light image taking into account the scattering angle θ, which depends on the inclination and PA. The phase function is given by the Henyey Greenstein (HG) function (Henyey & Greenstein 1941) for total intensity $$f_I(\theta )=\frac{1-g^2}{4\pi {(1+g^2-2g\cos \theta )}^{3/2}},$$(3)

where g is an anisotropic scattering factor that we leave as a free parameter. For 0 < g < 1, the scattering by dust particles is predominantly forward (isotropic if g = 0), and conversely backward for − 1 < g < 0.

The free parameters of our model are listed below. The initial guesses are obtained from previous work on this system (Currie et al. 2012; Boccaletti et al. 2012) and first-order estimations from the SPHERE images, totalling 23 040 models.

• inclination i (°): 87.5, 88.0, 88.5, 89.0;

• position of the ansae R₀ (au): 115, 120, 125, 130, 135, 140, 145, 150;

• power-law index α_in: 2, 5, 8, 10;

• power-law index α_out: −4, −5, −6, −7, −8 ;

• HG parameter g: 0.4, 0.5, 0.6, 0.7, 0.8, 0.9;

• disk aspect ratio h = H₀∕R₀: 0.01, 0.015, 0.02, 0.025, 0.03, 0.035.

The PA is kept constant at 47.6° as a result of the analysis presented in Sect. 3.2. Each GRaTer model is first injected into an empty data cube void of noise and convolved with the measured PSF. To account for the ADI photometric bias in the forward modeling, a given GRaTer model defined by a set of parameters is projected onto the KL basis which was formed and used for the data. Each model is then normalized to the maximum PSF intensity as formerly done for the data (Sect. 2). Looking for the best-fit model implies that we compare the science image with a series of model images in an appropriate region encompassing the pixels containing disk signal. A numerical mask ~ 0.15″ × 3″, aligned with the disk mid-plane, is applied, which is compared to observations in a χ² fashion. The central part with a stellocentric distance lower than ~0.15″ is removed from the mask. The northeast and the southwest parts are analyzed separately. A comparison between the science image and the noiseless model masked with an effective aperture is credible as the disk has high S/N assuming that there is no over-subtraction due to the combination of noise and disk. Also, in Boccaletti et al. (2019) a comparison between forward modeling with noiseless models and the injected model into the science image at a different PA with a certain flux resulted in similar χ² values. Therefore, we refrain from doing the latter due to reduced computational momentum.

The reduced χ² is calculated between the science image (S_i,j) and the models (M_i,j) at i, jth pixel and summed over the number of pixels in the mask (N_data). The calculation is described as follows: $${\chi }_{\nu }^2=\frac{1}{\nu }{\displaystyle \sum _{i,j=1}^{N_{\text{data}}}}{(\frac{S_{i,j}-a.M_{i,j}(p)}{{\sigma }_{i,j}})}^2,$$(4)

where ν is the degree of freedom equal to N_data − N_param, where N_param is the number of free parameters and p the parameter space explored in these models. The noise term (σ_i,j) is derived from the azimuthal standard deviation in the image masking the disk. In high-contrast imaging, the noise has a spatial structure related to the stellar residuals, varying with wavelengths, while χ² minimization applies to Gaussian errors and linear models. These two conditions are not rigorously met in our forward modeling technique, which imposes some limitations to this approach. The parameter a is the scaling factor between the data and the model.

To identify the best models, we select 1% of those with the lowest reduced χ² values. This approach is considered instead of taking models corresponding to the 1σ deviation of the reduced χ² distribution using the theoretical threshold of $\sqrt{2\nu }$. This is because there are very few models falling into this latter category, which puts overly strong restrictions on the measurements of error bars. Also, this threshold is theoretically applicable for Gaussian errors and linear models which, as discussed earlier, our models do not follow. As a final output we derive the histogram for each parameter, for the given set of best models. A Gaussian profile is fitted to the histograms and Table 2 provides the peak of the Gaussian and the 1σ deviation from the measured peak as errors. Albiet, in two scenarios a Gaussian cannot be fitted. The first is when the distribution of the histogram isflat. The second is when the number of values explored for the parameter are less than four. The best-fit model taken for creating the residual image as shown in Fig. 4 is the model which has the minimum value of the reduced χ² (χ² = 5.26).

The inclination we find, 88.3° ± 0.3°, is compatible with previous measurements (Boccaletti et al. 2012; Currie et al. 2012) but achieves a higher accuracy. The position of the ansae R₀ is 134.7 ± 9.3 au. The variation of the dust density inwards from the belt (α_in) is relatively difficult to constrain for inclined disks; only limited values are explored for this parameter and evidently the retrieved value of 6.0 ± 4.0 is not constrained. On the other hand, the values for α_out, g, and h are relatively well constrained with values of −6.17 ± 1.17, 0.55 ± 0.13, and 0.020 ± 0.006, respectively.

From Fig. 4, one can see that the residuals are still relatively large, which indicates that the adopted model does not perfectly explain the data. Using two HG phase functions instead of one has been proved effective to better fit inclined disks (Currie et al. 2012; Milli et al. 2017). Similarly, we also try to improve our model by using two HG phase functions to model the disk.

Fig. 3 Spine of the disk measured in total intensity image in BB_H (top) and in polarimetric image in BB_J (bottom). The spines are fitted with an ellipse.

Fig. 4 From top to bottom: IRDIS science image in BB_H, best GRaTer model with the parameters i = 88.5°, αout = −6, αin = 10, R0 = 130 au, g = 0.60, h = 0.020, KLIP processed image for the corresponding GRaTer model, and the residual image. The images are cropped at 3″ × 0.5″ and rotated at 90° – PA. The science image, KLIP processed image, and the residual image are scaled linearly in the range [ − 1 × 10−5, 1 × 10−5 ] and the GRaTer model is scaled linearly in the range [0.0, 0.5].

4.2 Modeling total intensity with two HG phase functions

The previous section uses a single HG function to model the scattering phase function but allowed to explore a large range if parameters. In the total intensity IRDIS image, a faint signature of the back side is visible (Fig. 1). In trying to better model the back scattered grains in the disk, we adopt a new phase function based on two HG functions as shown below. $$f_T(\theta )=w_1.f_I(g_1,\theta )+(1-w_1).f_I(g_2,\theta ),$$(5)

where g₂ is assumed to be negative as it models the backward scattering component of the disk.

We reduced the parameter space guided by our results from the previous section and thus the explored parameters for this case are:

• inclination i (°): 88.0, 88.5;

• position of the ansae R₀ (au): 120, 125, 130, 135, 140, 145;

• power-law index α_in: 2, 5, 8, 10;

• power-law index α_out: −5, −6, −7;

• first HG parameter g₁: 0.4, 0.5, 0.6, 0.7;

• second HG parameter g₂: −0.5, −0.4, −0.3, −0.2, −0.1;

• disk aspect ratio h: 0.02, 0.025, 0.03;

• weight w₁: 0.80, 0.83, 0.87, 0.90, 0.93.

The range of weights w₁ is chosen toencompass the values derived for (HD 32297 with KECK by Currie et al. 2012, w₁ = 0.90) and another inclined debris disk (HR 4796 Milli et al. 2017, w₁ = 0.83).

With the combination of the above listed parameters, 43 200 models are created, for which the reduced χ² is measured and the best models are identified following the same approach as described in Sect. 4.1. The parameters giving the best-fit model for this case are provided in Table 3.

The inclination is very well constrained to 88.2 ± 0.2°, but we highlight the fact that we considered only two values to explore the parameter space as this is relatively well-constrained from previous studies as well as in the previous section (Sect. 4.1). The position of the ansae and the inner power-law index are found to be R₀ = 134.4 ± 8.5 au and α_in = 6.0 ± 4.0, and α_out is −6.0 ± 1.0 which are within the error bars retrieved in the previous section. The value of h = 0.022 ± 0.002 is also consistent with results from the previous section.

We find g₁ = 0.69 ± 0.06, g₂ = −0.3 ± 0.2, and w₁ = 0.81 ± 0.05. These values do not agree with the best fits in K band in Currie et al. (2012), where a similar HG function to Eq. (5) was used to model the phase component with g₁ = 0.96, g₂ = −0.1, and w₁ = 0.9. It should however be considered that our data have better S/N compared to the observations of Currie et al. (2012) and we observe the back side of the disk in our data.

In Fig. 5, we plot our best-fit model together with the data and the residuals found after subtracting them both. The reduced χ² value for the best-fit model is 3.29. We find that our model provides a 31% better match to the data compared to models with single HG phase function. Irrespective of the improvement, some residuals remain, which may be due to limitations that are combinationof several factors. The simplicity of our model with only eight parameters cannot reproduce the complexity of the disk; the variation of the PSF during observation, the possibility of over-subtraction and indirect self-subtraction in the forward modeling (Pueyo 2016), along with nonlinear terms ignored in this study can all add to limitation of our best-fit model.

Independently fitting the disk observed in each spectral channel of IFS with the same set of parameters (R₀, α_in, α_out, h) but g₁, g₂, and w₁ being free parameters, we found that the variation of the anisotropic scattering factors is no larger than 8%. This gives confidence in selecting the same value for all spectral channels. We note that the disk is observed only partly in the IFS FOV, and therefore any degeneracy between morphological and phase parameters is not considered in this test.

Fig. 5 (a) From top to bottom: IRDIS science image in BB_H, best GRaTer model with the parameters i = 88°, αout = −6, αin = 2, R0 = 140 au, g1 = 0.70, g2 = −0.4, w1 = 0.80, h = 0.020, KLIP processed image for the corresponding GRaTer model, and the residual image. All the images are cropped at 3″ × 0.5″ and rotated to 90° – PA. The science image, KLIP processed image, and the residual image are scaled linearly in the range [ − 1 × 10−5, 1 × 10−5 ] and the GRaTer model is scaled linearly in the range [0.0, 0.25].

4.3 Modeling polarimetric images

Polarimetric observations complement total intensity data and depend differently on the morphological parameters. For instance, as shown in Engler et al. (2017), the phase function in polarimetry peaks close to ~ 90° phase angle as opposed to the total intensity phase function, which reaches a maximum at small phase angles. As a result, the ansae has stronger signatures in polarimetry than in total intensity. Modeling the polarimetric image is therefore crucial to derive the geometry of the debris disk (Olofsson et al. 2016).

We use GRaTer to create a grid of 32 256 geometrical models with a polarised phase function f_P (θ), given below: $$f_P(\theta )=f_I(\theta )\frac{1-{\cos }^2\theta }{1+{\cos }^2\theta },$$(6)

where θ is the scattering angle. The phase function used is the combination of the HG function (Eq. (5)) and Rayleigh scattering, to account for the angular dependence of linear polarization due to single scattering of an optically thin disk.

The back side of the disk is not visible in the polarimetric image, and therefore we used a single HG function to construct the phase function in the GRaTer models. The free parameters are as follows:

• inclination i (°): 87.5, 88.0, 88.5, 89.0;

• position of the ansae R₀ (au): 115, 120, 125, 130, 135, 140, 145, 150;

• power-law index α_in: 2, 5, 8, 10;

• power-law index α_out: −3, −4, −5, −6, −7, −8;

• HG parameter g: 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99;

• disk aspect ratio h: 0.01, 0.015, 0.02, 0.025, 0.03, 0.035.

We applied the same procedure as in Sect. 4.1 to measure the reduced χ² and derive the range of best models. The value of the reduced χ² is 3.52 for the best-fit model.

The parameters for the best-fit model of the polarimetric image are given in Table 4. We find that the inclination of the disk is 88.6° ± 0.3° and the position of the ansae (127.9 ± 8.0 au) is consistent with the one obtained in total intensity within error bars. Regarding the slopes of the surface density profile, polarimetry favors a steep inner edge (α_in = 8.1 ± 3.2) while this parameter was essentially unconstrained in total intensity. On the contrary, the outer slope α_out = −3.9 ± 0.8 is flatter.

The anisotropic scattering factor is significantly larger in polarimetry (0.84 ± 0.08) as compared to intensity irrespective of the number of g parameters we considered. Therefore, the dust grains which are probed in these data are more prone to forward scattering in polarimetry. It should be noted that in the case of strong forward scattering, the polarised phase function peaks at a smaller scattering angle as well as at the ansae (Milli et al. 2019), which is visible in our best-fit model as seen in Fig. 7. Figure 6 plots the scattering phase functions for total intensity and polarimetric best-fit models. The scattering angles probed for HD 32297 are between 6° and 175°. The phase function in polarimetry peaks at scattering angles smaller than ~ 90° as seen in Fig. 6 while in total intensity we observe an increase beyond 110° due to the back-scattering for a double HG function.

Finally, although the residuals displayed in Fig. 7 are much lower than for the modeling of total intensity, there is still some intensity left near the ansae, indicating that the model does not perfectly reproduce the disk. The residuals could also be an indication of our preference of larger g values over a possibility of an inner component at a separation of <40 AU. Modeling the polarimetric observation with consecutive inner and outer belts is beyond the scope of this paper.

Fig. 6 Scattering phase function adopted for the total intensity and polarimetric best-fit models.

Fig. 7 From top to bottom: DPI Qϕ science image, best GRaTer model with parameters i = 88.5°, R0 = 125 au, αin = 10.0, αout = −4.0, g = 0.8, h = 0.020, and the residual image. The images are cropped at 3″ × 0.5″ and rotated at 90° - PA. All images are scaled linearly in the range [ − 1 × 10−5, 1 × 10−5 ].

5 Photometry and analysis

5.1 Self-subtraction profile and photometry for the total intensity

There are several possible approaches to retrieving the spectrophotometry of the disk. We compare two methods which we use to retrieve the surface brightness of the disk and discuss the limitations associated with each of them.

Method 1

The first method is similar to the one used to measure the photometry of the HD 32297 disk in Boccaletti et al. (2012) and HD 15115 in Mazoyer et al. (2014). It was developed for highly inclined debris disks and relies on estimating the self-subtraction caused by the ADI process. The edge-on geometry allows simplification of this calculation to a 1D problem. Given the best fit model fitting the data identified in Sect. 4.2, we first extracted the radial profile of the GRaTer model (convolved with the PSF) and its associated KLIP image. The profiles for both the model and the KLIP image are measured in the same numerical mask as for the χ² minimization, and we average the flux in nonoverlapping concentric arcs of 0.15″ vertical width, and four-pixel (~0.05″) length.

The self-subtraction is derived from the ratio of these two profiles (Fig. 8, left) and can be significant at short separations when using KLIP (about 60 at 0.2″ in the H band). The unbiased surface brightness (in counts) of the disk is then obtained by multiplying the self-subtraction profile with the science (KLIP) image profile to compensate for the ADI effect, in each spectral channel of the IFS and each filter of IRDIS.

Self-subtraction could also be deduced directly from the ratio of the images of the GRaTer model (PSF convolved) and its KLIP version instead of using profiles. Milli et al. (2017) used this method in the case of a less inclined disk and using a specific ADI process (masked classical ADI), which minimizes the self-subtraction beforehand. Here, in the case of a highly inclined disk together with KLIP processing, this method would lead to the self-subtraction profile presented in the right panel of Fig. 8, which clearly cannot be used for photometric correction. The large variations are attributed to the strong positive to negative fluctuations resulting from the KLIP processing. Therefore, we conclude that this method of self-subtraction measurement should be avoided for KLIP-processed data.

Fig. 8 Self-subtraction measured for IRDIS in the BB_H filter for the method 1 (see text for details) using the ratio of radial profiles (left) and the ratio of images (right).

Fig. 9 Surface brightness of the disk measured for IRDIS in the BB_H filter for method 1 (left) and method 2 (right). The blue (resp. red) solid line represents the SW (resp. NE) side of the disk. The dotted lines show the errors in the measurements.

Method 2

The second solution does not require evaluation of the self-subtraction, but instead uses one step of the modeling when the reduced χ² is calculated. For the minimum value of reduced χ², the scaling parameter a in Eq. (4) directly provides the scaling factor between the KLIP image of the best-fit model and the data. Therefore, contrary to method 1, the surface brightness profile is evaluated in the scaled GRaTer model (PSF convolved) image instead of the real disk image.

For further calculation we first obtain the stellar flux by integrating over a masked PSF which contains 99.99% of its total flux. The radius obtained for the mask is 0.4″ for the PSF obtained by IRDIS and 0.3″ for IFS. Surface brightness profiles are then converted to magnitude arcsec⁻² taking into account the pixel size and normalizing with respect to the stellar flux. The error bars linearly combine two terms, the dispersion of the stellar+background residual intensity as measured in the mask rotated by 90° relative to the real disk for each radial bin, and the accuracy of the photometric extraction. The latter is estimatedwith a fake disk (same parameters as the best-fit model) injected into the raw data at 90° from the real disk at a contrast level of 5 × 10⁻⁴ and processed in the same way. The photometric extraction is found to be consistent with this initial contrast within ~ 3.3%.

The two methods provide very similar results as seen in Fig. 9. Hence, we derived surface brightness profiles for all spectral channels of IFS and IRDIS. Flowcharts representing the calculations using either method are shown in Fig 10.

At this stage, and contrary to Asensio-Torres et al. (2016), we do not confirm the presence of a dip or break in the surface brightness profiles near 0.75″ on the NE side or 0.65″ in the SW, and therefore we rule out the presence of a gap (Fig. 9).

5.2 Photometry for polarimetric images

In polarimetry, the disk photometry is not affected by self-subtraction and therefore can be directly measured from the image. Contrast is calculated directly from the science DPI BB_J image, instead of a model, using the same process as explained in method 2 of Sect. 5.1.

As a result of the scattering angle dependence (Sect. 4.3), the slope of the surface brightness profile is clearly different in polarimetry compared to total intensity, with a less steep decrease from 0.1″ to about 0.7″ (Fig. 11). At the location of the ansae (~0.9″), the surface brightness shows a peak instead of a break. As advocated in Engler et al. (2019) for the case of the inclined debris disk around HD 15115, this dependence of the phase function provides a better sensitivity for polarimetric data to pinpoint the inner edge or the ansae of a dust belt, or to reveal multiple structures. In the case of HD 32297, the polarimetric surface brightness does not show any particular signs of such multiple belts. The polarimetric surface brightness is about two magnitudes fainter than in total intensity for a stellocentric distance of 0.5″, which translates to a polarized fraction of about 15%.

Fig. 10 Flowchart representing surface brightness calculations using method 1 or method 2. In Method 1 the arrow anchored with “no” depicts that the calculation of SS with the previous process does not work, resulting in Fig. 8 (right). Next, we proceed with another process providing Fig. 8 (left) which is used further and therefore the arrow is anchored with a “yes”.

Fig. 11 Surface brightness profile of the disk measured in the IRDIS BB_J DPI Qϕ science image. The blue line represents the SW side of the disk and the red line represents its counterpart on the NE side. The dotted lines represent the error bars.

5.3 Average spectral reflectance

To derive the reflectance of the grains as a function of wavelength, we converted the surface brightness profiles into contrast with respect to the star and averaged the values between 0.2″ and 0.8″, for all spectral channels of IFS and IRDIS, as well as for polarimetric data. Even though, the intensity decreases by ~1 mag arcsec⁻² between 0.2 and 0.8″ for all spectral channels, this separation range is chosen as it corresponds to the minimal and maximal distances where the disk signal starts to dominate over the stellar halo (>0.2″) while encompassing the IFS field of view. The reflectance spectra obtained with the two methods of Sect. 5.1 are plotted in Fig. 12, in which the error bars are obtained from the averaged errors of the surface brightness between 0.2″ and 0.8″ for each wavelength. The main characteristic of the reflectance, irrespective of the method used, is a slow decreasing trend with wavelength which gives the disk a gray to blue color in the YJH spectral range (Fig. 12). Independently of any assumptions on the grain properties, by fitting a straight line to the spectra, we measured a contrast arcsec⁻² variation of − 0.013 ± 0.002 per μm.

Method 1 (Fig. 12, left) clearly produces more dispersion in the Y J band as a result of biases introduced by the self-subtraction estimation. Since method 2 (Fig. 12, right) is based on photometry extracted from models, it provides a smoother spectrum, but a slightly lower reflectance (~12–20%) than method 1. For further assurance, we checked the consistency of the two methods by injecting a fake disk 90° to the real disk PA and performing consecutive photometry for the real disk and the fake disk with both methods. We find that the uncertainty on the measurement between the real disk and the fake disk in Method 1 is twice thatfound when using Method 2. We therefore chose to use Method 2 for subsequent analyses.

Using previous measurements from the literature to confirm the trend of the reflectance on a larger spectral range wouldbe valuable, but photometric measurements are usually derived from various methods, and at various locations on the disk image, which makes the comparison difficult. However, in the case of HD 32297, Ks band data were obtained by some of us using similar (but not quite identical) methods (Boccaletti et al. 2012). In these data, the disk is only detected at stellocentric distances of 0.5–0.8″ for both the NE and SW sides and the surface brightness is ~ 15 ± 0.5 mag arcsec⁻² which translates to a contrast of (1.09 ± 0.5) × 10⁻³ arcsec⁻².

Considering the variation of the disk intensity in SPHERE images between the two ranges of separations 0.2″– 0.8″ and 0.5″– 0.8″ at H band, we can extrapolate this Ks band contrast to (1.99 ± 0.5) × 10⁻³. Therefore, at first order, the photometry in the Ks band confirms the spectral slope of the reflectance (Fig. 13). Observing the disk at the Ks band with SPHERE would further validate the current value of the slope found.

The globally blue behavior of the spectrum is a key element that can help to constrain the size distribution of the grains composing the HD 32297 disk.

Fig. 12 Average spectral reflectance of HD 32297 as measured with method 1 (left) and method 2 (right), for total intensity data in IFS YJ and IRDIS BB_H, as well as polarimetric data in IRDIS BB_J.

Fig. 13 Average spectral reflectance of HD 32297 as measured with method 2, for total intensity data in IFS YJ (in the 0.2–0.8″ range) and IRDIS BB_H (0.2–0.8″ and 0.5–0.8″), as well as polarimetric data in IRDIS BB_J (0.2–0.8″). The NACO K band measurements in the range 0.5–0.8″ (Boccaletti et al. 2012), and an extrapolated value in the range 0.2–0.8″ are over-plotted. The SPHERE Y JH band data are fitted with a straight line (in gray) to estimate a global slope which is (−0.013 ± 0.002) arcsec−2 per μm.

6 Grain modeling

Now that we have constrained the geometry and morphology of the disk, we return to the GRaTer radiative transfer code to constrain the particle size distribution (PSD) of the dust grains. We take as a reference morphology (R₀ ± dr) the best fit obtained for the “full disk” fit displayed in Table 3 and produce synthetic spectra of the disk as a function of grain-related parameters, which we then compare to the observed spectrum in the NIR.

We consider three different grain compositions: astro-silicates, porous astro-silicates with 80% porosity, and a mixture with 50% water ice and 50% astro-silicates. The two crucial free parameters are the minimum grain size s_min and the index κ of the power-law that the PSD is assumed to follow (dn(s) ∝ s^κds). We explore values of s_min between 0.1 and 10 μm with an increment of 0.1 μm and values of κ ranging from −5.0 to −3.0 with an increment of 0.1.

Each model spectrum is interpolated on the forty wavelength channels (39 for IFS and 1 for IRDIS) and globally scaled to the data using the same χ² minimization framework as in Sect. 4.1.

For all wavelength channels (i), the reduced χ² is measured between the average spectral reflectance D_i and each grain model spectrum M_i as given below $${\chi }_{\nu }^2=\frac{1}{\nu }{\displaystyle \sum _{i=1}^{40}}{(\frac{D_i-M_i(p)}{{\sigma }_i})}^2,$$(7)

where σ_i is given by the error bars for each data point in Fig. 9 and p is the parameter space. There are 40 wavelength channels (n_data) and two parameters (n_param; minimum grain size and distribution index), and therefore there are 38 degrees of freedom; ν = 38. The best models are selected as per 1σ deviation of the reduced χ² distribution, which is given by ${\chi }_{\nu ,\text{th}}^2={\chi }_{\nu ,\text{min}}^2+\Delta {\chi }_{\nu }^2$, with $\Delta {\chi }_{\nu }^2=\sqrt{2\nu }$.

The reflectance spectrum corresponding to the best fits obtained for the three considered compositions is displayed in Fig. 14 (right) and the best-fit parameters are provided in Table 5. An important result is that, for all considered compositions, s_min is well below the blow-out limit size s_blow. We indeed always have s_min∕s_blow ≤ 0.085, where we derive s_blow in gm cm⁻³ using the prescription by Wyatt (2008): $$s_{\text{blow}}=0.8\frac{L_{*}}{M_{*}}\frac{2.7}{\rho },$$(8)

where L_* and M_* are the stellar luminosity and mass expressed in solar values, and ρ is the bulk density of the material (given in Table 5). We take L_* = (8.4 ± 0.2) L_⊙ (Moór et al. 2017) and M_* = 1.8M_⊙ (Kalas 2005).

As for the best fit of the slope (Fig. 14, right) of the size distribution, we find κ = −3.79 ± 0.34 for astro-silicates, which is relatively close to the slope expected for collisional steady states (between −3.6 and −3.7; see Gáspár et al. 2012, and references therein). For the other two compositions, we find slightly steeper PSDs, with κ ~−4.5. For all three slopes, it is important to stress that the geometrical cross section, and thus the flux, is dominated by the smallest grains in the PSDs, that is, those close to s_min (Thebault & Kral 2019). Also, s_min = 2.2 μm, the minimum size found by Donaldson et al. (2013), also corresponds to s < s_blow grains given the very high 90% porosity they assumed. To further stress the crucial role of unbound s ≤ s_blow grains, we display in Fig. 14 (left) the best fits that would be obtained when considering a PSD stopping at s_min = s_blow and a canonical slope with κ = −3.5. As can be clearly seen, in the absence of the unbound grains no satisfying match can be obtained, in particular regarding the blue slope of the reflectance spectrum.

These results, especially those regarding the minimum grain size, only weakly depend on the morphology assumed for the disk. Taking a size for the planetesimal belt other than that of Table 3 indeed leads to relatively similar results. This is expected, as the slope of the reflectance spectrum is essentially imposed by the size-dependence of the scattering coefficient Q_sca, which does not depend on the location of the grains with respect to the star. The scattering anisotropy parameters g do depend on this location, but under the assumption made in Sect. 4 that there is no size dependence for g, this would not translate into a different slope of the reflectance spectra for a given PSD. Therefore, we note that at this stage it is safe to derive synthetic spectra without a dependence of anisotropic scattering parameter.

Fig. 14 Best fits of the reflectance spectrum (method 2) obtained for three different grain compositions, for a fiducial case where we force smin = sblow and κ = −3.5 (left), and when smin and κ are free parameters (right).

7 Discussion

7.1 Comparison of geometrical parameters to millimeter observations by ALMA

From the analysis presented in Sect. 4.3 we find that the disk is best described as a relatively narrow belt peaking at 132.3 ± 6.2 au (according to polarimetric data), which agrees with previous scattered light observations (Currie et al. 2012; Boccaletti et al. 2012,after the correction of the distance from the new Gaia measurements). HD 32297 was also observed recently with ALMA in the dust continuum at 1.3 mm by MG18. However, the beam size was about 0.76″ × 0.51″ which is more than ten times larger than the SPHERE angular resolution, meaning that the ALMA resolution is significantly poorer.

While considering a geometrical model that is similar to GRaTer, MG18 concluded that the disk is rather broad, extending from 78 ± 8 to 122 ± 3 au, with a surface density rising as r². The peak in density occurs slightly closer-in than in the SPHERE images, but this discrepancy could be the result of the angular resolution. In any case, the radial dependency of the surface density, even if relatively poorly constrained in our case, is considerablysteeper with SPHERE.

Beyond the planetesimal ring at 122 au, MG18 observed a halo extending out to 440 ± 32 au, where the surface density decreases as r⁻⁶. This is consistent with the SPHERE image in which we observe dust scattering as far as 3.3″ (equivalent to 440 au), with a comparable radial decrease in density. MG18 do not report any asymmetry in this outer part while the aforementioned concavity is obvious at shorter wavelengths. This could again be a resolution effect or because of the sensitivity of ALMA to bigger grains. Nevertheless, subtracting an axi-symmetrical disk model from the ALMA image leaves residuals co-located with the region where the concavity is detected by SPHERE and HST. The presence of millimeter grains at such stellocentric distances would, however, be at odds with the expected mechanism creating such concavity. For the same reason, our value for the disk inclination is significantly more precise (88.2 ± 0.3°) than the one derived by MG18 (83.6° ${}_{-0.4}^{+4.6}$).

In order to check the compatibility of the ALMA and SPHERE models, we used the MG18 parameters α_in, α_out, R_in, R_out, and took the other parameters (g₁, g₂, w₁, h) from our best-fit values (Table 6). As a result, we obtained a much larger χ² value of 8.5 for total intensity compared to our best fit. For polarimetry, χ² = 3.6, which is not far from the value we find for our best-fit model. This is because R₀ = 122 au, as used in MG18, is close to that of our polarimetric best-fit model R₀ = 125 au. We also explored models with g₁, g₂, w₁ as free parameters corresponding to the parameter space of Sect. 4.2 for total intensity and Sect. 4.3 for polarimetry in order to restrict possible degeneracies between dust density distribution corresponding to the ALMA values and phase function. The attempt resulted in larger χ² (8.3 for total intensity and 3.6 for polarimetry) values compared to our best fit indicating that the ALMA model does not match the SPHERE image. The disagreement between the models derived from these two instruments could stem from one of two factors, or both. The first is the angular resolution as mentioned above. The second is that ALMA and SPHERE probe different grain size of which the dynamics can be governed by different processes resulting in two distinct spatial distributions.

7.2 Disk color and sub-micron grains

7.2.1 Quantity of sub-micron grains produced naturally in debris discs

Our new observations and analysis confirm two striking characteristics of the HD 32297 disk: the blue color of the spectrum in the NIR and the significant presence of tiny grains much smaller than s_blow. Moreover, we confirm that there is an intrinsic coupling between these two characteristics, as was also inferred for the HD 15115 (Debes et al. 2008) or AU Mic systems (Augereau & Beust 2006; Fitzgerald et al. 2007b).

This link between a blue NIR spectrum and sub-micron grains has been quantitatively investigated in the recent study by Thebault &Kral (2019). This numerical exploration shows that for bright debris disks with high fractional luminosity f_d, a collisional cascade at steady-state can “naturally” produce a level of unbound sub-micron grains that is high enough to lead to a blue slope of the spectrum in the NIR. This is because for bright and dense disks, the drop in grain number density at the s = s_blow frontier, which is to a first order ∝ 1∕f_d, is much less pronounced than for fainter systems. For a very bright disk with f_d = 5 × 10⁻³ comparable to that of HD 32297, Thebault & Kral (2019) found a profile of the relative NIR L_d ∕L* spectrum that is qualitatively similar to the one obtained here (see Fig. 13 of that paper). However, the blue slope they obtained is not as steep as in the present case, with a flux ratio between the λ = 1 μm and λ = 1.6 μm fluxes that is ~ 1.1, as compared to ~ 1.6 here.

This could indicate an additional source of sub-micron grains, which cannot be explained by the steady-state collisional evolution of the system. One possible cause could be the so-called collisional “avalanche” mechanism (Grigorieva et al. 2007; Thebault & Kral 2018), initiated by the break-up of a large planetesimal closer to the star, which releases large amounts of unbound dust grains that then trigger a collisional chain-reaction as they sandblast at very high velocity through a dense outer disk. The ideal case for an avalanche-producing system is a double-belt configuration, with an inner belt (where the large planetesimal breaks up) at ~1–10 au with f_d ≳ 10⁻⁴, and a bright outer belt with f_d ≳ 10⁻³ (Thebault & Kral 2018). This could match the structure of HD 32297, for which an inner belt of brightness f_d ~ 6 × 10⁻⁴ has been inferred by Donaldson et al. (2013), even though the reality of this inner belt is still debated (e.g., Kennedy & Wyatt 2014). Observational confirmation of an inner belt would need to achieve contrasts significantly higher (a factor of 10) than those currently feasible with SPHERE. In this case, the level of s ≤ s_blow grains would vary stochastically, on a timescale t_av that is roughly a third of the typical dynamical timescale in the disk (Thebault & Kral 2018). This would, however, correspond here to t_av ~ 300–400 yrs, much too long to be observationally monitored. Moreover, it is not guaranteed that the rate at which large planetesimals break up in the inner regions is high enough for such an event to be likely to be witnessed (see discussion in Thebault & Kral 2018).

7.2.2 Effect due to gas on the presence of sub-micron grains

Another possibility is that the system is able to retain s ≤ s_blow grains significantly longer than the radiation pressure blow-out time. If there is enough gas in the system, gas drag could act to significantly increase the time for an unbound grain to leave the system. To check that, we first compute the stopping time for the case where grains are bound to their host star, equal to (Takeuchi & Artymowicz 2001) $$T_{\text{sb}}~2\left(\frac{\rho }{1.5\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-3}}\right)\left(\frac{s}{1\mu \text{m}}\right){\left(\frac{M_{\text{gas}}}{0.1M_{\oplus }}\right)}^{-1},$$(9)

where we assume M_⋆ = 1.8 M_⊙, R = 130 au, ΔR = 50 au to be consistent with results from Sect. 4. We also fix the gas temperature to 30 K and its mean molecular weight to 28 based on Cataldi et al. (2019), where they show that the gas mass is dominated by CO rather than carbon in HD 32297. Accounting for the observed neutral and ionized carbon (in addition to CO, Moór et al. 2019), the total gas mass barely goes above 0.1 M_⊕. If accounting for potentialCO₂ or water being released from planetesimals at the same time as CO (and thus producing extra oxygen not coming from CO and some extra hydrogen), and assuming a solar-system comet-like composition (e.g., Kral et al. 2016), the total gas mass could go up to 0.5 M_⊕ ².

Therefore, small bound grains close to the blow-out limit (i.e., 1–10 μm, see Table 5) will have a stopping time close to 1 and will be affected by gas drag over a few orbital periods before they have time to collisionally deplete. The smallest bound grains will likely (depending on the gas-pressure gradient) move outwards before being collisionally destroyed (Takeuchi & Artymowicz 2001) and will therefore be present for longer than usually assumed in a standard size distribution (e.g., Kral et al. 2013).

For unbound grains, Eq. (9) needs to be adjusted. The velocity of unbound grains can reach $\sqrt{2(\beta -1)}$ of the Keplerian velocity at which they are released initially, which would increase the velocity difference between gas and dust, hence entering the Stokes regime of drag (rather than the Epstein regime, Takeuchi & Artymowicz 2001). On top of that, the stopping time of Eq. (9) is calculated over an orbital period and unbound grains travel almost radially over ΔR, the width of the disk. Therefore, we calculate a new dimensionless Stokes stopping time for unbound grains T_su (based on Takeuchi & Artymowicz 2001), scaled by the crossing time over ΔR and therefore equal to T_sb[2c_s∕v_K][R∕ΔR], where v_K is the Keplerian velocity and c_s the sound speed: $$T_{\text{su}}~0.03\left(\frac{\rho }{1.5\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-3}}\right)\left(\frac{s}{0.1\mu \text{m}}\right){\left(\frac{M_{\text{gas}}}{0.1M_{\oplus }}\right)}^{-1},$$(10)

meaning that 0.1 μm grains will be slowed down significantly before they have time to leave the disk³. The exact orbit of unbound grains interacting with gas is time dependent and is not derived here as it goes beyond the scope of this paper, but generally speaking, an unbound grain would start on a very hyperbolic orbit and eventually be circularized around the star. These grains will then accumulate before being destroyed collisionally. Therefore, we find tantalizing evidence that gas observed in this system may be able to explain the blue color of the disk by allowing small unbound grains to be present for longer.

Fig. 15 Limits of detection in contrast (top) for IRDIS (solid line) and IFS (dashed line), and converted into Jovian masses (bottom, IRDIS only) for two age assumptions (10 and 30 Myr), using the COND evolutionary model.

7.3 Point-source detection limit

We do not detect any point source in the entire IRDIS field of view. The contrast curves at 5σ are measured with the SpeCal pipeline (Galicher et al. 2018) for KLIP-reduced data (Fig. 15). IRDIS provides contrasts of about 10⁻⁵ at 0.5″ and 10⁻⁶ at 1″. The IFS contrasts are similar or slightly better at <0.4″ and then degrade for larger separations. For this reason, we considered the IRDIS contrast curve only to derive the limit of detection in terms of mass. Figure 15 displays the limit of detection using the COND model (Allard et al. 2001), and assuming two possible ages of the system: 10 and 30 Myr. We note that if the system were 10 Myr old we should be able to detect a planet of 3 Jupiter masses at a separation of 0.5″ (and respectively 1.3″ for 30 Myr). The values below 1 M_J are unreliable in the framework of the COND model. Since the noise is estimated azimuthally in SpeCal, the disk itself contributes 40% to the contrast curve in the range 0.2″–1.2″.

8 Conclusions

The findings of this paper can be summarized as follows:

• We observed the debris disk of HD 32297 in the NIR in the Y, J, and H bands out to stellocentric distances of 3.3″, and for the first time as close as 0.15″. We obtained both total intensity and polarimetric images as well.

• At large separations, the disk is characterized by a concavity as reported by Schneider et al. (2014). At shorter separations (<1″), a bow-like shape is reminiscent of a very inclined belt of which we see mostly one side (northwest) due to the forward scattering by the grains. Noticeably, we were able to detect the back side of the disk which we modeled using two HG phase functions. This feature is not observed in polarimetric data possibly due to low S/N.

• Upon first inspection, the disk appears to be symmetrical in NE and SW sides and has no gapped structure in contrast to the claims of Asensio-Torres et al. (2016). This is confirmed unambiguously with our photometric study in which the surface brightness profiles do not show any significant brightness asymmetry between the two sides, or any gap.

• We present two methods for extracting the photometry of inclined disks observed in total intensity. The first method includes the estimation of the ADI self-subtraction in model images and accounting for this bias into the data. We find that this method can induce some irregularities depending on the accuracy of measurement of the ADI self-subtraction. In the second method we calculate a scaling factor between the ADI processed data and its corresponding model. The scaled model is used to measure the photometry instead of the data. As a drawback any departure from the model is not represented in the measurements.

• Comparing total intensity and polarimetry in the J band we derived a polarisation fraction of about 15% which is in accordance with other debris disks.

• From photometric measurements obtained in 40 spectral channels we obtained an average spectral reflectance and conclude that the disk is “gray to blue” color in the YJH spectral range. Using a radiative transfer module in GRaTer we were able to compare this measured reflectance with those expected for a variety of grain sizes and compositions. We found that irrespective of the composition, grains should be significantly smaller than the corresponding blowout size (sub-micron size for astrosilicates).

• Finally, we discussed that the presence of the small grains and the associated blue color of the disk can originate from a combination of several physical processes, including steady-state collisional evolution and the avalanche process. Given the amount of gas in this system (Greaves et al. 2016; Cataldi et al. 2019), we also found that the gas drag can retain smaller unbound grains over a longer timescale.

HD 32297 is amongst the very few known bright and extended debris disks with gas. The SPHERE observations are of unprecedented quality allowing the detection of this disk at high S/N in all spectral channels, and strong constraints to be derived on the grain properties. Confirming the trend of the spectral reflectance would require additional SPHERE observations in the K band in total intensity as well as polarimetric data in the H and K bands. It would be interesting to perform the very same type of observations and data analysis for other gas-rich debris disks and investigate if they share similarities with HD 32297 as an attempt to understand whether the presence of gas can fully explain the dust size distribution.

Acknowledgements

French co-authors acknowledge financial support from the Programme National de Planétologie (PNP) of CNRS-INSU in France. The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). Finally, this work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice) and Observatoire de Paris/LESIA (Paris). J.Ma. acknowledges support for this work was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51414.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.O. and C.P. acknowledge financial support from the ICM (Iniciativa Científica Milenio) via the Núcleo Milenio de Formación Planetaria grant. J.O. acknowledges financial support from the Universidad de Valparaíso, and from Fondecyt (grant 1180395). C.P. acknowledges financial support fron Fondecyt (grant 3190691).

References

All Tables

All Figures

	Fig. 1 Top: S/N map of total intensity with IRDIS in BB_H filter. Bottom: S/N map of polarimetric image observed with IRDIS in BB_J filter. The dashed line indicates approximate position of the ansa in both images and the arrows indicate front and back sides of the disk. Both images are rotated to 90° – PA and cropped at 6″ × 1.3″. The color bar shows the intensity in both the images.
In the text

	Fig. 2 Inner part of the disk with IRDIS in BB_H (KLIP, top), IRDIS in DPI (Qϕ, middle), and IFS (YJ combined, bottom). The top two images are cropped at 3″ × 0.8″ and the IFS image is cropped at 2″ × 0.5″. All images are rotated to 90° – PA. All images are scaled linearly in the range [ − 1 × 10−5, 1 × 10−5 ].
In the text

	Fig. 3 Spine of the disk measured in total intensity image in BB_H (top) and in polarimetric image in BB_J (bottom). The spines are fitted with an ellipse.
In the text

Fig. 4 From top to bottom: IRDIS science image in BB_H, best GRaTer model with the parameters i = 88.5°, αout = −6, αin = 10, R0 = 130 au, g = 0.60, h = 0.020, KLIP processed image for the corresponding GRaTer model, and the residual image. The images are cropped at 3″ × 0.5″ and rotated at 90° – PA. The science image, KLIP processed image, and the residual image are scaled linearly in the range [ − 1 × 10−5, 1 × 10−5 ] and the GRaTer model is scaled linearly in the range [0.0, 0.5].

In the text

Fig. 5 (a) From top to bottom: IRDIS science image in BB_H, best GRaTer model with the parameters i = 88°, αout = −6, αin = 2, R0 = 140 au, g1 = 0.70, g2 = −0.4, w1 = 0.80, h = 0.020, KLIP processed image for the corresponding GRaTer model, and the residual image. All the images are cropped at 3″ × 0.5″ and rotated to 90° – PA. The science image, KLIP processed image, and the residual image are scaled linearly in the range [ − 1 × 10−5, 1 × 10−5 ] and the GRaTer model is scaled linearly in the range [0.0, 0.25].

In the text

	Fig. 6 Scattering phase function adopted for the total intensity and polarimetric best-fit models.
In the text

	Fig. 7 From top to bottom: DPI Qϕ science image, best GRaTer model with parameters i = 88.5°, R0 = 125 au, αin = 10.0, αout = −4.0, g = 0.8, h = 0.020, and the residual image. The images are cropped at 3″ × 0.5″ and rotated at 90° - PA. All images are scaled linearly in the range [ − 1 × 10−5, 1 × 10−5 ].
In the text

	Fig. 8 Self-subtraction measured for IRDIS in the BB_H filter for the method 1 (see text for details) using the ratio of radial profiles (left) and the ratio of images (right).
In the text

	Fig. 9 Surface brightness of the disk measured for IRDIS in the BB_H filter for method 1 (left) and method 2 (right). The blue (resp. red) solid line represents the SW (resp. NE) side of the disk. The dotted lines show the errors in the measurements.
In the text

	Fig. 10 Flowchart representing surface brightness calculations using method 1 or method 2. In Method 1 the arrow anchored with “no” depicts that the calculation of SS with the previous process does not work, resulting in Fig. 8 (right). Next, we proceed with another process providing Fig. 8 (left) which is used further and therefore the arrow is anchored with a “yes”.
In the text

	Fig. 11 Surface brightness profile of the disk measured in the IRDIS BB_J DPI Qϕ science image. The blue line represents the SW side of the disk and the red line represents its counterpart on the NE side. The dotted lines represent the error bars.
In the text

	Fig. 12 Average spectral reflectance of HD 32297 as measured with method 1 (left) and method 2 (right), for total intensity data in IFS YJ and IRDIS BB_H, as well as polarimetric data in IRDIS BB_J.
In the text

Fig. 13 Average spectral reflectance of HD 32297 as measured with method 2, for total intensity data in IFS YJ (in the 0.2–0.8″ range) and IRDIS BB_H (0.2–0.8″ and 0.5–0.8″), as well as polarimetric data in IRDIS BB_J (0.2–0.8″). The NACO K band measurements in the range 0.5–0.8″ (Boccaletti et al. 2012), and an extrapolated value in the range 0.2–0.8″ are over-plotted. The SPHERE Y JH band data are fitted with a straight line (in gray) to estimate a global slope which is (−0.013 ± 0.002) arcsec−2 per μm.

In the text

	Fig. 14 Best fits of the reflectance spectrum (method 2) obtained for three different grain compositions, for a fiducial case where we force smin = sblow and κ = −3.5 (left), and when smin and κ are free parameters (right).
In the text

	Fig. 15 Limits of detection in contrast (top) for IRDIS (solid line) and IFS (dashed line), and converted into Jovian masses (bottom, IRDIS only) for two age assumptions (10 and 30 Myr), using the COND evolutionary model.
In the text