© The Authors 2025

1 Introduction

The direct imaging technique is able to constrain planetary system architectures by probing regions at large (>5 au) separations from the star. This approach can detect the near- and mid-infrared (NIR-MIR) emission of young self-luminous giant planets and resolve circumstellar disks, such as debris disks, made up of leftover material from planet formation processes.

In volume, they consist mainly of small, dust particles resulting from the collisional cascade of larger (i.e., ∼kilometer-sized) planetesimals (e.g., reviews from Wyatt 2008; Krivov 2010; Hughes et al. 2018; Pearce 2024).

To date, space telescopes and ground-based facilities equipped with extreme adaptive optics instruments have spatially resolved the starlight scattered off dust grains in debris disks down to separations of 0.1" and up to about 10". Such facilities include HST (Schneider et al. 1999, 2014; Augereau et al. 1999a; Choquet et al. 2014, 2018; Ren et al. 2017, 2023b), JWST/NIRCam (Rebollido et al. 2024; Lawson et al. 2024), VLT/SPHERE (Wahhaj et al. 2016; Boccaletti et al. 2018; Xie et al. 2022; Olofsson et al. 2022b), and Gemini/GPI (Esposito et al. 2020; Crotts et al. 2024). Such images require specific strategies of observation (e.g., adaptive optics and coronagraphy) and processing to block the stellar light and reveal the faint disk signal. Scattered-light images enable a better understanding of very small (≲ a few μm) dust particles by determining the degree of forward scattering and linear polarization of the dust grains and their color. This provides constraints on the properties of the dust grains, such as their size, shape, and composition.

On the other hand, images of debris disks from mid-infrared to millimeter provide a complementary view because they are sensitive mainly to the cold thermal emission of larger (sub-mm/mm) dust particles; for instance, with JWST/MIRI (Gáspár et al. 2023; Rebollido et al. 2024), Spitzer (Su et al. 2009), Herschel (Löhne et al. 2012; Wyatt et al. 2012; Kennedy et al. 2013, 2015, 2018), and ALMA (MacGregor et al. 2013, 2017; Su et al. 2017; Marino et al. 2017, 2019; Faramaz et al. 2019, 2021). For (sub)mm observations, the location of the large grains seen is thought to track more reliably the parent belt of large, colliding planetesimals, as small grains are subject to stellar radiation effects that make them drift from the parent belt where they are produced, thereby smearing any asymmetry gravitation-ally induced by massive perturbers (Thébault & Augereau 2007). The properties of debris disks, including their morphology and the properties (size, shape, and composition) of the constituent dust grains can therefore be studied more comprehensively by combining multi-wavelength observations.

In total, more than one hundred debris disks have been spatially resolved¹ so far. Among them, there exists a wide diversity of morphologies: narrow or broad rings, circular, eccentric or misaligned rings, arcs, warps, clumps, swept-back wings, needles, forks, or spiral arms (Hughes et al. 2018). A natural question is what causes this diversity (Wyatt et al. 1999; Lee & Chiang 2016), whether this is related mostly to the influence of planetary companions, or to other physical mechanisms such as disk self-gravity (which may or may not require the presence of planets; e.g., Ward & Hahn 1998; Hahn 2003; Jalali & Tremaine 2012; Sefilian et al. 2021, 2023), giant impacts (Jones et al. 2023), collisional dust avalanches (Grigorieva et al. 2007), stellar flybys (e.g., Kenyon & Bromley 2002; Reche et al. 2009; Lagrange et al. 2016), interstellar medium wind (Maness et al. 2009), or processes occurring during planet formation in protoplanetary disks (e.g., Najita et al. 2022).

In this context, the system HD 120326 (HIP 67497), located at 113.3 ± 0.4 pc (Gaia Collaboration 2021) belongs to the Upper Centaurus Lupus sub-group in the Sco-Cen association, and hosts a debris disk around the F0V-type young (16 Myr; Mamajek et al. 2002) star of mass of 1.6 M_⊙ (Chen et al. 2014). Using VLT/SPHERE near-infrared imaging in scattered light, Bonnefoy et al. (2017) reported the existence at intermediate scales (≲150 au) of a dust belt and a putative second dust structure, which could be another belt or a halo of dust grains. In addition, optical data seems to show scattered light of an extended halo, up to a projected separation of 700 au (Padgett & Stapelfeldt 2016; Bonnefoy et al. 2017, using HST/STIS).

Previously, the spectral energy distribution (SED) of HD120326 has been modeled with either a single-(Jang-Condell et al. 2015), or a double-belt (Chen et al. 2014) architecture, by using in particular measurements of the Infrared Spectrograph (IRS) data of Spitzer. From their SED modeling, Chen et al. (2014) found a bright (L_IR/L_⋆ = 1.1 × 10⁻³), relatively cold (127 ± 5 K) belt at 13.9 au, and a dimmer (L_IR/L_⋆ = 1.4 × 10⁻⁴), colder (63 ± 5 K) belt at 116.5 au, by assuming two blackbody models. Alternatively, Jang-Condell et al. (2015) modeled the data with a single belt made of a single grain population of a given size, temperature, and composition (a mixture of olivine and pyroxene). Their best fit was that of a belt with fractional luminosity of L_IR/L_* = 1.5 × 10⁻³ located at 8.8 ± 1.0 au, dust particles of temperature of 124 ± 5 K, size of 15.5 μm, and with a composition converging to pure olivine (instead of a mixture of olivine and pyroxene).

However, both locations of the dust belt are underestimated. This is common for representative radii of dust belts inferred via SED modeling (e.g., Pawellek & Krivov 2015) because small dust particles are inefficient emitters, resulting in higher temperatures than expected based on their separation from their host star. By using SPHERE total intensity observations, for the first time Bonnefoy et al. (2017) imaged the dust belt and derived a larger reference radius (r₀) of 59 ± 3 au, using a forward-modeling approach,. They also reported the possible existence of a second, outer disk component, which could be another belt (r₀ ∼ 130 ± 8 au) or a halo. However, Bonnefoy et al. (2017) were cautious about this second component, as it is faint, and could also be related to self-subtraction effects (Milli et al. 2012) caused by their data processing based on the angular differential imaging (ADI; Marois et al. 2006) technique. This purported second structure was not detected in polarized intensity light using the same instrument, SPHERE, at the same wavelength (1.6 μm; Olofsson et al. 2022b). This further complicates the interpretation of the dust distribution around HD120326.

In this study, we carried out a comprehensive, multi-wavelength study of the morphological and dust properties of the HD 120326 debris disk. We used a combination of optical, near infrared and millimeter data, including new observations using VLT/SPHERE and ALMA and archival data (HST/STIS and SPHERE).

In Sect. 2, we present the data used and our reduction process: SPHERE observations (Sect. 2.1; with SPHERE/IRDIS and SPHERE/ZIMPOL polarized intensity data in Sect. 2.1.1 and SPHERE/IRDIS and SPHERE/IFS total intensity data in Sect. 2.1.2), ALMA observations (Sect. 2.2), and STIS observations (Sect. 2.3). In Sect. 3, we present our joint modeling (Sect. 3.1) and analysis of the SPHERE total and polarized intensity observations. We derived new constraints on the properties of the inner dust belt, including one at 1.6 μm its scattering phase function (Sect. 3.2) and maximum fraction of polarization (Sect. 3.3), and between 1.0 μm and 1.8 μm its reflectance spectrum (Sect. 3.4). We also investigate the nature of the second, faint component at larger distances (Sect. 3.5) first mentioned by Bonnefoy et al. (2017), which we firmly detect in our new SPHERE total intensity observations. In Sect. 4, we report our study of the ALMA data, including the morphological analysis of the disk (Sect. 4.1), as well as the modeling of the SED (Sect. 4.2) and its dust mass (Sect. 4.3). In Sect. 5, we discuss the dust properties contextualized with other debris disks and the global morphology of the debris disk around HD 120326, based on the STIS, SPHERE, and ALMA images. We present our conclusions in Sect. 6.

2 Observations and processing of the SPHERE, ALMA, and STIS data

Here, we present the data used and our processing: VLT/SPHERE (Sect. 2.1), ALMA (Sect. 2.2), and HST/STIS (Sect. 2.3).

2.1 SPHERE NIR observations

The system HD 120326 (HIP 67497) was observed both in polarimetry and total intensity with the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE; Beuzit et al. 2019) using its dual imager IRDIS (Dohlen et al. 2008) and its Integral Field Spectrograph (IFS, Claudi et al. 2008) at the Very Large Telescope (VLT). The field of view of IRDIS is 11″ × 11″ and the plate scale of its detector is 12.25 mas (Maire et al. 2016); whereas for the IFS, these are 1.73″ × 1.73″ and 7.46 mas, respectively. The astrometric calibration of both IRDIS and IFS is based on past regular observations of the star crowded field 47 Tuc and includes measurements on the true north, distortion, and plate scale (Maire et al. 2021).

We summarize the log of the observations in Table A.1. The seeing and atmospheric coherence time (τ₀) were estimated by the Differential Image Motion Monitor (DIMM, Sarazin & Roddier 1990) and Multi-Aperture Scintillation Sensor (MASS, Kornilov et al. 2007) turbulence monitor at the Paranal Observatory. The Strehl ratio is computed by the real-time computer named Standard Platform for Adaptive optics Real Time Applications (SPARTA, Fedrigo et al. 2006) of the SPHERE extreme Adaptive Optics system (SAXO, Petit et al. 2014). In Sect. 2.1.1, we present the polarized intensity observations and our data processing. In Sect. 2.1.2, we report the total intensity observations and our data processing. We compare the polarized and total intensity observations in Sect. 2.1.3.

2.1.1 Polarimetric observations

As dust grains scatter light, they polarize it – and this process will depend on their intrinsic properties. Thus, observing circumstellar disks with polarimetry enables the study of dust grain properties. Despite the stellar light being very bright, this technique is ideally suited to reveal the much dimmer signal from the disk. Indeed, since stellar light is mostly unpolarized, polarimetric observations facilitate the subtraction of the stellar light via polarization differential imaging (PDI).

Two polarimetric datasets are available on HD 120326: NIR data using SPHERE/IRDIS and optical data using SPHERE/ZIMPOL. Below, we first describe the IRDIS data, where the disk signal has clearly been detected, and then the ZIMPOL data, where no disk signal has been detected.

By using the NIR instrument SPHERE/IRDIS in its dual-beam polarimetric imaging (DPI; De Boer et al. 2020) mode, the host star of HD 120326 was observed on the night of the June 1, 2018 (PI: Boccaletti) under well-suited observation conditions, with a seeing of 0.45" and a coherence time of the atmosphere of 4.2 ms. The observation was carried out in the field-stabilized mode in the broad band filter H (λ = 1.625 μm, Δλ = 0.29 μm).

The IRDIS-DPI mode allows us to image simultaneously the linearly polarized light in two orthogonal directions, by using polarizers with orthogonal transmission axes. A polarimetric observation with SPHERE/IRDIS is made of polarimetric cycles consisting of images acquired at four half-wave plate (HWP) angles to measure the Stokes parameters Q⁺, Q⁻, U⁺, and U⁻. These HWP angles are 0°, 45°, 22.5°, 67.5°, respectively, with the Q⁺ vector aligned with a preferred orientation, usually the meridian on-sky (see Fig. 1 from De Boer et al. 2020). This allows us to retrieve the two linear polarization contributions, Q and U expressed as $$Q=\frac{1}{2}(Q^{+}-Q^{-})$$(1) $$U=\frac{1}{2}(U^{+}-U^{-}).$$(2)

Expressing Q and U as such enables the removal of the instrumental polarization effects caused by the reflections in the instrument downstream of the HWP. From them, one can reconstruct the linearly polarized intensity via the equation $pI=\sqrt{Q^2+U^2}$, and also the degree of linear polarization following DoLP = pI/I, where I is the total intensity (De Boer et al. 2020). Tackling the instrumental polarization generated by the reflections upstream of the HWP require additional corrections, which are described in Canovas et al. (2011). To process the polarimetric observations, we used the public pipeline IRDAP² (van Holstein et al. 2020), which corrects for instrumental polarization both upstream and downstream the HWP.

Finally in practice, we work with the azimuthal Stokes parameters Q_ϕ and U_ϕ (Schmid et al. 2006), expressed as $$Q_{\varphi }=-Q\cos (2\varphi )-U\sin (2\varphi )$$(3) $$U_{\varphi }=+Q\sin (2\varphi )-U\cos (2\varphi )$$(4)

where ϕ is the angle between the north and the point of interest (from the north toward the easterly direction). A positive Q_ϕ signal corresponds to polarization oriented in azimuthal direction (with respect to the position of the star), while a negative Q_ϕ corresponds to radial polarization. The signal corresponds to polarization angles oriented at ±45° with respect to the radial (or azimuthal) direction.

In Fig. 1, we show the final Stokes vector components pI, Q_ϕ and U_ϕresulting from the polarimetric observation acquired with 10 HWP cycles, expressed in μJy/arcsec². Stellar light has been subtracted in this image, which was evaluated to 0.06 ± 0.07% at 1σ by the IRDAP pipeline. All the signal of the disk is seen positively in Q_ϕ, i.e., this corresponds to azimuthal polarization of the light scattered by the debris disk, but it is not seen in U_ϕ. Thus, we used this polarization to estimate the noise in the polarimetric data. To convert from detector units into flux unit, we normalized the coronagraphic images with the total flux of the star computed in a circular region of radius 25 pixels summed over the two channels and corrected from respective frame time exposure and transmission of the neutral density of the corona-graphic and non-coronagraphic observations. We then multiplied this by a stellar flux density of 0.977 ± 0.044 Jy (from Johnson H filter; Bourg6s et al. 2014) and divided by the plate scale area.

To derive the signal-to-noise (S/N) map of the disk seen in polarimetry, we divided the Q_ϕ component containing the signal of the disk by the azimuthal mean of the U_ϕ representing the noise (see Fig. 1). The disk is easily detected in the IRDIS polarimetric data, with a S/N of an average value of 4.6/pixel² for the eastern part, and 3.5/pixel² for the western part of the disk.

Our processing of the SPHERE/IRDIS polarized intensity dataset yields results similar to those of Olofsson et al. (2022b), if not identical. However, in contrast to their work, we provide additional outputs, such as the calibration of the reduced image in Jy/arcsec² and the signal-to-noise (S/N) map. In addition, we also processed the SPHERE/ZIMPOL polarimetric observations acquired on the night of August 1, 2015. The ZIMPOL data were reduced by the High-Contrast Data Centre³ (HC-DC; Delorme et al. 2017), which implements the pipeline developed at ETH Zürich and described in Hunziker et al. (2021). In particular, we measured and corrected the relative beam shift of the two orthogonal polarization states, subtracted the frame transfer smearing, and corrected for the residual telescope polarization and the intrinsic polarization of the star. The Stokes Q and U images were then transformed into the Q_ϕ and U_ϕ images following Eq. (4). Unfortunately, no circumstellar structure is visible in the reduced Q_ϕ image. Therefore, we do not discuss these data further in this paper.

Fig. 1 Polarimetric images of the system HD 120326 using VLT/SPHERE at 1.6 μm (broad band H; scattered light; epoch 2018-06-01). We processed the data by using the pipeline IRDAP. From left to right, we show the Stokes vector components pI (intensity), Qϕ (azimuthal or radial polarization), and Uϕ (polarization oriented at ±45° with respect to the azimuth) calibrated in μJy/arcsec2, and the signal-to-noise (S/N) map, which is derived on the basis of the Qϕ and Uϕ images. For these images, as all the ones in this paper: north is up and east is left, while the gray circular mask corresponds to the region within the inner working angle of the coronagraph.

2.1.2 Total intensity observations

In terms of the total intensity, the star HD 120326 was observed in several modes with VLT/SPHERE: IRDIFS H23+YJ (2015-04-09, 2016-04-05, 2016-06-03), IRDIS J23 (2016-06-13), and IRDIFS BBH+YJ (2019-06-26, 2019-07-09), see Table A.1. The IRDIFS mode enables simultaneous observations of the dualband imager IRDIS in the filter doublet H2H3 (λ_H2 = 1.593 ± 0.055 μm, λ_H3 = 1.667 ± 0.056 μm) and of the low-resolution integral field spectrograph IFS in the band YJ (0.95–1.35 μm). The filter doublet J2J3 corresponds to λ_J2 = 1.190 ± 0.042 μm, λ_J3 = 1.273 ± 0.046 μm. The filters J3 and H2 are supposed to be sensitive to the main emission peaks of cold companions, whereas the filters J2 and H3 should be sensitive to molecular absorptions. The broad band H (BBH) filter is the same filter as that used in polarimetric mode (Sect. 2.1.1). All observations were carried out in the pupil-tracking mode to enable the processing of the data with the ADI technique and achieve higher contrast at sub-arcsecond separations. In pupil-tracking observations, circumstellar signal (e.g., disks or exoplanets) rotates through the sequence of observation, while the stellar halo remains quasi-static as a function of the parallactic angles. The ADI technique therefore enables the removal of most of the stellar halo and the retrieval of faint circumstellar signal. For one epoch of observation, a dataset consists of the non-coronagraphic image of the star (to derive the stellar flux and have the point-spread function -PSF- of the instrument), the coronagraphic cube acquired at two (IRDIS) or thirty-nine (IFS) channels (x, y, t, λ), and the list of the parallactic angles corresponding to each temporal frame.

Among the total intensity data, we mainly focus on the unpublished high-quality epoch of observation 2019-07-09, which has a total time exposure of 119 min, with an average seeing of 0.54″, and a variation of parallactic angles of 58°, as displayed in Table A.1. In Fig. 2, we also present results on the epochs of 2016-04-05 (published in Bonnefoy et al. 2017) and 2016-06-03 (unpublished), but in which the circumstellar signal is fainter. Regarding the other datasets (see Table A.1), we do not consider the observation of 2015-04-09 as the focal plane mask was missing during the observation sequence. We processed the observation 2019-06-26, but we do not present the associated results in this paper as the data are of significantly lower quality than those of 2019-07-09, acquired in the same band, and in which the disk is hardly detectable. We processed the epoch 2016-06-03 acquired in the J23 band, but almost no disk signal was recovered.

The data are preprocessed by the HC-DC, using the SPHERE Data Reduction and Handling pipeline (Pavlov et al. 2008). This preprocessing corrects the bad pixels, dark current, flat non-uniformity, and sky background in both IRDIS and IFS data. In addition to the calibration of the IFS, the pre-processing corrects for the wavelength and cross-talk between the spectral channels. Coronagraphic images are centered via four satellite spots used to determine the accurate position of the star hidden behind the coronagraphic mask.

We homogeneously postprocessed the observations in total intensity with a principal component analysis (PCA; Amara & Quanz 2012; Soummer et al. 2012) algorithm as implemented in the open library VIP⁴ (Gomez Gonzalez et al. 2017; Christiaens et al. 2023). For the SPHERE/IRDIS data, we tested several methods to process the data. First, we considered either one of the two channels (combined before or after the PCA processing) and we applied between 1 and 20 components to model the stellar halo before its subtraction with the PCA algorithm. The ultimate goal was to have the cleanest image to retrieve the circumstellar signal. Although the circumstellar signal is drowned in the stellar light, the number of PCA components used to model the stellar halo has to be limited, because the model will capture circumstellar signal too and remove it. This effect is named self-subtraction and it is well known (Milli et al. 2012). By evaluating the S/N of the disk and examining the residuals, we concluded that applying PCA with ten components to each of the two channels, before averaging them into a final image, is a relatively satisfying compromise; however, we do note that some other options give comparable results.

We show the reduced SPHERE/IRDIS observations for the epochs 2016-04-05, 2016-06-03 and 2019-07-09 in Fig. 2, calibrated in contrast or in S/N. To convert to contrast, we normalized the coronagraphic observation by the stellar flux estimated in a circular aperture of 25-pixel radius, as for the polarimetric observations. To obtain the S/N map, we built the noise map by processing with PCA the coronagraphic cube for the opposite parallactic angles to retain a similar temporal dependence of the stellar residuals and average any circumstellar signals, as done in Pairet et al. (2019). We can see that the disk has a higher S/N in the last observation (2019-07-09, H band), reaching (on average) about 3.9/pixel² for the eastern part of the disk and 2.5/pixel² for the western part. Our postprocessed image of the epoch 2016-04-05 is similar to that of Bonnefoy et al. (2017). However, the flux scales of their images tend to be more saturated, in an attempt to highlight the disk signal (see their Figs. 1 and 3), and they do not show any S/N maps.

From Fig. 2, we can see two disk structures from the scattered total intensity light observations. These two structures will be referred in this work as (i) "the inner belt", which was already reported by Bonnefoy et al. (2017) based on their epoch 2016-04-05, and (ii) "the second, outer structure" (see Fig. 2), which is also reported by Bonnefoy et al. (2017), but which appeared particularly faint in their observation of 2016-04-05.

For the second, outer structure, the western part is brighter than the eastern part. Using the new observation epoch of 2019-07-09, we retrieved the western part more effectively than Bonnefoy et al. (2017), who used the epoch 2016-04-05. However, they retrieved the eastern part more effectively compared to our results. For the western part, we derived a median S/N of 3.4 ± 0.6 per pixel using the data of 2019-07-09, while using the data of 2016-04-05 we obtained a median S/N of 1.5 ± 1.0 per pixel.

In addition, we processed the epochs 2019-07-09 and 2016-04-05 with an iterative version of PCA (Pairet et al. 2021; Juillard et al. 2023), using ADI only (I-PCA ADI), reference differential imaging (RDI; Lafrenière et al. 2009) only (I-PCA RDI) or combined (I-PCA ARDI; Juillard et al. 2024). Reference differential imaging has proved to be efficient to remove the stellar halo and retrieve extended circumstellar signal of dozen of ground-based observations (Xie et al. 2022; Ren et al. 2023a). Using an iterative approach for PCA helps limit the effects of self-subtraction, over-subtraction and deformations of the disk. To apply I-PCA ARDI, we used a library of reference frames obtained from observations of other stars with SPHERE/IRDIS in the same filter. On HD 120326, contrary to I-PCA RDI, the processing I-PCA ADI and I-PCA ARDI retrieve the second, outer structure on the epochs 2016-04-05 and 2019-07-09, see Fig. B.1. In particular, the eastern and western parts of the second, outer structure are more visible than in the simple PCA ADI reductions shown in Fig. 2.

Concerning IFS observations covering Y and J bands, we processed them using PCA ADI and ten modes. Since the disk is faint, we consider only the best epoch, 2019-07-09, and we binned the spectral channels per six, to increase the signal-to-noise ratio. Over the 39 spectral channels of the IFS, we kept 30 channels, removing the first and last ones, and a few that were more impacted by artifacts. In these images (Fig. B.2, from top-left to bottom-middle), the first belt is visible, but not the second, outer disk structure. We recovered the western part of the second structure only by taking the average of the 30 spectral channels and applying a spatial binning (Fig. B.2, bottom-right). The second structure was previously also hardly detected by Bonnefoy et al. (2017), see their Fig. 1, using the epoch 2016-04-05.

In the rest of this work, we will mainly use for the SPHERE/IRDIS and SPHERE/IFS total intensity observations, the epoch 2019-07-09 because its signal-to-noise ratio is higher than the other exploitable observations (2016-04-05, 2016-06-03), and the IRDIS dataset is acquired in the exact same filter than the SPHERE/IRDIS polarized intensity data (see Table A.1). We will use the PCA ADI processing to have a fast-computing processing algorithm (compared to I-PCA ADI/ARDI) when modeling the data to derive the best-fitting disk parameters (see Sects. 3.1 and 3.4).

Fig. 2 Total intensity ("I") images at about 1.6 μm reduced using the processing PCA ADI on the system HD 120326, calibrated in contrast unit (top) or in S/N (bottom). Two structures can be seen: the inner belt and a second, outer structure (see the labels in the top-right image). All images have the same orientation, with north up and east on the left.

Fig. 3 Comparison of the disk structures in HD 120326 detected in scattered light (SPHERE, 1.6 μm) and in thermal emission (ALMA, 1.3 mm). On the left, we show the polarized intensity SPHERE image with in overlay the contours of the inner belt determined for a signal level of 2 mJy/arcsec2 on this image convolved with a gaussian to obtain smooth contours. The contours are added on the total intensity SPHERE image (middle) and on the zoom-out ALMA reconstructed image (right). In the middle image, we added a white arrow to point out the southwestern part of the second structure seen in total intensity, which does not overlap with the signal seen in polarized intensity.

2.1.3 Comparison of polarized versus total intensity data

The striking result from the previous Sects. 2.1.1 and 2.1.2 is the differences between the disk structures imaged in polarized and total intensity. Only one disk belt is seen in polarized intensity, while two components are detected in total intensity. In Fig. 3, we compare the location of both structures, and we show that the first structure in total intensity matches the location of the ring imaged in polarized intensity. The second structure is located further away.

2.2 ALMA millimeter observations

ALMA observations of HD 120326 were acquired as part of project 2022.1.00968.S (PI J. Miley). In total five separate executions were made on October 19 and 21, and December 17, 18, and 27, 2022. The observations use a spectral set-up using Band 6 receivers, which include two relatively coarse resolution spectral windows with 128 channels of width 15.6 MHz placed at central frequencies of 218.0 GHz (1.375 mm) and 233.0 GHz (1.287 mm). Two further spectral windows were placed to cover any potential bright molecular lines. The first was positioned with central frequency 230.5 GHz (1.301 mm) containing 3 840 channels of width 244 kHz (1.4 μm). The final spectral window was positioned at central frequency 220.0 GHz (1.363 mm) with 1920 channels of width 976 kHz (6.0 μm). We reduced the raw data by using the pipeline script provided by ALMA observatory using CASA version 6.4 (McMullin et al. 2007).

Although there was no emission from molecular lines detected, we flagged channels corresponding to ± 15 km/s of the rest frequency of the bright molecular line transitions covered by channels within the spectral windows; specifically: ¹²CO at 230.538 GHz, ¹³CO at 220.4 GHz, and ¹⁸CO at 219.6 GHz. This ensures that any potential contamination of the continuum image by line flux is avoided. The observed data was then averaged in frequency to produce consistent channel widths of 15.625 MHz before imaging and the five executions were then concatenated into a single dataset for imaging. Images were reconstructed using tclean task in CASA using the Hogbom deconvolver (Högbom 1974) and natural weighting⁵. This produced a continuum image with a synthesized beam of size of 1.321" × 0.950" (PA = –86.729°), and achieved a noise level of 17 μJy/beam. The clear detection of circumstellar material, which is unresolved, is shown in Fig. 3 (right image). There is no evidence of any offset from the expected stellar position.

Our measurement of the flux of the debris disk is achieved by fitting a Gaussian in the uv-plane, giving an integrated flux value of 561 ± 20 μJy in the continuum at 1.3 mm. In addition, we derived an upper limit of the ¹²CO value of 4.8 Jy km/s. This corresponds to three times the rms noise level measured over an area equal to the size of the synthesized beam.

2.3 HST/STIS optical observations

The archival HST/STIS observation of HD 120326 were already published in Padgett & Stapelfeldt (2016) and Bonnefoy et al. (2017). These coronagraphic observations were acquired on 2013-06-20 with the star centered behind the WedgeA1.0 at two different rolling angles of separation 30° and in the broad band filter centered at 0.59 μm (dλ = 0.44 μm). Both papers published reduced images featuring an extended halo around HD 120326, and strong stellar artifacts at close separations, without providing a flux calibration or deriving a S/N map.

We postprocessed the dataset again with a bespoke routine to remove most of the stellar light. We used as inputs the two calibrated images acquired at ORIENT1 (28^°) and ORIENT2 (58°), performed a roll-subtraction (IM1–IM2) and de-rotated the image, resulting in north being up and east on the left. The corresponding results are shown in Fig. 4. This image shows fewer instrumental artifacts than that of Bonnefoy et al. (2017), because we avoided duplicating the cross-pattern feature during postprocessing (see their Fig. 5). As seen here, we successfully recover the disk, which is impacted by self-subtraction (see black negatives feature at the bottom right of the wedge in Fig. 4). We derived the flux map (Fig. 4) by converting the ADU in Jy/arcsec², following the procedure described on the STScI website⁶. At large separations (>2.5″), the disk has a flux density of a few μJy/arcsec².

To better distinguish the disk features from the stellar residuals, we also computed the S/N map. We defined the signal map as the reduced STIS image. Regarding the noise map, we computed radially the noise as the standard deviation in annuli. We did not consider the full image, but only the top-right quadrant, in which we masked the (small) regions defined by the location of the WedgeB1.0, one of the stellar spikes, a few background sources or suspicious, bright pixels. We applied a spatial binning of a factor four (two in height and two in width) to increase the S/N. We show the S/N map in Fig. 5. A large asymmetric disk structure extends up to a projected separation of about 700 au. Its S/N per pixel is between 5 and 8. This spiral-like feature is discussed in Sect. 5.2. Some signal is still visible relatively close to the star, at the eastern edge of the WedgeA1.0 which has a width of 1". For scale comparison, in Fig. 5, we added the image of the inner structures detected with SPHERE in total intensity (top) and in polarized intensity (middle), and with ALMA (bottom). In Sect. 5.2 we examine the global architecture of the debris disk around HD 120326.

In the next two sections, we present our constraints on the properties of the inner disk structures of HD 120326, which are detected in the NIR with VLT/SPHERE (Sect. 3) and in the millimeter with ALMA (Sect. 4).

Fig. 4 Optical coronagraphic HST/STIS observations of HD 120326 centered behind the WedgeA1.0″. We applied a spatial binning 2 × 2 pixels2 → 1 × 1 pixel2. The flux unit is μJy/arcsec2. The areas hidden by the wedges (WedgeA and WedgeB) and the diffraction spikes (the cross-pattern) are masked in gray.

3 Morphology, photometry, and polarization of the disk structures from SPHERE NIR data

In Sects. 3.1–3.4, we describe how we constrained the morphology, polarization, and reflectance spectrum of the inner dust belt of HD 120326. To this aim, we modeled the VLT/SPHERE NIR polarimetric and total intensity data. In Sect. 3.5, we examine the leftover circumstellar signal in the residuals (indicated by the white arrow in Fig. 6, bottom-right) and investigate its origin.

In particular, to constrain the disk morphology and photometry, we developed the open-access Python-based tool Modeling Disks in Scattered light⁷ (MoDiSc) to constrain the morphology and photometry of a disk seen in scattered light. The tool MoDiSc was originally inspired by the code DiskFM (Mazoyer et al. 2020). The tool MoDiSc takes as input a configuration file setting the paths to the datasets of the observations and defining the variables used in the simulations. By exploring the parameter space with a MCMC algorithm (emcee package from Foreman-Mackey et al. 2013), MoDiSc looks for the disk model matching at best the observations. To initialize the values of the MCMC simulation, a first satisfying guess can be determined by running a Nelder-Mead minimization (Nelder & Mead 1965), which is implemented in MoDiSc. The disk model can be composed of one or several belts, each of them generated with the module fm.scattered_light_disk from VIP (Gomez Gonzalez et al. 2017; Christiaens et al. 2023), the module being based on the radiative transfer code GRaTeR (Augereau et al. 1999b). One or several observations in total and/or polarized intensity can be considered simultaneously, from one or several instruments.

3.1 Disk modeling with one belt

To derive the disk morphology, we modeled both jointly and independently the SPHERE/IRDIS total and polarized intensity observations (2019-07-09 and 2018-06-01, respectively) by using the open-access code MoDiSc. We first describe our methodology and results for our joint modeling, and will then briefly also mention the results obtained for the independent modeling following a similar methodology.

In our simulations, we constrained the seven following parameters related to the belt: its reference radius r₀, its position angle PA, its inclination i, its scattering anisotropy parameter g, its outer slope α_out and two scaling factors, one to match the flux of the total intensity observations (scaling_I) and the other one to match the flux of the polarized intensity observations (scaling_pI). To limit the number of free parameters, we fixed the eccentricity, ascending node and argument of periastron of the ring to zero, and a steep inner slope to α_in = 10. We used a gaussian vertical profile (γ = 2), a linearly flared ring (β = 1) and set the scale height ξ₀ to 1.5 au at the reference radius r₀ (aspect ratio of ∼3.75%), see Eqs. (C.1)–(C.3) and Augereau et al. (1999b) for a thorough description of the model for the dust density distribution. To model the scattering phase function (SPF) of the disk in total intensity, we used the Henyey–Greenstein (HG) function (Henyey & Greenstein 1941, see Eq. (C.5)) parametrized by the scattering anisotropy parameter g. The same goes for the polarized intensity, which is modeled by this HG function multiplied by the Rayleigh scattering function, see Eq. (C.6).

Our simulations aimed to minimize the reduced χ²_r,I+pI defined as $${\chi }_{r,I+pI}^2={\chi }_{r,I}^2+{\chi }_{r,pI}^2,$$(5) $$=\frac{{\chi }_I^2}{N_{\text{resol},I}-N_{\text{param},I}-1}+\frac{{\chi }_{pI}^2}{N_{\text{resol},pI}-N_{\text{param},pI}-1},$$(6)

where N_resol,I and N_resol,I are the number of pixels considered within the regions used to model the inner belt in the total and polarized intensity images (see Fig. C.1), and N_param,I and N_param,pI are the number of free parameters used to model the total and polarized intensity images, respectively. In our case, Nresol,I = Nresol,pI = 1968 and Nparam,I = Nparam,pI = 6.

The ${\chi }_I^2$ and ${\chi }_{pI}^2$ are expressed as follows $${\chi }_I^2=\sum _{\text{pixel}}{\left[\frac{\mathtt{\text{PCA(}}{\text{obs}}_I-M_I(r_0,PA,i,{\alpha }_{\text{out}},{\text{scaling}}_I)\mathtt{\text{)}}}{{\sigma }_I}\right]}^2,$$(7) $${\chi }_{pI}^2=\sum _{\text{pixel}}{\left[\frac{Q_{\varphi }-M_{pI}(r_0,PA,i,{\alpha }_{\text{out}},{\text{scaling}}_{pI})}{{\sigma }_{pI}}\right]}^2,$$(8)

where obs_I is the pre-processed total intensity science coronagraphic cube before applying PCA ADI, Q_ϕ is the IRDAP-processed polarized intensity image, M_I and M_pI are the convolved disk models in total (I) and polarized (pI) intensity, and σ_I and σ_pI are the noise maps of the observation in total and polarized intensity, respectively. Both the total and polarized intensity (I and pI, respectively) terms have similar values, they therefore give an equal contribution in this merit function.

To explore the parameter space, we first ran simulations using a Nelder-Mead minimization (Nelder & Mead 1965), to converge quickly to a minimum, before running a Markov Chain Monte Carlo (MCMC) algorithm with emcee (Foreman-Mackey et al. 2013) to explore the space parameter more extensively. The walkers were initialized by drawing random values in the interval set by the best parameters found with Nelder-Mead plus or minus 10%. The random draw was linearly uniform for the reference radius, the position angle, the scattering anisotropy parameter, the outer slope, the cosine of the inclination and the logarithm in base ten of the scaling flux factor. In total, we used 100 walkers, with 3000 iterations for each parameter. For each iteration, we generated a disk model with the radiative transfer code GRaTeR (Augereau et al. 1999b), available in the open-access library VIP (Gomez Gonzalez et al. 2017; Christiaens et al. 2023).

In the case of the total intensity data, we convolved directly the disk model M_I with the SPHERE/IRDIS non-coronagraphic observation of the star, corresponding to the PSF of the instrument. Then, we rotated the convolved synthetic disk model for each parallactic angle of the sequence of observation, and subtracted it from the science coronagraphic cube. We processed this new cube by applying PCA using ten modes. We divided the cube by the noise map σ_I (computed in Sect. 2.1.2), producing a normalized residual map. If the disk model matches the observations, no circumstellar signal should remain in this normalized residuals map, that should only represent noise.

On the other hand, in the case of polarimetric observations, we matched the disk model with the Q_ϕ image obtained by processing with the IRDAP pipeline (see Sect. 2.1.1). To convolve robustly the disk model, we constructed the synthetic Stokes Q and U images from the disk model generated with VIP. We convolved those Q and U images with the PSF of the star and reconstructed the Q_ϕ image. Applying robustly the convolution on the reconstructed Stokes Q and U images is important for disks located at small separations, to account for positive and negative polarization signals that cancel each other when observed with a limited angular resolution (e.g., Engler et al. 2018; Heikamp & Keller 2019). We obtained the residual map by subtracting the synthetic image M_pI from the processed image Q_ϕ, and normalized it by the noise map σ_pI (computed in Sect. 2.1.1).

By considering only pixels within the region of interest (where the disk signal may be located), we summed the square value of these pixels for the total and polarized intensity, to obtain ${\chi }_I^2$ and ${\chi }_{pI}^2$, respectively, and then ${\chi }_{r,I+pI}^2$ following Eq. (6). By minimizing ${\chi }_{r,I+pI}^2$, the MCMC exploration maximizes the likelihood directly related to it.

We list the best parameters derived for the dust belt based on our disk modeling in Table 1. The total surface brightness of the disk is determined to be (1.2 ± 0.2) × 10⁻⁴ in polarized intensity and (2.0 ± 0.2) × 10⁻³ in total intensity, corresponding to 0.12 ± 0.2 μJy and 2.0 ± 0.2 μJy, respectively, by using a stellar flux of 0.975 mJy at 1.6 μm (Ofek 2008). Figure 6 shows the synthetic best disk models and residuals for the observations in polarized and total intensity at 1.6 μm.

The posteriors on the disk parameters are reported in Fig. D.1, for which we removed the first 2000 iterations for each walker, considered as the burn-in phase. We can see that the reference radius is correlated with the inclination of the belt and the outer slope. The reference radius r₀ of ${39.9}_{-0.8}^{+0.4} \text{au}$ au corresponds to a semi-major axis at maximum dust surface density in the mid plane of ${43.7}_{-1.0}^{+0.6}$ au with a full-width at half maximum of 24.6 au and a shallow outer edge, thus decreasing slowly in density for higher radii.

As the inner and outer slopes used to model the disk in our work (α_in,fixed = 10 and ${\alpha }_{\text{out,fitted}}=-{2.17}_{-0.8}^{+1.0}$) differ from those used in Bonnefoy et al. (2017) and Olofsson et al. (2022b), respectively α_in = 10 and α_out = –5 and α_in = 12.0 ± 3.9 and α_out = –19 ± 0.1, it is more relevant to compare the maximum dust surface density between our studies than the reference radius. The disk parameters of Bonnefoy et al. (2017) and Olofsson et al. (2022b) correspond to a maximum dust surface density peaking at 64.0 ± 3.3 and 37.6 ± 4.0 au, respectively, using Eq. (1) from Augereau et al. (1999b). Therefore, our value (43.7 au) is much closer to that of Olofsson et al. (2022b). Potential factors that could explain the difference between the values are, first, the type of data used (polarized or total intensity), which may probe different dust grains (e.g., depending on their scattering angles, see Fig. 7). Second, for a given data type, the epoch of observation considered plays a role, as it may be more or less sensitive to the inner belt (Sect. 2.1.2). Third, the inclination and the distance of the belt to the star are correlated (see Fig. D.1). For instance, decreasing the inclination by 1° (for inclination values around 78^°) results in a reduction of about 2 au in both the reference radius and the maximum dust surface density radius. The inclinations found by Bonnefoy et al. (2017) and Olofsson et al. (2022a) are 80° ± 1° and 76.7° ± 0.6°, respectively. This may explain part of the differences, as the smallest inclination and smallest maximum dust surface density radius first correspond to the values from Olofsson et al. (2022b), then to those from our work, and, finally, to those from Bonnefoy et al. (2017).

For completeness, we also ran additional simulations to model independently the debris disk seen either in the SPHERE/IRDlS total intensity or polarized intensity. We display the best-fitting parameters obtained in each case in Tables D.1 and D.2, respectively. The results based on total intensity give higher values for the inclination, the reference radius and the anisotropic parameter compared to the results based on polarized intensity only: i_I = 79.4° ± 0.3° and $i_{pI}={77.1}_{-0.5}^{+0.6}\circ ,r_{0,I}={43.2}_{-0.9}^{+1.0} au$ and $r_{0,pI}={38.5}_{-0.7}^{+1.1} au$ (corresponding to maximum dust surface density at ${46.7}_{-0.8}^{+1.6}$ and ${42.3}_{-0.8}^{+1.7} au$, respectively), _gI = 0.81 ± 0.01 and $g_{pI}={0.74}_{-0.02}^{+0.03}$. Results are consistent with each other at 2σ. The uncertainties provided on the disk parameters are statistical uncertainties derived from the MCMC exploration process (see Fig. D.1). They are likely underestimated, which is common in morphological simulations of disks using MCMC exploration (see Mazoyer et al. 2020). We can see that the best parameters derived by jointly fitting the total and polarized intensity are roughly an average of the best parameters derived by independently fitting these observations (Tables 1, D.1 and D.2). Compared to the inclination i we determined from jointly fitting the total and polarized intensity data, the inclinations i_I and i_pI are closer to those found by Bonnefoy et al. (2017) and Olofsson et al. (2022b). To some extent the same applies to the maximum dust surface density radius. Nonetheless, the maximum dust surface density radius determined by our modeling of the total intensity data alone still differs significantly from that of determined by Bonnefoy et al. (2017). This could be related to the fact it is difficult to detect in total intensity the ansea of the belt around HD 120326, thus to constrain its semi-major axis and so its maximum dust surface density radius (see Fig. 2).

To summarize, our best model fitting jointly the observations acquired in total and polarized intensity seems to reproduce well the observed disk signals, excepting the signal in the residuals indicated by the white arrow (Fig. 6, bottom-right). This signal corresponds to the second, outer structure already mentioned in Sect. 2.1.2 (see also Fig. 2), and we will investigate its nature in Sect. 3.5. Since most of the signal is reproduced by modeling with one dust belt jointly the polarized and total intensity observations, the difference between both images can be mostly assigned to processing effects (Sect. 2.1) and the scattering phase function (Sect. 3.2).

Fig. 5 Overview of the young debris disk around HD 120326. The black-to-white images are the S/N map of the optical HST/STIS coronagraphic data, with a spatial binning 2 × 2 pixels2 → 1 × 1 pixel2. For reference, in the center of the images we display the inner structures detected with SPHERE in total intensity (top), in polarimetry (middle) and with ALMA (bottom). The areas hidden by the coronagraphs and the HST diffraction spikes are masked in dark gray.

Fig. 6 Results of our MCMC simulations, using the Pythonized version of the GRaTeR tool to generate the belt model. We show the best models in linear or logarithmic scales (top and middle), and the associated residuals normalized by the noise (“σ”, bottom). Scattered light images for the polarized intensity images are on the left and those for total intensity images on the right. Some signal is still visible in the total intensity residuals, and is indicated by the white arrow. Unlike the disk model in total intensity, the disk model in polarized intensity accounts for Rayleigh scattering, hence visual brightness differences between these disk models depending on the scattering angle. The best-fitting parameters are given in Table 1. As a reminder, since the disk inclination found is 78°, the scattering angle is 12° at the closest location of the belt, 90° in the anses, and 168° at its furthest location (see also Fig. 5 from Perrin et al. 2015, to have an insight of the scattering angles around a dust belt).

3.2 Scattering phase function and degree of polarization

In Fig. 7, we compare the scattering phase function of the joint and independent modeling of the polarized and total intensity observations derived from the disk modeling described above. Using the ad hoc formalism of Henyey-Greenstein (see Eqs. (C.5) and (C.6)), the scattering phase function is parameterized by one scattering anisotropy parameter g, which is the same for the polarized and total intensity data in the case of the joint modeling.

Both polarized SPF mostly agree with each other, with some differences at low scattering angles (10-30°) that correspond to a region of the disk close to or below the inner working angle of the coronagraph. In Fig. 7, we also compare our polarized SPF with the one derived by Olofsson et al. (2022b). They are consistent with each other at 2o% and mostly at 1cr, with differences for scattering angles between 105° and 125°, corresponding to backward scattering from the disk. We hide in the Fig. 7 (and Fig. 8) the parts of the SPFs for which we do not have access to the scattering angles, that is to say, below 12° (i.e., 90°- inclination) and above 178° (i.e., 90° + inclination) due to geometrical effects, and also above 130^°, because no disk signal is detected in either polarized or total intensity data.

As for total intensity, the scattering anisotropy parameter g increases when modeling independently the total intensity data, compared to modeling it jointly with polarimetric observations (see Tables 1 and D.1). This results in a disk which seems even more strongly forward scattering, in other words, an SPF that is steeper at low scattering angles (Fig. 7). Nonetheless, the height of the peak must be taken with caution, first because it cannot be constrained below 12^° (geometrical effects). Second, the inner belt of HD 120326 at low scattering angles (10-30°) corresponds to a projected separation below or at the edge of the inner working angle of the coronagraph. Moreover, at small separations, the self-subtraction effects caused by ADI removes most of the light (Milli et al. 2012). Estimating the SPF of the inner belt of HD 120326 is therefore particularly challenging at low scattering angles.

Figure 8 shows the comparison of the SPF of HD 120326 in total intensity at 1.6 μm with that of other extrasolar debris disks and dust populations in the Solar System. The dust populations in the Solar System include the zodiacal dust (Leinert et al. 1976), the ring D68 of Saturn (Hedman & Stark 2015), the comet 67/Churyumov-Gerasimenko (using the dataset MTP020 acquired on 2015-08-28; Bertini et al. 2017), and microasteroids of diameter 10 μm, 30 μm, and 100 μm (Min et al. 2010). The SPFs of these microasteroids, corresponding to dust grains covered by small regolith particles, were derived theoretically using Fraunhofer diffraction model and Hapke reflectance theory (Hapke 1981) by Min et al. (2010). All SPFs are normalized at 60^° for easier comparison.

Based on Fig. 8, one can see that the SPF of HD 120326 is somehow similar to the SPFs of HD 35841 (Esposito et al. 2018), HD 61005 (Olofsson et al. 2016) and HD 114082 (Engler et al. 2023), with a significant forward scattering peak and a continuous decline as the scattering angle increases. Nonetheless, we note that the derived SPF of HD 120326 is limited by the one-component HG function we used in our modeling, which prevented us from obtaining both forward and backward scattering contributions for the inner belt of HD 120326. This can be observed, for instance, in the cases of HR 4796 (Milli et al. 2017) and HD 11721 (Engler et al. 2020). Even though we did not observe any backward scattering for HD 120326 in the images, it may be present at a lower level that our observations might not be sensitive to. Regarding the dust populations in the Solar System, the SPF of HD 120326 is very similar to that of the zodiacal dust, as already reported for other debris disks in Hughes et al. (2018).

Fig. 7 Scattering phase function for our best model of the polarized and total intensity observation, fit independently or jointly. For comparison, we added the model retrieved by Olofsson et al. (2022b) for the polarized intensity data. We hide the part below 12° and above 130°, because we do not detect the disk for such scattering angles, neither in polarized or total intensity light (see text).

3.3 Maximum degree of linear polarization

We here derived the maximum polarization fraction of the inner belt, defined as the maximum of the linear degree of polarization over the scattering angles. We divided the scattering phase function of the polarized intensity data by that of the total intensity data for our best synthetic models, using the scaling flux found in our modeling to express them in the same unit. This led us to obtain the linear degree of polarization for our best disk model.

The maximum polarization fraction occurs at a scattering angle of 90°, because of the Rayleigh scattering assumed in our formalism to model the polarized SPF (Eq. (C.6)), and the fixed scattering anisotropy parameter g for both total and polarized intensity models. This sets the bell shape of the linear degree of polarization of the inner belt.

We found that the inner belt of HD 120326 has a maximum polarization fraction of 51 ± 6%. The 1cr uncertainty accounts for the error of the flux in total and polarized intensity (Table 1). This maximum polarization fraction is similar to HR4796 A (50% ± 3% at 1.9–2.19 tmi; Perrin et al. 2015), and higher than other debris disks (10-20% in the NIR, e.g., Beta Pictoris, HD 15115, HD 32297, HD 114082, see Tamura et al. 2006; Asensio-Torres et al. 2016; Bhowmik et al. 2019; Engler et al. 2023). A high maximum polarization fraction can for instance be associated with a low albedo of dust particles (Umov 1905), or also to small monomers in dust grain aggregates (Tobon Valencia et al. 2022).

3.4 Reflectance of the inner belt

We derived the reflectance of the inner ring based on the spectroscopic SPHERE/IFS (YJ bands) and photometric SPHERE/IRDIS (broad band H) data acquired simultaneous during the best epoch of observation (2019-07-09) in total intensity, and the additional photometric points acquired in 2016-04-05 in the narrow bands H2 and H3. The reflectance represents the surface brightness of the disk (SB_disk) divided by the total stellar flux (F_★).

In practice, we ran new MCMC simulations for the five IFS images (Fig. B.2) to estimate the surface brightness of the inner belt by using the open-access code MoDiSc. In the MCMC simulations, one free parameter (the flux scaling factor) was used, and the five others (r₀, PA, i, g, α_out) were fixed to the best fitting values (Table 1) found previously (Sect. 3.1) based on our modeling of the joint total and polarized intensity observations in the broad band H. We determined the surface brightness of the inner belt as the total flux of the best synthetic belt model. Then, we divided the surface brightness of the belt by the stellar flux estimated within a large circular region of diameter 12 FWHM, representing between 20 and 27 pixels depending on the wavelength and plate scale of the instrument (IFS or IRDIS).

We show the reflectance spectrum of the inner ring between 1.0 and 1.8 μm in Fig. 9. Based on the SPHERE data, we conclude that the inner belt seems to have a red color at 1σ between 1.0 and 1.3 μm, and a gray color between 1.3 and 1.8 μm. Although there are large uncertainties on the 1.6 μm measurement, there might be a tentative evidence of a break in the reflectance spectrum between 1.3 and 1.5 μm. Such a break could be caused by an absorption band around 1.5 tmi, possibly due to H₂0 in water ice (Ciarniello et al. 2021) or in minerals (e.g., phyllosilicates; Hu et al. 2012; Lane et al. 2024). The red spectral slope from 1.3 to 1.5 μm would be due to additional compound(s) mixed with H₂0 ice or H₂0-bearing minerals, possibly opaque minerals or organic matter (e.g., mixture of water ice and kerite in Ciarniello et al. 2021).

Regarding the uncertainties on SB_disk/F_*, we derived them by including the errors on the estimations of the surface brightness of the disk and the stellar flux. Concerning the surface brightness, we estimated the error via the posteriors of the MCMC simulations. The lower and upper errors are defined to cover an interval of 68% of the posteriors, centered on their median value (e.g., Fig. D.1). As for the stellar flux, we considered a conservative 5% variation during the sequence of observation, based on empirical experience.

Fig. 8 Scattering phase function of the inner belt of HD 120326 seen in total intensity compared to other extrasolar debris disks (left) and different dust populations in the Solar System (right).

Fig. 9 Reflectance spectrum of the inner belt of HD 120326. For each point, we derive the surface brightness of the disk (SBdisk) and the stellar flux (F⋆) over a spectral range shown by the horizontal line around each point, and also by the colored areas (IFS) or the filter transmission (IRDIS) at the bottom.

3.5 Nature of the second component of the disk (seen in the SPHERE total intensity data)

By examining the residual image in Fig. 6 (bottom-right), one can see that there is still some signal in the southwestern region of the disk, indicated by the white arrow. Such a signal is located beyond the inner belt seen in polarimetry (Fig. 1), and is unlikely to be caused by an effect of the ADI postprocessing on the inner belt. Indeed, the residual image is obtained by processing a inner-belt-free cube, since we subtracted the disk model to the cube of observation, before processing it with PCA ADI (Sect. 3.1). In addition, this second structure is imaged at different epochs and with different algorithms (Figs. 2 and B.1).

Therefore, we firmly confirmed this second disk structure, whose exact nature remains to be determined. Previously, Bonnefoy et al. (2017) considered whether it could be a halo of dust grains, or, a second belt. We revisit both options below, by using in particular the fact that the second structure is not detected in the polarimetric data (Olofsson et al. 2022b, and this work).

3.5.1 Hypothesis: Second belt

By assuming that the second structure is an additional belt, Bonnefoy et al. (2017) found that its reference radius would be 137 + 9 au⁸ and its flux ratio relative to the inner one would be 0.08 + 0.03 at the scattering angle for which the flux is maximum in total intensity. By using their disk models, we checked that this value of flux ratio is the same at scattering angles of ∼30°, namely, the maximum of the polarization intensity peak based on Fig. 7, and ∼70°, which still has a high level of polarization intensity and at which we do not detect the second structure in Fig. 1.

Thus, by assuming that the dust particles of the second belt have the same polarimetric properties as those of the inner disk, we expect the flux of the second belt to be 8 ± 3% of the flux of the inner belt, that is, 12.5 times smaller. However, such a flux is too faint to be recovered in our polarimetric data, because we only detect the inner belt with a S/N of 3-7 per pixel (see Fig. 1), so a flux 12.5 times smaller results in a S/N below one per pixel. Therefore, the second structure seen in total intensity could be a second ring, that remains undetected in our SPHERE polarimetric (2018-06-01) observations due to a lack of sensitivity.

3.5.2 Hypothesis: apocenter pile-up of small particles originating from the inner belt

Alternatively, the second structure could be a halo of small particles related to the inner belt. Indeed, these small particles could originate from the collisions of larger bodies within the inner belt ("birth ring", Strubbe & Chiang 2006, see also Lecavelier Des Etangs et al. 1996; Augereau et al. 2001) similarly to the debris disk around HD 129590 (Olofsson et al. 2023). This halo would correspond to the eccentric orbits of small particles, whose size sets their eccentricity. Instead of a smooth continuous halo, it could have a visual gap between the birth belt and the second disk structure. This could be explained by the apocenter pile-up model (Olofsson et al. 2023), that uses the fact that the small, eccentric particles spend most of the time at their apocenter. By considering many small particles in a narrow range of sizes, their eccentricities would be similar, and so too their apocenter distances (but with arguments of pericenter uniformly distributed between [0, 2π]). This results in a visual arc structure located at their apocenter distance for scattering angles giving a bright enough SPF. This arc structure could be enhanced by scattering efficiencies that depend on the size (and composition, porosity) of dust particles (Olofsson et al. 2023) and the wavelength of observation. We stress that the apocenter pile-up model introduces a visual pile-up (not a dynamical one), and that it is different to the apocenter glow, which requires an asymmetric disk (Pan et al. 2016, see also Wyatt et al. 1999).

To investigate this, we used the open-access betadisk tool⁹ (Olofsson et al. 2022a). This code can compute scattered light images for different grain sizes, parametrized by the parameter β This parameter is the ratio between the stellar radiation pressure and gravitational forces. It is inversely proportional to the particle size s in the case of stars bright enough, β < 0.5 and s ≳ 0.5 μm (e.g., Burns et al. 1979, Figs. 1 and 2 from Artymowicz 1988, and Fig. 7 from Lohne et al. 2012).

By assuming the disk parameters determined in Sect. 3.1 (disk inclination, position angle, reference radius, and the scattering anisotropy parameter), we modeled the inner belt with betadisk. This synthetic disk is a linear combination of 14 scattered light images computed for different grain sizes, corresponding to β from 0.01 to 0.4. We found out that a combination of large particles, mainly in the parent bodies corresponding to the inner belt, and eccentric, small particles could better model the disk structures than the previous model generated with the radiative transfer GRaTeR code, for which we used its fast-computing version. This version does not assume any size distribution, but an overall scattering efficiency for the dust without any underlying assumption on the grains. We can indeed see that the residual map corresponding to the betadisk model (Fig. 10) is to some extent cleaner than the one corresponding to the GRaTeR model (Fig. 6). The apocenter pile-up model may therefore explain some signal of the second structure.

Yet, this second structure is asymmetric, with a higher flux and larger extension in the west than in the east. In addition, this asymmetry is reversed for the inner belt, which is brighter and more extended in the east than in the west, as also seen in Fig. 2 here and Fig. 2 from Bonnefoy et al. (2017). This aspect still remains to be explained.

Fig. 10 Scattered light images of one dust belt generated using the code betadisk, with the associated residual map normalized by the noise ("σ"), in polarized and total intensity (left and right, respectively). The second structure seen with SPHERE (indicated by the white arrow) is somehow more removed, compared to the GRaTeR modeling (Fig. 6, bottom-right).

4 Morphology and photometry of the dust belt seen with ALMA millimeter observations

In Sect. 4.1, we constrain the morphology of the disk detected with ALMA (Fig. 3). Then, in Sect. 4.2, we describe our modeling of the SED, taking into account the new ALMA measurement. Last, in Sect. 4.3, we explain how we derived the dust mass ofHD120326.

4.1 Morphology of the disk

The ALMA observation does not resolve the disk (or only very marginally, see Fig. 3), so we used parametric modeling of the visibilities to obtain constraints on the disk parameters. We used an optically thin 3D model that takes line of sight effects into account, which was previously used to model debris disks (e.g., Kennedy 2020). To obtain constraints on parameters, we modeled the data using galario (Tazzari et al. 2018) to Fourier transform the model images and compare them to the observed visibilities, and emcee to perform Markov-Chain Monte-Carlo sampling and extract posterior distributions for parameters. We assumed a Gaussian torus model, for which the model parameters are x/y offsets from the observation phase center, the disk inclination and position angle, the disk flux, the ring radius and width. The scale height (H/r) was restricted to be small (< 0.05). We did not include the stellar flux as it is ∼ 100 times fainter than the disk at 1 mm. We ran the MCMC until all parameters converged, and then continued sampling to build up the posterior distributions.

Posterior distributions from the modeling are shown in Fig. D.2. The inclination (70 ± 5°) and position angle (−91.7 ± 4.5°) are consistent at 2σ with the scattered light results. While the resolution clearly limits the conclusions that can be drawn, the disk detected with ALMA peaks at radii of 35 ± 25 au at 1σ based on our modeling. At 1σ, the outer edge of the disk in the millimeter is 75 au at most, thus closer to the star than the second structure seen with SPHERE, which peaks at radii larger than 100 au (Sect. 3.5).

In particular, the millimeter-wave emission peaks at a location consistent with that of the inner belt seen in the NIR, for which scattered light models give disk radii between 30 and 65 au at lσ, depending on the modeling and the SPHERE data used (see Sect. 3.1). By assuming millimeter dust traces the presence of planetesimals, we would expect that the inner belt resolved with SPHERE is a planetesimal belt. To be more confident would require higher resolution ALMA observations.

4.2 Modeling of the spectral energy distribution

To construct a spectral energy distribution, we collected photometry from various sources, for instance: UBV, Stromgren, HlPPARCOS, Gaia, 2MASS, WISE, Spitzer, and Herschel. We also used the new ALMA photometric point derived in this work, and a Spitzer IRS spectrum. These data were fit with stellar and blackbody disk models as described previously by Yelverton et al. (2019). The best fitting model is shown in Fig. 11. The disk model used is a single component modified blackbody, as we found that this was sufficient to explain the data. We determined that the star has an effective temperature of 6940 ± 100 K, and a luminosity of 4.5 ± 0.1 L_⊙. We derived for the disk a temperature of 110 ± 1K and a fractional luminosity of 1.8 × 10⁻³. The ALMA photometry is somewhat lower than expected for a pure blackbody, but this is expected for debris disks where much of the emission comes from bodies that are significantly smaller than 1 mm.

Assuming blackbody absorption and emission, the dust temperature implies a dust population at 13 au from the star, which is clearly closer than it is observed. However, the small dust sizes also mean that their temperatures are higher than a blackbody at the same stellocentric distance, meaning that the 13 au distance is a lower limit. Empirically, Pawellek & Krivov (2015) found that the difference in disk radii at 4.5 L_⊙ is about a factor of four, which would make the SED estimate consistent with the semi-major axis derived for the SPHERE or ALMA data.

Fig. 11 Spectral energy distribution of HD 120326 based on archival data (see main text) and the new ALMA measurement at 1.3 mm.

4.3 Dust mass

Based on the results from our SED modeling, we derived the dust mass, which consists of the mass of dust particles of size below a few millimeters. This is a significant underestimation of the total mass budget of the debris disk, as most of the mass is represented by the largest bodies (e.g., planetary embryos or kilometer-size planetesimals). The dust mass is regularly derived in the literature, expressed as $$M_{\text{dust}}=\frac{F\nu d^2}{\kappa B_{\nu }(T_{\text{dust}})},$$(9)

where F_v is the flux density of the dust belt in W.m⁻².Hz⁻¹, d is the distance of the star hosting the dust belt, k_v is the dust opacity assumed to be 1.7 cm²g⁻¹ (850 μm/λ_obs) from Beckwith et al. (1990) (also assumed in Krivov & Wyatt 2020), B_v(T_dust) is the Planck function, and T_dust is the temperature of the dust grains chosen to be consistent with the observed SED (110 K, see Sect. 4.2).

We obtained a dust mass for HD 120326 of 0.067 M_⊕. In Fig. 12, we compare this dust mass to the sample from Krivov & Wyatt (2020), which is restricted to the best-quality ALMA data from Matra et al. (2018), observed either in band 6 (λ ∼ 1.3 mm) or 7 (λ ∼ 850 μm). From Fig. 12, HD 120326 has a fractional luminosity and dust mass very similar to β Pictoris (2.1 × 10⁻³ and 0.079 M_⊕; A6V-type star), and for the two next closest matches, to HD 61005 (2.3 × 10⁻³ and 0.13 M_⊕; G8Vk star) and HD 146181 (2.2 × 10⁻³ and 0.17 M_e; F6V star), based on the values derived in Krivov & Wyatt (2020) (dust mass) and Matrà et al. (2018) (fractional luminosity).

5 Discussion

We discuss our results on the inner belt (Sect. 5.1) below. We also investigate the global morphology of the debris disk around HD 120326 (Sect. 5.2).

Fig. 12 Dust mass and fractional luminosity of the debris disk around HD 120326 (diamond marker) compared to those of the sample used in Krivov & Wyatt (2020) (circle markers). The dust mass and fractional luminosity of HD 1203 26 is in particular very similar to β Pic. The color represents the location of the dust belt relative to its host star.

5.1 Inner belt

In Sects. 3.1–3.4, we discuss the constraints on different properties of the inner belt based on the VLT/SPHERE data, which we then identified as a planetesimal belt based on the ALMA data (Sect. 4.1). To summarize the results based on the SPHERE data, the inner belt has a semi-major axis of 43.7 au, with a FWHM of 24.6 au (fractional width of 0.56). Its dust grains are shown to be particularly forward scattering in terms of the total intensity. Assuming a Rayleigh scattering for the SPF in polarized intensity data does match the polarized intensity data well. Their relatively high maximum degree of polarization (51% ± 6%) could be related to various physical properties (e.g., low albedo or small monomers if particles are structured in dust aggregates; Sect. 3.3). The red (1.0–1.3 μm) and then gray (1.5–1.8μm) slopes of the dust particles can be interpreted as either related to the size or the composition of the dust grains. In the former case, dust particles would rather have a size (s) equal or larger than the one defined by the size parameter size (2πs/λ~1) at the turnover wavelength 1.3 μm (see Fig. 9). On the other hand, in the latter case, the slope and tentative break in the reflectance spectrum could be caused by absorption related to the presence of H₂O (see Sect. 3.4). Increasing the wavelength coverage and/or the spectral resolution could help leverage the degeneracies between different dust grain properties.

We compared the color 1.1–1.6 μm of the inner belt of HD 120326 with other debris disks (Ren et al. 2023b) in Fig. 13. The inner belt of HD 120326 has a color redder than 0.36 mag compared to what is predicted from Ren et al. (2023b) for a star of 4.5 L_⊙, lying at 3σ of their confidence region. Nonetheless, as they highlighted in their paper, the trend is based on a relatively small sample (ten disks), and some systematics could occur, as their data are from a different instrument, HST/NICMOS (but with a similar spectral coverage, with the filters F110W and F160W).

We point out that it would be particularly interesting to image the inner belt of HD 120326 with HST/STIS, which has a wavelength range centered at 0.59 μm (dλ = 0.44 μm), using a mask with a smaller inner working angle (e.g., the Cbar5). Indeed, Ren et al. (2023b) derived the color 0.59–1.1 μm for twenty debris disks, and reported that debris disks with stellar luminosity of 4.5 L_⊙ have a bluer color of 0.7 mag at these spectral ranges. If HD 120326 would also have a blue color 0.59 – 1.1 μm and assuming the color is mainly related to the size of dust grains rather than their composition, this could hint at two different dust populations (e.g., Fig. 2 from Thébault & Krai 2019). Very small grains would be responsible for a blue color between 0.59 and 1.1 μm, and larger ones would be responsible for a red color between 1.1 and 1.6 μm.

Fig. 13 Dust color at scattering angles between 80 and 100° for the debris disk around HD 120326 compared to other debris disks (Ren et al. 2023b). The trend and shaded areas corresponding to the uncertainties of the trend at 1σ, 2σ and 3σ are from Ren et al. (2023b) based on their analysis (see their Sect. 3.2).

5.2 Overall picture of the debris disk around HD 120326

The complex, global morphology of the debris disk around HD 120326 is illustrated in Fig. 5. Two key questions are, first, what is the nature of the extended asymmetric structure seen with HST/STIS at 100s of au and the second one seen with SPHERE in total intensity at smaller distances, and second, to what extent the inner structures detected with SPHERE and ALMA connect with the outer one seen with STIS.

Regarding the first question, the extended, asymmetric structure could be an eccentric or misaligned ring, or a spiral. The onset of an eccentric or misaligned disk is preceded by the appearance of spirals (Wyatt 2005b, assuming massless disks), and requires timescales that are an order of magnitude greater than the timescale governing the appearance of spirals, amounting to typically several 10-100 Myr at the edges (≳100 au) of planetary systems (e.g., Pearce & Wyatt 2014). Thus, dynamical timescales are shorter with respect to the creation of a spiral, as opposed to an eccentric or misaligned disk. Therefore, considering the young age (16 Myr; Mamajek et al. 2002) of the HD 120326 system, the asymmetric structure could is more likely to be a spiral than an eccentric or misaligned disk.

While spirals are regularly observed in gas rich proto-planetary disks, their detection is much rarer in debris disks. The reason for this is still unclear. The argument usually advanced is that they are only transient (Wyatt 2005b; Farhat et al. 2023). However, theoretical studies of planet debris disk interactions indicate that these features can be more persistent depending on the importance of disk self-gravity (Hahn 2003; Ward & Hahn 1998; Jalali & Tremaine 2012; Sefilian 2022). Notably, HD 120326 is one of the very rare debris disks with an imaged spiral-like feature. The other ones are the debris disks around TWA 7 (Olofsson et al. 2018), and possibly HD 106906 (Rodet et al. 2017; Farhat et al. 2023), and the gas-rich debris disk around HD 141569 (Konishi et al. 2016, but this one could be considered as a transition disk as with an age of 5 Myr). Among all of them, HD 120326 is the only one with a spiral-like feature spanning over 100s of au, up to 1000 au if co-planar with the inner belt. By assuming it is indeed a spiral, different mechanisms of dynamical interactions could be outlined to shape it. Stellar flybys can generate spirals (e.g., HD 141569; Reche et al. 2009), but this scenario remains very unlikely for HD 120326 based on the extensive search for such a dynamical perturber by Bertini et al. (2023).

As a result, the likely remaining possible mechanisms are, first, secular dynamical interactions with a planet (Wyatt 2005b; Farhat et al. 2023) possibly combined to disk self-gravity (Hahn 2003; Ward & Hahn 1998; Sefilian 2022; Sefilian et al. 2021, 2023), or considering disk self-gravity alone (Jalali & Tremaine 2012). It could also be collisional dust avalanches following the breakup of a large planetesimal-like object (Artymowicz 1997; Grigorieva et al. 2007), giant impacts (Jones et al. 2023), interstellar medium wind (as for HD 61005; Maness et al. 2009), or it could be potentially inherited from the protoplanetary disk phase (e.g., Najita et al. 2022).

To further investigate the hypothesis of dynamical interactions with a putative planet, we ran N-body simulations using the symplectic code SWIFT-RMVS3 (Levison & Duncan 1994), assuming a massless disk and a giant planet with various orbital configurations, to assess whether an extended asymmetric structure spanning 100 of au could be generated. We successfully reproduced a spiral structure extending over 100s of au within 16 Myr (i.e., the age of the system), by considering a giant exoplanet of 1 M_Jup on a highly eccentric (e = 0.9) orbit with a periastron of 750 au (see Fig. 14). A Jupiter-mass exoplanet remains below the detection limits achievable for this system using SPHERE (see Fig. B.1 from Bonnefoy et al. 2017). This highly eccentric, wide-orbit Jupiter-mass scenario is possible under the assumption of a massless disk. However, other scenarios involving exoplanet(s) – whether accounting for disk self-gravity or not – may also be possible (Sefilian 2022). A comprehensive investigation into the physical and orbital constraints required for an exoplanet to generate the extended spiral-like feature is beyond the scope of this paper. In addition, based on our N-body simulations, this giant exoplanet alone would not explain the asymmetric structure observed at 100-150 au, which likely has a separate origin.

Alternatively, we point out that the location of a bright pointlike source overlays with the southeastern end of the spiral-like feature (see Figs. 4 and 5). This bright source was determined to be a background star (see Appendix E and in particular Fig. E.1). We investigated whether the asymmetry of this extended disk structure could be caused by a fortuitous alignment of the background star with dust grains from the halo of dust grains around HD 120326. In such cases, the particularly forward-scattering dust grains could scatter some of the light from the background star. This could result in more scattered light on the southeastern grains than on the southwestern ones, because there is no bright background star there. This would mean that the spiral-like feature seen might actually be a fraction of an outer belt. Initially, we postulated that this candidate outer belt could be illuminated by the background star. Because of its location and strong forward scattering, the fact that the background star is located on the west side could explain the west-east asymmetry. However, as further described in Appendix E, the difference in scattering angles between the eastern and western sides is not sufficient to explain such a strong asymmetry.

To determine what is the true morphology of the extended debris disk around HD 120326, and ultimately identify the mechanism causing it, new HST/STIS data are required. Indeed, large areas where disk signal is expected are masked by the WedgeA1.0 or the cross-pattern of the spiders, or are affected by self-subtraction effects from the data processing (see Sect. 2.3 and Fig. 5). The new STIS data should be acquired with an optimized strategy of observation, using a color-matched reference star, more angular diversity (e.g., with three rolling angles instead of two), and an appropriate position of the wedges relatively to the disk, to probe a larger spatial coverage (e.g., Schneider et al. 2009).

Concerning the second, enigmatic structure seen with SPHERE in total intensity at smaller radii (e.g., Sect. 2.1.2), we investigate its nature in Sect. 3.5. Analysing this faint structure requires pushing SPHERE to its limits in terms of sensitivity, so we could not draw decisive conclusions about its true nature. We confirmed that the second structure could be a second ring, which was previously hypothesized by Bonnefoy et al. (2017), and we showed that due to its faintness, it is not expected to be detected in the SPHERE polarimetric observations (see Sect. 3.5.1). Alternatively, instead of a second ring, it could also be a halo of small particles that originate from the birth ring (see Sect. 3.5.2), identified as a planetesimal belt (see Sect. 4.1). This is supported by the apocenter pile-up model (Olofsson et al. 2023).

One way to distinguish between the belt and halo hypotheses is through a color analysis of the dust grains in the second structure. This would help determine the grain size and, consequently, the nature of the structure. If the grains have a blue color, this could indicate they are small, and smaller than those in the inner belt, given the red-gray reflectance spectra of the dust grains in the inner belt (Sect. 3.4). This would support a halo nature for the second structure. Unfortunately, we do not detect well the second structure at wavelengths shorter than those of the broad band H (1.48–1.78 μm) in the SPHERE/IFS data (see Fig. B.2). In these data, we only faintly detect the western part of the disk when combining 30 spectral channels (0.98–1.32 μm) and when applying a strong spatial binning (Fig. B.2, bottom-left image). In addition, this western structure is located very close to a bright, extended signal of unclear origin, which could bias the flux measurement. As a result, only a handful of pixels capture a small fraction of the second structure. Therefore, highly sensitive HST/STIS observations using the Cbar5 mask would be highly valuable for detecting the second structure at short wavelengths and analyzing its color.

Furthermore, we stress that the second structure is asymmetric, with a higher flux and larger extension in the south-west region. This could be explained by a giant impact at the opposite direction in the birth ring, generating small, eccentric dust particles pushed away by the stellar radiation pressure (Jones et al. 2023).

As for the second question, regarding whether both the outer (imaged with STIS) and inner (detected with SPHERE and ALMA) structures connect with each other, tackling it also suffers from the lack of spatial coverage of STIS, and in particular at small separations, hidden by the too large mask WedgeA1.0. In addition, some parts of the spiders coincide with the spirallike feature, in particular in the region where it seems it could connect with the inner ring, see Fig. 5.

Concerning the SPHERE and ALMA images, in both images we do not detect a disk structure above a projected separation of 1.0–1.5", where it could overlay with the disk signal detected with STIS data. However, focusing on the SPHERE NIR total intensity data, we cannot rule out the possibility that the second structure (or halo) extends further away. If such an extension exists, it may remain undetected in the SPHERE data due to sensitivity limitations. Indeed, the starlight reflected onto dust grains decreases with their distance from the star, following a 1/r² dependence. Alternatively, in the outer regions of the halo, the dust grains are likely smaller, resulting in scattered light reduced in the NIR and increased in the optical. In any case, the second structure detected with SPHERE in total intensity is intriguing, and lies in the direction where the spiral-like feature begins.

Therefore, the existing observations of HD 120326 cannot distinguish whether the intermediate-scales structures seen with SPHERE and ALMA are connected with the large-scale structure detected with STIS. However, new STIS data using masks with a smaller inner working angle,such as the Cbar5 or WedgeA0.6, may provide crucial constraints.

Fig. 14 Very eccentric, wide-separation Jupiter-mass exoplanet that could generate a spiral-like feature spanning 100 of au, during the timescale of the age of the system. The blue curve shows part of the planetary orbit used in our N-body simulations.

6 Conclusion

In this paper, we present a panchromatic study of the young (16 Myr) debris disk around HD 120326 surrounding a F5V-type star and located in the Sco-Cen OB association. We used optical (HST/STIS), NIR (VLT/SPHERE, both in total and polarized intensity), and millimeter (ALMA) data, with some of the data published for the first time in this work. These observations were sensitive either to the scattered starlight on dust grains (STIS, SPHERE) or their thermal emission (ALMA).

The debris disk around HD 120326 shows substructures both at intermediate and large distances, imaged with SPHERE and ALMA (30-150 au) and STIS (≤1000 au). From our analysis, we show that modeling the ALMA data with one belt results in a belt morphology and location that are consistent with that of the inner belt imaged with SPHERE. By assuming that millimeter dust particles trace the location of planetesimals, we identified the inner belt as a planetesimal belt. This belt has a semi-major axis at maximum dust volume density of ∼43.7 au and an inclination of ∼78°, based on our joint modeling of the SPHERE NIR polarized and total intensity data using our open-access code MoDiSc. Constraints on the belt morphology using millimeter data are looser and, thus, higher resolution ALMA observations are needed to be more accurate. Moreover, we derived from the ALMA data an integrated flux value of 561 ± 20 μJy in the continuum at 1.3 mm and an upper limit of ¹²CO emission of 4.8 Jy·km/s. Using the new photometric ALMA measurement, we modeled the SED and determined a fractional luminosity of 1.8 × 10⁻³ and a dust mass of 0.067 M_⊕, with both values being very similar to the debris disk around β Pictoris.

Based on the SPHERE data, still regarding the inner belt, we derived its scattering phase function both in total and polarized intensity and its maximum degree of polarization (51% ± 6%) at 1.6 μm. We derived its reflectance spectrum, which seems to show a red color between 1.0 and 1.3 μm and a gray color between 1.5 and 1.8 |mi. There might also be a tentative evidence of break in the reflectance between about 1.3 and 1.5 |mi. Such a break could be caused by various absorption features, for instance, water in hydrated minerals (e.g., phyllosilicates) or opaque mixtures of water ice and kerite.

Furthermore, we also explored the nature of the second faint structure (≤150 au), seen in several SPHERE epochs in total intensity, with various processing algorithms. We confirmed it could be a second dust belt or a halo of small particles. In the context of the halo hypothesis, the asymmetry structure, which is brighter and more extended in the southeast than in the southwestern part, may be explained by a giant impact that could have occurred in the opposite direction in the inner ring (northeast) and that would generate small, eccentric dust particles pushed away by the stellar radiation pressure (Jones et al. 2023). We show that the non-detection in the polarized intensity data of the second structure is not contradictory with the halo or second belt hypotheses. A further investigation of both hypotheses could be achieved through a color analysis of the dust grains in the second structure. This would require new observations, such as STIS imaging using the Cbar5 mask. Such observations would also help determine whether the second structure is connected to the extended structure previously detected in archival STIS data acquired using the WedgeA1.0 mask.

As for the extended asymmetric-like feature (≤1000 au), we investigated its morphology on the basis of a new processing of the STIS archival data and theoretical considerations. New, optimized STIS observations are required to precisely constrain its morphology and the physical mechanisms that could generate and shape it. If the spiral feature is confirmed, HD 120326 would be one of the very rare debris disks with an imaged spiral feature so far and the only one with a spiral extending over hundreds of au, up to 1000 au. Yet, it is interesting to note that spiral features are a common outcome of planet-debris disk interactions, both in massless (e.g., Wyatt 2005a; Pearce & Wyatt 2014; Farhat et al. 2023) and massive disks (with or without planets; e.g., Ward & Hahn 1998; Hahn 2003; Jalali & Tremaine 2012; Sefilian 2022).

To date, no planet has been discovered in HD 120326. Bonnefoy et al. (2017) ruled out the presence of >2 M_Jup companions beyond projected separations of 50 au (and below 750 au). Future observations with the upgrade of SPHERE (SPHERE+; Boccaletti et al. 2020), as well as high-contrast imaging instruments on extremely large telescopes (ELT and GMT) or on the JWST and Nancy Grace Roman Space Telescope would be sensitive to less massive exoplanets in HD 120326 and/or located closer to the host star. Such observations will provide major constraints on the global architecture of HD 120326, including any potential disk-planet interactions that are sculpting the observed disk structures. Ultimately, these advancements would help to determine how unique the architecture of the young system HD 120326 really is.

Acknowledgements

C.D. thanks for fruitful discussions Elodie Choquet regarding handling HST/STIS data, Gaspard Duchene concerning scattering phase function and disk modeling, Veronica Roccatagliata concerning any potential flybys, Natalia Engler about albedo derivation for dust grains, and Anne-Marie Lagrange, regarding the presence of putative planets in HD 120326. C.D. also thanks Karl Stapelfeldt, Elisabeth Matthews, Sophia Stasevic, Xie Chen, Bin Ren, and Oliver Absil for useful discussions. C.D. is also grateful to Anthony Boccaletti, Deborah Padgett, and Sasha Hinkley, who obtained as PI some of the datasets used in this papers. Last but not least, we thank the referee for their comments which improve the quality of this article. This work is based on observations collected at the European Southern Observatory at Paranal with SPHERE under ESO programmes 095.C-0607(A), 095.C-0487(A), 097.C-0060(A), 097.C-0949(A), 097.C-0865(F), 0101.C-0128(D), 0101.C-0016(A), with HST/STIS (program: 12998). In addition, this paper makes use of the following and also observations collected with ALMA (project ID: 2022.1.00968.S) andALMA data: ADS/JAO.ALMA#2022.1.00968.1. 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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work has made use of the High Contrast Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon/CRAL, and supported by a grant from Labex OSUG@2020 (Investissements d'avenir – ANR10 LABX56). C.D. is part of Labex OSUG (ANR10 LABX56). C.D. acknowledges support from the European Research Council under the European Union's Horizon 2020 research and innovation program under grant agreement No. 832428-Origins. G.M.K. is supported by the Royal Society as a Royal Society University Research Fellow. A.A.S. is supported by the Heising-Simons Foundation through a 51 Pegasi b Fellowship. T.D.P. is supported by a UKRI/EPSRC Stephen Hawking Fellowship. F.M. has received funding from the European Research Council (ERC) under the European Union's Horizon Europe research and innovation program (grant agreement No. 101053020, project Dust2Planets). J.M. acknowledges support from FONDECYT de Postdoctorado 2024 #3240612 V.F. acknowledges funding from the National Aeronautics and Space Administration through the Exoplanet Research Program under Grants No. 80NSSC21K0394 (PI: S. Ertel) and No 80NSSC23K0288 (PI: V. Faramaz). M.B. is supported by the European Union's Horizon 2020 research and innovation programme under grant agreement no. 951815 (AtLAST).

Appendix A Log of the VLT/SPHERE observations of the system HD 120326

In Table A.1, we list all the SPHERE observations on HD 120326 to date.

Appendix B NIR total intensity SPHERE observations

In the main text, we describe our PCA ADI processing of the SPHERE/IRDIS data of the epochs acquired in total intensity in Fig. 2. For the best epoch in the dual band H23 (2016-04-05) and the epoch in the broad band H (2019-07-09), we also processed the data with an iterative PCA, using ADI, RDI or ADI+RDI (ARDI) techniques. We only show the I-PCA ADI and ARDI reductions in Fig. B.1, because our I-PCA RDI processing is not able to retrieve the disk.

Concerning SPHERE/IFS observations, we show our PCA ADI processing of the best epoch 2019-07-09 in Fig. B.2. There are five IFS images with sub-coverages corresponding to binning of six spectral channels in the Y or J bands. The sixth IFS image is a mean of these thirty channels, thus covering YJ. We did not consider the first and last channels of the IFS, nor a few that are more impacted by artifacts.

Appendix C Morphological parameters of the belt based on SPHERE data

In this appendix, we include materials regarding our morphological analysis of the inner dust belt carried out in Sect. 3.1.

First, Fig. C.1 shows the areas considered to model the inner dust belt. Second, we describe below the equations setting the geometry of the disk, show the posteriors of our simulations jointly modeling the total and polarized intensity observations in Fig. D.1. Third, in Tables D.1 and D.2, we display the best-fitting parameters obtained when we independently modeled either the total or polarized intensity observations, respectively.

Following Augereau et al. (1999b), the dust density distribution may be expressed in cylindrical coordinates as, $$\rho (r,\theta ,z)={\rho }_0 \times R(r,\theta ) \times Z(r,z),$$(C.1)

where R(ρ, θ) and Z(r, z) describe the radial and vertical distribution profile, respectively, and ρ₀ is the normalized grain density.

This can be expressed as $$\rho \left(r,\theta ,z\right)={\rho }_0\times {\left[\frac{2}{{\left(\frac{r}{\tilde{r}\left(\theta \right)}\right)}^{-2{\alpha }_{in}}+{\left(\frac{r}{\tilde{r}\left(\theta \right)}\right)}^{-2{\alpha }_{out}}}\right]}^{1/2}\times \exp \left[-{\left(\frac{z}{\xi \left(r\right)}\right)}^{\gamma }\right],$$(C.2)

where γ accounts for the shape of the vertical profile, either gaussian (γ = 2), spread out (γ >> 2) or stepped (γ << 2), ξ(r) is the scale height of the ring. It is parameterized with β the flaring of the disk a the semi-major axis, and e the eccentricity as $$\xi (r)={\xi }_0{\left(\frac{r}{a(1-e^2)}\right)}^{\beta }.$$(C.3)

If the disk is circular, $\tilde{r}(\theta )$ is equal to the constant reference radius r₀, and otherwise parameterized following Milli et al. (2017) via $$\tilde{r}(\theta )=\frac{a(1-e^2)}{1+e\cos (\theta )},$$(C.4)

where θ is the polar angle.

In NIR observations, disk structures are seen via the scattering of the stellar light on the dust particles of the disk. The scattering phase function of a debris disk in total intensity light may be modeled for the sake of convenience with the Henyey-Greenstein (HG) function (Henyey & Greenstein 1941) as follows $$\text{HG}(g,\varphi ) = \frac{1}{4\pi }\times \frac{1-g^2}{{\left(1-2g\cos (\varphi )+g^2\right)}^{3/2}},$$(C.5)

which is parameterized by the HG scattering anisotropy parameter g. If there is only backward scattering, g is equal to -1, if only forward scattering, g is equal to 1, and if the scattering is isotropic, g is equal to 0. For debris disks, g is typically greater than 0 because small dust grains are known to be forward scattering.

For polarimetric observations, we multiply the SPF by the Rayleigh scattering function as follows $$\text{HG}(g,\varphi ) \to \text{HG}(g,\varphi )\times \frac{1-(\cos (\varphi ){)}^2}{1+(\cos (\varphi ){)}^2}.$$(C.6)

Appendix D Morphological parameters of the belt based on ALMA data

Figure D.2 shows some of the posteriors of the simulations constraining the morphology of the debris disk based on the ALMA data (Sect. 4.1).

Fig. B.1 Total intensity images in the broad band H and reduced using an iterative PCA ADI or ADI+RDI algorithm.

Appendix E The background contaminant

We show here that the bright point-like source on the extended spiral-like feature (see Fig. 4) is a background star and investigate whether it could explain the asymmetry of the extended feature seen with HST/STlS. This bright point-like source was detected in the five epochs HST/STIS and SPHERE/IRDIS from 2013 to 2019. By comparing the relative motion of this bright point source for different epochs we show in Fig. E.1 that this source follows the track of a background object. This track is determined by using the parallactic and proper motions of the primary star. HD 120326 has a proper motion of (−29.0 mas/yr in right-ascension, and -21.0 mas/yr in declination (Gaia Collaboration 2021).

Fig. B.2 Images depicting (from top-left to bottom-middle) the total intensity ("I") in different sub-coverages of the bands Y and J, reduced using PCA ADI. The bottom-right image corresponds to the average of the PCA ADI processed images over the full YJ bands, and that we spatially binned (3 × 3 pixels2 → 1 × 1 pixel2) to reveal the western part of the second, outer structure. North is up, east is left for all the images.

Fig. C.1 Region of the SPHERE polarized (left) and total (right) intensity images considered to model the inner belt in Sect. 3.1.

Regarding the photometry of the background object, we derived it based on SPHERE/IRDIS data in J2, J3, H2, H3, and BBH bands (see Table E.1). We extracted its consrast using a circular aperture of radius 2.5 FWHM on the SPHERE/IRDIS non-ADI¹⁰ postprocessed images by the High-Contrast Data Center (HC-DC; Delorme et al. 2017) produced using the SpeCal pipeline (Galicher et al. 2018). We subtracted the sky contribution in the contrast measurements. We computed the sky contribution in circular apertures of 2.5 FWHM, centered at the same separation than the background object, and with a close position angle. At 1σ, the background object does not show a significant color in J2-J3 (0.28 ± 0.35 mag) and H2-H3 (0.12 ± 0.16 mag), but it does between J3-H2 (1.12 ± 0.34 mag). This could hint to a M-type object (from M0 to M9) with a J3-H2 color reddened by the presence of dust along the line of sight, or to a L-type object (from L0 to L9), as described in Appendix C in Bonnefoy et al. (2018) and references therein. The absolute magnitude of a M-type object is about between 6 and 11 mag in H2, while for a L-type object it is between 11 and 15 mag in H2 (Bonnefoy et al. 2018). By deriving the distance of the background object to match its absolute magnitude in H2 with that of a M or L-type object, we would expect the object to have a distance between 42 and 420 pc if M type, or between 6.7 and 42 pc if L type. Since the object follows the track of a background object, without showing additional proper motion (Fig. 4), we would rather expect it to be a background M-type object than a foreground L/M-type object.

Below, we investigate whether the alignment of the background object with the southeastern part of the dust seen with HST/STIS could explain why the disk is brighter in the southeastern part than in the southwestern part, which would results in a visual (and not dynamical) asymmetry. We stress that such investigation is limited because the southwestern part of the disk lacks of spatial coverage in the HST/STIS image.

Fig. D.1 Posteriors of the MCMC exploration of the parameter space for the disk seen in NIR (1.6 μm) SPHERE observations in polarized and total intensity. The vertical dashed lines represents the values of the median and at ±1σ, which are defined as the 15.9% and 84.1% percentiles to encompass 68.2% of the solutions for a given parameter. The blue vertical line represent the best-fitting parameters, which are the values minimizing the ${\chi }_r^2$.

We consider two regions of interest where there are dust grains, one located at a close projected separation of the background star and another one located further away (see Fig. E.2). By deriving the mean average flux within each of this circular aperture (of radius five pixels, i.e., 254 mas), we obtained a flux of 0.94 ± 0.03 μJy and 0.52 ± 0.03 μJy, respectively.

The first goal is to express the scattering angle (ɸ) as a function of the distance between the system HD 120326 and the observer (d_HD120326), the distance between the background object and the observer (d_background), the angular separation between the background object and the dust grain of interest (θ, which also represents the angle between the line of sight observer-background object and the scattered wave on the dust grain of interest), see Fig. E.3. We assume that the distance between the dust grains around HD 120326 seen with HST/STIS (Fig. E.2) and the observer is approximately equal to the distance between the host star of HD 120326 and the observer, denoted as d_{HD 120326}.

Fig. D.2 Posteriors of the MCMC exploration of the parameter space for the disk seen in continuum millimeter (1.3 mm) interferometric ALMA observations. The parameters are: position angle (PA), inclination (i), radius (r), and width (dr).

Fig. E.1 Difference of the relative astrometry of the bright point-source lying on the spiral-like feature (see Fig. 4) between the HST/STIS epoch (2013-06-20, used as reference) and the SPHERE/IRDIS epochs (2016 to 2019, see Table A.1). The blue solid curve shows the expected track of a background object, based on the parallactic and proper motions of the primary star HD 120326. The point source is clearly identified as a background object.

Fig. E.2 Location of the two regions of interest of dust grains in the extended dust emission seen in the HST/STIS reduced image. The dotted blue circle corresponds to the region located at a projected distance close to the background star, while the solid pink circle corresponds to the region located a projected distance far to the background star. Unlike Fig. 4, no spatial binning is applied on this image.

By considering the triangle whose corners are the background star, the dust grain of interest and the observer, and the fact that the sum of angles in a triangle is equal to 180^°, one can express $$\theta +\psi +(180-\varphi )=180,$$(E.1) $$\psi =\varphi -\theta ,$$(E.2)

Fig. E.3 Sketch representing the parameters (angles and distances) used to express the scattering angle ɸ of the background stellar light scattered by a dust grain in the system HD 120326.

where ψ is the angle between the line of sight observer-background star and the incident wave from the background star on the dust grain of interest. As ψ is very small, tan(ψ) ∼ ψ, hence one can express $$\psi \sim \frac{d_{\text{HD 120326}}\times \theta }{d_{\text{background}}-d_{\text{HD 120326}}},$$(E.3)

and so by injecting Eq. E.2 in Eq. E.3, one gets $$\varphi \sim \theta +\frac{d_{\text{HD 120326}}\times \theta }{d_{\text{background}}-d_{\text{HD 120326}}}.$$(E.4)

This results in the following expression for the scattering angle $$\varphi \sim \frac{\theta }{1-d_{\text{HD 120326}}/d_{\text{background}}}.$$(E.5)

By assuming that the background object is located at 118.3 pc, namely, 5 pc further away than the system HD 120326 (d_HD120326 = 113.3 pc), we derived scattering angles ϕ = 0.0073° and ϕ = 0.065°, for dust grains at the close (θ = 22 pixels = 1.1″ = 0.00031°¹¹) and far (θ = 195 pixels = 9.9″ = 0.0028^°) projected separations in Fig. E.2, respectively. The distance of the background object is very likely underestimated, but this is done on purpose to maximize the difference of the scattering angles, and so to maximize the difference of the scattering phase function values, to be able to explain such a strong flux asymmetry.

By assuming Mie's theory, astrosilicate grains distributed from 0.1 to 100 μm with a power law in -3.5, we computed the scattering phase function using the open-access tool optool¹² (Dominik et al. 2021) at the central wavelength of the STIS filter, 0.59 μm. Due to sampling constraints, we obtained with the code optool values of the scattering phase function of 1.97 × 10⁶ and 1.54 × 10⁶ for scattering angles of 0° and 0.0625°, respectively. This corresponds to a ratio of 1.28, so a positive increase of 28%. By considering the scattering angles of interest, 0.0073° and 0.065°, the positive increase correspond to is therefore 28% × (0.065 – 0.0073)/(0.0625 – 0) = 26%. Consequently, we expect that the dust grains with a scattering angle of 0.065^° scatter 26% less starlight than those with a scattering angle of 0.0073°. This results in an expected flux of 0.94/1.26 = 0.75 μJy at 0.065°. This is higher than the flux (0.52 ± 0.03 μJy) measured in the STIS image at such a location.

One can formalize the expression of the flux expected at the far location F_far,expected as a function of ϕ_far – ϕ_close, and so as a function of d_background using the Eq. (E.5), as follows: $$\begin{gathered} F_{\text{far,predic}} \\ =\frac{F_{\text{close,predic}}}{\left({\displaystyle \frac{{\mathtt{\text{SPF}}}_{\mathtt{\text{optool}}}\mathtt{\text{(}}{\varphi }_1\mathtt{\text{)}}}{{\mathtt{\text{SPF}}}_{\mathtt{\text{optool}}}\mathtt{\text{(}}{\varphi }_2\mathtt{\text{)}}}-1}\right)\left({\displaystyle \frac{{\varphi }_{\text{far}}-{\varphi }_{\text{close}}}{{\varphi }_2-{\varphi }_1}}\right)+1} \end{gathered}$$(E.6) $$=\frac{F_{\text{close},\text{predic}}}{\left(\frac{{\text{SPF}}_{\text{optool}}\left({\varphi }_1\right)}{{\text{SPF}}_{\text{optool}}\left({\varphi }_2\right)}-1\right)\times \frac{\left({\theta }_{\text{far}}-{\theta }_{\text{close}}\right)}{\left(1-\frac{d_{\text{HD 120326}}}{d_{\text{background}}}\right)\left({\varphi }_2-{\varphi }_1\right)}}+1$$(E.7)

By isolating d_background, we get $$\begin{array}{l}{d}_{\text{background}}\\ =\frac{d_{\text{HD 120326}}\times \left(\frac{F_{\text{close,obs}}}{F_{\text{far,obs}}}-1\right)}{\left(\frac{F_{\text{close,obs}}}{F_{\text{far,obs}}}-1\right)\left(\frac{{\text{SPF}}_{\text{optool}}\left({\varphi }_1\right)}{{\text{SPF}}_{\text{optool}}\left({\varphi }_2\right)}-1\right)\left(\frac{{\theta }_{\text{far}}-{\theta }_{\text{close}}}{{\varphi }_2-{\varphi }_1}\right)}\end{array}$$(E.8)

Using the values given previously, we found that the distance of the background object is expected to be 114.9 pc in order to reproduce the observed flux asymmetry between the close and far projected locations shown in Fig. E.2.

We also tested different assumptions on the distribution of dust grains (0.1–1 μm, 0.1–10 μm, 0.1–50 μm, 1-10 μm, 150 μm, 1-100 μm). Regarding the distributions of dust grains 0.1–50 μm, 1-50 μm, and, 1-100 μm, the distance of the background object is predicted to be smaller, between 113.5 and 114.7 pc. In the case of the other distributions of dust grains (0.1–1 μm, 0.1–10 μm, and, 1-10 μm), the scattering phase function is actually slightly higher at ϕ₂ than ϕ₁, so the flux predicted at ϕ_far would be slightly higher than the flux at ϕ_close.

To conclude, based on our investigation presented in this appendix, we do not expect that the alignment of the background object with the southeastern part of the dust could explain why the flux asymmetry observed with HST/STIS is that strong between the southeastern part and the southwestern part of the debris disk around HD 120326.

References

All Tables

All Figures

Fig. 1 Polarimetric images of the system HD 120326 using VLT/SPHERE at 1.6 μm (broad band H; scattered light; epoch 2018-06-01). We processed the data by using the pipeline IRDAP. From left to right, we show the Stokes vector components pI (intensity), Qϕ (azimuthal or radial polarization), and Uϕ (polarization oriented at ±45° with respect to the azimuth) calibrated in μJy/arcsec2, and the signal-to-noise (S/N) map, which is derived on the basis of the Qϕ and Uϕ images. For these images, as all the ones in this paper: north is up and east is left, while the gray circular mask corresponds to the region within the inner working angle of the coronagraph.

In the text

	Fig. 2 Total intensity ("I") images at about 1.6 μm reduced using the processing PCA ADI on the system HD 120326, calibrated in contrast unit (top) or in S/N (bottom). Two structures can be seen: the inner belt and a second, outer structure (see the labels in the top-right image). All images have the same orientation, with north up and east on the left.
In the text

Fig. 3 Comparison of the disk structures in HD 120326 detected in scattered light (SPHERE, 1.6 μm) and in thermal emission (ALMA, 1.3 mm). On the left, we show the polarized intensity SPHERE image with in overlay the contours of the inner belt determined for a signal level of 2 mJy/arcsec2 on this image convolved with a gaussian to obtain smooth contours. The contours are added on the total intensity SPHERE image (middle) and on the zoom-out ALMA reconstructed image (right). In the middle image, we added a white arrow to point out the southwestern part of the second structure seen in total intensity, which does not overlap with the signal seen in polarized intensity.

In the text

	Fig. 4 Optical coronagraphic HST/STIS observations of HD 120326 centered behind the WedgeA1.0″. We applied a spatial binning 2 × 2 pixels2 → 1 × 1 pixel2. The flux unit is μJy/arcsec2. The areas hidden by the wedges (WedgeA and WedgeB) and the diffraction spikes (the cross-pattern) are masked in gray.
In the text

Fig. 5 Overview of the young debris disk around HD 120326. The black-to-white images are the S/N map of the optical HST/STIS coronagraphic data, with a spatial binning 2 × 2 pixels2 → 1 × 1 pixel2. For reference, in the center of the images we display the inner structures detected with SPHERE in total intensity (top), in polarimetry (middle) and with ALMA (bottom). The areas hidden by the coronagraphs and the HST diffraction spikes are masked in dark gray.

In the text

Fig. 6 Results of our MCMC simulations, using the Pythonized version of the GRaTeR tool to generate the belt model. We show the best models in linear or logarithmic scales (top and middle), and the associated residuals normalized by the noise (“σ”, bottom). Scattered light images for the polarized intensity images are on the left and those for total intensity images on the right. Some signal is still visible in the total intensity residuals, and is indicated by the white arrow. Unlike the disk model in total intensity, the disk model in polarized intensity accounts for Rayleigh scattering, hence visual brightness differences between these disk models depending on the scattering angle. The best-fitting parameters are given in Table 1. As a reminder, since the disk inclination found is 78°, the scattering angle is 12° at the closest location of the belt, 90° in the anses, and 168° at its furthest location (see also Fig. 5 from Perrin et al. 2015, to have an insight of the scattering angles around a dust belt).

In the text

Fig. 7 Scattering phase function for our best model of the polarized and total intensity observation, fit independently or jointly. For comparison, we added the model retrieved by Olofsson et al. (2022b) for the polarized intensity data. We hide the part below 12° and above 130°, because we do not detect the disk for such scattering angles, neither in polarized or total intensity light (see text).

In the text

	Fig. 8 Scattering phase function of the inner belt of HD 120326 seen in total intensity compared to other extrasolar debris disks (left) and different dust populations in the Solar System (right).
In the text

	Fig. 9 Reflectance spectrum of the inner belt of HD 120326. For each point, we derive the surface brightness of the disk (SBdisk) and the stellar flux (F⋆) over a spectral range shown by the horizontal line around each point, and also by the colored areas (IFS) or the filter transmission (IRDIS) at the bottom.
In the text

	Fig. 10 Scattered light images of one dust belt generated using the code betadisk, with the associated residual map normalized by the noise ("σ"), in polarized and total intensity (left and right, respectively). The second structure seen with SPHERE (indicated by the white arrow) is somehow more removed, compared to the GRaTeR modeling (Fig. 6, bottom-right).
In the text

	Fig. 11 Spectral energy distribution of HD 120326 based on archival data (see main text) and the new ALMA measurement at 1.3 mm.
In the text

	Fig. 12 Dust mass and fractional luminosity of the debris disk around HD 120326 (diamond marker) compared to those of the sample used in Krivov & Wyatt (2020) (circle markers). The dust mass and fractional luminosity of HD 1203 26 is in particular very similar to β Pic. The color represents the location of the dust belt relative to its host star.
In the text

	Fig. 13 Dust color at scattering angles between 80 and 100° for the debris disk around HD 120326 compared to other debris disks (Ren et al. 2023b). The trend and shaded areas corresponding to the uncertainties of the trend at 1σ, 2σ and 3σ are from Ren et al. (2023b) based on their analysis (see their Sect. 3.2).
In the text

	Fig. 14 Very eccentric, wide-separation Jupiter-mass exoplanet that could generate a spiral-like feature spanning 100 of au, during the timescale of the age of the system. The blue curve shows part of the planetary orbit used in our N-body simulations.
In the text

	Fig. B.1 Total intensity images in the broad band H and reduced using an iterative PCA ADI or ADI+RDI algorithm.
In the text

Fig. B.2 Images depicting (from top-left to bottom-middle) the total intensity ("I") in different sub-coverages of the bands Y and J, reduced using PCA ADI. The bottom-right image corresponds to the average of the PCA ADI processed images over the full YJ bands, and that we spatially binned (3 × 3 pixels2 → 1 × 1 pixel2) to reveal the western part of the second, outer structure. North is up, east is left for all the images.

In the text

	Fig. C.1 Region of the SPHERE polarized (left) and total (right) intensity images considered to model the inner belt in Sect. 3.1.
In the text

Fig. D.1 Posteriors of the MCMC exploration of the parameter space for the disk seen in NIR (1.6 μm) SPHERE observations in polarized and total intensity. The vertical dashed lines represents the values of the median and at ±1σ, which are defined as the 15.9% and 84.1% percentiles to encompass 68.2% of the solutions for a given parameter. The blue vertical line represent the best-fitting parameters, which are the values minimizing the ${\chi }_r^2$.

In the text

	Fig. D.2 Posteriors of the MCMC exploration of the parameter space for the disk seen in continuum millimeter (1.3 mm) interferometric ALMA observations. The parameters are: position angle (PA), inclination (i), radius (r), and width (dr).
In the text

Fig. E.1 Difference of the relative astrometry of the bright point-source lying on the spiral-like feature (see Fig. 4) between the HST/STIS epoch (2013-06-20, used as reference) and the SPHERE/IRDIS epochs (2016 to 2019, see Table A.1). The blue solid curve shows the expected track of a background object, based on the parallactic and proper motions of the primary star HD 120326. The point source is clearly identified as a background object.

In the text

Fig. E.2 Location of the two regions of interest of dust grains in the extended dust emission seen in the HST/STIS reduced image. The dotted blue circle corresponds to the region located at a projected distance close to the background star, while the solid pink circle corresponds to the region located a projected distance far to the background star. Unlike Fig. 4, no spatial binning is applied on this image.

In the text

	Fig. E.3 Sketch representing the parameters (angles and distances) used to express the scattering angle ɸ of the background stellar light scattered by a dust grain in the system HD 120326.
In the text