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 to even 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, $v\sin{i}$, 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\pm 0.7$ km s-1, $6.1^{+4.9}_{-3.8}$ km s-1, and $5.3^{+0.9}_{-1.0}$ km
s-1, 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) and taking 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\degr$ and $80\degr$ for DH Tau B and GSC
6214 B and between $0\degr$ and $45\degr$ 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 ${\sim}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 ${\sim}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\degr$ and $190\degr$ (see Fig. 9), whereas that of the disk of GSC 6214 B
is around $48\degr$ (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\degr$ and $2.5\degr$, 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\degr$
(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 ${\sim}30\degr$ (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) ${\leq}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-1. 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}\leq 10$ km s-1 whereas
all old brown dwarfs ($\geq 1$ Gyr) have $v\sin{i}\approx 30$ km s-1. 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\arcsec$ 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\arcmin$
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\text{\,}\mathrm{\SIUnitSymbolMicro 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\sigma$ polarimetric contrast of $3\cdot 10^{-8}$ at an angular separation
of $0.5\arcsec$ from the star and a contrast ${<}10^{-8}$ for separations
${>}2.0\arcsec$.

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, $\beta$ 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 atmospheric
conditions. 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.

RGvH thanks ESO for the studentship at ESO Santiago during which part of this
project was performed. TS acknowledges the support from the ETH Zurich
Postdoctoral Fellowship Program. The research of FS, JdB, and AJB leading to
these results has received funding from the European Research Council under
ERC Starting Grant agreement 678194 (FALCONER). CP 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 et al. 2013; Price-Whelan et al. 2018).

## References

  * Ackerman & Marley (2001) Ackerman, A. S. & Marley, M. S. 2001, ApJ, 556, 872
  * Alibert et al. (2005) Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343
  * Allard et al. (2001) Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357
  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765
  * Alonso-Martínez et al. (2017) Alonso-Martínez, M., Riviere-Marichalar, P., Meeus, G., et al. 2017, A&A, 603, A138
  * Alves & Franco (2006) Alves, F. O. & Franco, G. A. P. 2006, MNRAS, 366, 238
  * Amara & Quanz (2012) Amara, A. & Quanz, S. P. 2012, MNRAS, 427, 948
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Ayliffe & Bate (2009) Ayliffe, B. A. & Bate, M. R. 2009, MNRAS, 397, 657
  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58
  * Bailey et al. (2013) Bailey, V., Hinz, P. M., Currie, T., et al. 2013, ApJ, 767, 31
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42
  * Bate (2009) Bate, M. R. 2009, MNRAS, 392, 590
  * Bayo et al. (2017) Bayo, A., Joergens, V., Liu, Y., et al. 2017, ApJ, 841, L11
  * Beuzit et al. (2019) Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155
  * Boccaletti et al. (2003) Boccaletti, A., Chauvin, G., Lagrange, A. M., & Marchis, F. 2003, A&A, 410, 283
  * Boccaletti et al. (2020) Boccaletti, A., Chauvin, G., Mouillet, D., et al. 2020, arXiv e-prints, arXiv:2003.05714
  * Bohn et al. (2020) Bohn, A. J., Kenworthy, M. A., Ginski, C., et al. 2020, MNRAS, 492, 431
  * Bonnefoy et al. (2014) Bonnefoy, M., Chauvin, G., Lagrange, A. M., et al. 2014, A&A, 562, A127
  * Bonnefoy et al. (2010) Bonnefoy, M., Chauvin, G., Rojo, P., et al. 2010, A&A, 512, A52
  * Bonnefoy et al. (2016) Bonnefoy, M., Zurlo, A., Baudino, J. L., et al. 2016, A&A, 587, A58
  * Boss (1997) Boss, A. P. 1997, Science, 276, 1836
  * Bouvier et al. (1999) Bouvier, J., Chelli, A., Allain, S., et al. 1999, A&A, 349, 619
  * Bowler (2016) Bowler, B. P. 2016, PASP, 128, 102001
  * Bowler et al. (2015) Bowler, B. P., Andrews, S. M., Kraus, A. L., et al. 2015, ApJ, 805, L17
  * Bowler et al. (2011) Bowler, B. P., Liu, M. C., Kraus, A. L., Mann, A. W., & Ireland, M. J. 2011, ApJ, 743, 148
  * Bryan et al. (2018) Bryan, M. L., Benneke, B., Knutson, H. A., Batygin, K., & Bowler, B. P. 2018, Nature Astronomy, 2, 138
  * Bryan et al. (2020) Bryan, M. L., Chiang, E., Bowler, B. P., et al. 2020, AJ, 159, 181
  * Caceres et al. (2015) Caceres, C., Hardy, A., Schreiber, M. R., et al. 2015, ApJ, 806, L22
  * Cameron (1978) Cameron, A. G. W. 1978, Moon and Planets, 18, 5
  * Canovas et al. (2015) Canovas, H., Ménard, F., de Boer, J., et al. 2015, A&A, 582, L7
  * Carbillet et al. (2011) Carbillet, M., Bendjoya, P., Abe, L., et al. 2011, Experimental Astronomy, 30, 39
  * Chauvin et al. (2005a) Chauvin, G., Lagrange, A. M., Lacombe, F., et al. 2005a, A&A, 430, 1027
  * Chauvin et al. (2005b) Chauvin, G., Lagrange, A. M., Zuckerman, B., et al. 2005b, A&A, 438, L29
  * Chilcote et al. (2017) Chilcote, J., Pueyo, L., De Rosa, R. J., et al. 2017, AJ, 153, 182
  * Christiaens et al. (2019) Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJ, 877, L33
  * Crepp et al. (2012) Crepp, J. R., Johnson, J. A., Fischer, D. A., et al. 2012, ApJ, 751, 97
  * Crepp et al. (2014) Crepp, J. R., Johnson, J. A., Howard, A. W., et al. 2014, ApJ, 781, 29
  * Crepp et al. (2018) Crepp, J. R., Principe, D. A., Wolff, S., et al. 2018, ApJ, 853, 192
  * Crepp et al. (2015) Crepp, J. R., Rice, E. L., Veicht, A., et al. 2015, ApJ, 798, L43
  * Cuello et al. (2019) Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114
  * Cuello et al. (2020) Cuello, N., Louvet, F., Mentiplay, D., et al. 2020, MNRAS, 491, 504
  * Daemgen et al. (2017) Daemgen, S., Todorov, K., Quanz, S. P., et al. 2017, A&A, 608, A71
  * de Boer et al. (2020) de Boer, J., Langlois, M., van Holstein, R. G., et al. 2020, A&A, 633, A63
  * de Kok et al. (2011) de Kok, R. J., Stam, D. M., & Karalidi, T. 2011, ApJ, 741, 59
  * Delorme et al. (2017) Delorme, P., Schmidt, T., Bonnefoy, M., et al. 2017, A&A, 608, A79
  * Dohlen et al. (2008) Dohlen, K., Langlois, M., Saisse, M., et al. 2008, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, The infra-red dual imaging and spectrograph for SPHERE: design and performance, 70143L
  * Donati et al. (2012) Donati, J. F., Gregory, S. G., Alencar, S. H. P., et al. 2012, MNRAS, 425, 2948
  * Dong et al. (2016) Dong, R., Fung, J., & Chiang, E. 2016, ApJ, 826, 75
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
  * Dupuy et al. (2019) Dupuy, T. J., Brandt, T. D., Kratter, K. M., & Bowler, B. P. 2019, ApJ, 871, L4
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Galli et al. (2018) Galli, P. A. B., Loinard, L., Ortiz-Léon, G. N., et al. 2018, ApJ, 859, 33
  * Ginski et al. (2018) Ginski, C., Benisty, M., van Holstein, R. G., et al. 2018, A&A, 616, A79
  * Ginski et al. (2014) Ginski, C., Mugrauer, M., Neuhäuser, R., & Schmidt, T. O. B. 2014, MNRAS, 438, 1102
  * Goodman et al. (1992) Goodman, A. A., Jones, T. J., Lada, E. A., & Myers, P. C. 1992, ApJ, 399, 108
  * Gravity Collaboration et al. (2019) Gravity Collaboration, Lacour, S., Nowak, M., et al. 2019, A&A, 623, L11
  * Gray et al. (2006a) Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006a, AJ, 132, 161
  * Gray et al. (2006b) Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006b, AJ, 132, 161
  * Gray et al. (2003) Gray, R. O., Corbally, C. J., Garrison, R. F., McFadden, M. T., & Robinson, P. E. 2003, AJ, 126, 2048
  * Grinin (2000) Grinin, V. P. 2000, Astronomical Society of the Pacific Conference Series, Vol. 219, Photopolarimetric Activity of Pre-Main-Sequence Stars (Invited Review), ed. G. Garzón, C. Eiroa, D. de Winter, & T. J. Mahoney, 216
  * Guerri et al. (2011) Guerri, G., Daban, J.-B., Robbe-Dubois, S., et al. 2011, Experimental Astronomy, 30, 59
  * Haffert et al. (2019) Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nature Astronomy, 3, 749
  * Haffert et al. (2020) Haffert, S. Y., van Holstein, R. G., Ginski, C., et al. 2020, A&A, 640, L12
  * Hashimoto et al. (2012) Hashimoto, J., Dong, R., Kudo, T., et al. 2012, ApJ, 758, L19
  * Heiles (2000) Heiles, C. 2000, AJ, 119, 923
  * Herbig (1977) Herbig, G. H. 1977, ApJ, 214, 747
  * Heyer et al. (1987) Heyer, M. H., Vrba, F. J., Snell, R. L., et al. 1987, ApJ, 321, 855
  * Hiranaka et al. (2016) Hiranaka, K., Cruz, K. L., Douglas, S. T., Marley, M. S., & Baldassare, V. F. 2016, ApJ, 830, 96
  * Houk & Cowley (1975) Houk, N. & Cowley, A. P. 1975, University of Michigan Catalogue of two-dimensional spectral types for the HD stars. Volume I.
  * Isella et al. (2019) Isella, A., Benisty, M., Teague, R., et al. 2019, ApJ, 879, L25
  * Jensen-Clem et al. (2016) Jensen-Clem, R., Millar-Blanchaer, M., Mawet, D., et al. 2016, ApJ, 820, 111
  * Jensen-Clem et al. (2020) Jensen-Clem, R., Millar-Blanchaer, M. A., van Holstein, R. G., et al. 2020, AJ, 160, 286
  * Kalas et al. (2015) Kalas, P. G., Rajan, A., Wang, J. J., et al. 2015, ApJ, 814, 32
  * Keppler et al. (2018) Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44
  * Kharchenko & Roeser (2009) Kharchenko, N. V. & Roeser, S. 2009, VizieR Online Data Catalog, I/280B
  * Kratter et al. (2010) Kratter, K. M., Murray-Clay, R. A., & Youdin, A. N. 2010, ApJ, 710, 1375
  * Kraus et al. (2015) Kraus, A. L., Andrews, S. M., Bowler, B. P., et al. 2015, ApJ, 798, L23
  * Kraus et al. (2014) Kraus, A. L., Ireland, M. J., Cieza, L. A., et al. 2014, ApJ, 781, 20
  * Lagrange et al. (2010) Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57
  * Lagrange et al. (2016) Lagrange, A. M., Langlois, M., Gratton, R., et al. 2016, A&A, 586, L8
  * Langlois et al. (2018) Langlois, M., Pohl, A., Lagrange, A. M., et al. 2018, A&A, 614, A88
  * Liu et al. (2002) Liu, M. C., Fischer, D. A., Graham, J. R., et al. 2002, ApJ, 571, 519
  * Lombardi et al. (2008) Lombardi, M., Lada, C. J., & Alves, J. 2008, A&A, 480, 785
  * Luhman et al. (2010) Luhman, K. L., Allen, P. R., Espaillat, C., Hartmann, L., & Calvet, N. 2010, ApJS, 186, 111
  * MacGregor et al. (2017) MacGregor, M. A., Wilner, D. J., Czekala, I., et al. 2017, ApJ, 835, 17
  * Macintosh et al. (2014) Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proceedings of the National Academy of Science, 111, 12661
  * Maire et al. (2016a) Maire, A. L., Bonnefoy, M., Ginski, C., et al. 2016a, A&A, 587, A56
  * Maire et al. (2016b) Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016b, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, SPHERE IRDIS and IFS astrometric strategy and calibration, 990834
  * Mamajek (2008) Mamajek, E. E. 2008, Astronomische Nachrichten, 329, 10
  * Marley & Sengupta (2011) Marley, M. S. & Sengupta, S. 2011, MNRAS, 417, 2874
  * Marois et al. (2006) Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556
  * Martin & Lubow (2011) Martin, R. G. & Lubow, S. H. 2011, MNRAS, 413, 1447
  * Mawet et al. (2014) Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97
  * Ménard et al. (2003) Ménard, F., Bouvier, J., Dougados, C., Mel’nikov, S. Y., & Grankin, K. N. 2003, A&A, 409, 163
  * Ménard et al. (2002) Ménard, F., Delfosse, X., & Monin, J. L. 2002, A&A, 396, L35
  * Miles-Páez et al. (2013) Miles-Páez, P. A., Zapatero Osorio, M. R., Pallé, E., & Peña Ramírez, K. 2013, A&A, 556, A125
  * Miles-Páez et al. (2017) Miles-Páez, P. A., Zapatero Osorio, M. R., Pallé, E., & Peña Ramírez, K. 2017, MNRAS, 466, 3184
  * Millar-Blanchaer et al. (2020) Millar-Blanchaer, M. A., Girard, J. H., Karalidi, T., et al. 2020, ApJ, 894, 42
  * Millar-Blanchaer et al. (2015) Millar-Blanchaer, M. A., Graham, J. R., Pueyo, L., et al. 2015, ApJ, 811, 18
  * Millar-Blanchaer et al. (2016) Millar-Blanchaer, M. A., Perrin, M. D., Hung, L.-W., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Proc. SPIE, 990836
  * Milli et al. (2017) Milli, J., Hibon, P., Christiaens, V., et al. 2017, A&A, 597, L2
  * Min et al. (2009) Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155
  * Min et al. (2016) Min, M., Rab, C., Woitke, P., Dominik, C., & Ménard, F. 2016, A&A, 585, A13
  * Miret-Roig et al. (2020) Miret-Roig, N., Galli, P. A. B., Brandner, W., et al. 2020, A&A, 642, A179
  * Moneti et al. (1984) Moneti, A., Pipher, J. L., Helfer, H. L., McMillan, R. S., & Perry, M. L. 1984, ApJ, 282, 508
  * Montes et al. (2001) Montes, D., López-Santiago, J., Gálvez, M. C., et al. 2001, MNRAS, 328, 45
  * Mora et al. (2020) Mora, A., Wu, Y.-L., Bowler, B. P., & Sheehan, P. 2020, Research Notes of the American Astronomical Society, 4, 9
  * Müller et al. (2018) Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2
  * Neugebauer et al. (1984) Neugebauer, G., Habing, H. J., van Duinen, R., et al. 1984, ApJ, 278, L1
  * Nielsen et al. (2019) Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13
  * Ortiz-León et al. (2017) Ortiz-León, G. N., Loinard, L., Kounkel, M. A., et al. 2017, ApJ, 834, 141
  * Pearce et al. (2019) Pearce, L. A., Kraus, A. L., Dupuy, T. J., et al. 2019, AJ, 157, 71
  * Pecaut & Mamajek (2016) Pecaut, M. J. & Mamajek, E. E. 2016, MNRAS, 461, 794
  * Pérez et al. (2019) Pérez, S., Marino, S., Casassus, S., et al. 2019, MNRAS, 488, 1005
  * Pilbratt et al. (2010) Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123
  * Ricci et al. (2017) Ricci, L., Cazzoletti, P., Czekala, I., et al. 2017, AJ, 154, 24
  * Ricci et al. (2014) Ricci, L., Testi, L., Natta, A., et al. 2014, ApJ, 791, 20
  * Rizzo et al. (1998) Rizzo, J. R., Morras, R., & Arnal, E. M. 1998, MNRAS, 300, 497
  * Rizzuto et al. (2015) Rizzuto, A. C., Ireland, M. J., & Kraus, A. L. 2015, MNRAS, 448, 2737
  * Rodríguez-Ledesma et al. (2009) Rodríguez-Ledesma, M. V., Mundt, R., & Eislöffel, J. 2009, A&A, 502, 883
  * Rostopchina et al. (1997) Rostopchina, A. N., Grinin, V. P., Okazaki, A., et al. 1997, A&A, 327, 145
  * Schneider et al. (2011) Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79
  * Scholz et al. (2015) Scholz, A., Kostov, V., Jayawardhana, R., & Mužić, K. 2015, ApJ, 809, L29
  * Scholz et al. (2018) Scholz, A., Moore, K., Jayawardhana, R., et al. 2018, ApJ, 859, 153
  * Schwarz et al. (2016) Schwarz, H., Ginski, C., de Kok, R. J., et al. 2016, A&A, 593, A74
  * Scott (2015) Scott, D. W. 2015, Multivariate Density Estimation: Theory, Practice, and Visualization
  * Seifahrt et al. (2007) Seifahrt, A., Neuhäuser, R., & Hauschildt, P. H. 2007, A&A, 463, 309
  * Sengupta & Krishan (2001) Sengupta, S. & Krishan, V. 2001, ApJ, 561, L123
  * Sengupta & Marley (2009) Sengupta, S. & Marley, M. S. 2009, ApJ, 707, 716
  * Sengupta & Marley (2010) Sengupta, S. & Marley, M. S. 2010, ApJ, 722, L142
  * Sengupta & Marley (2016) Sengupta, S. & Marley, M. S. 2016, ApJ, 824, 76
  * Serkowski et al. (1975) Serkowski, K., Mathewson, D. S., & Ford, V. L. 1975, ApJ, 196, 261
  * Shabram & Boley (2013) Shabram, M. & Boley, A. C. 2013, ApJ, 767, 63
  * Sheehan et al. (2019) Sheehan, P. D., Wu, Y.-L., Eisner, J. A., & Tobin, J. J. 2019, ApJ, 874, 136
  * Snik et al. (2010) Snik, F., de Wijn, A. G., Ichimoto, K., et al. 2010, A&A, 519, A18
  * Soummer et al. (2012) Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28
  * Sparks & Axon (1999) Sparks, W. B. & Axon, D. J. 1999, PASP, 111, 1298
  * Stamatellos & Whitworth (2009) Stamatellos, D. & Whitworth, A. P. 2009, MNRAS, 392, 413
  * Stolker et al. (2019) Stolker, T., Bonse, M. J., Quanz, S. P., et al. 2019, A&A, 621, A59
  * Stolker et al. (2017) Stolker, T., Min, M., Stam, D. M., et al. 2017, A&A, 607, A42
  * Stolker et al. (2020) Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182
  * Szulágyi et al. (2017) Szulágyi, J., Mayer, L., & Quinn, T. 2017, MNRAS, 464, 3158
  * Takata & Stevenson (1996) Takata, T. & Stevenson, D. J. 1996, Icarus, 123, 404
  * Tata et al. (2009) Tata, R., Martín, E. L., Sengupta, S., et al. 2009, A&A, 508, 1423
  * Tinyanont et al. (2019) Tinyanont, S., Millar-Blanchaer, M. A., Nilsson, R., et al. 2019, PASP, 131, 025001
  * Tobin et al. (2016) Tobin, J. J., Kratter, K. M., Persson, M. V., et al. 2016, Nature, 538, 483
  * Torres et al. (2006) Torres, C. A. O., Quast, G. R., da Silva, L., et al. 2006, A&A, 460, 695
  * van der Plas et al. (2016) van der Plas, G., Ménard, F., Ward-Duong, K., et al. 2016, ApJ, 819, 102
  * van Holstein et al. (2020) van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020, A&A, 633, A64
  * van Holstein et al. (2017) van Holstein, R. G., Snik, F., Girard, J. H., et al. 2017, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10400, Proc. SPIE, 1040015
  * Veras et al. (2009) Veras, D., Crepp, J. R., & Ford, E. B. 2009, ApJ, 696, 1600
  * Vigan et al. (2020) Vigan, A., Fontanive, C., Meyer, M., et al. 2020, arXiv e-prints, arXiv:2007.06573
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Wahhaj et al. (2011) Wahhaj, Z., Liu, M. C., Biller, B. A., et al. 2011, ApJ, 729, 139
  * Wang et al. (2018) Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192
  * Wertz et al. (2017) Wertz, O., Absil, O., Gómez González, C. A., et al. 2017, A&A, 598, A83
  * Whittet et al. (1992) Whittet, D. C. B., Martin, P. G., Hough, J. H., et al. 1992, ApJ, 386, 562
  * Wiktorowicz et al. (2014) Wiktorowicz, S. J., Millar-Blanchaer, M., Perrin, M. D., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Proc. SPIE, 914783
  * Williams & Cieza (2011) Williams, J. P. & Cieza, L. A. 2011, ARA&A, 49, 67
  * Woitke et al. (2016) Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103
  * Wolff et al. (2017) Wolff, S. G., Ménard, F., Caceres, C., et al. 2017, AJ, 154, 26
  * Wood et al. (2019) Wood, C. M., Boyajian, T., von Braun, K., et al. 2019, ApJ, 873, 83
  * Wu et al. (2020) Wu, Y.-L., Bowler, B. P., Sheehan, P. D., et al. 2020, AJ, 159, 229
  * Wu et al. (2017a) Wu, Y.-L., Close, L. M., Eisner, J. A., & Sheehan, P. D. 2017a, AJ, 154, 234
  * Wu et al. (2015) Wu, Y.-L., Close, L. M., Males, J. R., et al. 2015, ApJ, 807, L13
  * Wu & Sheehan (2017) Wu, Y.-L. & Sheehan, P. D. 2017, ApJ, 846, L26
  * Wu et al. (2017b) Wu, Y.-L., Sheehan, P. D., Males, J. R., et al. 2017b, ApJ, 836, 223
  * Xuan et al. (2020) Xuan, J. W., Bryan, M. L., Knutson, H. A., et al. 2020, AJ, 159, 97
  * Yan et al. (2019) Yan, Q.-Z., Zhang, B., Xu, Y., et al. 2019, A&A, 624, A6
  * Zapatero Osorio et al. (2011) Zapatero Osorio, M. R., Béjar, V. J. S., Goldman, B., et al. 2011, ApJ, 740, 4
  * Zapatero Osorio et al. (2005) Zapatero Osorio, M. R., Caballero, J. A., & Béjar, V. J. S. 2005, ApJ, 621, 445
  * Zapatero Osorio et al. (2006) Zapatero Osorio, M. R., Martín, E. L., Bouy, H., et al. 2006, ApJ, 647, 1405
  * Zhou et al. (2014) Zhou, Y., Herczeg, G. J., Kraus, A. L., Metchev, S., & Cruz, K. L. 2014, ApJ, 783, L17

##  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. 15 (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.

Figure 15: 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. 15 (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 respect to 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 use
large apertures to perform the photometry with, such that the bad pixels
values (which are both positive and negative 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 $\beta$ Pic b

Companions at 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 the method 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 polarimetric data-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
$\beta$ Pic. In this section we demonstrate the method with the 2019-11-26
$H$-band data set of $\beta$ Pic and show how we calculate upper limits on the
degree of polarization of the companion $\beta$ Pic b. A total-intensity image
of the data after applying classical ADI with IRDAP is shown in Fig. 16
(left).

Figure 16: Reduced images of the 2019-11-26 data set of $\beta$ Pic. Left:
Total-intensity image after applying classical ADI with IRDAP to reveal the
companion $\beta$ 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
PynPoint666https://pynpoint.readthedocs.io 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^{6}$ 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. 17 (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.

(a)  (b)

Figure 17: Final probability distributions of the signal of $\beta$ Pic b from
the 2019-11-26 data set of $\beta$ 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.

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 $\beta$
Pic, the classical ADI 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 $\beta$ 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^{\circ}$. Indeed, as can be seen in Fig. 16 (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 following the 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 $\beta$ Pic we detect
no signals with an S/N higher than 0.9 in $Q$ and 1.9 in $U$ for any aperture
radius. Indeed, the reduced $Q$-image (see Fig. 16, 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
$\beta$ 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\sigma$ and $3\sigma$
confidence intervals, respectively. These intervals are calculated one-sided
and starting at zero because the degree of linear polarization is computed as
$P=\surd(q^{2}+u^{2})$ and therefore can only have positive values (see Sparks
& Axon 1999). Figure 17 (right) shows the distribution of the degree of
polarization for an aperture radius of 1.85 pixels with the two intervals
indicated. Figure 18 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
$\beta$ 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).

Figure 18: Upper limits on the degree of linear polarization of $\beta$ 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 $\beta$ 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 19 shows the $I_{Q}$-image of
this data set.

Figure 19: Reduced $I_{Q}$-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.

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. 19) and is described by a constant
(the $z$-intercept) and slopes in the $x$\- 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 ($\mathit{SSR}$):

$\mathit{SSR}=\sum\limits_{i=1}^{N}\left[\left(I_{Q,i}-\hat{I}_{Q,i}\right)^{2}+\left(I_{U,i}-\hat{I}_{U,i}\right)^{2}\right],$
(5)

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, $\hat{I}_{Q,i}$ and $\hat{I}_{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 now repeat 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-\dfrac{1}{2}\Bigg{[}N\ln\left(2\pi\sigma^{2}\right)+\dfrac{\mathit{SSR}}{\sigma^{2}}\Bigg{]},$
(6)

where $\mathit{SSR}$ is computed from Eq. (5) and $\sigma$ 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 $\sigma$, we include it among the parameters to be
fitted (i.e., we treat $\sigma$ as a nuisance parameter). We set the prior for
$\sigma$ proportional to $1/\sigma$, 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. 20. Figure 21 shows the
resulting 1D- and 2D-projections of the posterior distribution of the fitted
parameters. The distributions in Fig. 21 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.

Figure 20: Data, best-fit model and residual images of the MCMC fitting of the
stellar PSF at position of HD 19467 B to the reduced $I_{Q}$\- and
$I_{U}$-images of the 2018-08-07 data set of HD 19467. Figure 21: Posterior
distributions after using MCMC to fit the stellar PSF at the position of the
companion HD 19467 B to the reduced $I_{Q}$\- and $I_{U}$-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 $I_{Q}$ and
$I_{U}$, the background’s $z$-intercept and slopes in the $x$\- and
$y$-direction in $I_{Q}$ and $I_{U}$, and the standard deviation $\sigma$ 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 and
84th 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\sigma$, $2\sigma$ and $3\sigma$ 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 Contrast curve of $\beta$ Pic data

Figure 22 shows the $1\sigma$ and $5\sigma$ point-source contrast in $Q$ and
$U$ as a function of angular separation from the star for the mean-combined
data set of $\beta$ 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 ${\sim}0.2\arcsec$ and $2.0\arcsec$ the
polarimetric sensitivity is close to the photon-noise limit, with a
$1\sigma$-contrast of $7\cdot 10^{-8}$ to $1\cdot 10^{-8}$ and a
$5\sigma$-contrast of $5\cdot 10^{-7}$ to $5\cdot 10^{-8}$. At separations
larger than $2.0\arcsec$ the sensitivity is limited by read noise or
background noise and the $1\sigma$\- and $5\sigma$-contrast are ${<}1\cdot
10^{-8}$ and ${<}5\cdot 10^{-8}$, respectively.

Figure 22: $1\sigma$ and $5\sigma$ point-source contrast in $Q$ and $U$ as a
function of angular separation from the star for the mean-combined data set of
$\beta$ Pic. The azimuthally averaged flux in the total-intensity $I_{Q}$\-
and $I_{U}$-images and the corresponding photon noise are shown for
comparison.