Polarimetric differential imaging with VLT/NACO A comprehensive PDI pipeline for NACO data (PIPPIN)

Context. The observed diversity of exoplanets can possibly be traced back to the planet formation processes. Planet–disk interactions induce sub-structures in the circumstellar disk that can be revealed via scattered light observations. However, a high-contrast imaging technique such as polarimetric di ff erential imaging (PDI) must first be applied to suppress the stellar di ff raction halo. Aims. In this work we present the PDI PiPelIne for NACO data (PIPPIN), which reduces the archival polarimetric observations made with the NACO instrument at the Very Large Telescope. Prior to this work, such a comprehensive pipeline to reduce polarimetric NACO data did not exist. We identify a total of 243 datasets of 57 potentially young stellar objects observed before NACO’s decommissioning. Methods. The PIPPIN pipeline applies various levels of instrumental polarisation correction and is capable of reducing multiple observing setups, including half-wave plate or de-rotator usage and wire-grid observations. A novel template-matching method is applied to assess the detection significance of polarised signals in the reduced data. Results. In 22 of the 57 observed targets, we detect polarised light resulting from a scattering of circumstellar dust. The detections exhibit a collection of known sub-structures, including rings, gaps, spirals, shadows, and in-or outflows of material. Since NACO was equipped with a near-infrared wavefront sensor, it made unique polarimetric observations of a number of embedded protostars. This is the first time detections of the Class I objects Elia 2-21 and YLW 16A have been published. Alongside the outlined PIPPIN pipeline, we publish an archive of the reduced data products, thereby improving the accessibility of these data for future studies.


Introduction
Over 5 500 exoplanets 1 have been discovered to date, and they are extremely diverse in terms of their masses, compositions, and rings or cavities, are expected byproducts of planet formation and are indeed associated with the protoplanet-hosting PDS 70 (Keppler et al. 2018(Keppler et al. , 2019;;Haffert et al. 2019) and AB Aur systems (Currie et al. 2022), although the evidence for AB Aur b was recently disputed by Zhou et al. (2023).Multi-wavelength observations trace different disk regions, including the large, millimetre-sized dust grains near the midplane (e.g.ALMA Partnership et al. 2015) at longer wavelengths.Scattered light can be captured from the upper surfaces of the disk at optical and nearinfrared (NIR) wavelengths and provides information about the material through the measurements of phase functions and the degree of polarised light.Since central stars are observed close to the peak of blackbody emission, a high-contrast imaging technique is employed to reveal the faint structures in their immediate vicinities.Polarimetric differential imaging (PDI; Gledhill et al. 1991Gledhill et al. , 2001;;Kuhn et al. 2001) is especially well suited to observing the optical and NIR scattered light of a circumstellar disk.Unpolarised stellar light becomes polarised after being scattered by circumstellar dust grains, and PDI can be used to remove the stellar component, revealing the fainter polarised light structures below the diffraction halo of the star.
Several instruments, including the Subaru High-Contrast Coronagraphic Imager for Adaptive Optics (HiCIAO; Hodapp et al. 2008;Suzuki et al. 2010), the Gemini South Gemini Planet Imager (GPI; Macintosh et al. 2006Macintosh et al. , 2014)), the Very Large Telescope (VLT) Nasmyth Adaptive optics system COude (NACO) NIR camera (Lenzen et al. 2003;Rousset et al. 2003), and the VLT Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE; Beuzit et al. 2019), have exploited the PDI technique to observe a large number of young stellar objects (YSOs).These instruments utilise a polarised beam-splitter to separate the incoming light into two beams with orthogonal linear polarisations.The instrumental point spread function (PSF) is unchanged for both beams, as they are recorded simultaneously.The high contrast (∼ 10 −2 -10 −4 ; Avenhaus et al. 2018) between the faint scattered light disk and the bright stellar halo can be suppressed by subtracting measurements of the two orthogonal polarisation states.In particular, PDI is an effective imaging technique for circumstellar disks with low inclinations (e.g.HD 169142, i ≈ 13 • and TW Hya, i ≈ 7 • ; Hales et al. 2006;Apai et al. 2004;van Boekel et al. 2017), for which angular differential imaging (Marois et al. 2006) leads to the selfsubtraction of the face-on disk's signal.
Observations that use PDI have revealed a large number of disks with different sizes, surface brightnesses, and morphologies in scattered light.Scattered light observations trace the upper layers of a circumstellar disk since the micron-sized dust grains found there are optically thick at optical and NIR wavelengths.Hence, circumstellar disks must be flared to intercept the stellar radiation at large distances (Chiang & Goldreich 1997;de Boer et al. 2016;Ginski et al. 2016).Transition disks with large dust-depleted inner cavities are frequently detected (e.g.Mayama et al. 2012;Canovas et al. 2013;Keppler et al. 2018;Maucó et al. 2020), and the observed circumstellar disks commonly show rings at varying radii (e.g.Quanz et al. 2013;Muro-Arena et al. 2018;Avenhaus et al. 2018).Additionally, spiral features are frequently detected in scattered light (see Fig. 9 of Benisty et al. 2023).The gas perturbations, coupled to the small grains that are traced in scattered light, are suggested to emerge from interactions with a companion or with the environment.Furthermore, the combination with sub-millimetre observations can reveal dust filtering at pressure maxima (e.g.Garufi et al. 2013;Maucó et al. 2020) and help identify fragmentation, possibly resulting from gravitational instability (Weber et al. 2023).
In scattered light imaging, the misalignment of an (un-resolved) inner disk can cast a shadow onto the outer disk (Bohn et al. 2022).Depending on the magnitude of the misalignment, narrow shadow lanes (e.g.HD 100453; Benisty et al. 2017) or wideangle obscurations can appear (e.g.HD 143006; Benisty et al. 2018).In the case of stellar multiplicity, the geometry of the circumstellar environment can be assessed further by interpreting which stellar component is responsible for the dust illumination (Weber et al. 2023;Zurlo et al. 2023).Depending on the size, composition, and porosity of the small dust grains, different scattering phase functions can be measured (Shen et al. 2009;Tazaki et al. 2016Tazaki et al. , 2019)).By studying the dust properties in circumstellar disks, we can assess the efficiency of dust growth based on the size, composition, and porosity of the grains involved.PDI observations are not limited to Class II disks (Lada 1987): secondgeneration dust disks, or debris disks, are also observed with facilities such as VLT/SPHERE (e.g.HIP 79977; Engler et al. 2017, HR 4796A;Milli et al. 2019), Gemini South/GPI (e.g.HD 157587; Millar-Blanchaer et al. 2016), and Subaru/HiCIAO (e.g.HD 32297; Asensio- Torres et al. 2016).
For the most studied Class II disks (Lada 1987), the observed sub-structures are frequently explained by invoking the presence of planetary companions (e.g.ALMA Partnership et al. 2015;van der Marel et al. 2019;Long et al. 2019;Asensio-Torres et al. 2021).The existence of sub-structures suggests that planet formation is already underway and began when the YSOs were still embedded in their natal envelopes, during the Class 0 or I phases (t < 10 6 yr; Garufi et al. 2022a).Furthermore, measurements of the dust masses of Class II disks appear incompatible with predicted planet formation efficiencies and the masses of exoplanetary systems (Manara et al. 2018;Mulders et al. 2021).The higher dust masses of Class 0 and I disks determined by Tychoniec et al. (2020) could indicate that giant planet formation commences before the protostellar envelope has dissipated (Cridland et al. 2022;Miotello et al. 2023).Alternatively, the accretion of material from the surrounding cloud can continually replenish the mass of the protoplanetary disk.The total mass budget available for planet formation therefore exceeds the disk mass at any given time (Manara et al. 2018;Garufi et al. 2022a).The two explanations put forward to solve the missing mass problem demonstrate the important role of embedded Class 0 and I objects in the formation of planets.However, the earliest YSOs are particularly difficult to observe at optical wavelengths due to their embedded nature.As a consequence, the optical wavefront sensors (WFSs) of most modern extreme adaptive optics (AO) systems do not allow for an adequate AO correction of deeply embedded YSOs.The NIR AO188 system, part of SCExAO on the Subaru telescope, is an exception as it provides AO for polarimetric imaging in the northern hemisphere.However, embedded protostars in the south were only observable to some older, ground-based instruments that were equipped with an infrared WFS.For completeness, the retired NICMOS instrument on the Hubble Space Telescope measured the polarised light of the earliest YSOs (e.g.Silber et al. 2000;Kóspál et al. 2008;Perrin et al. 2009), though JWST is not equipped with polarimetric capabilities.
In this work, we present a re-reduction of polarimetric archival data from NACO, the Nasmyth Adaptive Optics System (NAOS), which together with the COude Near-Infrared CAmera (CONICA) forms the NACO instrument at the VLT (Lenzen et al. 2003;Rousset et al. 2003).Initially installed at the Nasmyth B focus of UT4 in 2001, NACO was reinstalled at the Nasmyth A focus of UT1 from 2014 until its decommissioning in 2019.NACO operated at wavelengths between 1 and 5 µm, and NAOS was equipped with a visible (0.45-1.0 µm) and infrared (0.8-2.5 µm) WFS, enabling observations of embedded YSOs despite their faint optical magnitudes.NACO was equipped with a Wollaston prism (and also wire grids; see Sect.2.2.4) to perform polarimetric observations, and a half-wave plate (HWP) was installed in 2003.In Sect. 2 we describe how our PDI PiPelIne for NACO data (PIPPIN)2 reduces the NACO polarimetric data with the PDI technique.Section 3 outlines a broad inspection of the PIPPIN-reduced data, and we present a novel method for assessing the detection significance of a polarised signal.Section 4 compares the reduced data with those of the SPHERE instrument.The conclusions are summarised in Sect.5, and the reduced data archive has been published online3 .

Selection of polarimetric observations
Since the polarimetric mode of NACO was not solely used to observe YSOs, we made a selection of observations of interest to this study.First, the European Southern Observatory (ESO) archive was searched for every polarimetric NACO observation classified as the SCIENCE data type.Using the object identifier and the astroquery Python package (Ginsburg et al. 2019), we searched the SIMBAD archive (Wenger et al. 2000) to select any object that was ever classified as one of the following categories: (candidate) Orion variable, (candidate) Herbig Ae/Be star, (candidate) T Tauri star, or a (candidate) YSO.However, a large number of observations have unclear object identifiers.In these instances, astroquery was utilised to locate the object closest to the target right ascension (RA) and declination (Dec.)coordinates.In total, we find 57 candidate Class 0 -III objects that potentially exhibit polarised light from circumstellar material.As these systems were observed in multiple filters, epochs, or with different instrument setups, we find a total of 243 datasets.Table A.1 lists the objects of interest and information on the observation setup for each dataset.

PDI PiPelIne for NACO data (PIPPIN)
A general pipeline to reduce NACO data is provided by ESO4 .However, this pipeline cannot reduce the polarimetric observations and thus previous works utilised custom, self-written pipelines (e.g.Apai et al. 2004;Quanz et al. 2011;Canovas et al. 2013).The different data reduction methods could lead to inconsistent scientific results.For instance, one of the rings of HD 97048 observed by Ginski et al. (2016) was not recovered from the same data in the earlier analysis of Quanz et al. (2012).Such discrepancies can be avoided by using a single, comprehensive pipeline.In this section, we describe the operation of our PIPPIN pipeline, which applies the PDI technique to polarimetric NACO observations.With the exception of an instrumental Mueller matrix model, PIPPIN largely follows the polarimetric data reduction outlined in de Boer et al. (2020).For a more detailed characterisation of the instrumental polarisation of NACO, we refer to de Boer et al. (2014) andMillar-Blanchaer et al. (2020).

FLATs, bad-pixel masks, and DARKs
To correct for any variations of the detector's gain, PIPPIN performs a FLAT-fielding of the SCIENCE images.In general, internal lamp FLATs were taken for each filter and detector (i.e.S13, S27, L27, S54, and L54) that were used during the night.The polarimetric mask, which prevents the ordinary and extraordinary beams from overlapping, is also inserted when measuring the FLAT fields.The FLATs are DARK-subtracted and subsequently normalised by being divided by the median counts.The bad-pixel masks are generated by assessing which pixels had a non-linear response in the FLAT fields.The linearity of the pixel response is determined by comparing the FLATs observed with the internal lamp switched on (FLAT on ) to FLATs made with the lamp turned off (FLAT off ).The factor by which the pixel-counts are expected to increase is computed by dividing the median of FLAT on by the median of FLAT off .Pixels were flagged when their response deviated by more than 2σ from the expected increase.Similar to the FLAT fields, the bad-pixel masks are computed for each filter and detector used throughout the night.

Pre-processing
The PIPPIN pipeline can be described in two parts: the preprocessing and the application of PDI.The pre-processing commences by reading parameters from a configuration file that allows users to customise the data reduction.The configuration file must be located in the same directory as the SCIENCE observations; otherwise, the pipeline creates a default file.Table B.1 outlines the parameter keywords in the configuration file along with the recognised values, descriptions, and default values.After reading the configuration parameters, PIPPIN groups observations by the utilised detector, window-size, observing ID (if requested), filter, exposure time, HWP usage, and whether the Wollaston prism or wire grids were used.Each observation is DARK-subtracted and FLAT-normalised by division.The pixels flagged in the bad-pixel mask are replaced by the median counts of the surrounding square of 5 × 5 pixels.
To retrieve the approximate positions of the ordinary and extra-ordinary beams, PIPPIN applies a minimum-filter with a specific kernel-shape to the images.The filter consists of two squares of 3 × 3 pixels that are offset by the approximate separation of the beams, which in turn depends on the pixel scale of the utilised detector.The maximum in the filtered image yields the approximate location of the ordinary and extra-ordinary beams.This method avoids any persisting bad pixels or image artefacts such as the polarimetric mask.Subsequently, the initial guesses are used to retrieve more accurate PSF locations via a userspecified fitting method.For each beam, PIPPIN can employ a single 2D Moffat function or subtract two Moffat functions from each other to reproduce the flat top of a saturated PSF.Alternatively, the pipeline can use a maximum-counts method for the asymmetric PSFs that are encountered in the case of deeply embedded stars.
The sky-subtraction is performed by subtracting two dithering positions or by subtracting the median per row of pixels.To avoid contamination from the target, a region around the fitted beam centres is excluded in the median sky-subtraction method.This region is defined with the sky_subtraction_min_offset parameter in the configuration file.Moreover, this parameter ensures that the two dithering positions are sufficiently offset to perform a sky-subtraction.In addition, horizontal gradients are removed by a linear fit that excludes the region around the beams.The linear fit is applied to the average of five rows of pixels and a 2D Gaussian filter with σ = 5 pixels is applied to smooth out the resulting background approximation.Some observations show a distinct horizontal pattern, which can be removed by fitting each row of pixels individually and without applying a Gaussian filter.
Next, the ordinary and extra-ordinary beams are cut out of the images by a user-specified crop-size.The maximum counts of the beams are evaluated with an iterative sigma-clipping to determine which observations suffered from a poor AO correction.

Polarimetric differential imaging
Polarised light can be described with the Stokes formalism and the Stokes vector: where I is the total intensity, Q and U are intensities of the linear polarisation components and V describes the circular polarisation.As NACO was not primarily designed for polarimetry, the observations suffer from instrumental polarisation (IP) and crosstalk effects.Reflections within the instrument can introduce polarised signal whose magnitude depends on the instrument configuration, altitude of the target object, etc.Furthermore, crosstalk between the linear and circular polarisation components reduces the polarimetric efficiency (Witzel et al. 2011).Hence, PIPPIN employs a multi-stage correction for these effects.A first-order correction for different transmission efficiencies is to impose that the stellar flux in the ordinary (I ord ) and extra-ordinary (I ext ) beams are the same, as described in Appendix C of Avenhaus et al. (2014a).Since the PSF core is often saturated in NACO observations, PIPPIN draws multiple, userspecified annuli and computes the total fluxes within them.For each annulus i, the ratio between the fluxes, is used to scale the ordinary and extra-ordinary images as I ord / X ord/ext,i and I ext X ord/ext,i , respectively.This method implicitly assumes that the total flux in annulus i is unpolarised, thereby ignoring any polarisation induced by the interstellar medium or any intrinsic polarisation originating from an unresolved inner disk, for example.We note that this correction could overcompensate for a true disk signal if the disk is not axisymmetric and if its scattered light comprises a considerable fraction of the stellar signal.
If the HWP has a rotation angle of θ = 0 • , the ordinary beam (I ord ) measures light polarised in the +Q direction and the extraordinary beam (I ext ) measures the perpendicularly polarised light in the −Q direction, both in the HWP reference frame.The I Q and Q components are found by addition and subtraction of the equalised beam intensities: Measurements of the U component are made by rotating the incoming beam by 45 • , which means that the HWP is rotated by θ = 22.5 • .The I U and U components are calculated with The top panels of Fig. 2 show the resulting median Stokes Q and U images for HD 135344B.The position angle (PA) is −35 • , so that the sky is rotated anti-clockwise to the axes of the detector as is indicated by the compasses in the figure.In the Q image, the positive signal aligns with the Y-axis and the negative signal aligns with the X-axis.The U image displays a similar butterfly pattern, but rotated by 45 • since it measures different components of the disk.
Instrumental polarisation introduced downstream of the HWP can be removed by recording the −Q and −U parameters at θ = 45 • and 67.5 • , respectively.The instrumental Q IP and U IP components are unaffected by this rotation of the HWP and contribute in the same manner as before: Using the double-difference method (Hinkley et al. 2009;Bagnulo et al. 2009), we can subtract the IP components: Similarly, the IP-corrected intensities are found with the doublesum: The total intensity is calculated as The second row of Fig. 2 shows the median Q and U images resulting from the double-difference method.Due to the IP removal, the butterfly patterns show more distinct features than the Q + and U + images and the recorded noise outside of the disk is reduced.
An additional correction is made for the IP introduced upstream of the HWP, following the method outlined in Canovas et al. ( 2011) and de Boer et al. ( 2020).The correction is performed for each HWP cycle to mitigate temporal differences in the IP as a result of changing angles of reflection.As before, it is assumed that the stellar light is unpolarised and polarised signal near the star is ascribed to instrumental polarisation (Quanz et al. 2011).The median Q/I signal is computed over the same annulus i from Eq. 2 to obtain a scalar c Q .To obtain c U , we calculate the median U/I signal over the same annulus.Per annulus, the IP-subtracted linear Stokes components are found by subtracting the product of these scalars and the respective I Q or I U image: By using multiple user-specified annuli, the pipeline retrieves various IP-subtracted results.The third row of panels in Fig. 2 displays the median Q IPS and U IPS images where the annulus was drawn between a radius of 3 and 6 pixels.As expected from the correction, the Q IPS measurement shows a decreased signal near the star compared to the Q image.In Fig. 2, the U IPS signal is lower than Q IPS as a result of crosstalk between the linear and circular Stokes components (Witzel et al. 2011).If a disk is unmistakably detected and approximately axisymmetric, this reduced efficiency of the Stokes U component relative to Q can be estimated following the method outlined by Avenhaus et al. (2014a).In an annulus with disk signal, the number of pixels where |Q IPS | > |U IPS | is expected to be equal to the number of pixels where We can multiply the U IPS image by a factor of 1/e U so that the above assumption holds.The crosstalk-corrected components are then where we assume an efficiency of 100% for Stokes Q.By modelling the NACO IP with standard star observations, Millar-Blanchaer et al. (2020) conclude that the Stokes Q has an efficiency of ∼ 90%.Since such a correction is not performed with PIPPIN, any quantitative polarimetry measurements on the reduced data products could be off by ∼ 10%.The efficiencycorrection should not be performed in instances with ambiguous signal and thus PIPPIN only makes the crosstalk-correction if requested.
Ord./ext.-equalisedIncomplete HWP cycles, with only measurements of Q ± (or U ± ), are removed.If only the Stokes Q + and U + (or only Q − and U − ) were recorded, PIPPIN will still be able to produce the final data products, but the double-difference method cannot be applied.At this point, the pipeline computes the median Q, U, I Q , I U , and I over all observations.The final polarisation images (PI, Q ϕ , and U ϕ ) are described below in terms of Q and U, but we note that these data products are also calculated with Q IPS /U IPS and Q CTC /U CTC , if possible.The total polarised inten-sity is calculated as This method of squaring Q and U can lead to the increase in noise in regions where the Q or U signal originating from the disk is low.A cleaner image can be found with the azimuthal Stokes parameters outlined in Monnier et al. (2019) where ϕ is the azimuthal angle and is calculated for each pixel with where (x star , y star ) are the pixel-coordinates of the central star.If the disk has a low inclination and the scattered light emerges from single scattering events, the polarisation is oriented azimuthally with respect to the star.Consequently, the Q ϕ image shows a positive signal as it measures polarisation angles of ±90 • .Simultaneously, the U ϕ image is expected to show a negligible signal as it measures polarisation angles of ±45 Otherwise, the offset angle ϕ 0 is set to 0 • .The median total intensity I, polarised intensity PI, and azimuthal Stokes parameters Q ϕ and U ϕ with different levels of IP correction are shown in Fig. 3 for HD 135344B.Once PIP-PIN has computed the final data products, these images are derotated using scipy.ndimage.rotate.Therefore, contrary to Fig. 2, the panels of Fig. 3 have north pointing up and east to the left.It is apparent from the total and polarised intensity images that the PDI technique applies an extremely effective suppression of the stellar signal, thus revealing the circumstellar disk and its spiral arms.In this example, we observe the U ϕ signal diminish as the IP corrections are performed.Since HD 135344B is observed at a low inclination and axisymmetric to a first order, we employed crosstalk-correction and U ϕ minimisation to produce the final Stokes images.For these Ks-band observations, we find a reduced efficiency of e U = 0.65, in agreement with Garufi et al. (2013) who find an efficiency of 0.61 (Avenhaus et al. 2014a).Similarly, the more extensive IP model presented by Millar-Blanchaer et al. (2020) results in an efficiency of e U = 0.7 ± 0.02 for Elia 2-25.Furthermore, we find an offset angle of ϕ 0 = 5.3 • , while 3.7 • is derived in the previous analysis of these data (Garufi et al. 2013;Avenhaus et al. 2014a).

Non-HWP and wire-grid observations
Prior to August 8, 2003, NACO was not equipped with a HWP.Rather than rotating the HWP to modulate the direction of polarisation, observers would alter the PA by rotating the instrument on its rotator ring.PIPPIN automatically diagnoses whether the de-rotator flange of the telescope support structure was used (Lenzen et al. 2003).For these data, the HWP angles θ = 0, 22.5, 45, and 67.5 • in Eqs. 7, 8, 9, and 10 are replaced by the PAs of the instrument: θ PA = 0, 45, 90, and 135 • .The Q ± , U ± , I Q ± , and I U ± images are also de-rotated to align the circumstellar structures before combining them with Eqs.11, 12, 13, and 14.For the rotator observations, the IP subtraction of Eqs.16 and 17 is also performed.
In the early stages of its operation, NACO was equipped with wire grids to carry out polarimetric observations, rather than the Wollaston prism.In our cross-validation of the ESO archive, we found four potentially young sources that were observed in this manner: V1647 Ori, NX Pup, Mon R2 IRS 3, and R Mon.PIP-PIN adopts the Pol_00, Pol_45, Pol_90, and Pol_135 wire grids as measurements of the Stokes Q + , U − , Q − , and U + components, respectively.The linear Stokes components are propagated in the presence of the HWP.The only beam that is present in the images is fit with a single Moffat function.Since the wire-grid observations are not limited by the height of the polarimetric mask, their final data products have a much larger field of view than those obtained with the Wollaston prism.

Supplemental data products
Since the disk is illuminated by the star, the scattered light brightness decreases by the inverse of the squared distance to the host star.To better visualise structures at larger separations from the star, PIPPIN also produces images that are multiplied by the squared, de-projected radius.The disk PA, PA disk , is used to calculate the offsets along the major axis, ∆x disk , and minor axis, ∆y disk , with where ∆(R.A.) and ∆(Dec.)are the right ascension and declination offsets with respect to the star.Subsequently, the deprojected radius r is computed with where i disk is the disk inclination.As is shown in Table B.1, the disk PA, PA disk and inclination i disk are specified in the configuration file for PIPPIN and are set to 0 • by default.In cases where the disk inclination and PAs are unknown, the default values ensure that the images are scaled by the projected separation from the host star.
The height of the final data products is limited to ∼ 3.0 arcsec (using the S27 detector) due to the polarimetric mask.Observations where the PA was rotated, rather than the HWP, cover a larger area of the sky.Since the sky rotates while the polarimetric mask remains stationary, the effective field of view is increased.Figure C.1 depicts this increased sky coverage.An eight-pointed star emerges where at least one Q and one U component are covered and thus the polarised intensity can be computed within this shape.An inner octagon appears where every positive and negative Stokes component is observed.We note that the signalto-noise decreases for areas outside of this octagon, due to the reduced number of observations.PIPPIN outputs the extended eight-pointed star images in addition to the data products resulting from the double-difference method, which are restricted to the inner octagon that has a complete coverage.

Identification of detections
The PIPPIN pipeline described above was used to reduce all observations listed in Table A.1.The table lists multi-epoch, multiwavelength observations as well as different exposure times and whether the wire grids were used or the Wollaston prism, with the HWP or PA.For each set of observations, we indicate the (non-)detection of circumstellar material in the final data products.The detection significance of polarised signal was assessed via a template-matching method, akin to cross-correlation, applied to the Stokes Q and U images.In the case of a detection, we expect that the signal is present in multiple, adjacent pixels and forms a specific butterfly pattern.Synthetic Q synth and U synth templates of the expected butterfly patterns were constructed with where ϕ is the azimuthal angle calculated with Eq. 23 and PA is the position angle of the observation, which is subtracted since PIPPIN de-rotates the final data products, including the Q IPS and U IPS images.Subsequently, the Q synth and U synth templates were divided into multiple annuli with increasing radius and width of 2 pixels, roughly corresponding to one resolution element in the H (42 mas) and Ks band (56 mas) at a pixel scale of 27 mas pixel −1 .Figure 4 shows an example of the Q synth and U synth templates and a single annulus for a PA of −35 • , corresponding to the observations of HD 135344B.The values in the templates range from −1 to +1 and pixels outside of the annulus are set to 0, thus ensuring that they do not contribute when calculating the cross-correlation coefficient.In annulus i, a crosscorrelation coefficient is calculated for the Q and U signals: where the sum is performed over every pixel within annulus i.
In this manner, a positive pixel increases the coefficient if the respective quadrant expects a positive signal.A negative signal in the negative quadrants of the template also contributes, whereas a discrepant signal reduces the cross-correlation coefficient.A cross-correlation function (CCF) is constructed by computing a coefficient for each annulus.For the narrowband NB_1.64 observations of HD 135344B, the rightmost panel of Fig. 5 displays the CCFs for the Q IPS and U IPS components in blue and red, respectively.The CCF was converted into a signal-to-noise (S/N) function by subtracting the mean coefficient between 35 and 50 pixels, and subsequently dividing by the standard deviation of coefficients within that same range, indicated by the grey shaded region in the right panel.The annulus-wise CCFs peak at a radius of 8 pixels with signal-to-noises of S/N ∼ 13 and ∼ 15, respectively for Q IPS and U IPS .These maxima surpass our 5σ detection threshold, thereby identifying this observation as a detection.Although the template-matching method generally worked well, it failed to flag two observations of HR 4796 as detections, despite the polarised signal evident from a visual inspection.These non-detections can be ascribed to the high inclination and narrow features of HR 4796, while the outlined template-matching analysis works optimally for face-on disks.
In this reduction of the NACO data, many of the nondetections are likely the result of small or faint disks, or the absence of polarised light.Notably, IM Lup, GQ Lup, and EX Lup do not show polarised light in the NACO data, despite their prominent detections with SPHERE/IRDIS (Avenhaus et al. 2018;van Holstein et al. 2021;Rigliaco et al. 2020).The data of IM Lup and GQ Lup were not previously published, whereas Kóspál et al. (2011)

Analysis of detected polarised light
As demonstrated in Table A.1, in 22 out of the 57 observed systems, we find at least one set of observations with polarised signal originating from circumstellar material.Figure 6 presents a gallery of these detections and highlights a diverse collection of morphologies.As mentioned in Sect.2, HD 135344B shows distinct spiral arms in its circumstellar disk while HD 142527 has spiral features in its eastern and western lobes.Furthermore, we detect rings in a large number of disks, including HD 169142, HD 163296, HD 97048, HR 4796, TW Hya, HD 142527, and Sz 91.HR 4796 is the only debris disk in our sample and its Q ϕ image shows a narrow ring.The disks around HD 163296, HD 97048, and HD 100546 are offset from the central star, suggesting that the scattering surface is above the disk midplane as confirmed by Monnier et al. (2017), Ginski et al. (2016), andSissa et al. (2018).The highly extended disk of HD 142527 shows a large inner cavity that is possibly cleared out by an inner companion (Biller et al. 2012;Close et al. 2014;Lacour et al. 2016;Claudi et al. 2019), undetected in polarised light.Moreover, we find narrow shadow lanes imprinted on the disks of HD 142527 and SU Aur, similar to Avenhaus et al. (2017) and Ginski et al. (2021).In these cases, misaligned inner disks prevent the stellar light from reaching certain areas of the outer disk.CR Cha, MP Mus, AK Sco, and Elia 2-25 show negligible structure due to their small sizes, but a significant butterfly pattern was detected, leading to their inclusion in Fig. 6.As a possible consequence of their small sizes, the polarimetric NACO data of CR Cha and MP Mus were previously unpublished.
Figure 6 displays a number of objects with extended features that appear inconsistent with circumstellar disks.SU Aur shows tail-like features, where Ginski et al. (2021)  presents the median I Q , which also includes unpolarised (scattered) light.Despite this, filamentary structure is unambiguously detected.As discussed in Sect.2.2.4, the two detections utilising polarimetric wire grids, R Mon and Mon R2 IRS 3, have much larger image sizes than the regular data products (i.e. compared to the upper rows of panels in Fig. 6).

Disk classification and brightness
Circumstellar polarised light is detected in nine out of the 14 Herbig stars in our sample.Grouping the Orion variables together with T Tauri stars, we detect polarised signal in eight out of 28 systems.For the YSOs, five of the 14 sources are flagged as detections with the template-matching analysis.The only highproper motion star in our sample, HR 4796, also constitutes the only detection of a debris disk.However, our selection of young stars based on the SIMBAD object type could have missed some of the older Class III disks.Since NACO frequently observed disks known to be extended and thus potentially observable in scattered light, the gathered sample is certainly not unbiased.For that reason, a statistical analysis of the disk occurrence per object type is somewhat arbitrary.
Here we examine the disk brightness of the sample of circumstellar disks detected by NACO.Since the disk inclination, disk extent, stellar brightness, and the distance to the source affect the total disk brightness, we made use of the polarised-tostellar light contrast δ pol (Garufi et al. 2017(Garufi et al. , 2022b;;Benisty et al. 2023).The polarised flux per unit area, F pol , was multiplied by the squared separation 4πr 2 to account for the reduced stellar illumination.Subsequently, we normalised it by the stellar flux F * and averaged radially along the disk's major axis.Thus, the polarised-to-stellar light contrast was computed as where r in and r out are the inner-and outermost radii, respectively.This measurement expresses the fraction of stellar photons that became polarised as a result of scattering by the resolved disk.
The contrast δ pol is set by the geometry and composition of the circumstellar disk and is sometimes referred to as the geometrical albedo.Following the method outlined in Appendix B of Garufi et al. (2017), we performed a 2-pixel-wide cut (one resolution element) of the Q ϕ image along the major axis of the disk.The photons measured along the major axis are scattered with angles close to 90 • .This cut reduces the impact of the disk inclination on its brightness due to any asymmetry in the scattering phase function.The PA of the major axis varies for the sources observed by NACO, but was set to 0 • when this axis could not be estimated.These ambiguous sources (e.g.CR Cha) were roughly azimuthally symmetric and therefore did not significantly affect the derived contrast.The inner-and outermost radii (r in , r out ) were determined by eye for each system detected in scattered light.The disk inclination and scale height were not taken into account when re-scaling the polarised flux by the projected separation.Similar to Garufi et al. (2017), we calculated the primary error of δ pol by deriving the standard deviation in a resolution element of the Q ϕ image.Subsequently, this noise estimate for each pixel was propagated through Eq. 31 to find σ δ pol .Fig. 7 displays the polarised-to-stellar light contrast δ pol against the apparent J-band magnitude m J , measured as part of the 2 Micron All Sky Survey (2MASS; Cutri et al. 2003).The source names are shown along the top axes.Blue, orange, and red markers indicate whether the observation was performed in the H, Ks, or L' band.Diamonds (T Tau), circles (Herbig), and squares (debris) mark the object types.We note that HD 135344B is shown as a Herbig star (circle), in line with Garufi et al. (2014), but contrary to the SIMBAD object type in Table A.1.Similarly, Parsamian 21 is depicted as an embedded YSO despite its Orion variable type reported in Table A.1.The δ pol values of saturated PSFs are presented as upper limits in Fig. 7 (triangles; 99.75-th percentile) because the stellar flux F * is underestimated.In some instances, the source was also observed with narrowband filters where the stellar PSF was not saturated due to the smaller filter width.Hence, we could estimate the broadband flux F BB * , using the narrowband flux F NB * , if the two filters had overlapping wavelength ranges.We calculated the stellar flux as where t BB exp and t NB exp are the exposure times of the broadband and narrowband observations, respectively.We integrated over the corresponding transmission curves T (λ) to estimate how many photons should be detected for each photon in the narrowband filter.For the H band, we used the NB_1.64 and NB_1.75 narrowband filters.For the Ks band, NB_2.17,IB_2.18,NB_2.12,NB_2.15, and IB_2.21 were employed to compute the broad-band flux and the NB_3.74 filter was used for saturated L'-band observations.
Since NACO was equipped with a NIR WFS, it could observe sources down to K ≈ 14 mag (Rousset et al. 2003).For comparison, SPHERE's optical WFS has a magnitude limit of R ≈ 14 mag (Beuzit et al. 2019), GPI has a limit of I ≈ 10 mag (Macintosh et al. 2014), and SCExAO/CHARIS on the Subaru telescope is limited by R ≈ 13 mag5 .For that reason, NACO could achieve a unique insight into embedded protostars, despite their faint optical magnitudes.The grey shaded region in Fig. 7 roughly indicates the sources that are inaccessible by the optical WFS installed on SPHERE.Since the estimated J-band magnitude limit of SPHERE depends on the spectral type of the observed source, the limit of m J ∼ 10 mag should be viewed as a crude assessment.From Fig. 7 we find that the four embedded protostars Parsamian 21, Elia 2-21, Elia 2-29, and YLW 16A, in addition to the low-mass star Sz 91, are likely not observable with modern PDI instruments.
The polarised-to-stellar light contrasts δ pol are listed in Table A.1.We note that the δ pol values are possibly underestimated by ∼ 10% (see Sect. 2.2.3) due to the absence of a correction for the reduced Q efficiency resulting from crosstalk.For the sources with available mass estimates (also included in Table A.1), we fail to detect a trend between δ pol and the stellar mass M * .Since the disk's dust mass is related to the stellar mass (e.g.Pascucci et al. 2016), the absence of a distinct trend reveals that the disk brightness in polarised light is not strongly correlated with the abundance of dust in the system.Instead, the scattered light brightness is affected by the amount of light that is intercepted.The geometry of the system primarily influences the polarised-to-stellar light contrast δ pol , with the dust composition acting as a secondary effect.Similarly, Garufi et al. (2022b)  For the extended systems SU Aur, R CrA, and Z CMa, we assessed the polarised-to-stellar light contrast ratio of a potential disk component at close separations, meaning that r out was limited to 25, 25, and 53 pixels, respectively.We find δ pol ∼ 1.5 • 10 −3 for SU Aur, ≲ 3 • 10 −2 for R CrA, and ≲ 6 • 10 −3 for Z CMa.The brightest disk is found around HR 4796 with δ pol ∼ 0.3 − 0.4.This finding is somewhat surprising, given that it is a flat debris disk and therefore should not intercept much stellar light.However, the high brightness is also reported in previous works where it is argued that the scattering phase functions are consistent with large (∼ 20 µm) aggregate dust particles composed of small monomers (Milli et al. 2017(Milli et al. , 2019)).For Sz 91, the lowest-mass star (M = 0.58 ± 0.07M ⊙ ; Maucó et al. 2020) where polarised light is detected, we determine upper limits of δ pol ≲ 8 • 10 −2 and ≲ 4 • 10 −2 in the Ks and H band, respectively.The estimated contrasts of multi-wavelength observations do not show any clear discrepancies between Hand Ks-band observations, owing to relatively large uncertainties.Hence, it is difficult to draw any conclusions about the dust composition by evaluating the disk colour.

Discussion
The morphologies shown in Fig. 6 are almost identical to the polarised intensity images presented in previous publications of these data (see Table A.1 for references).The data reduction performed by PIPPIN therefore appears to reproduce the results obtained by the custom pipelines in other works.
To study the performance between different instruments, Fig. 8 presents a comparison between the NACO and modern SPHERE observations of SU Aur (programme ID: 1104.C-0415(E), PI: Ginski), HD 142527 (programme ID: 099.C-0601(A), PI: Avenhaus), and TW Hya (programme ID: 095.C-0273(D), PI: Beuzit).In this comparison, we find the results of the different instrument characteristics.For instance, the NACO observations of TW Hya were made under better seeing conditions (∼ 0.5 arcsec) than those made by SPHERE (∼ 0.7 arcsec), but we find that the NACO polarised signal displays residual speckles over the circumstellar disk.The SPHERE Q ϕ image does not show similar artefacts, likely due to the superior AO instrument.As part of the NACO instrument, NAOS had fewer actuators (185 active actuators for NAOS against 1377 for SAXO; Blanco et al. 2022) shaping the deformable mirror and its optical WFS operated at a lower frequency (1200 Hz versus 444 Hz; Fusco et al. 2006;Rousset et al. 2003), thus resulting in typical H-band Strehl ratios of ∼ 10 -35% as opposed to ∼ 60 -80% for SPHERE observations (Fusco et al. 2014;van Boekel et al. 2017).Furthermore, the SPHERE NIR camera (IRDIS) has a pixel scale of ∼ 12 mas pixel −1 (Maire et al. 2018) while the most-used S27 detector on CONICA had ∼ 27 mas pixel −1 .The NACO instrument was not equipped with a coronagraph in its polarimetric mode and thus short exposure times were utilised to avoid saturation by the central star.Each of the NACO observations presented in Fig. 8 employed considerably shorter singleframe integration times than the respective SPHERE observations (SU Aur: 0.35 versus 32 s, HD 142527: 0.3454 versus 16 s, TW Hya: 5 versus 16 s), thereby inevitably reducing the signalto-noise.Ginski et al. (2021) trace the extended western structure of SU Aur out to ∼ 6 arcsec, whereas the NACO data only confidently show signal out to ∼ 4 arcsec.Moreover, the filamentary structure observed in the tails and disk of SU Aur (Ginski et al. 2021) are not resolved in the NACO data due to the reduced signal-to-noise.Lastly, the polarimetric mask of NACO limits the vertical extent of the final data products to ∼ 3 arcsec.Hence, the north-western spiral structure of HD 142527 is eventually obscured in the NACO data.

Conclusions
We have presented a complete catalogue of polarimetric NACO images for YSOs, reduced in a homogeneous manner with a new pipeline that employs the PDI technique.Via a crossexamination with the object types reported on SIMBAD, 57 targets were identified as potentially young systems with polarimetric NACO observations.As a result of multi-epoch and multiwavelength observations, a total of 243 datasets were reduced with the publicly available PIPPIN pipeline6 .PIPPIN can handle observations made with NACO's HWP as well as its de-rotator.In addition to the Wollaston prism, observations measured with wire grids can be reduced too.Various levels of corrections for instrumental polarisation are performed, depending on the type of observation.
The reduced data products were analysed with a templatematching method to evaluate the detection significance.This technique exploits the butterfly pattern in the Stokes Q and U images that should be present in the case of significant polarised light.We find that 22 out of the 57 observed systems revealed polarised light in at least one observation.These detections revealed a wide diversity of sub-structures, including rings, gaps, spirals, shadows, and in-or outflowing matter.Since NACO was equipped with a NIR WFS, unique polarimetric observations of embedded YSOs were made.To our knowledge, this is the first work to publish the reduced data products of the Class I protostars Elia 2-21 and YLW 16A.PIPPIN also revealed detections of polarised light in the L' band for HD 100546 (Avenhaus et al. 2014b) and Elia 2-21.This long-wavelength filter (3.8 µm) is not available on current, state-of-the-art PDI instruments such as SPHERE/IRDIS, SCExAO/CHARIS, or GPI.
Alongside this article, we publish an archive of the reduced data products generated with PIPPIN on Zenodo7 .As these observations were made in the past two decades, their combination with modern scattered light observations can be used to identify temporal changes in the sub-structures of planet-forming disks.
In turn, such morphological changes can be used to infer the presence of a perturbing companion (Ren et al. 2020).Recent studies of the NACO data of HD 97048 and SU Aur (Ginski et al. 2016(Ginski et al. , 2021) ) have led to the discovery of previously unidentified features.With this work, we hope to have outlined the utility of NACO observations reduced with the PDI technique.servatoire de Genève (Switzerland), ETH Zürich (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)(2014)(2015)(2016).This research was performed with the Python programming language.In particular, the SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), and astropy (Astropy Collaboration et al. 2013Collaboration et al. , 2018) ) packages were utilised.Herbig Ae/Be B5 (35)   3.02 ± 0. (c) Datasets indicated with † were observed after April 11, 2018, when the HWP rotation mechanism failed (Millar-Blanchaer et al. 2020).After its repair, the motor encoder position no longer corresponded to the same polarisation angle.PIPPIN is currently not equipped to correct for this systematic offset in the observed polarisation angle, and results from these datasets should therefore not be trusted.

Figure 1 Fig. 1 .
Fig. 1.Open AO-loop assessment of HD 135344B Ks-band observations.Left panel: Maximum counts of the ordinary (red) and extraordinary (blue) beams.The horizontal dashed blue and red lines are the 3σ bounds for the respective beam, indicating which observations have adequate AO corrections.The upper-right panel shows an example of an effective AO correction for the ordinary beam, and the lower-right panel shows the blurred result of an open AO loop.

Fig. 2 .
Fig. 2. Median Stokes Q and U images with different levels of IP corrections for HD 135344B Ks-band observations.From top to bottom: Q + and U + components after equalising the ordinary and extraordinary fluxes, Q and U resulting from the double-difference method, Q IPS and U IPS after subtracting the median IP within an annulus, and the crosstalk-corrected Q CTC and U CTC components where the reduced Stokes U efficiency is accounted for.The characteristic butterfly pattern is visible in each panel, and the compasses show the orientation of the detector and the sky.

Fig. 3 .
Fig. 3. Final PIPPIN data products with different levels of IP correction.From left to right: Median total intensity (I), polarised intensity (PI), and the azimuthal Stokes components Q ϕ and U ϕ of HD 135344B observed in the Ks band.From top to bottom: Equalised ordinary and extra-ordinary beams, IP-subtracted, crosstalk-corrected, and U ϕ -minimised results.The total intensity is shown with a logarithmic scale from 20 to 10 4 counts, whereas the other panels use a linear scale from −5 to +5 counts and a logarithmic scale up to ±90.The negative signal is depicted in orange, and in each image north points up and east to the left.
also report a non-detection of polarised light in the EX Lup observations.

Fig. 4 .
Fig. 4. Templates for observations of HD 135344B with a PA of −35 • .Top panels: Complete Q synth and U synth templates.Values range from −1 to +1, and pixels outside of the image are set to 0. Bottom panels: Example annuli used in computing the cross-correlation coefficients.

Fig. 5 .
Fig. 5. Detection analysis of HD 135344B observed in the narrowband NB_1.64.Left panels: Q IPS and U IPS images divided by the white lines into the four quadrants of the expected butterfly pattern.Right panel: Annulus-wise CCFs, with the S/N shown against the annulus radius in pixels.The results for the Q IPS and U IPS images are plotted in blue and red, respectively.The shaded region specifies the coefficients used in normalising and converting the CCF into a S/N function.The 5σ detection limit is indicated with a horizontal line.

Fig. 6 .
Fig. 6.Gallery of young systems detected with NACO and reduced with PIPPIN.Each panel shows the polarised light on a logarithmic scale ranging between different values to highlight sub-structures.The highest degree of IP correction is used where possible.Scale bars in the lowerleft corners of each panel indicate 100 AU at each object's distance.HD 169142, R CrA, and Parsamian 21 are shown in the H band, MP Mus is shown in the IB_2.06 filter, and the other panels use Ks-band observations.Mon R2 IRS 3 shows the median I Q image because the Stokes U component was not observed.The images of YLW 16A and Elia 2-21 present the first polarised light detections in the NACO observations.

Fig. 7 .
Fig.7.Polarised-to-stellar light contrast, δ pol , plotted against the apparent J-band magnitude.The right panel shows a zoomed-in view of the bright m J .The object names are listed along the top axes.The marker colours and symbols specify the observing filter and object type, respectively.Upper limits are shown when the stellar PSF was determined to be saturated.The error bars show the 3σ uncertainties.The grey shaded region shows the approximate magnitudes (m J ≳ 10) inaccessible by the SPHERE AO system.

Fig. 8 .
Fig. 8.Comparison between PIPPIN-reduced NACO Q ϕ observations (left panels) and the more recent SPHERE data (right panels).From top to bottom: SU Aur, HD 142527, and TW Hya observed in the H band with both instruments.The SPHERE observations were previously published by Ginski et al. (2021), Hunziker et al. (2021), and van Boekel et al. (2017) for SU Aur, HD 142527, and TW Hya, respectively.
The abbreviations of the SIMBAD object types are: Orion Var. for Orion variable stars; Herbig Ae/Be for Herbig Ae stars; T Tau for T Tauri stars; Herbig Ae/Be for Herbig Be stars; High-PM for high-proper motion stars; and YSO for young stellar objects.Abbreviations followed by a question mark are candidate object types and those listed in parentheses show previous identifications.(b)Datasets where the cross-correlation of Sect.3.1 could not be applied, due to incomplete coverage of both Q and U, present a hyphen (-) in the 'Detection' column.Instances where the cross-correlation analysis resulted in a non-detection despite clear signs of polarised light from a visual inspection are appended with '(Yes)'.
Avenhaus et al. (2014a)U ϕ signal can occur if there is crosstalk between Q and U, if the light is scattered multiple times(Canovas et al. 2015b), if the disk has a high inclination, and if an inadequate correction retains stellar or instrumental polarisation(Hunziker et al. 2021).If requested, PIPPIN can minimise the U ϕ signal in the same annulus used for the crosstalk-correction by fitting for the azimuth angle offset ϕ 0 , similar toAvenhaus et al. (2014a).