© ESO 2021

1 Introduction

Understanding the formation and evolution of young, self-luminous exoplanets and brown dwarf companions is one of the main goals of high-contrast imaging at near-infrared wavelengths (e.g., Nielsen et al. 2019; Vigan et al. 2021). Only a few of these directly imaged substellar companions have been detected close to the parent star and within a circumstellar disk (e.g., Lagrange et al. 2010; Keppler et al. 2018; Haffert et al. 2019); most companions are found at much larger separations (≳ 100 au; see e.g., Bowler 2016). Close-in planets and companions are generally believed to form through core accretion (Pollack et al. 1996; Alibert et al. 2005) or gravitational instabilities in the circumstellar disk (Cameron 1978; Boss 1997). Companionsat larger separations may form through direct collapse in the molecular cloud (Bate 2009) or disk gravitational instabilities at an early stage (Kratter et al. 2010). Alternatively, companions may form close to the star and subsequently scatter to wide orbits through dynamical encounters with other companions (e.g., Veras et al. 2009).

In all formation scenarios, the companion is generally expected to form its own circumsubstellar accretion disk (e.g., Stamatellos & Whitworth 2009; Szulágyi et al. 2017). Indeed, a handful of substellar companions show evidence for the presence of an accretion disk through hydrogen emission lines, red near-infrared colors, and excess emission at mid-infraredwavelengths (e.g., Seifahrt et al. 2007; Bowler et al. 2011; Bailey et al. 2013; Kraus et al. 2014; Zhou et al. 2014; Haffert et al. 2019). Interestingly, whereas ALMA and other radio interferometers have been successful at detectingthe dust and gas of disks around isolated substellar objects (e.g., Ricci et al. 2014; van der Plas et al. 2016; Bayo et al. 2017), attempts to detect such disks around substellar companions have almost exclusively yielded non-detections (Bowler et al. 2015; MacGregor et al. 2017; Wu et al. 2017a,b, 2020; Wolff et al. 2017; Ricci et al. 2017; Pérez et al. 2019). The only detection of a disk around a substellar companion at mm-wavelengths is that of PDS 70 c with ALMA by Isella et al. (2019). ALMA has also detected a disk around FW Tau C (Kraus et al. 2015; Caceres et al. 2015), but, from models of the Keplerian rotation of the gas, the companion appears to be a ~ 0.1M_⊙ star (Wu & Sheehan 2017; Mora et al. 2020). To explain their non-detections, Wu et al. (2017a, 2020) suggest that the disks around substellar companions must be very compact (≲1000 R_Jup or ≲ 0.5 au) and optically thick to be able to sustain several million years of accretion. Alternatively, there might be a dearth of large dust grains in circumsubstellar disks because the observed mid-infrared excess could also be explained by a gaseous disk with small micron-sized dust grains.

Although compact circumsubstellar disks cannot be spatially resolved with current 8-m class telescopes, they can create a measurable, integrated linear polarization at near-infrared wavelengths (Stolker et al. 2017). The polarization can be introduced through scattering of the companion’s thermal photons by dust within the disk, (partial) obscuration of the companion’s atmosphere by the disk, or self-scattering in the case of a high-temperature disk. In all cases, the disk must have a nonzero inclination because the polarization of a face-on viewed, rotationally symmetric disk integrates to zero and a low-inclination disk cannot obscure the companion’s atmosphere. Measuring polarization originating from circumsubstellar disks enables us to study the structure and physical properties of the disks.

Planets and brown dwarf companions without a disk can also be linearly polarized at near-infrared wavelengths. Late-M- to mid-L-type dwarfs are expected to have dusty atmospheres because their temperatures are sufficiently low for refractory material to condense (Allard et al. 2001; Ackerman & Marley 2001). This atmospheric dust scatters the thermal radiation emanating from within the object, linearly polarizing the light. Whereas the spatially integrated polarization signal of a spherical, horizontally homogeneous dusty atmosphere is zero, a net polarization remains when this symmetry is broken (Sengupta & Krishan 2001). Examples of these asymmetries are rotation-induced oblateness and an inhomogeneous distribution of atmospheric dust clouds (Sengupta & Marley 2010; de Kok et al. 2011; Marley & Sengupta 2011; Stolker et al. 2017), or even a large transiting moon (Sengupta & Marley 2016). Based on the models, the degree of linear polarization due to circumsubstellar disks and atmospheric asymmetries can be several tenths of a percent up to several percent in favorable cases.

Spatially unresolved polarimetric observations have already been used to study disks around pre-main sequence stars (e.g., Rostopchina et al. 1997; Bouvier et al. 1999; Grinin 2000; Ménard et al. 2003). In addition, optical and near-infrared polarization has been detected for dozens of field brown dwarfs (Ménard et al. 2002; Zapatero Osorio et al. 2005, 2011; Tata et al. 2009; Miles-Páez et al. 2013, 2017). In most cases, the polarization of these brown dwarfs is interpreted as being caused by rotation-induced oblateness or circumsubstellar disks, whereas an inhomogeneous cloud distribution has appeared harder to prove. However, Millar-Blanchaer et al. (2020) recently measured the near-infrared polarization of the two L/T transition dwarfs of the Luhman 16 system and found evidence for banded clouds on the hotter, late-L-type object.

With the adaptive-optics-fed high-contrast imaging instruments Gemini Planet Imager (GPI; Macintosh et al. 2014) and SPHERE-IRDIS (Beuzit et al. 2019; Dohlen et al. 2008) at the Very Large Telescope (VLT), we now have access to the spatial resolution and sensitivity required to measure the near-infrared polarization of substellar companions at small separations. After correction for instrumental polarization effects, the polarimetric modes of both instruments can reach absolute polarimetric accuracies of ≲ 0.1% in the degree of polarization (Wiktorowicz et al. 2014; Millar-Blanchaer et al. 2016; van Holstein et al. 2020). Early attempts to measure the polarization of substellar companions by Millar-Blanchaer et al. (2015) and Jensen-Clem et al. (2016) with GPI and by van Holstein et al. (2017) with SPHERE-IRDIS have been unsuccessful. Nevertheless, van Holstein et al. (2017) showed that SPHERE-IRDIS can achieve a polarimetric sensitivity close to the photon noise limit at angular separations > 0.5′′. Ginski et al. (2018) detected a companion to CS Cha using SPHERE-IRDIS and measured the companion’s polarization to be 14%, suggesting that it is surrounded by a highly inclined and vertically extended disk. However, recent optical spectroscopic observationswith MUSE show that the companion is not substellar in nature, but is a mid M-type star that is obscured by its disk (Haffert et al. 2020).

In this paper, we present the results of a survey of 20 planetary and brown dwarf companions with SPHERE-IRDIS, aiming to detect linear polarization originating from both circumsubstellar disks and atmospheric asymmetries. Our study is complemented by a similar survey of seven companions using GPI and SPHERE by Jensen-Clem et al. (2020). The outline of this paper isas follows. In Sect. 2 we present the sample of companions and the observations. Subsequently, we describe the data reduction in Sect. 3 and explain the extraction of the polarization signals in Sect. 4. In Sect. 5 we discuss our detections of polarization and the upper limits on the polarization for the non-detections. In the same section, we present images of five circumstellar disks that we detected in our survey. Because the most plausible explanation for the polarization of the companions is the presence of circumsubstellar disks, we perform radiative transfer modeling of a representative example of such a disk in Sect. 6. Finally, we discuss the implications of our measurements in Sect. 7 and present conclusions in Sect. 8.

2 Target sample and observations

2.1 Target sample

The sample of this study consists of 20 known directly imaged planetary and brown dwarf companions, out of the approximately 140 such companions that are currently known¹. Because the expected polarization of the companions is around a few tenths of a percent or less, our primary selection criterion was whether SPHERE-IRDIS can reach a high signal-to-noise ratio (S/N) in total intensity without requiring an excessive amount of observing time. Therefore, the selected companions are relatively bright, are at a moderate companion-to-star contrast, are at a large angular separation from the star, and/or have a bright star for good adaptive-optics (AO) performance (see van Holstein et al. 2017). Our sample contains the majority of the approximately two dozen known companions that match these requirements. Three of the remaining companions have been observed by Jensen-Clem et al. (2020) in their survey of seven companions. An overview of the properties of the companions of our sample is shown in Fig. 1, with the full details presented in Table 1². The sample is diverse, with the companions spanning spectral types from T5.5 to M7, masses between approximately 6 and 70 M_Jup, and ages between approximately 2 Myr and 11 Gyr. The companions orbit stars of spectral types A5 to M1. Six companions show evidence of hosting a circumsubstellar disk, mostly in the form of red near-infrared colors, excess emission at mid-infrared wavelengths, and hydrogen emission lines that reveal ongoing accretion. As can be seen particularly well from Fig. 1, the overall sample ranges from young, hot, accreting companions with spectral types between late M and early L, to old, cold, and massive companions of later spectral types. For the six companions that show evidence of hosting a circumsubstellar disk, we expect any polarization to be primarily due to this (spatially unresolved) disk, whereas for the other companions polarization would most likely be due to an inhomogeneous cloud distribution or rotation-induced oblateness.

Fig. 1 Properties of the companions of our sample showing the age, spectral type, mass (surface area of data points), effective temperature, and possible existence of a circumsubstellar disk (thick border). The data points of HR 8799 b and c, and of HR 8799 d and e, overlap.

2.2 Observations

All our observations were performed with the dual-beam polarimetric imaging (DPI) mode of SPHERE-IRDIS (de Boer et al. 2020; van Holstein et al. 2020). In this mode, linear polarizers are inserted in the left and right optical channels of IRDIS to simultaneously create images of the two orthogonal linear polarization states on the detector. A rotatable half-wave plate (HWP) modulates the incident linear polarization with switch angles 0°, 45°, 22.5°, and 65.5° (a HWP cycle) to measure Stokes Q and U. The observationswere carried out between October 10, 2016, and February 16, 2020, under generally good to excellent atmospheric conditions. An overview of the observations is shown in Table 2.

The observation strategy was as follows. We generally observed each target multiple times with typically over 30 min of on-source exposure time per visit. However, for some targets a single visit was enough to detect the companion with high S/N in total intensity. We mainly observed in broadband H, but sometimes used broadband J or K_s when we wanted to obtain data in an additional filter or in the case the companion was brighter in K_s than H. We used the apodized Lyot coronagraph with a mask diameter of 185 mas (for J and H) or 240 mas (for K_s) to suppress the starlight (Carbillet et al. 2011; Guerri et al. 2011). This allowed us to use longer integration times per frame to minimize the effects of read noise. However, we did not use integration times longer than 64 s to limit the effect of changing atmosphericconditions during a HWP cycle. In addition to the polarimetric science frames, we took star center frames to accurately determine the position of the star behind the coronagraph and star flux frames to measure the total stellar flux. We also took sky frames with the same instrument setup as the science and star flux frames to subtract the sky background from the respective frames.

For the majority of the observations, we used the pupil-tracking mode (van Holstein et al. 2017). In this mode the image derotator (K-mirror) rotates such that the telescope pupil is kept fixed with respect to the detector while the on-sky field of view rotates with the parallactic angle. The pupil-tracking mode has numerous advantages. With sufficient parallactic rotation we can apply angular differential imaging (ADI; Marois et al. 2006) to suppress speckle noise and accurately determine the total intensity of the companions located at small angular separations from the star. Furthermore, because the speckles are quasistatic, they are more effectively removed in the polarimetric data-reduction steps (and can be further suppressed by applying ADI to the polarimetricimages). In addition, the diffraction spikes created by the support structure of the telescope’s secondary mirror are suppressed by a mask added to the Lyot stop. Finally, the loss of signal due to the crosstalk produced by the image derotator is limited (see van Holstein et al. 2020). As a result, the polarimetric efficiency, that is, the fraction of the linearly polarized light incident on the telescope that is actually measured, is always high (typically ≳ 90%).

For a few targets, we used the field-tracking mode to be able to offset the derotator position angle and control the orientation of the image on the detector. For instance, the companions of AB Pic and HD 106906 are at such large angular separations (see Table 1) that we needed to place them in one of the corners of the 11′′ × 11′′ field of view to make them visible. In the case of 1RXS J1609 and DH Tau we switched to field-tracking mode after we discovered that both companions crossed a cluster of bad pixels during the pupil-tracking observations. In all cases, we chose the orientation of the image derotator such that the polarimetric efficiency was high (see de Boer et al. 2020).

3 Data reduction

We reduced the data with the publicly available and highly automated pipeline IRDAP³ (IRDIS Data reduction for Accurate Polarimetry), version 1.2.2 (van Holstein et al. 2020). IRDAP preprocesses the raw data by subtracting the sky background, flat fielding, correcting for bad pixels, extracting the images of IRDIS’ left and right optical channels, and centering using the star center frames. It then subtracts the right images from the left images (the single difference) for each of the measurements taken at HWP switch angles equal to 0°, 45°, 22.5°, and 67.5° to obtain the Q⁺-, Q⁻ -, U⁺ -, and U⁻ -images, respectively. IRDAP also adds these same left and right images (the single sum) to obtain the total-intensity $I_{Q^{+}}$-, $I_{Q^{-}}$-, $I_{U^{+}}$-, and $I_{U^{-}}$-images. Subsequently, IRDAP computes cubes of Q- and U-images from the double difference and the corresponding cubes of total-intensity I_Q - and I_U -images from the double sum, as: $$\begin{array}{ll}Q\hfill & =\frac{1}{2}\left(Q^{+}-Q^{-}\right),\hfill \\ I_Q\hfill & =\frac{1}{2}\left(I_{Q^{+}}+I_{Q^{-}}\right),\hfill \end{array}$$

and similar for U and I_U . For the two data sets of HD 4747 and the dataset of PZ Tel in J-band, strongly varying atmospheric seeing prevents the double difference from fully removing the signal created by transmission differences between the two orthogonal polarization directions downstream of the image derotator. To remove this spurious polarization, we used the normalized double difference (see van Holstein et al. 2020) instead of the conventional double difference for these three data sets.

After computing the double difference and double sum, IRDAP uses a fully validated Mueller matrix model to correct for the instrumental polarization (created upstream of the image derotator) and crosstalk of the telescope and instrument with an absolute polarimetric accuracy of ≲ 0.1% in the degree of polarization. IRDAP also derotates the images and corrects them for true north (see Maire et al. 2016b). This results in a total of four images: Q, U, I_Q , and I_U , that constitute our best estimate of the linear polarization state incident on the telescope. Finally, IRDAP computes images of the linearly polarized intensity PI = √(Q² + U²), and, following the definitions of de Boer et al. (2020), images of Q_ϕ and U_ϕ. Positive (negative) Q_ϕ indicates linear polarization in the azimuthal (radial) direction, and U_ϕ shows the linear polarization at ±45° from these directions. In Sect. 5.5 we use the polarized intensity and Q_ϕ - and U_ϕ -images to show the five circumstellar disks that we detected.

The model-corrected Q- and U-images often contain a halo of polarized light from the star. This polarization can originate from interstellar dust, (unresolved) circumstellar material, and spurious or uncorrected instrumental polarization. With IRDAP we can therefore determine the stellar polarization from the I_Q-, I_U -, and model-corrected Q- and U-images by measuring the flux in these images in a user-defined region that contains only starlight and no signal from a companion, background star, or circumstellar disk. For most data sets we measured the stellar polarization using a star-centered annulus placed over the AO residuals, or in the case that region contains little flux, a large aperture centered on the star. IRDAP then determines the corresponding uncertainty by measuring the stellar polarization for each HWP cycle individually and computing the standard error of the mean over the measurements. Finally, IRDAP creates an additional set of Q- and U-images with the stellar polarization subtracted. To this end, it scales the I_Q- and I_U-images with the measured fractional stellar polarization and subtracts the resulting images from the model-corrected Q- and U-images. Whenever discussing data in this paper, we always mean the reduction without the stellar polarization subtracted, unless explicitly stated.

For the observations taken in pupil-tracking mode, IRDAP additionally performs classical ADI and ADI with principal component analysis (PCA; Soummer et al. 2012; Amara & Quanz 2012) to suppress the stellar speckle halo and detect the companions in total intensity. IRDAP also processes the star flux frames by performing sky subtraction, flat fielding, bad-pixel correction, and registering through fitting the frames to a 2D Gaussian function. We obtained the final images of the stellar point spread function (PSF) by mean-combining the left and right processed star flux frames and scaling the pixel values to the integration time and system transmission (i.e., due to neutral-density filters) of the science frames. We separately reduced the data sets of targets that we observed multiple times and then used IRDAP to mean-combine the final images produced in each reduction.

The final Q- and U-images of most data sets still contain a small amount of speckle noise close to the star. For the data sets of HR 8799, HD 206893 and β Pic, which have companions at small separations from the star, we therefore performed additional reductions in which we apply classical ADI on the polarimetric images to further suppress these speckles (see van Holstein et al. 2017). To this end, we added a reduction step to IRDAP in which we median-combine the instrumental-polarization-subtracted Q-frames (and U-frames) and subtract the resulting median image from each of the frames before derotating them. In these reductions we skip the later step of determining and subtracting the stellar polarization because the ADI step has already removed the halo of polarized starlight.

4 Extraction of polarization of companions: detection of polarization of DH Tau B

With the data of all targets reduced, we can determine the polarization of the companions, or, in the case we do not detect significant polarization, place upper limits on the degree of polarization of the companions. For this we have developed a method similar to that employed by Jensen-Clem et al. (2020), which, in turn, is based on the method used by Jensen-Clem et al. (2016). In this method, we use aperture photometry to estimate the probability distributions of the companion signals in the I_Q-, I_U -, Q-, and U-images. We then use these distributions to calculate the probability distributions of the degree and angle of linear polarization, from which we retrieve the median values, uncertainties, and upper limits. We applied this method to the data sets of GSC 8047, CD-35 2722, AB Pic, HD 106906, GSC 6214, 1RXS J1609, DH Tau, and TYC 8998. In this section, we demonstrate the method using the 2019-10-24 H-band data set of DH Tau and exemplify the detection of the polarization of DH Tau B, a companion at a large angular separation from its star. For companions at close separations or with large star-to-companion contrasts, we have slightly adapted the method and determine the distributions in I_Q and I_U through ADI with negative PSF injection or fitting of the companion PSF. In Appendices C and D we demonstrate the two respective methods and show how we set upper limits on the polarization of β Pic b and HD 19467 B.

To start the analysis of the 2019-10-24 data set of DH Tau, we determine the center coordinates of the companion DH Tau B by mean-combining the I_Q- and I_U-images and fitting a 2D Moffat function to the resulting image at the position of the companion. We then make a cosmetic correction to the Q- and U-images (if necessary) to remove spurious structures that result from imperfect relative centering of the images, image motion, and parallactic rotation (see Appendix A). The I_Q -, Q-, and U-images (after the cosmetic correction) at the companion position are shown in Fig. 2. The signals in Q and (in particular) in U are clear indications that DH Tau B is polarized.

To determine the probability distributions of the companion signals in I_Q , I_U , Q, and U, we define a range of aperture radii from 1 to 10 pixels to be used for the photometry. For each aperture radius we perform the following five steps, after which we select the final aperture radius to be used for our results. Because at the end of this section we select a final aperture radius of 8 pixels, we use this radius in the examples of the five steps below.

As the first step, we place an aperture of the given radius at the position of the companion in each of the I_Q -, I_U -, Q-, and U-images (see Fig. 2) and sum the flux in the aperture. In the same images we then place a ring of comparison apertures around the star at the same separation as the companion to sample the background. We exclude those apertures that contain the first Airy ring of the companion, diffraction spikes from the star and the companion, and clusters of bad pixels. The resulting ring of apertures for an aperture radius of 8 pixels is shown superimposed on the I_Q -image in Fig. 3. In this figure the first Airy ring and the diffraction spikes created by the Lyot stop mask are clearly visible at the companion position, which is evidence of the extremely good atmospheric conditions during the observations (see Table 2). Finally, we sum the flux in each of the comparison apertures and compute the mean background as the mean of the aperture sums.

In step two, we calculate the probability density function (PDF) of the companion signal in I_Q , I_U , Q, and U, taking into account only the photon noise of the companion. To this end, we compute the companion signals in I_Q , I_U , Q, and U by subtracting the mean background from the summed flux of the companion aperture. We then compute the PDFs of I_Q and I_U from a Gaussian distribution with the mean and variance equal to the respective companion signals, while accounting for the conversion from counts to total number of detected photoelectrons and back to counts (using a detector gain of 1.75 e⁻/count). The resulting PDF of I_Q for an aperture radius of 8 pixels is shown in Fig. 4 (left). For large number of photons, the photon noise in Q and U is the same as that in I_Q and I_U . We therefore construct the PDFs of Q and U from a Gaussian distribution with the mean equal to the companion signals in Q and U, but the variance equal to that of the PDFs of I_Q and I_U . Figure 4 (right) shows the resulting PDF in Q.

For the third step, we estimate the PDF of the background in I_Q , I_U , Q, and U using the comparison aperture sums obtained in the first step. To not a priori assume a specific functional form of the PDF, we use kernel density estimation (KDE). In this method, the PDF is obtained by placing a Gaussian kernel of a given bandwidth (i.e., a Gaussian distribution with a given standard deviation) at each data point of the sample and summing the resulting kernels. We compute the bandwidth of the Gaussian kernel using Scott’s rule (Scott 2015), in this case yielding a bandwidth of ~ 84 counts for I_Q and I_U , and ~ 18 counts for Q and U. Histograms of the background samples and the PDFs as estimated via KDE for an aperture radius of 8 pixels are shown in Fig. 5. We note that for very close-in companions such as PDS 70 b, the number of comparison apertures is low enough that KDE does not produce accurate results. When there are fewer than 21 comparison apertures, we therefore account for the small-sample statistics by fitting the background samples with a Student’s t-distribution with the empirical standard deviation equal to s = s_bg√(1 + 1∕n), with s_bg the standard deviation of the comparison aperture sums and n the number of comparison apertures (see Mawet et al. 2014).

In stepfour, we compute the final probability distributions in I_Q , I_U , Q, and U that includeboth the photon noise of the companion and the uncertainty of the background. For this, we draw 10⁶ random samples from the previously constructed PDFs of the companion signal (step two) and the background (step three). Because we already subtracted the background when computing the PDF of the companion signal, we first subtract the mean background from the drawn background samples. We then compute the final distribution by subtracting the resulting background samples from the samples of the companion signal. Next, we compute the median values of the final distributions and determine the uncertainties from the two-sided 68.27% equal-tailed interval around the median, corresponding to the 1σ (one standard deviation) confidence interval of the Gaussian distribution. The resulting probability distributions for an aperture radius of 8 pixels, including the median values, uncertainties, and S/Ns (i.e., the median value divided by the largest uncertainty), are shown in Fig. 6 (top row). The data are clearly photon-noise limited in Q and U because the distributions are nearly Gaussian and the uncertainties are close to the standard deviation shown in Fig. 4 (right). It follows that we detect DH Tau B with a very high S/N in total intensity and also have significant detections of polarization, especially in Stokes U.

As the fifth and final step, we use the I_Q-, I_U -, Q-, and U-samples to compute the distributions of normalized Stokes q = Q∕I_Q, normalized Stokes u = U∕I_U, the degree of linear polarization P = √(q² + u²), and the angle of linear polarization χ = 1∕2arctan(u∕q). We compute the median values and uncertainties in the same way as we did for I_Q , I_U , Q, and U. The results of these computations for an aperture radius of 8 pixels are shown in Fig. 6 (bottom row).

After performing the five steps above for each defined aperture radius, we plot the median values and uncertainties of q, u, the degree and angle of polarization, and the S/N in q, u, and the degree of polarization as a function of aperture radius in Fig. 7. From this figure we see that, within the uncertainties, the polarization of the companion is constant with changing aperture radius. We select a final aperture radius of 8 pixels, as indicated by the vertical dashed lines in Fig. 7, because at this radius the S/N in q and u is maximized and the aperture is sufficiently large to suppress (average out) the spurious signals resulting from incompletely removed bad pixels (see Appendix B). We conclude that for this 2019-10-24 data set, we measure DH Tau B to have a degree of polarization of 0.51 ± 0.04% and an angle of polarization of 56 ± 3° (east of north) in H-band.

Fig. 2 Reduced IQ-, Q-, and U-images (after applying the cosmetic correction described in Appendix A) at the position of the companion DH Tau B of the 2019-10-24 data set of DH Tau, showing an aperture of radius 8 pixels centered on the companion. The IU -image, which is not shown, is very similar to the IQ-image.

Fig. 3 Reduced IQ-image of the 2019-10-24 data set of DH Tau, showing an aperture of radius 8 pixels at the position of the companion DH Tau B (red) and the ring of comparison apertures of the same radius around the star (white).

Fig. 4 PDF of the signal of DH Tau B in IQ (left) and Q (right) fromthe 2019-10-24 data set of DH Tau, using an aperture radius of 8 pixels and taking into account only the photon noise of the companion. The mean and standard deviation of the distributions are shown above the graphs, with the latter also indicated by the light-blue shaded area.

Fig. 5 Histograms of the background in IQ, IU , Q, and U of the 2019-10-24 data set of DH Tau, as obtained through summing the flux in the 8-pixel-radius comparison apertures of Fig. 3. The mean and standard deviation of the samples are shown above the histograms. The blue curves show the PDFs as estimated through KDE and the red curves show the best-fit Gaussian distributions for comparison.

Fig. 6 Final probability distributions of the signals of DH Tau B in IQ, IU , Q, and U (top row), and in normalized Stokes q and u, degree of linear polarization, and angle of linear polarization (bottom row) from the 2019-10-24 data set of DH Tau, using an aperture radius of 8 pixels. The median values of the distributions, as well as the uncertainties computed from the two-sided 68.27% equal-tailed interval around the median, are shown above the graphs. The S/N, i.e., the median value divided by the largest uncertainty, is shown within parentheses. The 68.27% intervals are also indicated by the light-blue and light-red shaded areas.

Fig. 7 Normalized Stokes parameters q and u and degree of linear polarization (top), angle of linear polarization (center), and S/N in q, u, and the degree of linear polarization (bottom) of DH Tau B as a function of aperture radius for the 2019-10-24 data set of DH Tau.The uncertainties of the measured values are shown with error bars. The final selected aperture radius of 8 pixels is indicated with the dashed vertical lines.

5 Results

After careful analysis of our data with the methods as described in Sect. 4 and Appendices C and D, we detected unresolved polarization originating from DH Tau B and GSC 6214 B. We consider these measurements detections because the measured polarization signals are significant (i.e., have an S/N of at least 5 in q or u) and are very likely intrinsic to the companions (i.e., are not due to interstellar dust). We present these results in Sects. 5.1 and 5.2. We also marginally detected polarization from 1RXS J1609 B, but we show in Sect. 5.3 that this polarization is best explained by interstellar dust. For the other 17 companions we do not detect significant polarization. In Sect. 5.4, we place upper limits on the degree of polarization of 1RXS J1609 B and these other companions. Finally, in Sect. 5.5, we briefly describe five circumstellar disks that we detected in our survey and of which two had not been imaged in polarized scattered light before.

5.1 Detection of intrinsic polarization of DH Tau B

In this section we present the detection of polarization originating from DH Tau B. Table 3 shows the measured H-band degree and angle of polarization of DH Tau B, including the uncertainties and the attained S/Ns, for each of the various data sets and the data set created by mean-combining the final images of the three data sets taken at favorable atmospheric conditions (i.e., the 2019-10-24, 2019-10-25A, and 2019-10-25B data sets; see Table 2). For each data set the measured q- and u-signals are within the uncertainties constant with aperture radius. We determined the final values of the polarization signals using apertures of radius 8 pixels, which is at, or close to, the radius where the S/N in q and u is maximized for the various data sets (see Sect. 4). As shown in Table 3, we detect significant polarization from DH Tau B, reaching S/Ns of around 10 for the three data sets taken at favorable atmospheric conditions. The measured degree and angle of polarization for the different data sets are overall consistent. From visual inspection of the images, we find that the small differences among the data sets are primarily due to small biases caused by incompletely removed bad pixels (see Appendix B). These differences can additionally be caused by time-varying atmospheric conditions and AO performance, the limited accuracy of the Mueller matrix model with which the data have been corrected (see van Holstein et al. 2020), and other unknown systematic effects. From the mean-combined images, we measure DH Tau B to have a degree and angle of polarization of 0.48 ± 0.03% and 58 ± 2° (east of north), respectively, with an S/N of 7.7 in q and 16.1 in u.

Table 3 also lists the stellar degrees and angles of polarization as measured with an annulus at the location of the AO residuals (see Sect. 3). For the mean-combined data set we determined the uncertainty on the stellar polarization by propagating the uncertainties from the individual data sets using a Monte Carlo calculation and assuming Gaussian statistics. The measurements of the stellar polarization are very likely affected by some systematic effects because the signals are less consistent than those of the companion and show differences among the data sets that are much larger than the calculated (statistical) uncertainties. The most likely explanation for these differences is that time-varying atmospheric conditions and AO performance cause the effective coronagraphic extinction to vary from frame to frame. Because the companion is not affected by the coronagraph, this can also explain why the polarization measured for the companion is more consistent among the data sets. The stellar polarization measurements show that the star could be truly polarized because the angles of polarization for the three data sets taken at favorable conditions (2019-10-24, 2019-10-25A, and 2019-10-25B) are quite similar. Importantly, the measured polarization of the companion differs significantly from that of the star in all data sets, with the companion having a significantly larger degree of polarization and a very different angle of polarization.

DH Tau, at a distance of 135 pc⁴, is located atthe front side of the Taurus molecular cloud complex that extends from at least 126 to 163 pc (Galli et al. 2018). To determine whether DH Tau B is intrinsically polarized, we therefore need to determine the contribution of interstellar dust to the measured polarization. The interstellar polarization is a result of dichroism by elongated dust grains that are aligned with the local (galactic) magnetic field. Because interstellar dust creates the same polarization for the companion and the star, this contribution can often be determined from the measured stellar polarization (e.g., for 1RXS J1609, see Sect. 5.3, and ROXs 42B, see Jensen-Clem et al. 2020). However, we cannot do that in this case because the star hosts a disk that we spatially resolve in our images (see Sect. 5.5 and Fig. 12, top left) and therefore the stellar polarization is likely a combination of intrinsic and interstellar polarization.

To investigate the contribution of interstellar dust to the polarization of DH Tau B, we show in Fig. 8 a map of the polarization of DH Tau A and B and a few dozen nearby stars. The map is superimposed on a Herschel-SPIRE (Pilbratt et al. 2010) image at 350 μm that shows the concentrations of interstellar dust in the region. White lines show optical measurements of stars at the periphery of the B216-B217 dark cloud from Heyer et al. (1987). Yellow lines display measurements from Moneti et al. (1984) of the three nearest bright stars to DH Tau. Of these stars, HD 283704 (58 pc) is unpolarized as it is located in front of the clouds, whereas HD 283705 (170 pc) and HD 283643 (396 pc) are located behind the clouds and are both polarized with an angle of polarization of 26 ± 1°. Because the stars from Heyer et al. (1987) and Moneti et al. (1984) are generally much older than DH Tau and are therefore not expected to have a circumstellar disk that significantly polarizes their light, their polarization must primarily originate from interstellar dust. Comparing the angles of polarization of DH Tau A (128 ± 1°) and B (58 ± 2°) with those of the reference stars in Fig. 8, we conclude that the polarization of both DH Tau A and B must include an intrinsic component.

We now set limits on the interstellar degree of polarization of the DH Tau system. To this end, we convert the optical measurements of the degree of polarization of the nearby stars HD 283705 and HD 283643 (2.48 and 1.27%) from Moneti et al. (1984) to H-band. For this conversion we use Serkowski’s law of interstellar polarization (Serkowski et al. 1975): $$P=P_{\text{max}}\exp \left[-K{\ln }^2\left({\lambda }_{\text{max}}/\lambda \right)\right],$$(3)

where λ is the wavelength of the light, P_max is the maximum degree of polarization, and λ_max is the wavelength at which this maximum occurs. The parameter K is computedfollowing Whittet et al. (1992): $$K=0.01+1.66{\lambda }_{\text{max}},$$(4)

with λ_max in micrometers. Because the observations were taken without color filter, we retrieve the spectral response of a Ga-As photomultiplier tube similar to that used for the measurements⁵ and multiply it with the transmission of the Earth’s atmosphere. With the resulting spectral transmission, we can compute the degree of polarization that the instrument measures from the transmission-weighted average over the curve from Serkowski’s law. Assuming λ_max = 0.55 μm, which is the average value for the 16 bright stars in Taurus observed by Whittet et al. (1992), we fit P_max for both stars. From the fitted curves we then compute the degree of polarization at H-band, yielding 0.9% for HD 283705 and 0.5% for HD 283643. Because DH Tau is located at the front side of the clouds (rather than behind the clouds as are the comparison stars), the interstellar polarization of DH Tau is most likely below 0.9%, probably below 0.5%. This is in agreement with the H-band degrees of polarization of three nearby T Tauri stars whose archival SPHERE-IRDIS polarimetric data we analyzed (see Fig. 8). Of these stars, DF Tau (125 pc) is unpolarized, and DK Tau A (128 pc), which does not have a disk, and IQ Tau (131 pc), which has a very faint disk, are 0.33% and 0.34% polarized, respectively, both with an angle of polarization of ~ 30°.

Although we do not know the exact interstellar degree of polarization for DH Tau, the angle of polarization is likely close to 26°, which is the angle of both HD 283705 and HD 283643. To see whether DH Tau B is intrinsically polarized, we take the polarization signal that we measured in the mean-combined images (0.48 ± 0.03% at 58 ± 2°; see Table 3) and subtract interstellar polarization signals with an angle of polarization of 26° and a range of degrees of polarization. The resulting intrinsic degree and angle of polarization of DH Tau B versus the interstellar degree of polarization is shown in Fig. 9 (top). We see that the intrinsic polarization decreases for interstellar degrees of polarization between 0% and 0.2% and increases for larger interstellar polarizations. The intrinsic polarization increases because an ever larger interstellar polarization needs to be canceled to produce the measured polarization. For the range plotted, the intrinsic angle of polarization increases from 60° to 100°. Most importantly, the intrinsic degree of polarization is always higher than 0.4%, showing that DH Tau B should be intrinsically polarized if the interstellar polarization indeed has an angle of polarization of 26°.

From the measurements by Heyer et al. (1987) (white lines in Fig. 8), we see that there are slight variations in the angle of polarization of the stars in the region. Goodman et al. (1992) determined that the angles of polarization of these stars are Gaussian distributed with a mean of 27° and a standard deviation of 15°. Using this distribution of angles, we take a more probabilistic approach and perform a Monte Carlo simulation in which we compute for a range of interstellar degrees of polarization the probability distribution of the intrinsic polarization. The histograms of the resulting distributions for each value of the interstellar degree of polarization are displayed in Fig. 9 (bottom). In this figure we have normalized the distribution of each column to one. It follows that the curves of Fig. 9 (top) are in fact among the most probable scenarios. We also see that DH Tau B must be at least 0.2% intrinsically polarized for interstellar degrees of polarization between 0 and 0.3% or higher than 0.7%, regardless of the interstellar angle of polarization. Only for interstellar degrees of polarization between 0.3 and 0.7% there is a small possibility (~ 8%) that DH Tau B is not intrinsically polarized. Based on these findings, we conclude that DH Tau B is very likely intrinsically polarized.

Fig. 8 Map of the linear polarization of DH Tau B and nearby stars superimposed on a Herschel-SPIRE map at 350 μm. The length and orientation of the lines indicate the degree and angle of linear polarization, respectively. The black line shows the H-band polarization we measure for DH Tau B, and the orange lines display the SPHERE-IRDIS H-band measurements of DH Tau A and three other nearby T Tauri stars whose archival data we analyzed. White lines show optical measurements by Heyer et al. (1987). Yellow lines indicate the H-band polarization of three bright stars closest to DH Tau as derived from optical measurements by Moneti et al. (1984). The length of the H-band vectors are scaled by a factor of four with respect to the optical vectors.

Fig. 9 Intrinsic polarization of DH Tau B after subtracting interstellar polarization signals from the measured polarization of the companion. Top: intrinsic degree and angle of linear polarization of DH Tau B as a function of the degree of polarization due to interstellar dust, assuming an angle of 26° for the interstellar polarization. The bands around the curves show the uncertainties of our measurements. Bottom: probability distributions of the intrinsic degree of polarization of DH Tau B for a range of degrees of polarization due to interstellar dust, assuming the angle of the interstellar polarization to have the same distribution as that determined by Goodman et al. (1992) for the B216-B217 dark cloud adjacent to DH Tau. The probability distribution of each column is normalized to one.

5.2 Likely detection of intrinsic polarization of GSC 6214 B

In this section we present the likely detection of intrinsic polarization originating from GSC 6214 B. Figure 10 shows the reduced I_Q-, Q-, and U-images in H-band at the position of the companion of the data set created by mean-combining the final images of the three data sets. Table 4 shows the measured polarization of GSC 6214 B for the three individual data sets and the mean-combined one. Similar to the DH Tau data, the measured polarization signals of each data set are within the uncertainties constant with aperture radius. We select a final aperture radius of 4 pixels, corresponding to the (approximate) radius where the S/N in q and u is maximized in each of the data sets. Overall the measured degree and angle of polarization of the data sets are consistent within the uncertainties. The slightly different results of the 2019-02-22 data set compared to the other two data sets could be caused by the relatively strong time-varying atmospheric conditions that the observations were taken under (see Table 2). Whereas the q- and u-measurements of the three data sets individually do not reach the required 5σ-limit for a detection, the mean-combined measurement does, reaching an S/N of 5.2 in u. From the mean-combined data we therefore conclude that we detect significant polarization from GSC 6214 B, with a degree and angle of polarization of 0.23 ± 0.04% and 138 ± 5°, respectively.

Table 4 also shows the stellar degrees and angles of polarization. Because we do not spatially resolve a disk around GSC 6214 A, we used a star-centered aperture extending up to and including the AO residuals to maximize the S/N. The measured signals show significant differences and are overall inconsistent among the data sets. The signals average to a degree of polarization of only 0.10%. The measurements of the stellar polarization are therefore most likely dominated by spurious signals. To determine whether the companion is truly polarized, we need to investigate the potential origins of these spurious signals and the effect they have on the measurement of the companion polarization.

If the stellar polarization primarily results from uncorrected instrumental polarization, which to first order equally affects the star and thecompanion, we would need to subtract these signals from the images. Using the mean-combined images with the stellar polarization subtracted, we measure for the companion a degree and angle of polarization of 0.32± 0.04% and 141 ± 4°, respectively, with an S/Nof 1.4 in q and 7.2 in u. This polarization signal is larger and more significant than that measured from the images without the stellar polarization subtracted (see Table 4). However, the measured signals are less consistent among the data sets, suggesting that uncorrected instrumental polarization may not be the principal cause of the stellar polarization.

A more likely scenario seems that the stellar polarization signals are dominated by systematic effects due to time-varying atmospheric conditions and AO performance in combination with the coronagraph, similar to the case of DH Tau (see Sect. 5.1). Also in the case of GSC 6214, the systematic effects do not affect (as much) the companion measurements because those measurements are overall consistent among the data sets. This suggests that the measurements of the companion are more reliable than those of the star. Because the companion polarization is significantly different from the stellar polarization in all data sets, particularly in the angle of polarization (see Table 4), and we measure significant polarization from the companion for both the reduction with and without the stellar polarization subtracted (reaching S/Ns of 7.2 and 5.2 in u, respectively), we conclude that the companion is most likely truly polarized.

To determine whether the polarization of GSC 6214 B is intrinsic to the companion or caused by interstellar dust, we show in Fig. 11 a map of the angles of polarization of nearby bright stars from the catalog by Heiles (2000). The map is displayed over an IRAS (Neugebauer et al. 1984) 100 μm map that shows the dust concentrations in the region of the Ophiuchus molecular cloud complex where GSC 6214 is located. Comparing the angle of polarization of GSC 6214 B and the nearby stars, it may seem that the companion is polarized by interstellar dust. However, GSC 6214 is located at 109 pc, whereas estimates for the distance of the Ophiuchus molecular cloud complex range from approximately 120 to 150 pc (e.g., Mamajek 2008; Lombardi et al. 2008; Ortiz-León et al. 2017; Yan et al. 2019). Indeed, the three stars closest to GSC 6214 in Fig. 11 are located at 128 to 131 pc. We therefore consider it more likely that GSC 6214 is located in front of the main concentrations of dust. In addition, if the companion were polarized by interstellar dust, we would expect to measure in all data sets a stellar polarization with the same angle of polarization as the companion (which is the case for 1RXS J1609; see Sect. 5.3). In principle it is possible that GSC 6214 A is not significantly polarized because the interstellar polarization is canceled by intrinsic polarization due to an unresolved circumstellar disk. However, this scenario seems very unlikely because Bowler et al. (2015) do not detect a disk with ALMA and put an upper limit on the disk’s mass as low as 0.0015% of the mass of the star. Taking into account all considerations, we conclude that it is likely that the polarization we measure for GSC 6214 B is intrinsic to the companion, but we stress that we are less confident than for DH Tau B.

Fig. 10 Reduced mean-combined IQ-, Q-, and U-images (after applying the cosmetic correction described in Appendix A) at the position of the companion GSC 6214 B, showing an aperture of radius 4 pixels centered on the companion. The IU -image, which is not shown, is very similar to the IQ-image.

Fig. 11 Map of the angle of linear polarization of the companion GSC 6214 B, the star 1RXS J1609 A, and other nearby bright stars superimposed on an IRAS map at 100 μm. The angles of GSC 6214 B and 1RXS J1609 A are from the SPHERE-IRDIS H-band measurements from this work, whereas for the other stars the angles are taken from the catalog of optical measurements by Heiles (2000). The length of the lines is arbitrary and contrary to Fig. 8 does not indicate the degree of polarization. HD 147700 in unpolarized. White lines indicate stars at a distance between 128 pc and 142 pc, and blue and orange lines show objects closer or farther away, respectively. We note that the region shown is much larger than that ofFig. 8, and therefore the angular separation among the stars is much larger as well.

5.3 Detection of interstellar polarization from 1RXS J1609 B

In this section we present the detection of polarization in the 1RXS J1609 system. In all six data sets of 1RXS J1609, we consistently measure within the uncertainties the same degree and angle of polarization for the central star 1RXS J1609 A. In the mean-combined data set, which uses the four highest-quality data sets (2019-08-06, 2019-08-29, 2019-09-17B, and 2019-09-23), we measure for the star a degree and angle of polarization of 0.21± 0.01% and 97 ± 2°, respectively. In the same data set, we measure for the companion 0.2 ± 0.1% and 95 ± 26°, using an aperture radius of 5 pixels. Although the measurement of the companion polarization is not a significant detection, itis striking that it is within the uncertainties the same as the measured stellar polarization. In the six individual data sets we also measure the polarization of the companion to be consistent with that of the mean-combined data, although with higher uncertainties. Finally, using the mean-combined data, we measure for the relativelybright background object that is also visible in the field of view a degree and angle of polarization of 0.3 ± 0.1% and 103 ± 50°, respectively.Because all three objects have within the uncertainties the same degree and angle of polarization, their polarization likely originates from the same source, that is, from interstellar dust.

To confirm this scenario, we turn to Fig. 11, which shows that 1RXS J1609 is located in the Ophiuchus molecular cloud complex a few degrees west from GSC 6214. Contrary to GSC 6214, 1RXS J1609, at a distance of 139 pc, is definitely located within the dust clouds that are located at a distance of approximately 120 to 150 pc (see Sect. 5.2). Indeed, the measured angles of polarization of 1RXS J1609 A, 1RXS J1609 B and the background object agree well with those of the nearby bright stars located at a similar distance (see Fig. 11). Serkowski et al. (1975) have fitted their multiwavelength optical measurements of the stars HD 144470 and HD 144217 (142 and 129 pc; see Fig. 11) to Serkowski’s law of interstellar polarization (see Eq. (3)) and determined the values of P_max and λ_max for both stars. Using these values, we find that in H-band the degrees of polarization are equal to approximately 0.4 and 0.3%, respectively. These values are similar to the degree of polarization we measure for the star, the companion, and the background object in our images of 1RXS J1609, where the slight differences are likely due to the inhomogeneous spatial distribution of the interstellar dust. We conclude that the polarization we measure for 1RXS J1609 B originates from interstellar dust and therefore set an upper limit on the degree of polarization in Sect. 5.4.

5.4 Upper limits on polarization of other companions

In this section we present upper limits on the degree of polarization of the 18 companions for which we do not reach the 5σ-limit in q or u to claim a detection. For the majority of the companions, the S/N in q and u is typically ≲ 2–3 for any aperture radius. For four companions the maximum S/N in q or u reaches a value of almost 4. However, in these four cases the signals in the Q- and U-images (after the cosmetic correction described in Appendix A) do not resemble scaled-down positive or negative versions of the total-intensity PSF as one would expect for real signals, but show strong pixel-to-pixel variations caused by incompletely removed bad pixels (see Appendix B).

Table 5 shows for each target the upper limits determined from the 68.27 and 99.73% intervals, as described in Appendix C. For targets for which we obtained multiple data sets, we computed the upper limits from the mean-combined images. For the majority of the companions, which are generally the fainter ones, we determined the upper limits using an aperture radius equal to half times the full width at half maximum (FWHM) of the stellar PSF. This aperture radius is on average 1.9 pixels in H-band and 2.6 pixels in K_s-band, and is at, or close to, the radius at which the upper limit is minimized. For seven, generally brighter companions (CD-35 2722 B, AB Pic b, HD 106906 b, GQ Lup B, GSC 8047 B, PZ Tel B in J-band, and 1RXS J1609 B) we used an aperture radius of 5 pixels to average out and suppress the spurious signals created by incompletely removed bad pixels (see Appendix B). However, the bad pixels generally still create a bias in the q- and u-signals, and so we have to accept that this increases the upper limits. For the data sets where this bias is really strong (i.e., CD-35 2722, PZ Tel in J-band, and TYC 8998), we excluded from the data reduction those frames that contribute strong bad pixels at the position of the companion in the final images. Because HD 106906 b is located at an angular separation of 7.1′′ from the central star, which is larger than the isoplanatic angle during the observations, its PSF is strongly elongated in the radial direction from the star. To account for this, we used an elliptically shaped aperture. Finally, for the companions of HR 8799, HD 206893, and β Pic, we computed the upper limits using the polarimetric images from the reduction with the added classical ADI step (see Appendix C).

Table 5 also shows for each target the stellar degree and angle of polarization. For the majority of the stars the degree of polarization is around 0.1%. To be conservative and because we generally do not know the origin of these low polarization signals (intrinsic, interstellar dust or spurious), we interpret the signals as biases. For these targets we therefore computed the upper limits on the companion polarization from both the reductions with and without the stellar polarization subtracted, and show the highest values in Table 5. For three targets we measure a stellar polarization higher than approximately 0.1%. In the case of GQ Lup and PDS 70 this stellar polarization is caused by a circumstellar disk (see Keppler et al. 2018 and Sect. 5.5). Although GQ Lup is located in the Lupus I cloud, the contribution of interstellar dust is likely small because HD 141294, the nearest bright star to GQ Lup (at 14.3′ and a distance of 153 pc compared to 151 pc for GQ Lup), is unpolarized at optical wavelengths (Rizzo et al. 1998; Alves & Franco 2006). For PDS 70 and GQ Lup we therefore determined the upper limits using only the images without the stellar polarization subtracted. For 1RXS J1609 on the other hand, the stellar polarization is caused by interstellar dust (see Sect. 5.3), and we therefore used the reduction where the stellar polarization is subtracted.

Examining the upper limits in Table 5, we see that for 11 companions the 68.27% upper limits are ≤ 0.3%, with the lowest upper limit equal to 0.06% for PZ Tel B in H-band. These low upper limits are in almost all cases dominated by the photon noise from the companion in the Q- and U-images or the bias dueto incompletely removed bad pixels. The upper limits are still larger than the (minimum) polarimetric accuracy of the Mueller matrix model with which the data have been corrected (see van Holstein et al. 2020). For the companions of HR 8799, HD 19467, and HD 206893, which are fainter or located at a much smaller separation than the other companions, the 68.27% upper limits are dominated by the uncertainty of the background in Q and U and have values between 0.4 and 0.8%. For the very close-in planet PDS 70 b we reach upper limits of 5.0% in K_s -band and 9.2% in H-band. These upper limits are so high because the comparison apertures contain signal from the inner circumstellar disk of PDS 70 A (see Fig. 12) and the Student’s t-distribution imposes a large statistical penalty for the low number of available comparison apertures (see Appendix C). We note that for PDS 70 c (Haffert et al. 2019), the circumstellar disk prevents us from measuring the polarization altogether. Finally, we reach the highest polarimetric point-source contrast in the mean-combined data set of β Pic, with a 1σ-contrast of 3 × 10⁻⁸ at a separations of 0.5′′ and a contrast below 10⁻⁸ for separations > 2.0′′ (see Appendix E). Overall, it follows that our measurements are sensitive to polarization signals of around a few tenths of a percent.

5.5 Detection of circumstellar disks of DH Tau, GQ Lup, PDS 70, β Pic, and HD 106906

In our survey we also detected the five circumstellar disks displayed in Fig. 12. Although the disks of DH Tau and GQ Lup have already been detected at mm-wavelengths (Wolff et al. 2017; MacGregor et al. 2017; Wu et al. 2017b), here we present the first images in polarized scattered light, revealing various interesting features. For PDS 70, HD 106906, and β Pic near-infrared polarimetric images already exist (Keppler et al. 2018; Hashimoto et al. 2012; Kalas et al. 2015; Millar-Blanchaer et al. 2015), but our images are generally deeper, reveal new features, or confirm features that were previously observed. In this section, we therefore briefly discuss these disks, although we consider a detailed analysis beyond the scope of this paper.

Figure 12 (top left) shows the polarized intensity image of the DH Tau system, with the circumstellar disk visible in the top right corner of the panel. The relatively small disk has a diameter of approximately 0.50′′ or 67 au at 135 pc. From ALMA measurements of the Keplerian rotation of the disk, Sheehan et al. (2019) have determined an inclination of 48° and a position angle of 2.5° (east of north), with the northern side of the disk rotating toward us (i.e., blue shifted). In our images the disk has a smooth intensity profile with no visible disk gap, rings, or spirals. A strong brightness asymmetry is visible between the eastern and western sides of the disk, which could be related to the viewing angle of the disk and the dust scattering properties. This asymmetry is compatible with the position angle derived from ALMA: if the side inclined toward the Earth appears brightest due to enhancedforward scattering, then the eastern side is the forward-scattering near side of the disk. Alternatively, this brightness asymmetry could result from shadowing by an unresolved inner disk component because the brightness changes quite abruptly with azimuth. The brightness asymmetry might extend toward the inner (coronagraphically masked) parts of the disk because the angle of polarization that we measure for the average stellar polarization (128°, see Table 3) agrees well with the angle of polarization one obtains when integrating over the non-obscured parts of the disk. In the bottom left corner of the panel the polarization signal of DH Tau B is visible, where the angle of polarization is indicated with the two lines protruding from the circle around the companion.

Figure 12 (top row, second column) shows the polarized intensity image of the circumstellar disk and companion of GQ Lup. From ALMA images (MacGregor et al. 2017; Wu et al. 2017b), which show a rather featureless disk, the disk inclination and position angle are known to be 60° and 346°, respectively. Our scattered light images show a remarkable north-south asymmetry in the circumstellar disk, with the southern part of the disk extending out to 0.84′′ (127 au at 151 pc) and the northern part only out to 0.49′′ (74 au). Two spiral-like features can be seen protruding eastward from the southern part of the disk. The disk asymmetry and spiral-like features are reminiscent of those of the disk around RY Lup (Langlois et al. 2018) and could be the result of periodic close passes of GQ Lup B (see e.g., Dong et al. 2016; Cuello et al. 2019, 2020). The orbital analyses presented by Schwarz et al. (2016) and Wu et al. (2017b) indeed show that the orbit of GQ Lup B is almost certainly eccentric and that it is quite likely that the inclinations of the orbit and the disk are similar. However, Wu et al. (2017b) argue that although the inclinations may be similar, the disk and companion orbit are likely not in the same plane. We also find that the starlight of GQ Lup is polarized due to the unresolved part of the circumstellar disk, with an angle of polarization (83 ± 1°; see Table 5) approximately perpendicular to the position angle of the disk. GQ Lup B appears to be polarized in Fig. 12 (top row, second column), but this polarization is spuriously created by subtracting the stellar polarization from the image. We will present a dynamical analysis of the complete system and detailed radiativetransfer and hydrodynamical modeling of the circumstellar disk in a future paper (van Holstein et al. in prep.).

Figure 12 (top row, third and fourth columns) show the H-band Q_ϕ - and U_ϕ -images of the circumstellar disk around PDS 70. The disk is seen at a position angle of 159° and an inclination of 50°, with the southwestern side being the near side (Hashimoto et al. 2012). The Q_ϕ -image clearly shows the known azimuthal brightness variations of the outer disk ring, as well as bright features close to the coronagraph’s inner edge that most likely originate from the inner disk (see Keppler et al. 2018). The U_ϕ -image contains significant signal, with the maximum value equal to ~49% of the maximum in the Q_ϕ-image, revealing the presence of non-azimuthal polarization. The pattern in U_ϕ agrees well with the radiative transfer models by Canovas et al. (2015), indicating that part of the photons are scattered more than once. The Q_ϕ -image also shows a weak spiral-like feature extending toward the east from the northern ansa of the disk and perhaps a similar feature at the southern ansa. With these features the disk resembles the model images by Dong et al. (2016) for the inclination and position angle of the PDS 70 disk. We may therefore be seeing the effect of two spiral arms in the outer disk ring, potentially induced by PDS 70 b.

Figure 12 (top right) shows the Q_ϕ -image of the debris disk of HD 106906, which is viewed close to edge-on. The forward-scattering near side of the disk can be seen passing slightly to the north of the star. The image clearly shows the known east-west brightness asymmetry of the disk, which had until now only been detected in total intensity (Kalas et al. 2015; Lagrange et al. 2016). Because our data are particularly deep (i.e., 120 min total on-source exposure time), we detect the backward-scattering far side of the disk to the west of the star, just south of the brighter near side of the disk (see Kalas et al. 2015).

Finally, Fig. 12 (center) shows the Q_ϕ -image of the nearly edge-on-viewed debris disk of β Pic. The disk extends from one side of the 11′′ × 11′′ IRDIS field of view to the other. Earlier near-infrared scattered light images reported by Millar-Blanchaer et al. (2015) show the disk only to ~ 1.7′′ or 33 au at 20 pc due to the smaller field of view of GPI. In our images we see the disk extending to at least 5.8′′ or 115 au on both sides of the star. The disk midplane is seen slightly offset to the northwest of the star (up in Fig. 12 center) due to the disk’s small inclination away from edge-on. Our image also shows the apparent warp in the disk (see Millar-Blanchaer et al. 2015, and references therein) that extends eastward (to the bottom left) in the northeastern (left) part of the disk and westward (to the upper right) in the southwestern (right) part of the disk. This warp is particularly well visible in Fig. 12 (bottom), which shows a total-intensity image after applying ADI with PCA using IRDAP.

Fig. 12 H-band images of the five circumstellar disks detected in our sample, showing the linearly polarized intensity (PI) for DH Tau and GQ Lup, Qϕ and Uϕ for PDS 70, and Qϕ for HD 106906 and β Pic. A total-intensity image after applying ADI with PCA (subtracting four principal components) is additionally shown for β Pic. The Uϕ -image of PDS 70 is shown on a linear scale, whereas all other images are shown on a logarithmic scale. All polarimetric images have the polarized stellar halo subtracted. The angular scale and sky orientation are indicated in each image. Gray circles maskregions obscured by the coronagraph.

6 Modeling of polarization from circumsubstellar disks

As discussed in Sects. 5.1 and 5.2, we (very) likely detected intrinsic polarization from DH Tau B and GSC 6214 B, with a degree of polarization of several tenths of a percent in H-band. The host stars of these two companions are among the youngest in our sample ( ≲20 Myr) and the companions have indicators for the presence of circumsubstellar disks through hydrogen emission lines, red near-infrared colors, and excess emission at mid-infrared wavelengths (see Fig. 1 and Table 1). Therefore, the most plausible explanation for the polarization in these cases is scattering of the companion’s thermal emission by dust within a spatially unresolved circumsubstellar disk. However, we note that the late M to early L spectral types of these low-mass companions (see Fig. 1 and Table 1) suggest their atmospherescould be dusty. As a result, the polarization could also originate from rotation-induced oblateness, an inhomogeneous cloud distribution, or a combination of these atmospheric asymmetries and a disk (see Stolker et al. 2017). Still, it seems reasonable to assume that the polarization is solely caused by a disk because the companions have low projected rotational velocities (Bryan et al. 2018; Xuan et al. 2020), and out of the 20 companions observed, we only detect intrinsic polarization for the companions that have hydrogen emission lines.

In this section we perform (spatially resolved) radiative transfer modeling of a representative example of a circumsubstellar disk to investigate whether our detections of polarization of several tenths of a percent can really be explained by such disks. To this end, we first describe the setup of the radiative transfer model in Sect. 6.1. We then examine the generation of an integrated (i.e., spatially unresolved) polarization signal in Sect. 6.2 and thedependence of the polarization on the properties of the disk in Sect. 6.3. We stress that we consider an isolated circumsubstellar disk (i.e., it is not embedded in a circumstellar disk) and that our models are general and not tailored to either DH Tau B or GSC 6214 B. Because we only study the degree and angle of polarization produced by the disk, the exact spectrum of the companion has little effect on the results. In Sects. 7.1 and 7.2 we use the results of our modeling to interpret and discuss our measurements.

6.1 Setup of the radiative transfer model

To quantify the expected near-infrared polarization from a self-luminous atmosphere with a circumsubstellar disk, we computed a radiative transfer model with MCMax (Min et al. 2009), which is a Monte Carlo radiative transfer code for axisymmetric disks that is optimizedfor the high optical depths in protoplanetary disks. The model considers a passive, irradiated disk around a self-luminous substellar atmosphere (the contribution from the light of the central star is negligible). We selected a synthetic spectrum from the BT-Settl atmospheric models (Allard et al. 2012) at an effective temperature T_eff= 2000 K and surface gravity log g = 4.0 dex. We then scaled the spectrum to a luminosity of 10⁻⁴ L_⊙ by assuming a radius for the atmosphere of 2 R_Jup at an age of ~10 Myr (e.g., Baraffe et al. 2015). We then modeled the circumsubstellar disk as a scaled down version of a circumstellar disk (see e.g., Williams & Cieza 2011). We parametrized the structure of the circumsubstellar disk with a profile for the dust surface density that is inversely proportional to the radius, Σ ∝ r⁻¹. Using a surface density at the inner radius of Σ_in = 0.07 g cm⁻² and an inner and outer disk radius of R_in = 0.003 au and R_out = 0.01 au, we computed the total mass residing in the solids. For the pressure scale height, we used a linear dependence with the disk radius, h ∝ r, with a (constant) aspect ratio of h∕r = 0.1. The dust opacities contain by volume 60% silicates, 15% amorphous carbon, and 25% porosity (Woitke et al. 2016). Furthermore, we used a maximum hollow volume ratio of 0.8 for the distribution of hollow spheres, which approximates the irregularity of the dust grains (Min et al. 2016). The size distribution of the grains was chosen in the range of 0.05–3000 μm with a power-law exponent of − 3.5. Dust settling is included with the prescription from Dubrulle et al. (1995), which assumes an equilibrium between turbulent mixing and gravitational settling. In this way, the dust scale height is a function of disk radius and grain size, which is controlled by the viscosity parameter α = 10⁻⁴.

6.2 Origin of the spatially integrated polarization

After setting up the disk and dust properties, we can now perform the radiative transfer computations to study the generation of a spatially integrated polarization signal from the disk. We propagate the Monte Carlo photons through the disk to compute the thermal structure and the local source function. We then run a monochromatic ray tracing at 1.62 μm (the central wavelength of the IRDIS H-band filter) to compute the synthetic total-intensity and Stokes Q- and U-images. Figure 13 displays an example image of the total-intensity surface brightness for a disk inclination of 70°. In this figure, the length and orientation of the lines indicate the local degree and angle of polarization, respectively. Finally, we compute the spatially integrated polarization using the sum of the pixel values in each of the Stokes images. In Fig. 13, this results in an integrated degree and angle of polarization of 0.24% and 0°, respectively. Indeed, the polarized flux is largest along the major axis of the disk, where scattering angles are closest to 90°, yielding a net polarization that is oriented perpendicular to the major axis of the disk. In fact, the angle of polarization is always perpendicular to the position angle of the disk, independent of the disk inclination.

For the interpretation of a nonzero integrated polarization, we need to consider various effects that are visible in the spatially resolved image of the disk in Fig. 13. To produce a measurable degree of polarization, the linearly polarized intensity should have a nonzero value while lowering the total intensity will further enhance the degree of polarization. In the example of Fig. 13, most of the polarized flux comes from the inner edge of the disk, where the flux in total intensity is about 10 to 100 times lower than the atmospheric emission. Part of the polarization signal is canceled because there is both horizontally and vertically polarized flux, but a net vertically polarized flux remains. The local degree of polarization increases along the major axis of the disk toward larger separations because of reduced multiple scattering. However, the total intensity is also lower in these regions such that the polarized intensity is also low there. This means that the integrated polarization depends primarily on the inner radius and the surface density, whereas the outer radius, and therefore the total disk mass, are much less relevant. Because the inner radius of the disk is at ~6 R_Jup and the inclination is 70°, part of the photosphere of the central object is obscured by the near side of the disk. This reduces the total intensity of the system such that the net degree of polarization is enhanced compared to a situation in which the full atmosphere would be visible.

Fig. 13 Synthetic image of a self-luminous companion (Teff = 2000 K) with a compact circumsubstellar disk at a distance of 150 pc. The total intensity surface brightness is shownon a logarithmic scale and the lines indicate the local degree and angle of linear polarization. The spatially integrated degree of polarization is 0.24% and the angle of polarization is 0°, that is, the spatially unresolved light is linearly polarized along the minor axis of the disk.

6.3 Dependence on the inner radius and surface density

We now investigate the dependence of the spatially integrated degree of polarization on the inner radius and the surface density at the inner radius. To this end, we run a grid of 10 × 10 radiative transfer models with a varying inner radius (5–41 R_Jup) and dust surface density at the inner radius (10⁻²–10² g cm⁻²). All other parameters are the same as in Sect. 6.1, except for the outer radius which we changed from 0.01 to 0.4 au. In this way, the disk remains radially sufficiently extended even though the outer radius of the disk has a negligible impact on integrated degree of polarization because most of the polarized flux comes from the inner edge. Because thetotal disk mass depends on the inner radius and surface density, it is different for each model. We note that the estimated polarization may rely on additional properties of both the disk structure and the dust grains.

As discussed in Sect. 6.2, the integrated degree of polarization depends strongly on the fractional occultation of the substellar atmosphere by the disk. This effect occurs at a high enough inclination if the projected disk reaches close to the atmosphere and/or the vertical extend of the disk (which scales with the dust surface density) is sufficiently large. To resolve with a high precision the obscuration of the atmosphere, we perform the ray tracing at sufficient spatial resolution. We set the disk inclination i to 70° and 80° because for geometry reasons detections are biased toward highly inclined disks and, more importantly, a nonzero polarization from a circumsubstellar disk is only to be expected if the disk is sufficiently inclined. For example, we find that for i < 45° and i < 20°, the degree of polarization is <0.15% and <0.03%, respectively. We calculate the integrated degree of polarization as before and present the results for each combination of disk inner radius and the surface density at the inner disk radius in Fig. 14.

At an inclination of 70° (see Fig. 14, top), the polarization reaches a maximum value of 0.4–0.5% when the inner radius is 5–10 R_Jup and the surface density is ≳1 g cm⁻². At small inner radii, there is a correlation with the surface density because increasing the inner radius can be counteracted by an increase in surface density in order to maintain the same integrated degree of polarization. At a given surface density, the degree of polarization converges to a constant value at larger inner radii because the atmosphere is no longer obscured and most of the polarized flux originates from the cavity edge. For higher surface densities, this turnover point occurs at larger disk radii because the scattering surface is higher.

A more extreme picture emerges when the inclination is increased to 80° (see Fig. 14, bottom). Whereas for surface densities ≲0.1 g cm⁻¹ the correlation with the inner radius is comparable to the i = 70° case, at higher surface densities the substellar atmosphere becomes fully obscured by the disk. There is a peak in the degree of polarization when the vertically extended disk obscures the atmosphere along the minor axis of the disk while there is still some disk surface visible close to the major axis. As a result, the total intensity of the atmosphere is strongly reduced while the polarized flux at scattering angles close to 90° is less attenuated, leading to a degree of polarization as high as ~15%. For even higher surface densities, the degree of polarization remains approximately constant at a value of ~8% because both the substellar atmosphere and the cavity edge are obscured by the vertical extent of the disk. In this case, only light that scatters through the surface layer of the disk will reach the observer, which is therefore no longer dependent on the surface density and the inner radius.

Fig. 14 Dependence of the integrated degree of linear polarization on the inner radius of the circumsubstellar disk and the surface density of the dust at the inner radius. The grid of radiative transfer models is shown for a disk inclination of 70° (top) and 80° (bottom).

7 Discussion

In Sect. 6 we performed radiative transfer modeling of a generic circumsubstellar disk to study the origin of the integrated polarization and the dependence of this polarization on the disk properties. We use the results of our modeling in Sect. 7.1 to interpret our likely detections of polarization from DH Tau B and GSC 6214 B and our non-detection for GQ Lup B. In Sect. 7.2 we then briefly examine the non-detections of polarization for 1RXS J1609 B, HD 106906 b, and PDS 70 b, which also have evidence for the existence of a circumsubstellar disk. Subsequently, we outline the implications of our upper limits on the polarization of the other companions with respect to the presence of atmospheric asymmetries in Sect. 7.3. Finally, in Sect. 7.4, we discuss potential measurements with various instruments to confirm and further characterize the circumsubstellar disks of DH Tau B and GSC 6214 B.

7.1 Circumsubstellar disks, rotational periods, and formation of DH Tau B, GSC 6214 B, and GQ Lup B

As discussed in Sect. 6, the most plausible explanation for the polarization of DH Tau B and GSC 6214 B is the presence of a circumsubstellar disk. From the radiative transfer modeling in that same section, we see that the integrated degree of polarization of such a disk depends on many parameters and that estimating disk properties is therefore a degenerate problem. Nevertheless, we can still put constraints on the dust grain sizes and the disk’s inclination and position angle, and through that constrain the rotational periods and formation mechanisms of the companions.

Whereas we most likely detect polarization from the disks of DH Tau B and GSC 6214 B, no emission has been detected from these companions at mm-wavelengths (Bowler et al. 2015; Wu et al. 2017a, 2020; Wolff et al. 2017). It is therefore possible that these disks contain mostly micrometer-sized dust grains and only little mm-sized grains, or, as suggested by Wu et al. (2017a), that the disks are compact and optically thick at mm-wavelengths. From our polarimetric measurements we cannot determine whether the disks are really compact because most of the polarized flux originates from the inner edge of the disk (see Sect. 6.2). Because we do not spatially resolve the disks, we can put a limit on the disk size from the measured FWHM of the PSF. The FWHM corresponds to a maximum disk radius of ~ 3 au for both companions. This radius is much smaller than one-third of the Hill radius (14–20 au, assuming the companions are on a circular orbit), which is the distance at which the disks are suspected to be truncated due to tidal interactions with the primary star (Ayliffe & Bate 2009; Martin & Lubow 2011; Shabram & Boley 2013). However, it is possible that the disks extend beyond 3 au, but that we do not reach the sensitivity and contrast to detect the flux at the outer regions.

From the measured degree of polarization we can put constraints on the inclination of the disks. With degrees of polarization of a few to several tenths of a percent, the disks of DH Tau B and GSC 6214 B must have a high inclination because a low-inclination disk will have a very low, nearly zero degree of polarization below the sensitivity that we reach with our measurements (see Sect. 6.3). In fact, it could be that GQ Lup B hosts such a low-inclination disk because we do not detect significant polarization although the measured hydrogen emission lines are stronger than those of DH Tau B and GSC 6214 B (Zhou et al. 2014). We also see that the inclination of the disks of DH Tau B and GSC 6214 B cannot be close to edge-on so that disk completely obscures the companion’s atmosphere because in that case we would measure polarization degrees of several toeven ten percent. Such a high degree of polarization of 14% has been measured for CS Cha B in J- and H-band by Ginski et al. (2018), which the authors indeed interpret as being caused by a highly inclined and vertically extended disk. This interpretation was recently confirmed by Haffert et al. (2020) using medium-resolution optical spectroscopy with MUSE.

The projected rotational velocity, vsini, has been measured for DH Tau B, GSC 6214 B and GQ Lup B through high-resolution spectroscopic observations (Xuan et al. 2020; Bryan et al. 2018; Schwarz et al. 2016), finding values of 9.6 ± 0.7, ${6.1}_{-3.8}^{+4.9}$, and ${5.3}_{-1.0}^{+0.9}$ km s⁻¹, respectively.Assuming that the spin axes of the companions are perpendicular to the plane of their disks (the regular moons of our solar system’s giant planets, which are believed to have formed in circumsubstellar disks, orbit near the equatorial plane of the planets) andtaking the companion radii and uncertainties from Xuan et al. (2020) and Schwarz et al. (2016), we constrain the rotational period of the companions using a Monte Carlo calculation. We assume the inclination to be uniformly distributed in cos i, with values between 60° and 80° for DH Tau B and GSC 6214 B and between 0° and 45° for GQ Lup B. We find rotational periods equal to 29–37 h for DH Tau B, 22–77 h for GSC 6214 B, and 19–64 h for GQ Lup B,within the 68% confidence interval. These estimates of the rotational periods are roughly an order of magnitude larger than the average periods expected from the period-mass relation as determined from observations of free-floating low-mass brown dwarfs of similar ages (e.g., Rodríguez-Ledesma et al. 2009; Scholz et al. 2015, 2018). This discrepancy can be explained by the companions hosting circumsubstellar accretion disks. The estimated slow rotation of the companions, which is at ~ 0.1% of their break-up velocities (see Xuan et al. 2020), is consistent with a scenario in which the companions lose angular momentum to their disks during accretion and should still spin up as they contract (see Takata & Stevenson 1996; Bryan et al. 2018; Xuan et al. 2020). The long rotational periods we find also show that rotation-induced oblateness does not significantly contribute to the measured polarization because polarization > 0.1% is generally expected only for rotational periods of ~6 h or less (Sengupta & Marley 2010; Marley & Sengupta 2011).

As we discussed in Sect. 6.2, the angle of polarization we measure from an unresolved disk is always oriented perpendicular to the position angle of the disk. Therefore, the position angle of the disk of DH Tau B is likely between 150° and 190° (see Fig. 9), whereas that of the disk of GSC 6214 B is around 48° (see Table 4). Because we already found that both disks likely have large inclinations, we have strong constraints on the 3D orientation of the disks. The disk of DH Tau B is most likely misaligned with the circumstellar disk of DH Tau A, which has an inclination and position angle of 48° and 2.5°, respectively (see Sect. 5.5 and Fig. 12), although the position angles could possibly be aligned. Such a misalignment of disks is also found for CS Cha A and B (Ginski et al. 2018). Although GSC 6214 A is not known to host a circumstellar disk, orbital motion has been detected for GSC 6214 B (Pearce et al. 2019). However, the orbital elements are not sufficiently constrained to conclude on possible (mis)alignments of the disk and the orbit. If a low-inclination disk exists around GQ Lup B, it would be misaligned with the circumstellar disk that has an inclination of 60° (see Sect. 5.5 and Fig. 12) and the orbit of GQ Lup B that likely has a similar inclination. However, the circumsubstellar disk could be aligned with the spin axis of GQ Lup A that has an inclination of ~ 30° (Donati et al. 2012).

The misalignment of the disks of DH Tau A and B, and possibly also of GQ Lup A and B, suggests that the companions may have formed in situ through direct collapse in the molecular cloud, akin to the formation mechanism of binary stars. Indeed, CS Cha B, with its misaligned disk, was initially thought to be of substellar nature (Ginski et al. 2018), but was recently found to actually be a low-mass star (Haffert et al. 2020). However, formation through gravitational instabilities in the circumstellar disk is also possible because this mechanism can form companions at separations of up to at least 300 au (Tobin et al. 2016). Although one might expect the circumstellar and circumsubstellar disks to be coplanar in this scenario, misalignments can result if the companion formed away from the midplane of the original disk, the original disk was asymmetric, or the circumstellar disk or other objects disturbed the circumsubstellar disk (Stamatellos & Whitworth 2009; Bryan et al. 2020). It seems unlikely, however, that DH Tau B, GSC 6214 B, and GQ Lup B formed close to their stars and were subsequently scattered to a higher orbit through dynamical encounters with massive inner bodies. This is because tidal interactions would most likely severely disturb or even destroy the circumsubstellar disks (see Stamatellos & Whitworth 2009; Bailey et al. 2013; Bowler et al. 2011), and no massive objects at small separations, nor the gaps they would create in the circumstellar disks, have been detected (see Fig. 12; Pearce et al. 2019; Wu et al. 2017b).

7.2 Circumsubstellar disks of 1RXS J1609 B, HD 106906 b, and PDS 70 b

There is also evidence for disks around 1RXS J1609 B, HD 106906 b, and PDS 70 b (see Table 1), but we do not detect intrinsic polarization from these companions (see Sect. 5.4). It could be that these companions host a disk but that the properties and geometry of these disks is such that they do not produce a measurable degree of polarization. However, for 1RXS J1609 B no hydrogen emission lines are detected. Instead, the evidence for the existence of a disk is based on red near-infrared colors, weak mid-infrared excess that is spatially unresolved between the primary star and the companion, and a moderate extinction (Bailey et al. 2013; Wu et al. 2015). Because we find that the companion is polarized by interstellar dust (see Sect. 5.3), it seems more likely that these properties are caused by interstellar dust rather than a circumsubstellar disk.

As discussed in Sect. 5.4, we placed an upper limit of 0.2% on the degree of polarization of HD 106906 b, with a 68.27% confidence level. Because also no hydrogen emission lines are detected for this companion, a possible explanation for the non-detection is that the companion simply does not host a circumsubstellar disk. In the case of PDS 70 b we do not reach a very high sensitivity and placed a 68.27% upper limit of 5.0% on the degree of polarization in K_s -band. Therefore, we can conclude that if PDS 70 b hosts a disk, the inclination is probably not so high that it completely obscures the planet’s atmosphere. Because we only detected polarization for companions with hydrogen emission lines, it seems that these lines are the best non-polarimetric indicators for the existence of a circumsubstellar disk.

7.3 Atmospheric asymmetries of the companions

Of the 18 companions for which we do not detect significant polarization, 14 show no clear evidence of hosting a circumsubstellar disk (see Fig. 1 and Table 1). Because the majority of the companions have late-M to mid-L spectral types and are therefore expected to have dusty atmospheres, we could expect to detect polarization due to rotation-induced oblateness or an inhomogeneous cloud distribution. Indeed, polarization between several tenths of a percent to a percent has been detected at near-infrared wavelengths (in particular in J-, Z-, and I-band) for more than a dozen late-M to mid-L field brown dwarfs (Miles-Páez et al. 2013, 2017). In our survey, we reached sensitivities (upper limits) ≤0.3% for 11 companions (see Sect. 5.4 and Table 5), and so we might have expected to detect polarization for a few of the companions. Because we do not detect polarization due to atmospheric asymmetries for any of the companions, these asymmetries either do not exists for the companions observed or they produce polarization below the sensitivity reached.

In the majority of cases, the polarization of field brown dwarfs is interpreted to be caused by rotation-induced oblateness. In that sense our non-detections are particularly surprising because the companions observed are generally young and have a low surface gravity (see Table 1), which should result in a more oblate atmosphere for a given rotation rate and therefore more polarization (Sengupta & Marley 2010; Marley & Sengupta 2011). It is important to note, however, that the field brown dwarfs observed by Miles-Páez et al. (2013) are old (ages 0.5–5 Gyr) and have measured projected rotational velocities v sin i > 30 km s⁻¹. Indeed, in their sample of several dozen field brown dwarfs, Zapatero Osorio et al. (2006) found that about half of the very young field brown dwarfs (1–10 Myr) have v sin i ≤ 10 km s⁻¹ whereas all old brown dwarfs (≥1 Gyr) have v sin i ≈ 30 km s⁻¹. Very young brown dwarfs rotate slowly because they are still in the process of spinning up as they cool and contract. Looking at Fig. 1, we can divide our sample roughly into a large group of young, high-temperature companions with late-M to mid-L spectral types, and a smaller group of older, lower-temperature companions of mid-L to T spectral types. A possible explanation for our non-detections is that while the companions of the first group may have dusty atmospheres, they rotate too slowly to produce a measurable level of polarization. And on the other hand the second group may rotate faster, but due to their later spectral types their upper atmospheres may lack the scattering dust to produce polarization (Allard et al. 2001; Sengupta & Marley 2009). A more in-depth analysis of the degrees of polarization produced due to rotation-induced oblateness is presented in Jensen-Clem et al. (2020).

There could also be other explanations for our non-detections. It could be that the dust grains in the upper atmosphere are submicron sized, as also suggested by studies of the emission spectra of (field) brown dwarfs and planets (Hiranaka et al. 2016; Bonnefoy et al. 2016), and that we therefore need to observe at shorter wavelengths than H-band (i.e., in Y - or J-band). Miles-Páez et al. (2017) observed one of the field brown dwarfs in Z-, J-, and H-band and found that the degree of polarization decreases strongly with increasing wavelength, with the maximum polarization in Z-band and no detection in H-band. Another possibility, as suggested by Miles-Páez et al. (2017), is that the low-gravity atmospheres of young objects have thicker dust clouds, resulting in strong multiple scattering and therefore a low integrated degree of polarization. Finally, our non-detections may also indicate that the atmospheric dust clouds are homogeneously distributed, or that the inhomogeneities do not produce a measurable degree of polarization. Indeed, Millar-Blanchaer et al. (2020) recently detected polarization that is likely due to cloud banding on Luhman 16 A, but the measured degree of polarization is only 0.03% in H-band.

7.4 Confirmation and further characterization of the disks of DH Tau B and GSC 6214 B

To confirm that the polarization from DH Tau B and GSC 6214 B is truly intrinsic and not caused by interstellar dust, we should perform follow-up observations. For this we can use the recently implemented star-hopping technique for SPHERE-IRDIS to quasi-simultaneously measure the stellar polarization from nearby diskless reference stars. As a reference for DH Tau we can observe DI Tau, a very close neighbor to DH Tau located at a separation of only 15.3′′ and at the same distance from Earth, and whose spectral-energy distribution (classified as class III; Luhman et al. 2010) and very low mass-accretion rate (Alonso-Martínez et al. 2017) indicate it very likely does not host a circumstellar disk that creates significant intrinsic stellar polarization. For GSC 6214 we can observe BD-20 4481, which is of similar spectral type as GSC 6214 A and is located at a separation of 13.3′ and at a distance of 113 pc (compared to 109 pc for GSC 6214). We can use the measurements of the stellar polarization of the reference stars to subtract the interstellar component of the companions’ polarization, and with that accurately determine the intrinsic polarization of the companions. We can also measure the polarization of DH Tau A and GSC 6214 A with a different instrument than SPHERE, for example with the WIRC+Pol near-infrared spectropolarimeter (Tinyanont et al. 2019) on the Hale Telescope at Palomar Observatory. Using WIRC+Pol we can determine the stellar polarization as a function of wavelength, enabling us to quantify the interstellar polarization by comparing the measurements with Serkowski’s law of interstellar polarization.

We can further characterize the circumsubstellar disks of DH Tau B and GSC 6214 B, as well as the companions themselves, by performing measurements with various current and future instruments. We can perform IRDIS polarimetric measurements in multiple filters to constrain the distribution of particle sizes in the disks. By combining these measurements with optical and near-infrared medium-resolution spectroscopy using MUSE and ERIS on the VLT, we can constrain the fundamental parameters of the companions. If we are able to detect the disks with ALMA, we can derive their dust mass from the continuum emission, their gas mass from the CO spectral-line emission, and their effective temperature from the emission in two different wavelength bands. Similar to ALMA observations, the dust mass and effective temperature of the disk can be determined from mid-infrared photometric and spectroscopic observations, for example with MIRI on board the James Webb Space Telescope (JWST), METIS on the Extremely Large Telescope (ELT), or even VISIR on the VLT. Finally, with the sensitivity of MIRI and METIS we could detect silicate emission features at 10 and 18 μm.

8 Summary and conclusions

We measured the near-infrared linear polarization of 20 young planets and brown dwarf companions using SPHERE-IRDIS. We reduced the data using the IRDAP pipeline to correct for the instrumental polarization and crosstalk of the optical system with an absolute polarimetric accuracy < 0.1% in the degree of polarization.To retrieve the polarization of the companions, we employed a combination of aperture photometry, ADI, and PSF fitting. We achieved a best 1σ polarimetric contrast of 3 × 10⁻⁸ at an angular separation of 0.5′′ from the star and a contrast <10⁻⁸ for separations > 2.0′′.

We report the first detection of polarization originating from substellar companions, with a measured degree of polarization of several tenths of a percent for DH Tau B and GSC 6214 B in H-band. By comparing the measured polarization with that of nearby stars, we find that this polarization is unlikely to be caused by interstellar dust. Because the companions have previously measured hydrogen emission lines and red colors, we conclude that the polarization most likely originates from circumsubstellar accretion disks. Through radiative transfer modeling we constrain the position angles of the disks and find that the disks must have high inclinations to produce these measurable levels of polarization. For GQ Lup B, which has stronger hydrogen emission lines than DH Tau B and GSC 6214 B, we do not measure significant polarization. This implies that if GQ Lup B hosts a disk, this disk has a low inclination. Assuming that the spin axes of the companions are perpendicular to the plane of their disks, we use previously measured projected rotational velocities to constrain the rotational periods of DH Tau B, GSC 6214 B, and GQ Lup B to be 29–37 h, 22–77 h, and 19–64 h, respectively, within the 68% confidence interval. Finally, we find 1RXS J1609 B to be marginally polarized by interstellar dust, which suggests that the red colors and extinction that are thought to indicate the presence of a disk are more likely caused by interstellar dust.

The disk of DH Tau B, and possibly that of GQ Lup B, are misaligned with the disks around the primary stars. These misalignments suggest that these wide-separation companions have formed in situ through direct collapse in the molecular cloud, although formation through gravitational instabilities in the circumstellar disk cannot be excluded. Formation at close separations from the star followed by scattering to a higher orbit is unlikely because dynamical encounters with other bodies would most likely severely disturb or even destroy the circumsubstellar disks.

For 18 companions we do not detect significant polarization and place upper limits of typically < 0.3% on their degree of polarization. These non-detections may indicate that young companions rotate too slowly to produce measurable polarization due to rotation-induced oblateness, or that any inhomogeneities in the atmospheric clouds are limited. Another possibility is that the upper atmospheres of the companions contain primarily submicron-sized dust grains. This implies that we should perform future measurements in Y - or J-band, although these bands are more challenging in terms of companion-to-star contrast and contrast performance of the instrument.

In our survey, we also detected the circumstellar disks of DH Tau, GQ Lup, PDS 70, β Pic, and HD 106906, which for DH Tau and GQ Lup are the first disk detections in scattered light. The disk of DH Tau is compact and has a strong brightness asymmetry that may reveal the forward- and backward-scattering sides of the disk or may be caused by shadowing by an unresolved inner disk component. The disk of GQ Lup shows a pronounced asymmetry and two spiral-like features that could be the result of periodic close passes of GQ Lup B. The PDS 70 disk shows significant non-azimuthal polarization indicating multiple scattering. We also detect one or two weak spiral-like features protruding from the ansae of the disk that may be the result of two spiral arms in the outer disk ring, potentially induced by PDS 70 b.

Our measurements of the polarization of companions are reaching the limits of the instrument and the data-processing techniques. We find that incompletely corrected bad pixels can cause systematic errors of several tenths of a percent in the measured polarization. To minimize this effect, we recommend to use the field-tracking mode without dithering for future observations that aim to measure the polarization of companions. However, for companions at close separations or with large star-to-companion contrasts, pupil-tracking observations are still preferred to retrieve the companions’ total intensity with ADI. These close-in companions can alternatively be observed in field-tracking mode when using the recently implemented star-hopping technique to enable reference star differential imaging. We also find that the measurements of the stellar polarization are affected by systematic errors related to the use of the coronagraph in combination with time-varying atmosphericconditions. We therefore recommend to take additional stellar polarization measurements without coronagraph.

To further characterize the circumsubstellar disks of DH Tau B and GSC 6214 B, as well as the companions themselves, we can perform follow-up observations with SPHERE-IRDIS, ALMA, JWST-MIRI, MUSE and ERIS on the VLT, and METIS on the ELT. Our polarimetric detections of the disks of DH Tau B and GSC 6214 B are a first step in building a complete picture of the companions, their formation, and evolution, and pave the way to detecting polarization of young planets with for example SPHERE+ (Boccaletti et al. 2020) and the future planet-characterization instrument EPICS (or PCS) on the ELT.

Acknowledgements

R.Gv.H. thanks ESO for the studentship at ESO Santiago during which part of this project was performed. T.S. acknowledges the support from the ETH Zurich Postdoctoral Fellowship Program. The research of F.S., J.dB., and A.J.B. leading to these results has received funding from the European Research Council under ERC Starting Grant agreement 678194 (FALCONER). C.P. acknowledges financial support from Fondecyt (grant 3190691) and from the ICM (Iniciativa Científica Milenio) via the Núcleo Milenio de Formación Planetaria grant, from the Universidad de Valparaíso. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France; NASA’s Astrophysics Data System Bibliographic Services; Scipy, a free and open-source Python library used for scientific computing and technical computing (Virtanen et al. 2020); and Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

Appendix A Cosmetic correction of spurious structure in Q- and U-images

If a companion is polarized, we would expect the polarization signals in the Q- and U-images to resemble scaled-down positive or negative versions of the corresponding total-intensity images I_Q and I_U . However, for many data sets the Q- and U-images show spurious structure with adjacent positive and negative signals. For example, for the 2019-10-24 data set of DH Tau, as shown in Fig. A.1 (first column), we see that the Q-image contains positive and negative signals and that the signal in U is offset from the center coordinates of the companion. These spurious structures result from imperfect relative centering of the images of IRDIS’ left and right optical channels and image motion during the observations. For pupil-tracking observations the spurious structures can additionally originate from image rotation between the two measurements of the double difference.

Fig. A.1 Reduced Q- and U-images (first column), model images of the spurious structures in Q and U (second column), and Q- and U-images corrected for the spurious structures (third column) at the position of the companion DH Tau B of the 2019-10-24 data set of DH Tau. An aperture of radius 8 pixels centered on the companion is shown superimposed on the images.

In the case these spurious structures are visible in the Q- and U-images of a data set, we make to each image individually a (cosmetic) correction similar to that described in Snik et al. (2010). For this we retrieve a positive and negative copy of the I_Q- or I_U-image at the companion position and create a model image in which the two copies are symmetrically shifted in opposite directions from the center coordinates of the companion. We then subtract this model image from the Q- or U-image and fit the shifts in the x- and y-directions by minimizing the sum of squared residuals in an aperture of radius 8 pixels in the resulting image. Because the aperture sum in the model images is zero, subtracting these images only suppresses the spurious structures and does not alter the net polarization signals in Q and U.

For the data set of DH Tau, we find a total relative shift equal to 0.015 pixels for Q and 0.013 pixels for U. Only small relative shifts are needed because the maximum values of the total-intensity PSFs are more than 100 times larger than the maximum values of the positive and negative signals of the spurious structure. The model images and the corrected Q and U-images are shown in Fig. A.1 (second and third column). The spurious structure has clearly disappeared in the corrected images, with the Q-image only having negative signal and the signal in U being positioned at the companion’s center coordinates.

Appendix B Systematic errors due to bad pixels

A few percent of the pixels of the IRDIS Hawaii 2RG detector are bad, that is, they are dead, nonlinearly responding, or hot pixels. When preprocessing the raw frames with IRDAP, bad pixels are identified with a bad pixel map followed by sigma-filtering and then replaced by the median value of the surrounding pixels. These data-reduction steps correct the majority of the bad pixels, but some bad pixels remain uncorrected or are replaced by a value that is not accurate. This results in systematic errors of the pixel values. Whereas these small errors are not a real problem for photometry of point sources in total intensity or imaging of circumstellar disks in polarized light, they become quite problematic when trying to measure the polarization of point sources at a level as low as a few tenths of a percent of the total intensity.

Incomplete correction of bad pixels only marginally affects data taken in field-tracking mode. This is because in field-tracking mode the PSF of the companion is approximately stationary on the detector and only moves very slightly due to variations in AO performance. The bad pixels, which are at a fixed position on the detector, are therefore replaced by approximately the same (median) value in consecutive frames, and so are strongly suppressed when computing the double difference. In addition, any uncorrected or incompletely corrected bad pixels that remain are further averaged out over the various HWP cycles. However, this averaging over HWP cycles only partially applies for our data because we generally observed in field-tracking mode with dithering in which case the detector moves by one to a few pixels each HWP cycle. We note that for total-intensity imaging, for which we compute the median over many exposures rather than differences of exposures, dithering does help suppress bad pixels.

Data taken in pupil-tracking mode are typically more strongly affected by incomplete correction of bad pixels. In pupil-tracking observations the companion moves over the detector, and so in each frame the bad pixels are at a different location with respectto the companion PSF. Therefore, the bad pixels are replaced by very different median values, and relatively large systematic errors remain after the double difference. For data sets with fast parallactic rotation (e.g., the data sets of GSC 8047 and TYC 8998), the bad pixels are more problematic than for data sets with only little rotation (e.g., GSC 6214). For data sets with many HWP cycles the bad pixel effect averages out somewhat, but the systematic errors are still much larger than for field-tracking data.

We attempted to remove the systematic errors due to bad pixels by creating a more aggressive bad pixel map from the dark and flat frames, performing aggressive sigma-filtering, locally replacing the bad-pixel values with cubic spline interpolation rather than with the median filter, and computing the median over the Mueller-matrix-model-corrected Q- and U-images of each HWP cycle. Unfortunately, we were not able to identify all bad pixels in the data and completely remove the effect. This is primarily because part of the bad pixels cause systematic errors of only several percent or less of the total intensity. Such small deviations from the true value are practically impossible to detect in the images and only become evident when computing differences of images as we do in polarimetry.

Although we were not able to completely correct for the bad pixels, we can mitigate their effect by excluding those frames that contribute strong bad pixels to the final images or that show bad pixels at the position of the companion in the bad pixel map. In addition we can average out the systematic error due the bad pixels by mean-combining several data sets. We can also uselarge apertures to perform the photometry with, such that the bad pixels values (which are both positive andnegative in polarimetric images) average out and sum to a lower spurious signal. Future observations aimed at measuring the polarization of companions should preferably be performed in field-tracking mode without dithering.

Appendix C Retrieval of total intensity through ADI: upper limit on polarization of β Pic b

Companionsat small separations or at large star-to-companion contrasts are swamped in the halo of starlight in total intensity. For data sets that were taken in pupil-tracking mode and have sufficient parallactic rotation, we have therefore slightly adapted themethod described in Sect. 4 and determine the probability distribution of the total intensity of the companion by performing ADI with negative PSF injection. We still determine the distributions of Q and U using aperture photometry because the stellar speckle halo is almost completely removed in the polarimetricdata-reduction steps, in particular for the reductions with the added classical ADI step (see Sect. 3). We applied this adapted method to the data sets of HR 8799, HD 206893, PDS 70, and β Pic. In this section we demonstrate the method with the 2019-11-26 H-band data set of β Pic and show how we calculate upper limits on the degree of polarization of the companion β Pic b. A total-intensity image of the data after applying classical ADI with IRDAP is shown in Fig. C.1 (left).

Fig. C.1 Reduced images of the 2019-11-26 data set of β Pic. Left: total-intensity image after applying classical ADI with IRDAP to reveal the companion β Pic b. Right: Q-image after combining polarimetry with classical ADI, showing the aperture of radius 1.85 pixels at the position of the companion (yellow) and the ring of comparison apertures of the same radius around the star (black).

To perform the ADI with negative PSF injection, we use the PynPoint⁶ pipeline, version 0.8.2 (Amara & Quanz 2012; Stolker et al. 2019), and closely follow the steps described in Stolker et al. (2020). In short, we fetch the preprocessed science frames and the stellar PSF image from the reduction with IRDAP in Sect. 3. Subsequently, we iteratively subtract scaled copies of the stellar PSF from the preprocessed science frames at the approximate position of the companion and apply ADI with PCA (in this case subtracting three principal components) to minimize the residuals at that same location. Using Markov chain Monte Carlo (MCMC) we then sample the posterior distributions of the companion’s angular separation, position angle and contrast with respect to the star. We take the median of the posterior distribution of the contrast as the final contrast value and determine the corresponding statistical uncertainties from the 16th and 84th percentiles. We also estimate the systematic uncertainty on the contrast by injecting fake companions at various position angles (but the same separation and contrast as the real companion), retrieving them, and computing the distribution of the difference between the retrieved and injected contrasts. This systematic uncertainty accounts for the azimuthal variations of the noise around the central star and is generally 1 to 5 times larger than the statistical uncertainty (similar to the results of Wertz et al. 2017). Finally, we compute the overall uncertainty as the quadratic sum of the statistical and systematic uncertainties. After these steps, we determine the probability distribution of the companion’s total intensity (expressed in counts) for a range of aperture radii from 1 to 10 pixels. To this end, we draw 10⁶ samples from a Gaussian distribution with the mean and standard deviation equal to the companion-to-star contrast and uncertainty we retrieved with PynPoint. We then sum the flux in the stellar PSF image using an aperture of the given radius and multiply the Gaussian samples by this summed flux. The resulting total-intensity distribution of the companion, which we use for both I_Q and I_U , is shown in Fig. C.2 (left) for an aperture radius of 1.85 pixels. This radius is the final aperture radius we select at the end of this section and corresponds to half times the full width at half maximum (FWHM) we measure on the stellar PSF.

To determine the probability distributions in Q and U, we use the images from the reduction with the added classical ADI step. In the case of β Pic, the classicalADI step does not only further suppress the speckle noise, but also removes most of the signal from the nearly edge-on-viewed circumstellar disk that crosses the position of β Pic b (see Fig. 12, center and bottom). The classical ADI step suppresses the disk signal because the disk is broad and the parallactic rotation of the observations in only 19.8°. Indeed, as can be seen in Fig. C.1 (right) for the Q-image, the disk is almost completely removed and there are only few speckles left at the separation of the companion. Any polarization originating from the companion should still be visible in the Q- and U-images because point sources are much less affected by ADI-induced self-subtraction.

We analyze the Q- and U-images by followingthe exact steps as described in Sect. 4, but with one exception. Before performing the aperture photometry, we quantify the throughput of the ADI procedure by performing a simulation in which we inject and retrieve an artificial source at the separation of the actual companion. We then correct the Q- and U-images for the self-subtraction by dividing them by the calculated throughput, which for this data set is 49%. After performing all the steps, we determine the companion’s polarization for each of the defined aperture radii (as in Fig. 7 for DH Tau B). Contrary to the data of DH Tau, for this data set of β Pic we detect no signals with an S/Nhigher than 0.9 in Q and 1.9 in U for any aperture radius. Indeed, the reduced Q-image (see Fig. C.1, right) and U-image contain only noise at the position of the companion. We therefore conclude that we do not detect significant polarization originating from the companion β Pic b.

We now proceed to set limits on the degree of polarization of the companion. To this end, we determine, for each defined aperture radius, two upper limits from the probability distribution of the degree of polarization. We compute these upper limits from the 68.27 and 99.73% intervals, which for a Gaussian distribution would correspond to the 1σ and 3σ confidence intervals, respectively. These intervals are calculated one-sided and starting at zero because the degree of linear polarization is computed as P = √(q² + u²) and therefore can only have positive values (see Sparks & Axon 1999). Figure C.2 (right) shows the distribution of the degree of polarization for an aperture radius of 1.85 pixels with the two intervals indicated. Figure C.3 shows the two upper limits as a function of aperture radius. From this figure we see that the upper limits are relatively constant for an aperture radius smaller than approximately 3.5 pixels. For larger apertures, the upper limits increase as more noise is included in the companion aperture and the uncertainty of the background due to the low number of background samples increases. We select our final aperture radius to be 1.85 pixels, equal to half times the FWHM of the stellar PSF, and conclude that the 68.27 and 99.73% upper limits on the degree of polarization of β Pic b are equal to 0.2 and 0.4%, respectively. We note that while for this particular data set the upper limits monotonically increase with aperture radius, for several other data sets this is not the case due to incompletely removed bad pixels (see Sect. 5.4 and Appendix B).

Fig. C.2 Final probability distributions of the signal of β Pic b from the 2019-11-26 data set of β Pic, using an aperture radius of 1.85 pixels. Left: probability distribution of the total intensity. The mean and 68.27% uncertainties of the distribution are shown above the graph, with the latter also indicated by the light-blue shaded area. The S/N is shown within parentheses. Right: probability distribution of the degree of linear polarization. The upper limits computed from the one-sided 68.27 and 99.73% intervals are indicated by the light-red and darker red shaded areas, respectively.

Fig. C.3 Upper limits on the degree of linear polarization of β Pic b computed from the one-sided 68.27 and 99.73% intervals as a function of aperture radius from the 2019-11-26 data set of β Pic. The final selected aperture radius of 1.85 pixels, equal to half times the FWHM of the stellar PSF, is indicated with the dashed vertical line.

Appendix D Retrieval of total intensity through PSF fitting: upper limit on polarization of HD 19467 B

Several observations of faint or close-in companions were not executed in pupil-tracking mode (i.e., they were executed in field-tracking mode) or have little parallactic rotation. For these observations we cannot retrieve the probability distribution of the companion’s total intensity through ADI with negative PSF injection as described in Appendix C. We also cannot use aperture photometry with comparison apertures as outlined in Sect. 4 because the spatially varying stellar halo at the separations of these companions prevents accurate determination of the background. We therefore use MCMC to fit the stellar PSF image to the I_Q - and I_U -images at the companion position and determine the corresponding probability distributions. We applied this method to the data sets of PZ Tel, HR 7682, HD 19467, GQ Lup, and HD 4747. To confirm that the PSF fitting is accurate, we also used the method to retrieve the total intensities of HR 8799 b, c, and d, and find that the results differ only 0.03 to 0.07 mag with those obtained with PynPoint (see Appendix C). In this section we demonstrate the PSF fitting method with the 2018-08-07 H-band data set of HD 19467 and set upper limits on the degree of polarization of the companion HD 19467 B. Figure D.1 shows the I_Q -image of this data set. As the first step in our analysis, we obtain a rough estimate of the companion’s contrast and x- and y-coordinates. To this end,we fit a model consisting of a 2D Moffat function and an inclined plane to the reduced I_Q -image at the companion position. The inclined plane accounts for the (approximately) linearly varying local background due to the stellar PSF and the stellar diffraction spikes (see Fig. D.1) and is described by a constant (the z-intercept) and slopes in thex- and y-direction. We then fit a model containing the stellar PSF and an inclined plane to cropped versions of the I_Q and I_U -images, using the results from the Moffat fit for the initial estimates of the fit parameters. We use the Nelder-Mead method as implemented in the Python function scipy.optimize.minimize to minimize the sum of squared residuals (SSR): $$\text{SSR}={\displaystyle \sum _{i=1}^N\left[{\left(I_{Q,i}-{\widehat{I}}_{Q,i}\right)}^2+{\left(I_{U,i}-{\widehat{I}}_{U,i}\right)}^2\right]},$$(D.1)

Fig. D.1 Reduced IQ-image of the 2018-08-07 data set of HD 19467, showing the position of the companion HD 19467 B in the white circle. The asymmetric wind-driven halo and the stellar diffraction spikes are clearly visible.

where I_Q,i and I_U,i are the flux values in the i-th pixel of the cropped I_Q- and I_U-images, Î _Q,i and Î _U,i are the corresponding modeled flux values, and N is the total number of pixels in each of the cropped images. We minimize the residuals in the I_Q - and I_U -images simultaneously to obtain a single set of x- and y-coordinates for the companion position. For the other parameters we fit separate values for I_Q and I_U .

We nowrepeat the PSF fitting using MCMC to obtain the final values of the fit parameters and the corresponding posterior distributions. We use the MCMC sampler from the Python package emcee (Foreman-Mackey et al. 2013) and let 32 walkers explore the probability space with 20 000 steps each (resulting in a total of 640 000 samples). We randomly generate the starting values of the walkers from Gaussian distributions centered on the best-fit parameter values from our previous fit. We use a Gaussian distribution for the log-likelihood function: $$\ln \mathcal{L}\propto -\frac{1}{2}[N\ln \left(2\pi {\sigma }^2\right)+\frac{\text{SSR}}{{\sigma }^2}],$$(D.2)

where SSR iscomputed from Eq. (D.1) and σ is the standard deviation that accounts for the noise in the images. Because there is no region in the I_Q and I_U -images from which we can determine a representative value for σ, we include it among the parameters to be fitted (i.e., we treat σ as a nuisance parameter). We set the prior for σ proportionalto 1∕σ, that is, Jeffrey’s prior, to make sure it is non-informative. For the other parameters we use uniform priors. We remove the first 822 steps of each walker, equal to five times the maximum autocorrelation time, and check by visual inspection that the chains of all parameters have converged. The cropped I_Q - and I_U -images and the best-fit model and residual images are shown in Fig. D.2. Figure D.3 shows the resulting 1D- and 2D-projections of the posterior distribution of the fitted parameters. The distributions in Fig. D.3 are visually very close to being Gaussian and show correlations only between the companion’s contrast in I_Q or I_U and the corresponding z-intercept of the background.

Fig. D.2 Data, best-fit model and residual images of the MCMC fitting of the stellar PSF at position of HD 19467 B to the reduced IQ- and IU -images of the 2018-08-07 data set of HD 19467.

Fig. D.3 Posterior distributions after using MCMC to fit the stellar PSF at the position of the companion HD 19467 B to the reduced IQ- and IU -images of the 2018-08-07 data set of HD 19467. The fitted parameters are the companion position in x and y, the companion-to-star contrast in IQ and IU , the background’s z-intercept and slopes inthe x- and y-direction in IQ and IU , and the standard deviation σ that accounts for the noise in the images. The diagonal panels show the marginalized 1D distributions of the fitted parameters and the off-diagonal panels show the 2D projections of the posterior, revealing the covariance of the parameter pairs. The median and uncertainties (computed as the 18th and84th percentiles) of the distributions are shown above the histograms and are indicated with the dashed vertical lines. The contours superimposed on the off-diagonal panels indicate the 1σ, 2σ and 3σ confidence levels assuming Gaussian statistics. The figure is created using the Python package corner (Foreman-Mackey 2016).

We now determine the companion’s probability distributions in I_Q and I_U (expressed in counts) for a range of aperture radii from 1 to 10 pixels. Similarly to the method described in Appendix C, we sum the flux in the stellar PSF image using an aperture of the given radius and multiply the MCMC contrast samples in I_Q and I_U by this flux. For the remainder of the analysis we follow the steps described in Sects. 4 and Appendix C, with the only exception that we sample the PDFs in Q and U with the same number of samples as used for the MCMC analysis. After performing the complete analysis, we detect no signals with an S/N higher than 1.4 in Q and 2 in U for any aperture radius. Finally, using an aperture radius of 1.86 pixels, equal to half times the FWHM of the stellar PSF, we determine the 68.27 and 99.73% upper limits on the degree of polarization of HD 19467 B to be equal to 1.0 and 2.0%, respectively.

Appendix E Contrastcurve of β Pic data

Figure E.1 shows the 1σ and 5σ point-source contrast in Q and U as a function of angular separation from the star for the mean-combined data set of β Pic as constructed with IRDAP. The curves are computed by summing the flux in rings of apertures around the star, computing the standard deviation over the aperture sums, and normalizing the result with the total stellar flux retrieved from the star flux frames. At small separations the correction for small-sample statistics is applied (see Mawet et al. 2014). For comparison the figure also shows the azimuthally averaged flux in the total-intensity I_Q - and I_U -images and the corresponding photon noise. At angular separations between ~ 0.2′′ and 2.0′′ the polarimetric sensitivity is close to the photon-noise limit, with a 1σ-contrast of 7 × 10⁻⁸ to 1 × 10⁻⁸ and a 5σ-contrast of 5 × 10⁻⁷ to 5 × 10⁻⁸. At separations larger than 2.0′′ the sensitivity is limited by read noise or background noise and the 1σ- and 5σ-contrast are <1 × 10⁻⁸ and <5 × 10⁻⁸, respectively.

Fig. E.1 1σ and 5σ point-source contrast in Q and U as a function of angular separation from the star for the mean-combined data set of β Pic. The azimuthally averaged flux in the total-intensity IQ- and IU -images and the corresponding photon noise are shown for comparison.

References

All Tables

All Figures

	Fig. 1 Properties of the companions of our sample showing the age, spectral type, mass (surface area of data points), effective temperature, and possible existence of a circumsubstellar disk (thick border). The data points of HR 8799 b and c, and of HR 8799 d and e, overlap.
In the text

	Fig. 2 Reduced IQ-, Q-, and U-images (after applying the cosmetic correction described in Appendix A) at the position of the companion DH Tau B of the 2019-10-24 data set of DH Tau, showing an aperture of radius 8 pixels centered on the companion. The IU -image, which is not shown, is very similar to the IQ-image.
In the text

	Fig. 3 Reduced IQ-image of the 2019-10-24 data set of DH Tau, showing an aperture of radius 8 pixels at the position of the companion DH Tau B (red) and the ring of comparison apertures of the same radius around the star (white).
In the text

	Fig. 4 PDF of the signal of DH Tau B in IQ (left) and Q (right) fromthe 2019-10-24 data set of DH Tau, using an aperture radius of 8 pixels and taking into account only the photon noise of the companion. The mean and standard deviation of the distributions are shown above the graphs, with the latter also indicated by the light-blue shaded area.
In the text

	Fig. 5 Histograms of the background in IQ, IU , Q, and U of the 2019-10-24 data set of DH Tau, as obtained through summing the flux in the 8-pixel-radius comparison apertures of Fig. 3. The mean and standard deviation of the samples are shown above the histograms. The blue curves show the PDFs as estimated through KDE and the red curves show the best-fit Gaussian distributions for comparison.
In the text

Fig. 6 Final probability distributions of the signals of DH Tau B in IQ, IU , Q, and U (top row), and in normalized Stokes q and u, degree of linear polarization, and angle of linear polarization (bottom row) from the 2019-10-24 data set of DH Tau, using an aperture radius of 8 pixels. The median values of the distributions, as well as the uncertainties computed from the two-sided 68.27% equal-tailed interval around the median, are shown above the graphs. The S/N, i.e., the median value divided by the largest uncertainty, is shown within parentheses. The 68.27% intervals are also indicated by the light-blue and light-red shaded areas.

In the text

Fig. 7 Normalized Stokes parameters q and u and degree of linear polarization (top), angle of linear polarization (center), and S/N in q, u, and the degree of linear polarization (bottom) of DH Tau B as a function of aperture radius for the 2019-10-24 data set of DH Tau.The uncertainties of the measured values are shown with error bars. The final selected aperture radius of 8 pixels is indicated with the dashed vertical lines.

In the text

Fig. 8 Map of the linear polarization of DH Tau B and nearby stars superimposed on a Herschel-SPIRE map at 350 μm. The length and orientation of the lines indicate the degree and angle of linear polarization, respectively. The black line shows the H-band polarization we measure for DH Tau B, and the orange lines display the SPHERE-IRDIS H-band measurements of DH Tau A and three other nearby T Tauri stars whose archival data we analyzed. White lines show optical measurements by Heyer et al. (1987). Yellow lines indicate the H-band polarization of three bright stars closest to DH Tau as derived from optical measurements by Moneti et al. (1984). The length of the H-band vectors are scaled by a factor of four with respect to the optical vectors.

In the text

Fig. 9 Intrinsic polarization of DH Tau B after subtracting interstellar polarization signals from the measured polarization of the companion. Top: intrinsic degree and angle of linear polarization of DH Tau B as a function of the degree of polarization due to interstellar dust, assuming an angle of 26° for the interstellar polarization. The bands around the curves show the uncertainties of our measurements. Bottom: probability distributions of the intrinsic degree of polarization of DH Tau B for a range of degrees of polarization due to interstellar dust, assuming the angle of the interstellar polarization to have the same distribution as that determined by Goodman et al. (1992) for the B216-B217 dark cloud adjacent to DH Tau. The probability distribution of each column is normalized to one.

In the text

	Fig. 10 Reduced mean-combined IQ-, Q-, and U-images (after applying the cosmetic correction described in Appendix A) at the position of the companion GSC 6214 B, showing an aperture of radius 4 pixels centered on the companion. The IU -image, which is not shown, is very similar to the IQ-image.
In the text

Fig. 11 Map of the angle of linear polarization of the companion GSC 6214 B, the star 1RXS J1609 A, and other nearby bright stars superimposed on an IRAS map at 100 μm. The angles of GSC 6214 B and 1RXS J1609 A are from the SPHERE-IRDIS H-band measurements from this work, whereas for the other stars the angles are taken from the catalog of optical measurements by Heiles (2000). The length of the lines is arbitrary and contrary to Fig. 8 does not indicate the degree of polarization. HD 147700 in unpolarized. White lines indicate stars at a distance between 128 pc and 142 pc, and blue and orange lines show objects closer or farther away, respectively. We note that the region shown is much larger than that ofFig. 8, and therefore the angular separation among the stars is much larger as well.

In the text

Fig. 12 H-band images of the five circumstellar disks detected in our sample, showing the linearly polarized intensity (PI) for DH Tau and GQ Lup, Qϕ and Uϕ for PDS 70, and Qϕ for HD 106906 and β Pic. A total-intensity image after applying ADI with PCA (subtracting four principal components) is additionally shown for β Pic. The Uϕ -image of PDS 70 is shown on a linear scale, whereas all other images are shown on a logarithmic scale. All polarimetric images have the polarized stellar halo subtracted. The angular scale and sky orientation are indicated in each image. Gray circles maskregions obscured by the coronagraph.

In the text

Fig. 13 Synthetic image of a self-luminous companion (Teff = 2000 K) with a compact circumsubstellar disk at a distance of 150 pc. The total intensity surface brightness is shownon a logarithmic scale and the lines indicate the local degree and angle of linear polarization. The spatially integrated degree of polarization is 0.24% and the angle of polarization is 0°, that is, the spatially unresolved light is linearly polarized along the minor axis of the disk.

In the text

	Fig. 14 Dependence of the integrated degree of linear polarization on the inner radius of the circumsubstellar disk and the surface density of the dust at the inner radius. The grid of radiative transfer models is shown for a disk inclination of 70° (top) and 80° (bottom).
In the text

	Fig. A.1 Reduced Q- and U-images (first column), model images of the spurious structures in Q and U (second column), and Q- and U-images corrected for the spurious structures (third column) at the position of the companion DH Tau B of the 2019-10-24 data set of DH Tau. An aperture of radius 8 pixels centered on the companion is shown superimposed on the images.
In the text

	Fig. C.1 Reduced images of the 2019-11-26 data set of β Pic. Left: total-intensity image after applying classical ADI with IRDAP to reveal the companion β Pic b. Right: Q-image after combining polarimetry with classical ADI, showing the aperture of radius 1.85 pixels at the position of the companion (yellow) and the ring of comparison apertures of the same radius around the star (black).
In the text

Fig. C.2 Final probability distributions of the signal of β Pic b from the 2019-11-26 data set of β Pic, using an aperture radius of 1.85 pixels. Left: probability distribution of the total intensity. The mean and 68.27% uncertainties of the distribution are shown above the graph, with the latter also indicated by the light-blue shaded area. The S/N is shown within parentheses. Right: probability distribution of the degree of linear polarization. The upper limits computed from the one-sided 68.27 and 99.73% intervals are indicated by the light-red and darker red shaded areas, respectively.

In the text

	Fig. C.3 Upper limits on the degree of linear polarization of β Pic b computed from the one-sided 68.27 and 99.73% intervals as a function of aperture radius from the 2019-11-26 data set of β Pic. The final selected aperture radius of 1.85 pixels, equal to half times the FWHM of the stellar PSF, is indicated with the dashed vertical line.
In the text

	Fig. D.1 Reduced IQ-image of the 2018-08-07 data set of HD 19467, showing the position of the companion HD 19467 B in the white circle. The asymmetric wind-driven halo and the stellar diffraction spikes are clearly visible.
In the text

	Fig. D.2 Data, best-fit model and residual images of the MCMC fitting of the stellar PSF at position of HD 19467 B to the reduced IQ- and IU -images of the 2018-08-07 data set of HD 19467.
In the text

Fig. D.3 Posterior distributions after using MCMC to fit the stellar PSF at the position of the companion HD 19467 B to the reduced IQ- and IU -images of the 2018-08-07 data set of HD 19467. The fitted parameters are the companion position in x and y, the companion-to-star contrast in IQ and IU , the background’s z-intercept and slopes inthe x- and y-direction in IQ and IU , and the standard deviation σ that accounts for the noise in the images. The diagonal panels show the marginalized 1D distributions of the fitted parameters and the off-diagonal panels show the 2D projections of the posterior, revealing the covariance of the parameter pairs. The median and uncertainties (computed as the 18th and84th percentiles) of the distributions are shown above the histograms and are indicated with the dashed vertical lines. The contours superimposed on the off-diagonal panels indicate the 1σ, 2σ and 3σ confidence levels assuming Gaussian statistics. The figure is created using the Python package corner (Foreman-Mackey 2016).

In the text

	Fig. E.1 1σ and 5σ point-source contrast in Q and U as a function of angular separation from the star for the mean-combined data set of β Pic. The azimuthally averaged flux in the total-intensity IQ- and IU -images and the corresponding photon noise are shown for comparison.
In the text