© A. Boccaletti et al. 2019

1 Introduction

Debris disks correspond to a late stage in the evolution of planetary systems when the primordial material has been expelled out of the systems or incorporated into planets and other bodies. The observed dust results from collisions among the rocky planetesimals or is deposited by comets. Planets, if already formed, are expected to produce indirect signatures in the form of a departure from a pure axisymmetrical disk morphology. Recently, the advance of high-contrast imaging, in particular with the installation of Spectro-Polarimetic High contrast imager for Exoplanets REsearch (SPHERE, Beuzit et al. 2019) and Gemini Planet Imager (GPI, Macintosh et al. 2008), has yielded a significant number of new discoveries in this field either revealing new disks (Lagrange et al. 2016; Kalas et al. 2015; Currie et al. 2015; Bonnefoy et al. 2017; Sissa et al. 2018) or new structures in known disks (Boccaletti et al. 2015; Perrin et al. 2015; Garufi et al. 2016; Perrot et al. 2016; Milli et al. 2017). These observations are definitely pointing to the presence of planets.

Of particular interest is the ever-growing number of disks in which multiple belts are observed due to significant gain in angular resolution and contrast, both in the thermal emission and scattered light regimes. On the one hand, the sub-millimeter interferometer Atacama Large Millimeter Array (ALMA) has revealed obvious cases of gaps in several protoplanetary disks likely sculpted by sub-Jupiter-like planets (Dipierro et al. 2015; Nomura et al. 2016), as well as in one debris disk (Marino et al. 2018). On the other hand, at shorter wavelengths in scattered light, gaps were also found in the protoplanetary disks of TW Hya (Rapson et al. 2015a; van Boekel et al. 2017) and V4046 Sgr (Rapson et al. 2015b), for instance. A few cases featuring an alternance of gaps and belts were also observed in some debris disks such as HD 131835 (Feldt et al. 2017) and HD 141569 (Perrot et al. 2016). However, these disks contain gas (Kral et al. 2017), which might also be responsible for developing belts (Takeuchi & Artymowicz 2001; Kral et al. 2018) or arcs (Lyra & Kuchner 2013; Richert et al. 2018). Despite the fact that a double-belt structure has been inferred for several systems from analyses of spectral energy distribution (SED; Pawellek et al. 2014; Pawellek & Krivov 2015), so far a single gasless debris disk featuring two belts has been unambiguously imaged around HIP 67497 by Bonnefoy et al. (2017). This potentially new class of debris disks is different from the systems with inner (warm) and outer (cold) components as inferred from photometric measurements at mid- and far-IR (Chen et al. 2014). Instead the multiple belts we see in scattered light with high-contrast imaging are all located at several tens of astronomical units, hence rather cold. It is tempting to hypothesize that multiple-belt systems could be one particular stage in the history of debris disks, later evolving as single-belt systems with a broad inner depleted cavity once all the planets have cleared out their orbits. Therefore, the moment when these substructures form and when they disappear is crucial to understand the formation and architecture of planetary systems. Observing one specific system provides a single snapshot in a disk lifetime. To capture the “big picture”, several systems with various morphologies and at different stages of evolution must be found and studied.

The star NZ Lup (HD 141943, TYC 7846-1538-1, G2, V = 7.97, H = 6.41) is known to harbor a debris disk first inferred from Spitzer photometry (Hillenbrand et al. 2008). Gaia Collaboration (2018) measured a distance of ${60.34}_{-0.18}^{+0.19}$ pc. The IR excess in the SED is modeled as two blackbodies peaking at equivalent temperatures of 197 and 60 K and corresponding to physical separations of respectively 4 and 122 au (assuming d = 67 pc, Chen et al. 2014). The cold component is in fact poorly constrained with just a single photometric point measured by Spitzer/MIPS at 70 μm for which F_{70 μm}∕F_* = 15.66 ± 2.3 (Hillenbrand et al. 2008). As a young star (see Sect. 2) it has been a target for exoplanet searches by direct imaging (Chauvin et al. 2010; Galicher et al. 2016; Vigan et al. 2017) but none of these studies reported hints of a disk in scattered light. While the inner belt is presumably too close to the star (~67 mas) to be resolved, the outer belt was finally detected in HST/NICMOS (Hubble Space Telescope/Near Infrared Camera and Multi-Object Spectrometer) archival data in which it appears nearly edge-on (i = 85°) with a maximum intensity of ~0.25 mJy.arcsec⁻² (Soummer et al. 2014). The angular resolution of NICMOS did not permit confirmation of the physical size of the belt inferred from the SED, but the signal of the scattered light is detected from 0.7″ to 2.5″. Using the high-precision spectrometer HARPS, Lagrange et al. (2013) found no planet more massive than 1–5 M_J for periods shorter than ~100 days (~0.4 au).

In this paper we present the discovery of two cold belts in the debris disk of NZ Lup. The characteristics of the star are presented inSect. 2, while the observations and data reductions are provided in Sect. 3. A general description of the disk morphology is presented in Sect. 4, and this geometry is studied in more detail using modeling in Sect. 5. The SED is revisited in Sect. 6. We provide the astrometric characterization of the point sources contained in the field of view as well as the estimation of the limits of detection in Sect. 7. Finally, we discuss the implications of double belt structure with respect to the presence of planets (Sect. 8).

2 Host-star properties

The object NZ Lup is a bona fide young star (age < 50 Myr), as revealed consistently by a variety of indicators (e.g., strong lithium line, fast rotation, strong magnetic activity). It lies in front of the Upper Centaurus Lupus (UCL) group and shares similar kinematic parameters. For this reason, it was proposed as a UCL member (age 17 Myr) by Song et al. (2012).

There are indications of a modest amount of reddening (E(B–V) = 0.03–0.05), comparing the expected colors for a G2 star (spectral type from Torres et al. 2006) to the Pecaut & Mamajek (2013) sequence for young stars.

A robust isochrone age determination was prevented recently by a lack of trigonometric parallax. Exploiting Gaia DR2 (Gaia Collaboration 2018) parallax and coupling it with effective temperature from spectral type (5870 K) and observed magnitudes in different bands shows that the star has not yet settled on the main sequence and is therefore very young (Table 1 and Fig. 1). Comparison of the stellar parameters from Table 1 (including the small amount of reddening we derived) with the models of Bressan et al. (2012) yields an age of 16 ± 2 Myr. Marsden et al. (2011) mentioned the possibility of binarity, because of a marginal radial-velocity variability. Unresolved binarity could explain the off-sequence position on the color-magnitude diagram (CMD). However, higher-precision measurements by Lagrange et al. (2013) with High Accuracy Radial velocity Planet Searcher (HARPS) indicate only a moderately large scatter linked to magnetic activity, as demonstrated by the RV-line bisector correlation. An analysis of SPHERE images presented in the following section (including the noncoronagraphic ones taken during target acquisition) is also missing any evidence of multiplicity down to a very small separation of about 40 mas (2.4 au). We therefore dismiss the possibility of binarity and consider the derived isochrone age to be reliable.

The isochrone age is further supported by indirect indicators, such as lithium and magnetic/coronal activity, that are broadly compatible with ages of between 10 and 50 Myr. The rotation period is slightly longer than that of stars of similar color in the β Pic moving group (Messina et al. 2017), consistent with the slightly younger age derived with isochrone fitting.

Finally, the kinematic parameters are very similar to the UCL ones, although the star lies at a shorter distance than the bulk of UCL members (about 140 pc). The Banyan Σ online tool (Gagné et al. 2018) yields a probability of 55% of UCL membership, and no significant membership probability for other groups. We therefore conclude that NZ Lup is a star with UCL age and kinematics but that it is in front of the main body of the group. A link between this target and the UCL seems probable but the evaluation of the actual extension of the Sco–Cen groups at distances much smaller than 100 pc is beyond the scope of the present paper. We therefore adopt an age of 16 Myr.

Fig. 1 Color-magnitude diagram of NZ Lup. Overplotted isochrones for 5-12-20-30 and 70 Myr from Bressan et al. (2012).

3 Observations and data reduction

SPHERE (Beuzit et al. 2019) is the extreme AO (Fusco et al. 2014) instrument installed at the Very Large Telescope (VLT, ESO-Chile) equipped with coronagraphs (Boccaletti et al. 2008), which routinely delivers high-contrast imaging data of a large survey of young stars that began in 2015 (SpHere INfrared survey for Exoplanets, SHINE, Chauvin et al. 2017). We observed NZ Lup during guaranteed time on May 25, 2016, May 31, 2016 and April 30, 2017, using the IRDIFS mode of SPHERE, in which both the Infra-Red Dual-band Imager and Spectrograph (IRDIS) (Dohlen et al. 2008), and the Integral Field Spectrograph (IFS) (Claudi et al. 2008), are operated simultaneously. The observing log is displayed in Table 2. The first epoch was set up with a broad band filter (BB_H) for disk-detection purposes, while the other two epochs were using narrow bands (H2 = 1.593 μm, H3 = 1.667 μm, R ~ 30) for exoplanet-detection purposes (Vigan et al. 2010). The IFS was configured in YJ mode (0.95–1.35 μm, R ~ 54).

The sequence of observations is as follows: (1) the target acquisition optimizes the position of the star onto the coronagraphic mask to set up reference slopes for the wavefront sensor; (2) then the star is offset by ~ 0.5″ from the mask to obtain a flux calibration (PSF); (3) the star image is sent back onto the mask and a waffle pattern is applied on the deformable mirror to create four satellite spots at a separation of 14 λ∕D for centering the target in detector coordinates; (4) the waffle pattern is removed and the science exposures start while the detector is dithered on a 4 × 4-pixel grid to further reject bad pixels; (5) both centering frames and flux calibration are repeated; (6) finally the AO loop is opened and the telescope moves to the sky background. All coronagraphic images were acquired with an APodized Lyot Coronagraph (APLC, Soummer 2005), the focal mask of which is 185 mas in diameter combined to an apodizer which transmits 67% of the light (Carbillet et al. 2011).

The IRDIS and IFS data are reduced at the SPHERE Data Center¹ (Delorme et al. 2017) using the SPHERE pipeline (Pavlov et al. 2008) and following a standard cosmetic reduction (sky subtraction, flat field correction, bad pixel removal). Raw frames are corrected for distortion (Maire et al. 2016). The star position in the image is determined from the satellite spots as detailed in Boccaletti et al. (2018) and no further recentering is performed. More details are provided in Boccaletti et al. (2018). The north orientation is calibrated with astrometric reference fields (Maire et al. 2016). The IRDIS and IFS pixel scales are 12.25 and 7.46 mas, respectively.

Starting from the output of this reduction, the four-dimensional data cubes (spatial, spectral and temporal dimensions) were processed with SpeCal, the differential imaging implementation at the SPHERE data center (Galicher et al. 2018). Several types of angular differential imaging (ADI) techniques were considered (Marois et al. 2006, 2014; Lafrenière et al. 2007; Soummer et al. 2012). The study presented here is based on the principal component analysis (the KLIP algorithm; Soummer et al. 2012), but other algorithms provide similar images. In practice, the resulting KLIP image depends on the number of the lowest-order modes that are kept to build a reference frame, which acts as a spatial filtering of low-frequencies to reject stellar residuals.

4 General description

The images processed with ADI using the KLIP algorithm are displayed in Fig. 2. The top panel presents the two IRDIS observations from May 2016, while those from April 2017 are shown for both IRDIS and IFS at the bottom. The best-quality data are achieved for the last epoch, as is obvious from the images; these correspond to the best seeing and coherence time parameters (Table 2) as well as the largest contrasts (see Sect. 7). The reference frames that are subtracted out from the data set in the KLIP procedure were obtained with 10 modes for all reductions, except the one from May 25, 2016 (IRDIS), which used 40 modes. This larger number of modes (compared to the amount of frames; see Table 2) is the direct consequence of lower-quality data and a broader stellar halo. The disk analysis below is based on the April 2017 data.

The disk is oriented southeast to northwest at a position angle (PA) of 146.53 ± 0.15°, measured following the procedure described in Boccaletti et al. (2018), in which the error bar includes the measurement uncertainty (~ 0.1°) together with the TN uncertainty (~ 0.1°). The general aspect of the disk corresponds to a very inclined ring (~85°) of which the ansae are located at about 1.3″ southeast and 1.4″ northwest, while only one single side is visible. This bright side should presumably be the front side if forward scattering is predominant as expected for small dust grains. Although the main disk stops at ~ 1.5″ on both sides, we see a signal of dust scattering in the disk direction out to ~ 2.2″ on both sides. This global morphology agrees well with the NICMOS image although a higher angular resolution is achieved with SPHERE (Soummer et al. 2014).

A closer examination of the disk reveals two unusual characteristics. First the southeastern ansae features a break, coincident with the location of a point source (indicated with an arrow in Fig. 2). This object was removed in each frame of the data cube by subtraction of a scaled PSF and the break is still observed, indicating that it is not induced by an ADI artifact caused by the overlap of a point source on top of the disk image. In any case, this point source is not related to the system but is flagged as a background star (Sect. 7), and so cannot be responsible for a dynamical effect on the disk.

Secondly, and most importantly, the disk splits in two parts at radii closer than ~ 0.6″ as if it were an inner ring. This disk splitting is detected unambiguously in both IRDIS and IFS images from April 2017 (arrows/ellipse in Figs. 2 and A.1, bottom) but is also identified in lower-quality data from May 2016 (Figs. 2 and A.1, top right). We measured the spine of the disk (Fig. 3) by fitting a Lorentzian profile on the vertical cross section for each stellocentric distance, and considering either a single component or two components (to account for the disk splitting mentioned above). The single-component fit (black line) has two minima located at about 1.5″ on both sides, which correspond to the edges of the main disk, the one at the southeast being steeper than in the northwest. The maximum elevation with respect to the midplane is about 0.13″. The fit with two components (blue and red lines) was ordered according to the intensity of the component (the brightest being the closest to the midplane). The brightest component (red line) agrees well with the single component especially in the southeast. The faintest (blue line) deviates from the main spine starting at 1″ from the star and inwards while it culminates at ~0.17″. This departure comes in fact with decreased elevation of the brightest component inwards of ~ 0.6″ (in accordance with the ring-like shape in the image). The spine is noisier in the northwest, making it more difficult to distinguish the two components in Fig. 3. Overall, the measurement of the spine quantifies the main pattern that is seen in the image, where the inner ring appears superimposed on the main disk.

While the images convey the idea that a second ring would be sitting on top (higher elevation from the star) of the main ring, one should consider that the actual distribution of dust is altered by the ADI process, and cross talks are to be expected for intricate geometry. In any case, such a configuration would be difficult, if not impossible, to explain dynamically. Instead, we posit that the disk of NZ Lup is composed of two belts of different sizes, separated by a gap, and mutually inclined by a few degrees. Such an assumption was already suggested for the very inclined debris disk of HD 15115 (Engler et al. 2019). The following section is dedicated to the modeling of this structure in order to test our hypothesis.

Fig. 2 Images of the disk obtained with IRDIS in the H band (top and bottom left), and with the IFS in the YJ band in a narrower field of view (bottom right). The intensity scale is arbitrary for each epoch to optimize the disk visibility. The disk midplane is aligned with the horizontal direction and the displayed field of view is 6″ × 3″ for IRDIS and 2.4″ × 1.2″ for the IFS. The main features reported in the text are annotated with arrows and ellipse.

Fig. 3 Spine of the disk measured with a single component assumption (black line) and for two components (red and blue lines).

5 Modeling of the belts

Given that the structures observed in the NZ Lup disk are rather fine and relatively faint, we proceed in several steps for the modeling. We first consider a single belt scenario in Sect. 5.1, where we present the global assumptions to produce scattered light images of synthetic disks. We analyze the residuals between this one-belt model and the actual image to motivate a more refined analysis including two belts. In Sect. 5.2.1, considering that the two-belt model has a rather large number of parameters in regards of the signal-to-noise ratio (S/N), we define density functions to explore the distribution of the dust in between the two belts more specifically. We firstly assume that the outer belt is more inclined than the inner belt with respect to the line of sight as it appears more relevant from the image. We then check that this hypothesis is valid for one single type of model in Sect. 5.2.2.

5.1 One-belt scenario

We used a simplified version of GRaTer (Augereau et al. 1999) to produce synthetic images of debris disks with no particular assumption about the grain composition. We assume that dust is collisionally produced from parent bodies located in an axisymmetric narrow birth ring, at a separation r₀ from the star at which we assume adensity n₀, the maximum density in the disk. The edge of the ring observed at about 1.5″ corresponds to approximately 90 au for a distance of 60 pc. We further assume that the parent belt is small with respect to the angular resolution, and that the dust seen outside the parent belt corresponds to very small grains placed there either by PR-drag (in the r < r₀ region) or by high-eccentricity orbits induced by radiation pressure (in the r > r₀ region; Strubbe & Chiang 2006; Thébault & Wu 2008).

The density function is modeled by three components as follows. $$n(r,z)\propto n_0.R(r).Z(r,z)$$(1)

The radial profile function R(r), with r being the radial dimension in the disk plane, is described by power laws with a maximum at the location of the belt and decreasing inward (α_in > 0) and outward (α_out < 0). $$R(r)={\left({\left(\frac{r}{r_0}\right)}^{2{\alpha }_{\text{in}}}+{\left(\frac{r}{r_0}\right)}^{2{\alpha }_{\text{out}}}\right)}^{-1/2}.$$(2)

The vertical profile function (Z(r, z), with z being the vertical dimension perpendicular to the midplane, is assumed to be Gaussian, and the disk scale height (h) varies linearly with the radius (h = H₀∕r₀, H₀ being the height of the disk at the position r₀). $$Z(r,z)=exp\left({\left(-\frac{\left|z\right|}{h\times r}\right)}^2\right).$$(3)

Finally, the scattering phase function (θ being the scattering angle) is approximated analytically with the Henyey–Greenstein function controlled with the anisotropic scattering factor g. $$f(\theta )=\frac{1-g^2}{4\pi {(1+g^2+2g\times \text{cos}(\theta ))}^{3/2}}.$$(4)

Assuming the scattering of dust particles is preferentially directed forward, g is positive and greater than zero but smaller than one (isotropic scattering). The model is inclined with respect to the line of sight (i) and rotated in the sky plane (PA). We used a forward modeling approach which consists in generating a grid of models, processing them the same way as the data, and calculating a χ² metrics. As a first guess, we started with a single belt geometry and generated a grid of 9600 models with the following parameters (with some priors regarding the acceptable range).

• position angle: PA (°) = 146.5

• inclination: i (°) = 84, 85, 86

• radius of the belt: r₀ (au) = 65, 70, 75, 80, 85, 90, 95, 100

• anisotropic scattering factor: g = 0.4, 0.5, 0.6, 0.7

• slope of the surface density inward: α_in = 2, 4, 6, 8, 10

• slope of the surface density outward: α_out = −2, −3, −4, −5, −6

• scale height: h = 0.01, 0.02, 0.03, 0.04.

To compare the models with the data and to determine the best model parameters, we proceed as in Engler et al. (2019). Each model is convolved with the PSF, and projected onto the same KLIP basis as used to reduce the data in order to reproduce a similar level of self-subtraction on synthetic images. This new grid of models is the input of the minimization procedure. We define a rectangular aperture globally aligned with the disk axis to encompass the disk image, the length of which is 6″ in total and is 0.4″ in width. Because half of the disk (the northeast) is visible, the aperture is offset by 0.1″ with respect to the star along the minor axis to avoid including too much noise in the aperture. The reason the aperture extends out to 3″ although the main ring is located at ~1.5″ is to take into account the scattering signal observed out to at least 2.2″ on both sides, which is crucial to constrain the α_out parameter accounting for the halo of small grains beyond r₀. The central part of the image is also removed numerically to hide the strongest stellar residuals. We tested three different cases for which the central masked region is 1.0″, 0.5″, and 0.3″ in radius (mask1, mask05, mask03). However, since we are approximating a two-belt disk with a single belt model, we are mostly interested in the disk parameters at large separations in this first step. A 2 × 2-pixel binning is applied to the data and the models. For each model in the considered aperture, we derived the intensity scale which minimizes the quadratic difference between the data ($\mathcal{O}$) and the model ($\mathcal{M}$). We then obtained the reduced χ² with the standard relation: $${\chi }^2=\frac{1}{\nu }{\displaystyle \sum _{i,j}^{N\text{data}}}\frac{{(\mathcal{O}(i,j)-\mathcal{M}(i,j))}^2}{\sigma {(i,j)}^2}.$$(5)

The degree of freedom (ν) is taken as N_data − N_params with N_params = 6 in this first case, with N_data being the number of data points in the aperture after binning. The minimum value of χ² provides the parameters for the best model. To determine the range of models which best matches the data in a conservative way, we considered a threshold at 1% of the lowest χ² values instead of adopting the standard $\sqrt{2/\nu }$ threshold. The latter theoretically applies to Gaussian noise and linear models, two conditions not necessarily satisfied in the present case. The parameter values and dispersions of the best-fitting models are provided in Table 3, together with parameter frequencies of the 1% best models in Fig. B.1. The surface density slopes are consistent for all three cases suggesting a steeper density inwards than outwards (α_in > 5, while α_out = −4 to −5). The latter is also steeper than the canonical value of −1.5 expected for a collisionally produced halo of small grains placed by radiation pressure beyond a main birth ring (Strubbe & Chiang 2006; Thébault & Wu 2008). However, the acceptable values of α_in can be very broad when the mask size decreases, so it is difficult make conclusions from this first approach. The best inclination is about approximately 85− 86°, the position of the planetesimal belt is at approximately 90−95 au, and the asymmetric scattering factor is g ~ 0.6 indicating arather high value of forward scattering (but similar to other debris disks: Lagrange et al. 2016; Olofsson et al. 2016; Bonnefoy et al. 2017). Regarding the scale height, both h = 0.01 and h = 0.02 fit equally well at least in the case of model mask1. Images of the model compared to the data, together with residuals, are shown in Fig. 4. The single belt is clearly unable to reproduce the disk splitting seen at stellocentric distances shorter than 1″.

Fig. 4 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the best forward model (i = 85°, r = 95 au, αin = 10, αout = − 5, g = 0.6, h = 0.02) and the residuals. The field of view is 6″ × 0.8″.

Fig. 5 Examples of the density (blue) and intensity (red) functions corresponding to the four model families.

5.2 Two-belt scenario

We now consider a more realistic configuration of two separate belts at stellocentric distances r₁ and r₂, and with independentparameters except for their relative inclinations. Due to the time-consuming nature of the forward modeling approach, we restrained the range of solutions by fixing some of the parameters. Following the single-belt modeling, we fix g = 0.6 for the two belts. The ADI process is known to significantly bias the images of disks, and in particular can emphasize the darkness of a central cavity for ring-like geometries. Here, we are particularly interested in the density distribution in between the two belts, and therefore to investigate whether this region is filled or empty we defined four families of radial densities by changing the surface density slopes (α_in-1, α_out-1, α_in-2, α_out-2²): fill, ring, gap, 2gaps. The values of the density slopes are provided in Table 3 and an illustrative sketch is depicted in Fig. 5. For the gap and 2gaps cases the density slopes between r₁ and r₂ were chosen to be steep enough (+ ∕−20) to avoid mutual overlapping. In addition, the 2gaps models feature an additional inner gap (r < r₁). For the fill and ring cases the outer belt is truncated inwards of r₂. In the fill models the outer density slope of the first belt is fixed at −3.5 to avoid discontinuities between the belts. However, this slope is an average value since it does not match for any separations between r₁ and r₂. Finally, the ring case features an increasing density function to allow a flat intensity variation. The outer disk density slope is the same for any model family (− 4 or − 5). In the case of gap and fill models, the inner slope is symmetrical with the outer one (α_in-1 = − α_out-2). Similarly, the scale height is allowed to take two values but be identical for the two belts (h = 0.01, 0.02).

5.2.1 Edge-on outer belt

Once density functions are defined, the next parameters are the inclinations and the locations of the belts. Intuitively, from the visual inspection of the image we assume that of these two belts the most inclined should be the outer one. In that case, we defined the following ranges of parameters, again with some priors regarding the acceptable range:

• inclination of the outer belt: i₂ (°) = 85, 86, 87, 88

• relative inclination i₂−i₁: (°) = 4, 5, 6, 7

• radius of the inner belt: r₁ (au) = 60, 65, 70, 75, 80, 85, 90

• radius of the outer belt: r₂ (au) = 90, 95, 100, 105, 110, 115, 120.

In the parameter space, we intentionally matched the largest value explored for r₁ with the minimal value explored for r₂ to check if unrealistic situations (same radii but different inclinations) could come out from the analysis. Here, we used GRaTer twice to model each belt individually and then the two synthetic images were co-added to produce a two-belt model. This is not exactly similar to defining a global analytical expression of the radial density but it allows more flexibility in adapting the relative weight and inclination of the belts (for optically thin disks the intensity is proportional to the density). Taking the phase function into account, we expect very strong forward scattering for very inclined disks, so a disk image in scattered light should be very strongly peaked at small phase angles. Furthermore, this peak quickly decreases in intensity as the scattering angles increase (or conversely the inclination decreases). In the two-belt configuration the less inclined (inner) belt therefore appears much fainter for a similar dust density to the more inclined (outer) belt, although it receives more stellar flux. This behavior is illustrated in Fig. C.1. However, we observe in the image a nearly identical intensity for these two belts, so implicitly the inner belt should have a larger density to match the images. Qualitative tests led us to multiply the scattered light image of the inner belt by a factor two.

The aperture in which the χ² is evaluatedis similar to that of the single-belt case but reduced to 3″ instead of 6″, the outer slope of the surface density being already constrained at the former stage. However, the choice of the aperture size has clearly some impact on the model fitting. A larger aperture tends to increase the χ² and lead to similar solutions for all four of the families of models. We tested two approaches: either a forward modeling, where, as before, a perfect model (no noise) is compared directly to the data, or alternatively, the same model is injected into the data cube with a brightness ratio of 6 × 10⁻⁶ (slightly brighter than the real disk), and processed the same way, but reversing the parallactic angle sequence to cancel out the disk while keeping the noise structure. In both cases, the eigenvectors of the PCA are determined on the data without any fake disk; these are then stored and reapplied on the noiseless model, or when the model is added to the data. These two approaches yield mostly the same outcome and so the current results are based here on the latter for practical reasons.

The χ² for all tested models ranges from 4.15 to 11.75, globally indicating poor fits. The smallest χ² are obtained for the 2gaps family (${\chi }_{\text{min}}^2=4.15$), then for ring and gap (${\chi }_{\text{min}}^2=4.32$ and 4.43 respectively), while fill is clearly worse (${\chi }_{\text{min}}^2=4.64$). Overall, these values do not differ significantly, which means that the χ² metric has some limitations to identify fine structures at low S/N in this particular case, while visually it can be straightforward to make a distinction between some models. As in the one-belt scenario, we selected the best 1% models for each of the four model families. The frequency of a parameter is defined as the occurrence of a value in the 1% best models. Table 3 provides the most likely values for each parameter together with the dispersion estimated from a Gaussian fit of the frequency distribution when relevant (the parameter sampling is large enough and/or there is more than a single possible value). We do not provide a dispersion for α_in-1, α_out-2, and h since only two values were tested. The images and the frequencies of each parameter for the best model family, 2gaps, are displayed in Figs. 7 and 6, respectively. Similar figures are shown in the Appendix for the other model families (Figs. D.1–D.3, and B.2–B.4). It is rather clear that only the best models for 2gaps and gap feature a disk splitting similar to the one observed in the data. However, there are some parameter combinations in the 1% lowest χ² for fill and ring which can produce this pattern too.

For the inclination of the belts, the results are highly significant. The outer belt inclination (i₂) has a very high peak with frequencies as large as 80–100% (2gaps, gap, fill), while for the inner belt (i₁) the significance is lower but still achieves frequencies larger than 50%. The ring model family provides flatter histograms. Based on the minimum χ² values, we obtained i₁ = 82° and i₂ = 87°. We observed more dispersion on the radii of the belts, which is a direct consequence of the disk being seen at high inclination. If we exclude the ring family, which converges to r₁ = r₂ = 90 au, a likely unrealistic situation, then the external belt has its most likely position at about r₂ ~ 105−115 au, while the inner belt is located at r₁ ~ 85−90 au. The frequency peak for r₁ and r₂ can be as small as ~30%, hence rather broad uncertainties of 5–10 au. The best model family, 2gaps, provides the largest separation (30 au ≈ 0.5″) between r₁ and r₂, while the opposite is true for the fill case (only 15 au ≈ 0.25″). As for the scale height, there is no strong prevalence of one value with respect to another except for the 2gaps model which favors h = 0.02. This parameter is therefore not very well constrained in the two-belt analysis, but this is not very surprising since h = 0.02 corresponds to a height of 2 au (33 mas) at a distance of 100 au (about the position of the outer belt), which is also roughly the angular resolution of the images (40 mas at λ = 1.6 μm). Contrary to the one-belt analysis, considering two belts in the system favors a shallower inner/outer slope of 4/−4 instead of 5/−5 although this is still compatible within error bars. We note that we observed degeneracies between the inclination of the inner belt (i₁) and its radius (r₁), in the sense that for the KLIP image to show a disk splitting, this belt could be smaller (approximately 60− 70 au) than what is obtained from the lowest χ², but with a higher relative inclination with respect to the outer belt. However, such models are not part of any of the 1% best models selected. We would need to increase the threshold to ~5% to start observing this degeneracy. Finally, the differences between the 2gaps and gap cases are small enough to consider both of them as reliable descriptions of the disk image.

The intensity of the inner ring was multiplied by a factor of two as an initial guess and so this parameter was not included in the minimization. To provide further constraints on this intensity we now consider the best model (2gaps case, i₁ = 82°, i₂ = 87°, r₁ = 85 au and r₂ = 115 au, h = 0.02, α_out-2 = −4), and we vary the intensityweight of the inner belt from 1 to 3 (hence neglecting possible degeneracies between parameters). The optimal flux ratio is 1.8 ± 0.1 although it corresponds to a very small variation of χ² compared tothe former solution (4.13 instead of 4.15). In fact, the difference in the residual map is very difficult to appreciate by eye.

As a concluding note, there are strong limitations to the modeling of the images of the disk of NZ Lup because of the high inclinations of the belts. The residual maps also show signs of asymmetries along the major axis which were not accounted for. Still, wecan reasonably conclude that the disk is made up of two belts with different sizes, inclinations, and densities.

Fig. 6 Frequencies of parameters obtained for the 1% best models in the two-belt 2gaps case.

Fig. 7 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the corresponding forward model, the best model injected in the data, and the residuals (2gaps case: i1 = 82°, i2 = 87°, r1 = 85 au, i2 = 115 au). The field of view is 6″ × 0.8″.

5.2.2 Edge-on inner belt

Another possible configuration not considered above corresponds to the opposite geometry, where the inner belt is at higher inclination (closer to edge-on) than the outer belt. We first considered the same grid as before, but changing i₂ with i₁ and i₂ − i₁ with i₁ − i₂. The relative intensity is also adapted for the same reasons as provided in the previous section, so the outer belt synthetic image is multiplied by a factor of two. For simplicity, we only focused on the gap model and restrain the parameters α_in-1, α_out-2 and h to a single value (respectively 5, −5, and 0.02). The minimum χ² for this grid of models is 9.81 (corresponding to i₁ = 88°, i₂ = 83°, r₁ = 60 au, r₂ = 90 au), much higher than in the previous configuration. In fact, the best models have a much thicker outer belt than in the data, and the apparent projected separation between the two belts is too large. Therefore, we modified the range of inclination for the outer belt with i₁ − i₂ = 1°−4°. Even then however, the minimum χ² is 6.22 (corresponding to i₁ = 86°, i₂ = 85°, r₁ = 80 au, r₂ = 95 au) showing again that these models do not provide a good match to the data. In that case, the relative inclination is too small to generate an observable cavity between the two belts as in the data. As a consequence, we can confidently rule out this geometry.

5.2.3 Summary

From the modeling work we were able to establish that the one-belt case does not correctly match the data, leaving a significant residual in the inner part of the disk image. Therefore, the derived inclination and radius of the belt are biased and correspond to averaged values of a two-belt geometry. Using a two-belt scenario, we found that the presence of a gap in between the belts is a better match to the data. However, the differences between models with gap(s), a smooth transition, or no gap at all are small because the mutual inclinations and the distances of the belts could mimic a gap in the scattered light image. Some cases can still be rejected nevertheless, like for instance when the model converges to the same radii for the two belts when we consider a broad ring. Finally, models where the inner belt is more inclined than the outer belt can be confidently ruled out. In summary, the modeling favors a geometry with two belts at 85 and 115 au, separated with a sharp gap, and at inclinations of 85° and 87°, respectively.

6 SED modeling

6.1 Stellar parameters and model setup

Considering the identification of two belts in SPHERE observations, we decided to revisit the SED of the system to check for signs of such a structure. The contribution of the stellar photosphere to the total flux density is taken from a PHOENIX model (Hauschildt et al. 1999) assuming a stellar luminosity of 2.9 L_⊙ and a temperature of 6000 K (Chen et al. 2014) for the host star. The photometric data were collected from several catalogs and from different papers summarized in Table 4. The SED of the debris disk was fitted using the SONATA-code (Müller et al. 2010; Pawellek et al. 2014) using all data points longward of 10 μm where we expect the dust excess emission to start, but the lack of data points in the far-IR (apart from those at 70 and 250 μm) is clearly a limitation, preventing us from providing reliable constraints on the dust properties.

As inferred from the scattered light image we assume the disk to lie between 85 and 115 au. The SED is fitted using a grain size and radial distribution model with a power law of the following shape. $$N(r,s)~s^{-q}{r}^{-p},$$(6)

where N(r, s)drds is the number of grains with a size ranging from s to s + ds at a disk radius between r and r + dr. The parameter q is the size distribution and p the radial distribution index. The fitting method is similar to the study of Pawellek et al. (2014). We fix the maximum grain size to 1000 μm and fit the minimum grain size and the size distribution index. The radial distribution index is fixed to 1.5 resembling the profile in a small grain halo beyond a collision-dominated ring (Strubbe & Chiang 2006). We assume pure astronomical silicate (Draine 2003) as dust composition.

6.2 Single broad disk fit

Given the scarcity of the photometric data in the spectral range where the disk emission dominates over the photosphere, and the relative proximity of the two belts as inferred from direct imaging, the fitting of the two components in the SED is an underdetermined problem. For the sake of simplicity, we assume a single broad disk between 85 and 115 au as a first approach. The SED fitis shown in Fig. 8.

We find a minimum grain size of 2.27 ± 0.40 m and a size distribution index of 3.45 ± 0.22. The reducedχ² is found to be 11.2. This high value is caused by the relatively small error bars given in the literature. The corresponding mass is 1.60× 10⁻² M_⊕.

Fig. 8 Flux density as a function of wavelength. The disk lies between 85 and 115 AU.

6.3 Double belt fit

In the following approach, to avoid fitting an overly small number of data points with too many model parameters, we assume an inner component to be a pure blackbody ring without any size distribution and an outer component between 105 and 115 au. We fitted the blackbody radius for the inner component, while for the outer component we fitted the minimum grain size and the size distribution index. Comparable to the SED fitting of AU Mic (Pawellek et al. 2014) the code finds the best fit by setting the mass of the inner component to zero. This means that the best fit is a single-component model, and that a pure SED fit cannot discriminatebetween one single disk extending from 85 to 115 au and two rings in the same region. This is not completely unexpected, given the relatively narrow radial extent of the whole region and the relatively limited differences in terms of grain temperatures and thermal emissions.

However, our detailed SED fit rules out the need for an additional warm belt at ~4 au that was suggested by the SED fit of Chen et al. (2014) assuming pure blackbodies. This is because while a broad SED can only be explained by different spatial locations when assuming blackbodies, a model that takes into account size distribution and size-dependent radiative properties can result in an extended range of temperatures, and thus a broader SED, for a single spatial location. We note that we do not rule out the presence of a warm belt at 4 au, which would be undetectable with the angular resolution of SPHERE, but rather its signature in the global SED.

7 Point sources and limit of detection

We identified 16 point sources in the IRDIS field of view (Fig. E.1). We used SpeCal to determine the position of these sources in May 2016 and April 2017. A full description of astrometric errors is provided in Galicher et al. (2018). We used the model of planet image procedure that fits an estimation of the image of a point-source in the TLOCI reduced image. The accuracy of the fitting in the 2016 and 2017 images is 5 to 10 mas. We accounted for the proper motion and the parallax of NZ Lup to predict the positions of the point-sources in the 2017 image from the measured positions in 2016 as if they were background sources. We then compared these predictions to the measured positions in the 2017 image. None of the detected sources share the motion of NZ Lup. They are consistent with background stars (see Fig. E.2). Nevertheless, some dispersions are observed between the measured and estimated positions for 2017, possibly indicating that some of these background stars also have detectable proper motions.

The limits of the detection of point sources were estimated from the KLIP reductions both with IRDIS and IFS following the procedure described in Galicher et al. (2018). The contrast at a particular radius is calculated from the standard deviation of the pixel contained in an annulus of 0.5 ×FWHM in width (about 2 pixels) centered on the star. The self subtraction inherent to ADI is estimated with fake planetsinjected into the data (along a spiral pattern to cover a range of separations and azimuths), and is compensated to produce the contrast plot in Fig. 9. Inside the control radius (~ 0.8″ in the H band), we reached a contrast of about 10⁶ in April 2017. At a separation of 2.5″ the contrast is as large as 5 × 10⁷. The IFS outperforms IRDIS at stellocentric distances shorter than 0.4″ when using ADI. In this particular case, the gain when using the spectral diversity of the IFS in addition to ADI to further improve the contrast was found to be negligible and is therefore not presented.

These contrast values are converted into masses considering several atmosphere models calculated for the SPHERE filters, and assuming a stellar magnitude of H = 6.41 for IRDIS and J = 6.74 for the IFS spectral range. The limits of detection for BT-Settl (Allard 2014) are shown in Fig. 9 and the two other cases, DUSTY and COND (Allard et al. 2001), are provided in Fig. F.1. The models are considered unreliable for masses lower than 0.5 M_Jup and therefore lower values are not plotted; although the achieved contrast would in principle allow for the detection of smaller/lighter planets. A planet of 1 M_Jup would have been detected at a separation of 0.5″ (corresponding to about 30 au in projection) according to the BT-Settl model.

Fig. 9 Limits of detection in contrast (top) and converted into Jovian masses (bottom) for two age assumptions (20 and 30 Myr), using the BT-Settl atmosphere model.

8 Discussion

Despite some difficulties to clearly distinguish between the gap and 2gaps families of solutions, our analysis hasestablished that the most likely configuration is that the NZ Lup disk extends from ~85 au to ~115 au and displays both a discontinuity in density (be it with a gap or a sharp transition) and a tilt in inclination between the inner and the outer parts.

We note that HIP 67497 is the only debris disk in which two cold but coplanar belts are observed at distances of about 60 and 130 au (Bonnefoy et al. 2017). The obvious reference for a disk with an inclination tilt is the beta Pictoris system seen edge-on, with its inner ≲ 50 au region tilted by a ~5° angle with respect to the outer disk (e.g., Lagrange et al. 2012). This tilt is likely created by an outward-propagating warp induced by an off-plane inner planet (Mouillet et al. 1997; Augereau et al. 2001), later identified as being the imaged β Pic b planet (Lagrange et al. 2010, 2012). HST images of the beta Pic disk Golimowski et al. (2006) give the impression of two separate disks due to the fortuitous alignment of the line of nodes and the line of sight. The NZ Lup system could have a morphology similar to β Pictoris, but with the line of nodes nearly perpendicular to the line of sight. The mutual inclinations of the belts are of the same order, approximately 4− 5°. However, in the present case, we believe that a single undetected off-plane planet cannot account for all the observed characteristics of the disk. A planet located between the inner and outer parts of the disk (at around 90− 100 au) would for example carve out a gap or a density drop (e.g., Lazzoni et al. 2018) and launch both inward- and outward-propagating waves that would warp the disk inside and outside its orbit. However, we would expect the inner and outer warp to propagate at roughly the same speed, which would tend to align both of them towards the same plane. One way to alleviate this problem would be for the potential perturbing planet to be placed inside the inner edge of the disk (≤ 85 au). In this case, the warp would propagate outwards and would first reach the inner regions of the disk before affecting its outer parts, thus creating a de facto tilt between these two regions. However, such a planet would not be able to create a density gap or discontinuity in the middle of the disk.

An additional characteristic of the NZ Lup disk that the planet-induced warp scenario cannot explain is its very small vertical thickness, with h ≲ 0.02. This is at odds with the theoretical prediction that the warped regions should be puffed up to an opening angle roughly equal to twice the inclination of the planet (Nesvold & Kuchner 2015), meaning that this opening angle should be comparable to the off-plane tilt of the warp. Indeed, in the present case, this would lead to an opening angle of ~ 5°, which is at least a factor four higher than the value constrained by our parametric fitting.

Our modeling results suggest that the region inwards of 80 au is cleared of dust. If one assumes that planetesimals exist at all radii in the initial stage of planet formation and that gaps or cavities in debris disks are created by the emerging planets, one can assess the number of planets and their minimum masses necessary to create such a gap between 10 and 80 au during the lifetime of the system (i.e., 16 Myr). To do that, we use the numerical work of Shannon et al. (2016) to constrain the lower limit in mass for planets to have depleted the inner regions. To clear a gap between 10 and 80 au over 16 Myr, we find that two planets of > 0.65 M_Jup are needed. Shannon et al. (2016) assume equal-mass planets separated by 20 mutual Hill radii (i.e., the typical separation between planets in Kepler multi-planet systems). Using the planet mass upper limits from this paper together with the lower limit that we have just derived, we find that the potential planet able to warp the inner cold disc would have a mass of between 0.65 and 2 M_Jup.

Lee & Chiang (2016) investigated a variety of debris disk morphologies considering only a few parameters, the viewing angles, and a planet eccentricity orbiting inside a parent-belt of planetesimals. In particular, their model is able to produce a double wing geometry (also referred to as a moth geometry), if the planet eccentricity is large enough (~0.7), the apoastron of dust orbits is towards the observer, and the system is nearly edge-on (5–10°). Depending on whether the tail of dust particles is directed towards or away from the observer, the double wing geometry can morph into a bar geometry. This picture can describe rather well the case of HD 32297 where a faint bar with an offset from the main ring is seen from the HST STIS image (Schneider et al. 2014). However, the bar is nearly parallel to the midplane, while the two wings deviate from each other with the stellocentric distance. These behaviors do not match the image of NZ Lup. The disk spine in Fig. 3 shows that the second component is converging towards the midplane, which compounds the hypothesis of two separated rings.

Some multiple belts observed in a few debris disks (HD 131835 or HD 141569) are also suspected to be the result of the photoelectric effect generated by a gas–disk interaction (Lyra & Kuchner 2013), although inclinations of such belts are expected to be uniform. However, no gas detection has been reported so far for NZ Lup and therefore a planet scenario would be more likely. In any case, the two aforementioned debris disks feature several thin belts as expected for scenarios involving gas, while NZ Lup has only two. It is therefore not easy to find a straightforward explanation for the peculiar characteristics of this system. More sophisticated dynamical and numerical explorations should be carried out in the future in order to address this issue as well as a systematic search for planets.

9 Summary

Here we summarize the results of our analysis on the NZ Lup debris disk.

• The angular resolution and contrast achieved with SPHERE reveal a globally very inclined debris disk (~85°) with a disk splitting attributed to the presence of two noncoplanar planetesimal belts. The sizes of these belts are derived from modeling with GRaTer. Assuming some variations of the density between the belts, we determine the plausible range of stellocentric distances: 80–95 au for the inner belt, and 95–120 au for the outer belt.

• The modeling favors a configuration where the outer belt is more inclined than the inner belt by ~5°. The reversed geometry (inner belt is more inclined) is significantly worse in terms of matching the data and is therefore ruled out.

• The relative intensity between these two belts necessarily implies a variation of density with a discontinuity, which naturally generates a gap in scattered light. Models with a true dust depletion (or density gap) between the belts more closely match the images, although the difference is marginal with other model families where the region in between the belts is not completely depleted.

• The fine structure revealed by direct imaging is not measurable in the SED. However, the scarcity of photometric data in the far-IR supports the argument for a broad component, which is compatible with the position of the belts inferred from direct imaging. A closer warm dust component as proposed from a former analysis of the SED cannot be confirmed either from the reanalysis of the SED or from SPHERE observations.

• No co-moving candidates could be identified in the SPHERE data for all three available epochs, but the limit of detection reaches a lower limit of 1 M_Jup at a separation of 0.6″ (projected separation of ~4 au).

• The scattered light of the disk is also identified for the first time in the visible using HST archival data. Some characteristics can be recovered as compared to near-IR imaging with SPHERE but the data lacks contrast at short angular separations to confirm the presence of the double belt.

• Explaining the distribution of the dust in the form of two belts with a single planet has some theoretical shortcomings. Multiple planets could be required but their masses would be much lower than the limit of detection.

The debris disk around the G-type star NZ Lup features a rare case of multiple belts observed in a planetary system. Only four such systems have been identified with direct imaging, at least two of which are gas-rich. Therefore, this system is crucial in the context of understanding the last stages of planetary system evolution. It should be a prime target for future high-contrast facilities. The identification of the double belt pattern in SPHERE data was made difficult due to the high inclination of the disk and the ADI-induced artifacts. We are planning a follow-up investigation of this target at shorter wavelengths, with ZIMPOL for instance, to take advantage of a higher angular resolution, as well as in polarimetry (both near-IR and visible) to provide diversity in the phase function.

Acknowledgements

SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF-Osservatorio di Padova (Italy), Observatoire de Genéve (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004-2008), grant number 226604 for FP7 (2009-2012) and grant number 312430 for FP7 (2013-2016). French co-authors also acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU in France. This work has also been supported by a grant from the French Labex OSUG@2020 (Investissements d’avenir - ANR10 LABX56). The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). Italian co-authors acknowledge support from the “Progetti Premiali” funding scheme of the Italian Ministry of Education, University, and Research. This work has been supported by the project PRIN-INAF 2016 The Cradle of Life - GENESIS- SKA (General Conditions in Early Planetary Systems for the rise of life with SKA). C.P. acknowledges support from the ICM (Iniciativa Científica Milenio) via the Núcleo Milenio de Formación Planetaria grant. Finally, this work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice) and Observatoire de Paris/LESIA (Paris). We thank P. Delorme and E. Lagadec (SPHERE Data Centre) for their efficient help during the data reduction process.

Appendix A Signal-to-noise-ratio maps

Fig. A.1 Signal-to-noise-ratio map per resolution element (indicated in the color bar). The disk midplane is aligned with the horizontal direction and the displayed field of view is 6″ × 3″ for IRDIS and 2.4″ × 1.2″ for the IFS.

Appendix B Frequencies of model parameters

Fig. B.1 Frequencies of parameters obtained for the 1% best models in the single-belt mask1 case.

Fig. B.2 Frequencies of parameters obtained for the 1% best models in the two-belt gap case.

Fig. B.3 Frequencies of parameters obtained for the 1% best models in the two-belt fill case.

Fig. B.4 Frequencies of parameters obtained for the 1% best models in the two-belt ring case.

Appendix C Inclination effect on the second belt

Fig. C.1 Two-belt models assuming the outer component is inclined by 87° while the inner belt is (from top to bottom) inclined by 82, 80, 78 and 76°. The relative intensity evolves according to the mutual inclination.

Appendix D Comparison of data and best model images

Fig. D.1 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the corresponding forward model, and the best model injected in the data, and the residuals (gap case: i1 = 82°, i2 = 87°, r1 = 90 au, r2 = 115 au). The field of view is 6″ × 0.8″.

Fig. D.2 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the corresponding forward model, and the best model injected in the data, and the residuals (fill case: i1 = 83°, i2 = 87°, r1 = 90 au, r2 = 105 au). The field of view is 6″ × 0.8″.

Fig. D.3 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the best model injected in the data, and the residuals (ring case: i1 = 86°, i2 = 82°, r1 = 900 au, r2 = 90 au). The field of view is 6″ × 0.8″.

Appendix E Characterization of point sources in the IRDIS field of view

Fig. E.1 Positions of all identified points sources in the IRDIS field of view.

Fig. E.2 Positions of all points sources in the IRDIS field of view as measured at two epochs (2016-05-31 and 2017-04-30) compared to the expected positions if they were background stars.

Appendix F Limits of detection for DUSTY and COND atmosphere models

Fig. F.1 As in Fig. 9 but for the DUSTY (left) and COND (right) model.

References

All Tables

All Figures

	Fig. 1 Color-magnitude diagram of NZ Lup. Overplotted isochrones for 5-12-20-30 and 70 Myr from Bressan et al. (2012).
In the text

Fig. 2 Images of the disk obtained with IRDIS in the H band (top and bottom left), and with the IFS in the YJ band in a narrower field of view (bottom right). The intensity scale is arbitrary for each epoch to optimize the disk visibility. The disk midplane is aligned with the horizontal direction and the displayed field of view is 6″ × 3″ for IRDIS and 2.4″ × 1.2″ for the IFS. The main features reported in the text are annotated with arrows and ellipse.

In the text

	Fig. 3 Spine of the disk measured with a single component assumption (black line) and for two components (red and blue lines).
In the text

	Fig. 4 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the best forward model (i = 85°, r = 95 au, αin = 10, αout = − 5, g = 0.6, h = 0.02) and the residuals. The field of view is 6″ × 0.8″.
In the text

	Fig. 5 Examples of the density (blue) and intensity (red) functions corresponding to the four model families.
In the text

	Fig. 6 Frequencies of parameters obtained for the 1% best models in the two-belt 2gaps case.
In the text

	Fig. 7 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the corresponding forward model, the best model injected in the data, and the residuals (2gaps case: i1 = 82°, i2 = 87°, r1 = 85 au, i2 = 115 au). The field of view is 6″ × 0.8″.
In the text

	Fig. 8 Flux density as a function of wavelength. The disk lies between 85 and 115 AU.
In the text

	Fig. 9 Limits of detection in contrast (top) and converted into Jovian masses (bottom) for two age assumptions (20 and 30 Myr), using the BT-Settl atmosphere model.
In the text

	Fig. A.1 Signal-to-noise-ratio map per resolution element (indicated in the color bar). The disk midplane is aligned with the horizontal direction and the displayed field of view is 6″ × 3″ for IRDIS and 2.4″ × 1.2″ for the IFS.
In the text

	Fig. B.1 Frequencies of parameters obtained for the 1% best models in the single-belt mask1 case.
In the text

	Fig. B.2 Frequencies of parameters obtained for the 1% best models in the two-belt gap case.
In the text

	Fig. B.3 Frequencies of parameters obtained for the 1% best models in the two-belt fill case.
In the text

	Fig. B.4 Frequencies of parameters obtained for the 1% best models in the two-belt ring case.
In the text

	Fig. C.1 Two-belt models assuming the outer component is inclined by 87° while the inner belt is (from top to bottom) inclined by 82, 80, 78 and 76°. The relative intensity evolves according to the mutual inclination.
In the text

	Fig. D.1 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the corresponding forward model, and the best model injected in the data, and the residuals (gap case: i1 = 82°, i2 = 87°, r1 = 90 au, r2 = 115 au). The field of view is 6″ × 0.8″.
In the text

	Fig. D.2 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the corresponding forward model, and the best model injected in the data, and the residuals (fill case: i1 = 83°, i2 = 87°, r1 = 90 au, r2 = 105 au). The field of view is 6″ × 0.8″.
In the text

	Fig. D.3 From top to bottom: original data from April 2017 (same as in Fig. 2), the raw model before PSF convolution, the best model injected in the data, and the residuals (ring case: i1 = 86°, i2 = 82°, r1 = 900 au, r2 = 90 au). The field of view is 6″ × 0.8″.
In the text

	Fig. E.1 Positions of all identified points sources in the IRDIS field of view.
In the text

	Fig. E.2 Positions of all points sources in the IRDIS field of view as measured at two epochs (2016-05-31 and 2017-04-30) compared to the expected positions if they were background stars.
In the text

	Fig. F.1 As in Fig. 9 but for the DUSTY (left) and COND (right) model.
In the text