<figure><img src="image_1.png" /></figure>

# Detection of Polarization due to Cloud Bands in the Nearby Luhman 16 Brown Dwarf Binary

Maxwell A. Millar-Blanchaer<fnref n="1" />,<fnref n="2" />,<fnref n="12" /> ID, Julien H. Girard<fnref n="3" />,<fnref n="4" /> ID, Theodora Karalidi<fnref n="5" /> ID, Mark S. Marley<fnref n="6" /> ID, Rob G. van Holstein<fnref n="7" />, Sujan Sengupta<fnref n="8" /> ID, Dimitri Mawet<fnref n="2" /> ID, Tiffany Kataria<fnref n="1" /> ID, Frans Snik<fnref n="7" />, Jos de Boer<fnref n="7" /> ID, Rebecca Jensen-Clem<fnref n="9" /> ID, Arthur Vigan<fnref n="10" />, and Sasha Hinkley<fnref n="11" /> ID

## Abstract

Brown dwarfs exhibit patchy or spatially varying banded cloud structures that are inferred through photometric and spectroscopic variability modeling techniques. However, these methods are insensitive to rotationally invariant structures, such as the bands seen in Jupiter. Here, we present H-band Very Large Telescope/NaCo linear polarization measurements of the nearby Luhman 16 L/T transition binary, which suggest that Luhman 16A exhibits constant longitudinal cloud bands. The instrument was operated in pupil tracking mode, allowing us to unambiguously distinguish between a small astrophysical polarization and the ~2% instrumental linear polarization. We measure the degree and angle of linear polarization of Luhman 16A and B to be $p_A = 0.031 \pm 0.004$% and $\psi_A = -32^\circ \pm 4^\circ$, and $p_B = 0.010 \pm 0.004$% and $\psi_B = 73_{-11}^{+13}^\circ$, respectively. Using known physical parameters of the system, we demonstrate that an oblate homogeneous atmosphere cannot account for the polarization measured in Luhman 16A, but could be responsible for that of the B component. Through a nonexhaustive search of banded cloud morphologies, we demonstrate a two-banded scenario that can achieve a degree of linear polarization of $p = 0.03\%$ and conclude that the measured polarization of the A component must be predominantly due to cloud banding. For Luhman 16B, either oblateness or cloud banding could be the dominant source of the measured polarization. The misaligned polarization angles of the two binary components tentatively suggest spin-orbit misalignment. These measurements provide new evidence for the prevalence of cloud banding in brown dwarfs while at the same time demonstrating a new method—complementary to photometric and spectroscopic variability methods—for characterizing the cloud morphologies of substellar objects without signs of variability.

*Unified Astronomy Thesaurus concepts: Near infrared astronomy (1093); Very Large Telescope (1767); Polarimetry (1278)*

## 1. Introduction

Brown dwarfs occupy a unique parameter space, with effective temperatures ($T_{eff}$), masses, and radii in between those of giant exoplanets and stars. After their initial formation, they radiatively cool over time, moving from late-M through L, T, then Y spectral types, experiencing both chemical and atmospheric evolution. Brown dwarfs at the L/T spectral-type transition are believed to undergo an evolution from extremely dusty/cloudy atmospheres, where the clouds are mostly made of corundum, iron, and silicates, to nearly clear atmospheres that eventually begin to form clouds from other families of condensates such as Cr, MnS, Na₂S, ZnS, and KCl (Burgasser et al. 2002; Marley et al. 2010; Morley et al. 2012). This theory

is bolstered by photometric and spectroscopic variability studies that have revealed increased variability across the transition (e.g., Radigan et al. 2012; Crossfield et al. 2014; Biller 2017; Artigau 2018), suggestive of patchy clouds (Karalidi et al. 2016) or longitudinally varying cloud bands (Apai et al. 2017). Understanding cloud morphology in brown dwarfs is important as clouds affect their disk-integrated spectra and colors, and directly relate to the radiative, advective, and chemical processes taking place within their atmospheres (Showman & Kaspi 2013). Studies of brown dwarf clouds can also serve as probes of cloud formation and transport on directly imaged gas giants, which can have similar effective temperatures and surface gravities (e.g., Bowler 2016).

Polarimetry is a useful tool for studying clouds and hazes in brown dwarfs and is highly complementary to photometry and spectroscopy. As the emitted light of a brown dwarf is scattered by clouds and hazes in its atmosphere, it can locally acquire a preferred linear polarization as it gets redirected toward the observer (Sengupta & Krishan 2001; Sengupta & Marley 2009, 2010). This preferred polarization will cancel itself out in an

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹ Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91109, USA; max.a.millar-blanchaer@jpl.nasa.gov</footnote>

<footnote n="2">² Department of Astronomy, California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA</footnote>

<footnote n="3">³ Space Telescope Science Institute, Baltimore, MD 21218, USA</footnote>

<footnote n="4">⁴ Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="5">⁵ Department of Astronomy, UC Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="6">⁶ NASA Ames Research Center, Mountain View, CA 94035, USA</footnote>

<footnote n="7">⁷ Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="8">⁸ Indian Institute of Astrophysics, Koramangala 2nd Block, Bangalore 560 034, India</footnote>

<footnote n="9">⁹ Astronomy Department, University of California, Berkeley, Berkeley, CA 94720, USA</footnote>

<footnote n="10">¹⁰ Aix-Marseille Univ., CNRS, LAM, Laboratoire d'Astrophysique de Marseille, F-13013 Marseille, France</footnote>

<footnote n="11">¹¹ School of Physics, College of Engineering, Mathematics & Physical Sciences, University of Exeter, Stocker Road, Exeter, EX4 4QL, UK

Received 2019 May 21; revised 2020 January 8; accepted 2020 January 13; published 2020 May 5</footnote>

<footnote n="12">¹² NASA Hubble Fellow.</footnote>
<table><caption>**Table 1**
Selected Literature Properties for Luhman 16</caption><thead><tr><th rowspan="2">Property</th><th colspan="2">Component</th><th rowspan="2">References</th></tr><tr><th>A</th><th>B</th></tr></thead><tbody><tr><td>Spectral Type</td><td>L7.5 ± 1</td><td>T0.5 ± 1</td><td>Burgasser et al. (2013), Kniazev et al. (2013)</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>1310 ± 30</td><td>1280 ± 74</td><td>Faherty et al. (2014)</td></tr><tr><td>Mass (M<sub>☉</sub>)</td><td>33.5 ± 0.3</td><td>28.6 ± 0.3</td><td>Lazorenko & Sahlmann (2018)</td></tr><tr><td>Rotation Period (hr)</td><td>4.5–5.5<br>or 8</td><td>4.87 ± 0.01<br>or 5.05 ± 0.10</td><td>Buenzli et al. (2015a), Gillon et al. (2013)<br>Mancini et al. (2015), Burgasser et al. (2014)</td></tr><tr><td>Inclination (°)</td><td>56 ± 5 (for 5 hr period)<br>18 ± 8 (for 8 hr period)</td><td>26 ± 8</td><td>Karalidi et al. (2016), Crossfield (2014)<br>Karalidi et al. (2016)</td></tr><tr><td>Distance (pc)</td><td colspan="2">1.9937 ± 0.0003</td><td>Lazorenko & Sahlmann (2018)</td></tr></tbody></table>

unresolved measurement (which is always the case for brown dwarfs) and result in a net zero polarization unless there is some type of asymmetry in the atmosphere, such as rotationally induced oblateness, or patchy/banded clouds (de Kok et al. 2011). In addition to its oblateness, the measured degree of linear polarization for a given brown dwarf depends on the cloud particles in two ways: their scattering properties (determined by their size, shape, and composition) and their spatial distribution (Stolker et al. 2017). In this respect, polarization can act as a very effective diagnostic of cloud properties in brown dwarfs. These same effects are also present for giant extrasolar planets, where lower surface gravities can result in higher levels of oblateness for the same rotational periods (Marley & Sengupta 2011).

Polarization has been measured in over two dozen brown dwarfs (Ménard et al. 2002; Zapatero Osorio et al. 2005, 2011; Goldman et al. 2009; Tata et al. 2009; Miles-Páez et al. 2013). These measurements have revealed linear degrees of polarization between ~0.1% and 2.5%, spanning *R* to *J* bands. While many of these detections could possibly be explained by oblateness (e.g., Sengupta & Marley 2010), polarimetric monitoring of a handful of sources has revealed both short-term and long-term variability, suggesting that the time-varying morphology of the clouds also plays a significant role (Miles-Páez et al. 2015, 2017). Although in some circumstances polarization can be attributed to the presence of a circum-brown dwarf disk (Hashimoto et al. 2009; Zapatero Osorio et al. 2011; Miles-Páez et al. 2013; Ginski et al. 2018), for the vast majority of polarimetric measurements, the true origin of the net polarization remains unknown, due in part to the degeneracies in the atmospheric model parameters and the limited amount of physical characterization available from other measurements.

Here, we present near-infrared (NIR) linear polarimetric observations of the Luhman 16 brown dwarf binary system, obtained with the NaCo imager (Lenzen et al. 2003; Rousset et al. 2003) at the ESO Very Large Telescope (VLT). We measure a linear degree of polarization of $p_A = 0.031\% \pm 0.004\%$ and $p_B = 0.010\% \pm 0.04\%$ in Luhman 16A and B, respectively, corresponding to detection significances of $8\sigma$ and $2.5\sigma$. The existing extensive characterization of the system allows us to explore possible sources of polarization in greater detail than for all previous polarimetric measurements of brown dwarfs to date.

### 1.1. Luhman 16

The closest brown dwarf system to Earth is Luhman 16 (WISE J104915.57–531906.1AB; Luhman 2013), a brown

dwarf binary at a distance of ~1.99 pc (Lazorenko & Sahlmann 2018). The binary is of particular interest because the two components span the L/T transition, with spectral types of L7.5 ± 1 and T0.5 ± 1 for components A and B, respectively (Burgasser et al. 2013; Kniazev et al. 2013). The system's relative brightness compared to more distant brown dwarfs has resulted in many detailed studies that have been able to constrain the mass, rotation period, and inclination of both components (Table 1). Previous linear polarization measurements have put an upper limit of 0.07% in the *I* band on the unresolved binary (Kniazev et al. 2013).

Numerous photometric and spectroscopic variability studies have revealed that both components are variable, with variability amplitudes that change significantly from epoch to epoch (Biller et al. 2013; Gillon et al. 2013; Burgasser et al. 2014; Buenzli et al. 2015a, 2015b; Karalidi et al. 2016). In all cases, Luhman 16B is found to be the more variable of the two components, with peak to peak amplitudes up to 11% (*I* + *z* band; Gillon et al. 2013) and Luhman 16A having a maximum measured variability of only ~4.5% (between 0.8 and 1.15 µm; Buenzli et al. 2015a). Further, high-resolution spectroscopic monitoring of Luhman 16B has produced the first two-dimensional Doppler-imaging cloud map of a brown dwarf, revealing patchy variations in the cloud cover (Crossfield et al. 2014).

In general, the variability of both components of Luhman 16 has been attributed to patchy clouds, but recently longitudinally varying cloud bands with planetary-scale brightness variations have been used to explain the photometric variability of three other L/T transition objects, suggesting that this phenomenon may provide a more complete explanation (Apai et al. 2017). In the Luhman 16 system, preliminary evidence already hints at the possibility of cloud bands. Both the near-exact repetition of a light curve feature seen in two *Hubble Space Telescope* (HST) data sets separated by over a year (Karalidi et al. 2016), as well as the change in state of Luhman 16A from low variability to high variability (Buenzli et al. 2015a) could both be explained with variable cloud bands with slightly different periods, beating over time.

## 2. Observations and Data Reduction

Observations of Luhman 16 were obtained with VLT/NaCo in visitor mode, starting at 2018 April 11 23:22 UT and lasting until 2018 April 12 06:26 UT, covering a total of 7 hr and 4 minutes (Table 2). The observations span ~1.4 or ~0.9 rotation periods for Luhman 16A, assuming a rotation period of 5 or 8 hr, respectively, and ~1.4 rotations for Luhman 16B. The
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top row: the parallactic angle, (relative) Y-detector position of the Luhman 16A component in the extraordinary beam, air mass, and seeing as a function of time from the start of our observations. Bottom row: left—the parallactic angle vs. Y-detector position of Luhman 16A in the extraordinary beam. Center left, center right, and right—the flux ratio of A to B as a function of y-detector position, air mass, and seeing. Although slight correlations are measured between the flux ratio and these three parameters (Section 2.2), it is difficult to explain the variations in the ratio seen in Figure 3 as being significantly influenced by the y-detector position, air mass, or seeing.</figcaption></figure>

<table><caption>**Table 2**
Summary of Observations</caption><thead><tr><th>Object</th><th>UT Date</th><th>T<sub>int(s)</sub></th><th>N<sub>coadds</sub></th><th>N<sub>obs</sub></th><th>t<sub>total(s)</sub> =<br>T<sub>int</sub>N<sub>coadds</sub>N<sub>obs</sub></th><th>θ<sub>HWP</sub> Sequence (°)</th></tr></thead><tbody><tr><td>Luhman 16</td><td>2018 Apr 11</td><td>3.0</td><td>20</td><td>276</td><td>16,560.0</td><td>0, 45, 67.5, 135</td></tr><tr><td>Twilight</td><td>2018 Apr 11</td><td>8.0</td><td>4</td><td>20</td><td>640.0</td><td>0, 45, 67.5, 135</td></tr><tr><td rowspan="2">Elias 2-25 (Polarized Standard)</td><td>2018 Jun 01</td><td>0.345</td><td>60</td><td>8</td><td>165.6</td><td>0, 45, 67.5, 135, 22.5</td></tr><tr><td>2018 Jun 13</td><td>0.345</td><td>60</td><td>8</td><td>165.6</td><td>0, 45, 67.5, 135, 22.5</td></tr><tr><td rowspan="2">HD 162973 (Unpolarized Standard)</td><td>2018 May 28</td><td>30.0</td><td>1</td><td>10</td><td>300.0</td><td>0, 45, 67.5, 135, 22.5</td></tr><tr><td>2018 May 30</td><td>30.0</td><td>1</td><td>10</td><td>300.0</td><td>0, 45, 67.5, 135, 22.5</td></tr></tbody></table>

instrument was operated in NaCo's dual-channel polarimetry mode, using the SL27 camera, which provides a pixel scale of 27 mas/pixel and a 27 × 27" field of view. At the time of the observations, NaCo was mounted on UT1 on the Nasmyth A platform. NaCo's polarimetry mode includes a rotatable half-wave plate (HWP), a focal plane mask, and a Wollaston prism that splits the incident beam along the detector y-axis into ordinary ($I_\parallel$) and extraordinary ($I_\perp$) beams with a splitting angle that shifts the beams on the detector by the on-sky equivalent of 3"/5. We have labeled the two beams "ordinary" and "extraordinary" to be consistent with previous NaCo documentation; however, the modulation of the HWP means each beam effectively acts in both capacities. To prevent beam overlap, the focal plane mask blocks strips of 27 × 3"/5, alternating between blocked regions and transmitted regions with a width of 3"/5 along the detector y-axis.

Observations were obtained in four-image groups with the HWP cycling between position angles of 0°, 45°, 67°/5, and 135°. Note that this was unintentionally different from the standard 0°, 45°, 22°/5, and 67°/5, and the result of a user error. This sequence of HWP angles is not recommended. At each HWP position, we obtained one image with an exposure time of 3 s (DIT = 3 s) and 20 coadds (NDIT = 20) in NaCo's cube mode. After each HWP cycle, the telescope was dithered along the detector's x-direction by ±5", such that the observations were taken in an ABAB dither pattern.

The NAOS adaptive-optics (AO) system was operated using the K-band dichroic, which transmits 1.9–2.5 µm light to the wavefront sensor with 90% efficiency, while sending

0.45–1.8 µm light to the CONICA imaging system, also with 90% efficiency.

To minimize time-varying instrumental polarization effects, the instrument was operated in pupil tracking mode in a "crossed configuration." In this configuration, the entrance fold mirror of NaCo is oriented such that any instrumental polarization it introduces will be of opposite sign to that of M3, thus minimizing the cumulative instrumental polarization reaching the HWP (e.g., Witzel et al. 2011; de Juan Ovelar 2013). The combined use of the pupil tracking mode with NaCo's polarimetry mode was first employed by de Juan Ovelar (2013), who attempted to measure the polarization of the directly imaged planets of HR 8799. In NaCo's pupil tracking mode, the entire instrument rotates as the telescope altitude changes, such that the orientation of M3 relative to the instrument remains constant throughout the observations. Thus, under this configuration, the instrumental polarization is both stable and minimized. Over the course of our observations, the A and B components rotated with parallactic angle relative to their center of light due to the pupil tracking configuration. Our one-minute exposures and four-exposure modulation cycles are significantly faster than the parallactic rotation.

Figure 1 provides a summary of the parallactic angle, seeing, and air mass as a function of time. The seeing across the observations ranged from 0"/24 to 1"/09. During the observations, the AO loops opened several times, causing two minor interruptions (i.e., less than 10 minutes) and one ~45 minute interruption. After just over seven hours of observing, a fatal mechanical malfunction in the HWP unit forced us to stop
observing. We obtained a total of 276 individual exposures, amounting to a total integration time of 4 hr and 36 minutes.

Observations of the evening twilight sky were obtained on 2018 April 11 UT, with the telescope pointed with an altitude of 90° and an azimuth angle of 45° east of north. Although predicting the exact degree of linear polarization is difficult, the angle of linear polarization is expected to be oriented 90° away from the vector connecting the telescope pointing location and the Sun (e.g., Harrington et al. 2011; deBoer et al. 2014). We obtained five sets of four HWP positions. In order to mimic the Luhman 16 observing configuration, the K-band dichroic was inserted into the beam and the instrument was oriented in the same “crossed configuration” as above (verified by checking that the ESO ADA PUPILPOS header keyword was equal to the expected position of 90° for the crossed configuration). See Table 2 for a summary of the observations.

Due to the failure of the HWP rotation mechanism on 2018 April 11, observations of a polarized and unpolarized standard were not possible until after an intervention to repair the mechanism, which occurred in late 2018 May. The failure occurred because an axle that drives the rotation of the HWP snapped and part of the mechanism had to be replaced. The replacement process lost the known calibration between the motor encoder position and angle of the HWP. As a result, the observations after the intervention display a systematic offset in the measured polarization angle from the Luhman 16 data. This offset was fit for as part of our data analysis (see Section 3).

After the repair, observations of the polarized standard Elias 2-25 ($p = 4.13\% \pm 0.02\%$, $\psi_{PA} = 24^\circ \pm 1^\circ$ in the H band; Whittet et al. 1992) and the unpolarized standard HD 162973 ($p = 0.09\% \pm 0.055\%$, $\psi_{PA} = 104^\circ \pm 17^\circ$ in the B-band; Mathewson & Ford 1970) were obtained in queue mode, with each target observed both before and after meridian passage for one HWP cycle in each of the same two dither positions obtained for Luhman 16. For both of these targets, one HWP cycle consisted of HWP positions of 0°, 45°, 67°, 5, 135°, and 22°.5. The observations were obtained in exactly the same observing mode as Luhman 16: pupil tracking, H-band filter, field mask inserted, and using the K-dichroic. See Table 2 for a summary of the observations.

## 2.1. Luhman 16 Data Reduction

A master dark and a master flat field (using lamp flats) were created using the NaCo img dark and NaCo img lamp flat recipes, respectively, from the ESO Gasgano pipeline with the default calibration files provided by ESO. Each individual exposure was then dark subtracted and divided by the master flat field. Although the data was taken in NaCo's cube mode, we opted to carry out our analysis on the “single” frame images (i.e., the average of all the NDITS in each image). Background subtraction was carried out by subtracting from each frame an image at the opposite dither position but with the same HWP angle, from either the following (for the “A” dither position) or the preceding (for the “B” dither position) HWP cycle.

In each exposure, both Luhman 16A and B components were detected with a signal-to-noise ratio (S/N) between ~1000 and ~1500, in both the ordinary and extraordinary beams (Figure 2). The two objects were clearly resolved in all images. However, the two point spread functions (PSFs) can be seen to overlap with each other slightly.

For each image, the photutils (Bradley et al. 2017) Python package was used to measure the location of each of the

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Images of Luhman 16A and B from a single 60 s observation with NaCo. The pixel scale is $27.053 \pm 0.019$ mas/pixel. The data have been dark subtracted, flat-fielded, and background subtracted. NaCo's Wollaston prism creates two images of the binary, the extraordinary (left) and ordinary (right) beam. Luhman 16A is the northern point source and B is to the south. On the ordinary beam, apertures sizes of 6, 12, and 17 pixels have been overlaid. The original data have been rotated to orient the field north up and east left.</figcaption></figure>

two binary components in both the ordinary and extraordinary beams and extract their photometry. Each of the four sources was first found using the DAOStarFinder routine and then photometry was extracted using a circular aperture. Several different radii for the aperture were tested (6, 10, 12, 13, 16, and 17 pixels), and ultimately, we found a radius of 17 pixels ($0.46''$) to provide a balance between the formal S/N measured by photutils and the true scatter of the data points. This radius was the largest possible radius without the circular apertures from A and B overlapping. Uncertainties on each photometric measurement were provided by photutils, assuming a read noise of 4.2 ADU and a gain of 11 electrons/ADU.

## 2.2. Total Intensity Data Reduction

To obtain the total intensity for each component at each time step, the aperture sums for the ordinary and extraordinary beams were added together (see Table A1). The absolute variability of the individual components appears to be highly correlated and is likely due to changing atmospheric conditions and imperfect AO correction. Around the 4 hr mark, the measured flux dropped by 30%–40%, due in part to poor AO correction and an enlarged PSF. The apertures sizes could not be increased further to compensate because of the separation between the two components and the potential for overlap.

It was not possible to derive accurate absolute photometry due to changing atmospheric properties and adaptive-optics systematics. The flux ratio between the two components displays quasi-periodic variability on timescales of less than 2 hr, with the values of the peaks, troughs, and the peak-to-valley distance changing over the length of the observations (Figure 3). In order to explore the fidelity of this signal, we measured the Spearman-r and Kendall-Tau correlation coefficients (and the corresponding null hypothesis *p*-values) of the flux ratio against the detector *x* and *y* positions, the FWHM of each component, air mass, and seeing (Table 3). We measure no significant correlations, but some low-level correlations with *y*-position, air mass, and seeing. Although these small correlations may be real, when the flux ratio is plotted against these quantities, it is clear that the main quasi-periodic photometric signal cannot be explained by these correlations (see Figure 1).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top: the relative photometry of Luhman 16A and B as a function of time. Quasi-periodic variation can be seen throughout the sequence. The noisier region near the 4 hr mark is associated with poor adaptive-optics correction. Bottom: a Lomb–Scargle periodogram of the relative photometry. The data obtained between 4 and 4.75 hr after the start of the sequence have been excluded. A prominent peak can be seen at 1.64 hr.</figcaption></figure>

<table><caption>**Table 3**
Correlation Coefficients Calculated between Our *q* Measurements and Other Data Parameters</caption><thead><tr><th rowspan="2">Data Set</th><th rowspan="2">Comparison</th><th colspan="2">Spearman-r</th><th colspan="2">Kendall-Tau</th></tr><tr><th>Value</th><th>p-value</th><th>Value</th><th>p-value</th></tr></thead><tbody><tr><td>F<sub>A</sub>/F<sub>B</sub></td><td>x-position</td><td>-0.01</td><td>0.84</td><td>-0.01</td><td>0.75</td></tr><tr><td>F<sub>A</sub>/F<sub>B</sub></td><td>y-position</td><td>-0.18</td><td>0.004</td><td>-0.13</td><td>0.002</td></tr><tr><td>F<sub>A</sub>/F<sub>B</sub></td><td>FWHM<sub>A</sub></td><td>-0.02</td><td>0.70</td><td>-0.02</td><td>0.67</td></tr><tr><td>F<sub>A</sub>/F<sub>B</sub></td><td>FWHM<sub>B</sub></td><td>-0.03</td><td>0.68</td><td>-0.02</td><td>0.72</td></tr><tr><td>F<sub>A</sub>/F<sub>B</sub></td><td>Air mass</td><td>-0.19</td><td>0.002</td><td>-0.11</td><td>0.008</td></tr><tr><td>F<sub>A</sub>/F<sub>B</sub></td><td>Seeing</td><td>0.15</td><td>0.02</td><td>0.09</td><td>0.024</td></tr><tr><td>q<sub>A</sub></td><td>x-position</td><td>0.73</td><td>5.7 &times; 10<sup>-7</sup></td><td>0.52</td><td>1.0 &times; 10<sup>-5</sup></td></tr><tr><td>q<sub>A</sub></td><td>y-position</td><td>0.67</td><td>9.2 &times; 10<sup>-6</sup></td><td>0.45</td><td>1.5 &times; 10<sup>-4</sup></td></tr><tr><td>q<sub>B</sub></td><td>x-position</td><td>-0.26</td><td>0.12</td><td>-0.18</td><td>0.12</td></tr><tr><td>q<sub>B</sub></td><td>y-position</td><td>-0.08</td><td>0.62</td><td>-0.06</td><td>0.60</td></tr><tr><td>q<sub>A</sub></td><td>FWHM<sub>A</sub></td><td>0.03</td><td>0.86</td><td>0.02</td><td>0.85</td></tr><tr><td>q<sub>B</sub></td><td>FWHM<sub>B</sub></td><td>0.001</td><td>0.99</td><td>0.008</td><td>0.94</td></tr></tbody></table>

Figure 3 displays a Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) of the flux ratio data, excluding the measurements between 4 and 4.85 hr where the AO correction was poor. A strong peak can be seen at 1.64 hr, roughly a factor of 3 less than the previously measured 4.87 ± 0.01 hr rotation period of Luhman 16B (Gillon et al. 2013). In Figure 3, we also explore whether our selection of aperture size affected the peak in the periodogram. For all apertures sizes, we find the same significant peak at 1.64 hr.

Using the relative positions of each source measured by the DAOStarFinder, we measured a separation of 934 ± 2 mas (34.52 ± 0.07 pixels) by averaging over all measurements in our observing sequence, including both dither positions and both the ordinary and extraordinary beams. The relative positions in pixels were converted to an on-sky separation using a pixel scale for the S27 camera of 27.053 ± 0.019 mas/pixel.<fnref n="13" /> The quoted errors represent the standard deviation of

the measurements added in quadrature to the error due to the plate scale uncertainty. For each image, we also measured a relative angle between the two sources in detector coordinates, which was then converted to an angle on sky by adding the parallactic angle from the header for each observation (taken as the average of the ESO TEL PARANG START and ESO TEL PARANG END header keywords). By averaging all measurements, we obtained a relative position angle between A and B of 147.0° ± 0.1°, measured east of north. Neither north angle correction nor distortion correction has been applied.

## 2.3. Polarimetry Data Reduction

In order to extract the measured polarization from the photometric measurements of each component, we built a Mueller matrix model of the system to describe how on-sky polarization relates to the intensities measured on the detector. The on-sky Stokes vector, $S_{\text{sky}} \equiv [I_{\text{sky}}, Q_{\text{sky}}, U_{\text{sky}}, V_{\text{sky}}]$, is related to the measured Stokes vector at the detector, $S_{\text{det}}$, through a system Mueller matrix, $M_{\text{sys}}$:

$$S_{\text{det}} = M_{\text{sys}}(\theta_{\text{HWP}}, \theta_{\text{PA}})S_{\text{sky}}, \quad (1)$$

where

$$M_{\text{sys}} = M_{\text{Woll}}^{o,e} M_{\text{HWP}}(\theta_{\text{HWP}}) M_{\text{inst}} M_{\text{tel}} M_{\text{rot}}(\theta_{\text{PA}}) S_{\text{sky}}, \quad (2)$$

and $M_{\text{Woll}}$, $M_{\text{HWP}}$, $M_{\text{inst}}$, and $M_{\text{tel}}$ represent the Mueller matrix for the Wollaston prism, the half-wave plate, the optics of NaCo in front of the wave plate, and the telescope, respectively. A rotation matrix, $M_{\text{rot}}$, that depends on the parallactic angle, $\theta_{\text{PA}}$, compensates for the parallactic rotation of the sky relative to the instrument frame. The *o* and *e* superscripts for $M_{\text{Woll}}$ represent the ordinary and extraordinary beams. The Mueller matrix for the HWP depends on its position angle, $\theta_{\text{HWP}}$, obtained from the ESO INS RETA2 ROT header keyword. The Mueller matrices for the Wollaston, HWP, and rotation follow the conventions presented in Goldstein (2003). We follow the IAU definitions of *Q*, *U*, and *V* for our on-sky frame of reference.

We define the Stokes vector seen by the waveplate, $S_{\text{NaCo}}$, as

$$S_{\text{NaCo}} = M_{\text{inst}} M_{\text{tel}} M_{\text{rot}}(\theta_{\text{PA}}) S_{\text{sky}}, \quad (3)$$

such that

$$S_{\text{det}} = M_{\text{Woll}}^{o,e} M_{\text{HWP}}(\theta_{\text{HWP}}) S_{\text{NaCo}}, \quad (4)$$

Most detectors, including NaCo's, are only sensitive to the Stokes *I* term. By evaluating the Mueller matrices for the Wollaston prism and the HWP, it can be shown that $q_{\text{NaCo}} \equiv Q_{\text{NaCo}}/I_{\text{NaCo}}$ can be measured from the Stokes $I_{\text{det}}$ values of the ordinary and extraordinary beams by taking the normalized difference of the two beams (i.e., a “single difference”):

$$q_{\text{NaCo}} = \frac{I_{\text{det}}^{o}(0^{\circ}) - I_{\text{det}}^{e}(0^{\circ})}{I_{\text{det}}^{o}(0^{\circ}) + I_{\text{det}}^{e}(0^{\circ})}, \quad (5)$$

where the angle in parentheses refers to the angle of the HWP. However, noncommon path errors and different detector systematics between the ordinary and extraordinary beams can introduce an error/bias term, $\epsilon$, rewriting Equation (5) as

$$q_{\text{NaCo}} + \epsilon = \frac{I_{\text{det}}^{o}(0^{\circ}) - I_{\text{det}}^{e}(0^{\circ})}{I_{\text{det}}^{o}(0^{\circ}) + I_{\text{det}}^{e}(0^{\circ})}. \quad (6)$$

<footnote n="13">¹³ VLT NaCo User Manual Issue 102.</footnote>
Rotating the HWP by 45° swaps the sign of $q_{NaCo}$ in Equation (6):

$$q_{NaCo} - \epsilon = - \frac{I_{det}^{o}(45^{\circ}) - I_{det}^{e}(45^{\circ})}{I_{det}^{o}(45^{\circ}) + I_{det}^{e}(45^{\circ})}, \quad (7)$$

enabling the measurement of an unbiased $q_{NaCo}$ via a “double difference”:

$$q_{NaCo} = 0.5 * \left( \frac{I_{det}^{o}(0^{\circ}) - I_{det}^{e}(0^{\circ})}{I_{det}^{o}(0^{\circ}) + I_{det}^{e}(0^{\circ})} - \frac{I_{det}^{o}(45^{\circ}) - I_{det}^{e}(45^{\circ})}{I_{det}^{o}(45^{\circ}) + I_{det}^{e}(45^{\circ})} \right). \quad (8)$$

The error term $\epsilon$ can be calculated as

$$\epsilon = 0.5 * \left( \frac{I_{det}^{o}(0^{\circ}) - I_{det}^{e}(0^{\circ})}{I_{det}^{o}(0^{\circ}) + I_{det}^{e}(0^{\circ})} + \frac{I_{det}^{o}(45^{\circ}) - I_{det}^{e}(45^{\circ})}{I_{det}^{o}(45^{\circ}) + I_{det}^{e}(45^{\circ})} \right). \quad (9)$$

Likewise, a single differenced $u_{NaCo} \equiv U_{NaCo}/I_{NaCo}$ value can be obtained with the HWP at an angle of 22°:5:

$$u_{NaCo} + \epsilon = \frac{I_{det}^{o}(22^{\circ}5) - I_{det}^{e}(22^{\circ}5)}{I_{det}^{o}(22^{\circ}5) + I_{det}^{e}(22^{\circ}5)}, \quad (10)$$

or with the HWP at an angle of 67°:5:

$$u_{NaCo} - \epsilon = - \frac{I_{det}^{o}(67^{\circ}5) - I_{det}^{e}(67^{\circ}5)}{I_{det}^{o}(67^{\circ}5) + I_{det}^{e}(67^{\circ}5)}, \quad (11)$$

and the double-differenced $u_{NaCo}$ can be calculated in a similar fashion to Equation (8).

In general, the bias term $\epsilon$ represents the noncommon path errors between the ordinary and extraordinary beams that do not depend on the incident polarization state (and hence the HWP position). This error changes with time due to time-varying instrument and detector drifts, motivating the standard $\theta_{HWP} = [0^{\circ}, 45^{\circ}, 22^{\circ}5, 67^{\circ}5]$ HWP sequence, where $\epsilon$ is calculated for every pair of $q$ and $u$ measurements. In our case, for the twilight and Luhman 16 data, without measurements with $\theta_{HWP} = 22.5^{\circ}$, each $u_{NaCo}$ measurement (obtained with $\theta_{HWP} = 67.5^{\circ}$) is corrected for the $\epsilon$ term using the bias measured from the preceding measurements with $\theta_{HWP} = 0^{\circ}$ and $45^{\circ}$ (Equation (9); i.e., not standard double-differencing). The time-varying nature of this systematic means that this nonstandard correction is not perfect and is expected to result in additional systematic noise in $u_{NaCo}$ relative to $q_{NaCo}$. We attempt to compensate for this extra added noise in our modeling (the $\sigma_{SD}$ term in Section 3).

A 90° degeneracy in HWP angles means that the measurements at 45° and 135° should be nearly identical, modulo changes in $\epsilon$ over time. Before calculating $q_{NaCo}$ or the $\epsilon$ term for a given HWP cycle, we average together the two measurements at 45° and 135° to increase the S/N.

For the Luhman 16A and B data, using the photometry measured for the $o$ and $e$ beams, we calculated $q_{NaCo}$ and $u_{NaCo}$ for each HWP cycle from Equations (8) and (10) (hereafter $q_A$ and $u_A$ for Luhman 16A and $q_B$ and $u_B$ for Luhman 16B). A visual inspection of the resultant values revealed a systematic offset in $q_{NaCo}$ and $u_{NaCo}$ between the two dither positions in both components, possibly suggesting a spatial dependence in the instrumental polarization along the $x$ spatial direction (on the order of ~0.02%). To compensate for this spatial dependence, we average the $q_{NaCo}$ and $u_{NaCo}$ values for each AB dither pair. Figure 4 (and Figures A1 and A2) shows the

$q_{NaCo}$ and $u_{NaCo}$ for Luhman 16A and B after averaging the two dither positions (also see Table A2). Our correlation analysis (Section 2.4.2) suggests that after this averaging, there is no longer any significant spatial dependence in the signal. Errors on $q_{NaCo}$ and $u_{NaCo}$ were propagated from the original photometric errors through Equations (6)–(11).

A detection of polarization in Luhman 16A can be seen as a near-sinusoid in the $q_{NaCo}$ frame. In the instrument frame, the instrumental polarization holds a constant value over time (because of the cross-configuration), and any astrophysical polarization modulates between $q_{NaCo}$ and $u_{NaCo}$ as the sky (and therefore the angle of polarization) rotates relative to the instrument according to the parallactic angle. For Luhman 16A, this can most easily be seen by comparing the first half of the data to the second half. Although it cannot as easily be identified by immediate inspection, polarization is also detected for Luhman 16B, but at a lower significance. In Figure A3, we show the same data but combined to calculate the linear degree of polarization ($p_{NaCo} = \sqrt{q_{NaCo}^2 + u_{NaCo}^2}$) and the angle of linear polarization ($\theta_{NaCo} = 0.5 \arctan 2(u_{NaCo}/q_{NaCo})$) in the instrument frame as a function of parallactic angle. Before calculating $p_{NaCo}$ and $\theta_{NaCo}$, we subtracted the mean $q_{NaCo}$ and $u_{NaCo}$ values from each measurement. This step acts as a preliminary instrumental polarization subtraction; because any real signal presents itself as a modulation on top of the instrumental polarization and our observations are relatively well centered on meridian passage, there should be nearly equal signal above and below the instrumental polarization values of both $q_{NaCo}$ and $u_{NaCo}$. The $p_{NaCo}$ measurements suffer from the squaring of the noise of $q_{NaCo}$ and $u_{NaCo}$, and can be difficult to interpret given the S/N of each measurement. However, for Luhman 16A, the angle of linear polarization (which does not suffer from the same increase in noise) shows a clear linear trend with the parallactic angle, indicating the presence of an astrophysical signal that is rotating with the sky relative to the fixed instrument frame. The same signal is not seen for Luhman 16B. Translating the instrument-frame measurements into calibrated sky-frame polarization measurements (i.e., $S_{sky}$ from Equation (1)) requires further instrument modeling, as described below (Section 3).

## 2.4. Real Detection or Residual Instrumental Polarization?

The signal presented in Figures 4 and 15 for Luhman 16A represents a significant jump in accuracy compared to previous polarimetry measurements obtained with NaCo, and it is therefore prudent to question whether this signal could be due to other systematic effects. Here we consider previous efforts to model the instrumental polarization of NaCo to put our data set in context, and we also discuss potential systematic effects that could affect the interpretation of our data.

### 2.4.1. Previous Work

Two papers have previously developed detailed Mueller matrix models for NaCo’s Wollaston/HWP polarimetry mode, both using a field-stabilized observing mode: Witzel et al. (2011) developed a Ks-band model and deBoer et al. (2014) developed an H-band model. Although de Juan Ovelar (2013) also operate in pupil tracking mode, they did not present uncertainties on their measurements of instrumental polarization, and they assume a simplified Mueller matrix model that assumes negligible $Q \leftrightarrow U$ crosstalks. As a result, we have
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Measured *q* (left) and *u* (right) in NaCo's instrument frame as a function of time for Luhman 16A (top) and B (bottom). The degrees of polarization of A and B are related to the peak-to-trough distance, and the angle of linear polarization is related to the parallactic angle at which the peaks and troughs occur. The displayed error bars include both the original read-noise and photometric errors (black), as well as the extra systematic errors (red; see Methods) included in our MCMC fitting (Section 3.2). Increased scatter and larger systematic error bars are seen for the *u* measurements, due to the HWP angles that were inadvertently used in the observation set. Overlaid on top of the data are accepted MCMC fits to the instrumental polarization (orange; $IP_Q$ and $IP_U$), as well as the best-fit models for A and B (blue). We have also included the results of our simple six-parameter fit (see Section 3) as a purple line.</figcaption></figure>

opted to focus our discussion on the work of Witzel et al. (2011) and deBoer et al. (2014).

Witzel et al. (2011) developed a Fresnel reflection-based model of the system upstream of the HWP composed of Mueller matrices for the tertiary mirror (M3), NAOS, and associated rotations between their relative frames of reference. The Mueller matrices of M3 and NAOS were constructed as linear diattenuating retarders, where the linear diattenuation and retardation were calculated using assumed material properties for aluminum (for M3) and Silflex (for two NAOS 45° fold mirrors), and the known angles of incidence upon the different mirrors in the system. This system model was then compared to measurements of four polarized sources in the IRS 16 cluster of the Galactic center. In order to match the observations, the authors tweak the imaginary index of refraction for aluminum and the retardance of the Silflex while simultaneously fitting for linear degrees and angles of polarization for the IRS 16 sources until a minimum $\chi^2$ value is obtained. They quote a final error on their polarization estimates of the IRS 16 sources to be 0.8%. Unfortunately, the details of the fitting procedure and the derivation of this final error are not included in the paper and thus make the error difficult to compare to our data. They also

compare their tweaked model to observations of three unpolarized standards (obtained under a different instrument configuration, with the HWP out of beam) and one polarized standard, and find that their model is accurate to ~0.3% ± 0.2%, where the uncertainty represents the median deviation of their residuals.

deBoer et al. (2014) carry out a strategy of using twilight polarization measurements, and different telescope and instrument orientations in order to constrain the Mueller matrices of M3 and NaCo in the *H* band. By rotating the telescope to different azimuthal positions and NaCo to different derotator angles, they take advantage of the known twilight sky polarization angle to partially recover the Mueller matrices of both M3 and NaCo. Using this strategy, to fully describe the Mueller matrices of both components assumptions on their analytical forms must be made (e.g., Fresnel reflections with imposed material properties). Nonetheless, the authors were able to obtain errors on the instrumental polarization of 0.4%. Although this method was able to provide a good first step at developing an accurate instrument model, the authors note several discrepancies in their fitting that suggest their model is not fully self-consistent and suggest future avenues for a more
in-depth characterization. In particular, their estimate of the
instrumental polarization appears to be systematically different
from that of de Juan Ovelar (2013), obtained with observations
of an unpolarized standard, and may therefore be called into
question.

The apparent signal in Figure 4 is of order ~0.03%, an order
of magnitude less than the model accuracy presented in Witzel
et al. (2011) and deBoer et al. (2014). In comparison to Witzel
et al. (2011), the Luhman 16 data set presented herein is of
significantly higher quality than their IRS 16 data (with
individual error bars of 0.3% in both Q and U) and their
unpolarized standard star data (with a median deviation relative
to the instrumental polarization model of 0.2%); our formal
photometric errors on each Q and U measurement are on the
order of 0.005% and 0.008%, respectively, nearly two orders of
magnitude better than the IRS 16 data. Further, the standard
deviation of our Luhman 16A/B measurements are 0.03%/0.02%
and 0.03%/0.03% (even in the presence of a potential
astrophysical signal) for Q and U, respectively, demonstrating
the increased accuracy and stability of our data set relative to
0.2% measured by Witzel et al. (2011). Given the preliminary
nature of the deBoer et al. (2014) system model and their final
error bars on the instrumental polarization, we also consider our
data to also be of superior quality given our error bars and the
scatter of our data points relative to their 0.4% errors. Overall,
the stability of our data set is well below the errors quoted in
both these previous works and emphasizes the superior quality
of our data set.

2.4.2. *Instrumental and Data Extraction Systematics*

While our data set may have the statistical power to detect a
signal of ~0.03%, before considering the signal real, instru-
mental systematics must be ruled out as a cause of the observed
change in polarization. Because the detection of a signal in our
data relies on a time-variable signal, here we focus on potential
systematics that change on the same timescale as the signal: on
the order of one or more hours. Instrumental or data extraction
biases occurring on faster timescales are fit for in our analysis
as an additional noise source (see Section 3), and the best-fit
values for this noise are visualized in Figure 4 as the red error
bars superimposed upon the photometric error bars in black.
These “fast” varying biases do not affect what appears to be the
signal in Luhman 16A. Such fast-varying biases could include,
for example, changing atmospheric conditions (i.e., seeing)
between two subsequent frames, which has the potential to
affect double-difference measurements, among other effects.

One major strength of this data set is the nearly equal-
magnitude binary nature of our targets. This allows us to rule
out many potential sources of bias, because they would affect
the measured polarization of both sources in the same way. For
example, this includes slowly varying misalignments in the
optical train, second-order flexure effects on the telescope and/
or instrument, and possible polarization induced by thin layers
of clouds. In particular, this also includes changes in the
systematic alignment between M3 and NaCo's rotation ring.

Spatially varying instrumental polarization is another
potential effect to be explored. In Section 2.3, we noticed a
spatially dependent instrumental polarization offset between
the two dither positions that we corrected for by averaging
together the *q* and *u* values from back-to-back dithers. We
searched for residual spatial dependence after this dither-
averaging by comparing *q* measurements against the measured

x and y detector positions of each source (using the x and y
positions in the ordinary beam) using the Spearman-r and
Kendall-Tau correlation coefficients. Because of the dither-
averaging, we calculated the correlation coefficients using the
mean x position for each dither pair. The results are
summarized in Table 3. The table also includes p-values for
the null hypothesis that each coefficient (or greater) would be
obtained with a random data set. However, the individual q
data sets for each source only contain 35 data points, and so
these values may not be fully reliable. Here we limit ourselves
to searching for correlations only in q because that is where the
more significant (potential) signal is seen.

The correlation coefficients immediately suggest that there is
a significant spatial correlation for the *q* measurements of
Luhman 16A. However, because we were operating in pupil
tracking and the two objects were rotating about their common
center, the source's spatial location is correlated with the
parallactic angle—as is the expected signal—and so it is hard to
draw any conclusions from this correlation. On the other hand,
there does not appear to be a significant correlation between the
*q* measurements of Luhman 16B and its detector position. If
there were a significant correlation, it would be difficult to
distinguish between an astrophysical signal and a spatially
dependent instrumental polarization, but the lack of correlation
strongly suggests that there is no significant residual spatially
dependent instrumental polarization after our dither-averaging.
As a secondary check for the lack of spatial dependence in B's
*q* values, we randomly shuffled the *q* measurements 1000 times
and measured the Spearman-r and Kendall-Tau values for each
shuffle. The correlation values for B presented in Table 3 fall
within 1σ of the mean of the sample for both Spearman-r and
Kendall-Tau metrics for both *x* and *y* detector positions,
confirming the lack of significant spatial dependence in B and
implying that the signal seen for A is also not due to any
spatially dependent instrumental polarization. We also searched
for correlations as a function of distance and angle from the
detector center, as well as distance and angle from the mean
position of the binary. In all cases, we measured a significant
spatial correlation for the Luhman 16A, but not for B, again
suggesting that there is no remaining spatial dependence in the
instrumental polarization. Finally, we revisit the spatial
dependence in Section 3.2, where we demonstrate that our
best-fit models are sufficient to explain the observed spatial
correlations.

In our data extraction procedure, the only tunable parameter
is the aperture size used to extract the flux of each star. Figure 8
displays the *q* and *u* measurements for A and B, extracted from
a range of aperture sizes. Although the individual data points
appear to vary slightly for different aperture sizes, the large-
scale trends appear to be consistent for all of the aperture sizes
explored. This also suggests that the perceived trends in *q*
cannot be attributed to a varying PSF overlap between the two
stars, as the amount of intensity that leaks into one object's
aperture from the other object changes with aperture size,
although the signal appears relatively constant.

Finally, we searched for correlations between our measure-
ments and the FWHM of each source. We first fit each PSF in
each image to a 2D Moffat profile and then extracted the
FWHM of the Moffat profile. We then took the average
FWHM across eight frames (the number of frames needed to
get one set of dither-average *q* and *u* measurements). Table 3
displays the Spearman-r and Kendall-Tau correlation
<figure><img src="image_7.png" /><figcaption>Figure 5. The measured twilight polarization (black), as a function of MJD on UT 2018 April 11. Shown in light green is the expected twilight polarization for randomly chosen system models from the accepted MCMC chains, given the known solar azimuth.</figcaption></figure>

coefficients for our *q* measurements against the FWHM in the ordinary beam. Here, we have compared the FWHM of A against the *q* measurements of A, and similarly for B. No significant correlations are found.

We conclude that the long-term variations seen in our data cannot be attributed to changing instrument systematics, observing conditions, or our own data extraction. Further, when compared to previous efforts to characterize the NaCo system, our data set represents a new standard in terms of depth and stability. Without any evidence to the contrary, hereafter we consider the the long-term trends seen in the data to be astrophysical in nature.

## 2.5. Twilight Data Reduction

Flat-fielding and dark correction for the twilight data were carried out as described for Luhman 16. For each twilight sky observation, the intensity in the ordinary and extraordinary beam was measured in the same region of the detector as the Luhman 16 observations by summing the counts in a rectangular aperture. The apertures covered the [x, y] regions of [50:450, 821:920] and [50:450, 693:792] for the ordinary and extraordinary beams, respectively. The $q_{NaCo}$ and $u_{NaCo}$ values were obtained for each HWP cycle, and the angle of linear polarization was calculated as $0.5 \arctan(2(u_{NaCo}/q_{NaCo}) + 45^{\circ})$, where the extra $45^{\circ}$ was added to compensate for the telescope azimuth position of $45^{\circ}$ east of north. The measured angle of linear polarization can be seen in Figure 5, along with the expected position angle for a range of accepted system models from the Markov Chain Monte Carlo (MCMC) fitting (Section 3.2). Errors on each angle measurement were estimated by propagating the standard deviation in each rectangular aperture through Equations (6)–(11) using standard error propagation techniques.

## 2.6. Elias 2-25 and HD 162973 Data Reduction

Flat-fielding, dark correction, and background subtraction were carried out on the Elias 2-25 (polarized standard) and HD 162973 (unpolarized standard) data as described for Luhman 16. Photometry was extracted from all observations of Elias 2-25 and HD 162973 using a similar method as for

Luhman 16 (Section 2.2). Measurements of $q_{NaCo}$ and $u_{NaCo}$ were then calculated using a similar method to Luhman 16, except that both $q_{NaCo}$ and $u_{NaCo}$ were both calculated using double differencing because observations with the HWP at 22.5° were obtained for these two targets. As with Luhman 16, the two dither positions were averaged. For HD 162973, we average together the measurements from the two different dates, as we expect the instrumental polarization to strongly dominate over any potential stellar polarization. The $q_{NaCo}$ and $u_{NaCo}$ measurements for Elias 2-25 and HD 162973 can be seen in Figures 6 and 7, respectively. Errors on $q_{NaCo}$ and $u_{NaCo}$ were propagated from the original photometric errors through Equations (6)–(11).

## 3. Analysis

The ultimate goal of our analysis was to determine the degree and angle of linear polarization of both Luhman 16A ($p_A$ and $\psi_A$) and B ($p_B$ and $\psi_B$). We began our analysis with a simple instrument model to characterize the degree of polarization measured by the instrument for each object, where we were not (yet) concerned with the true on-sky degree of polarization or angle of polarization. We assumed an instrument and telescope (i.e., $M_{inst}$ and $M_{tel}$ from Equation (3)) model that included instrumental polarization, but assumed perfect efficiencies and no $Q \leftrightarrow U$ crosstalks. With this model, the measured $q_{NaCo}$ and $u_{NaCo}$ at each parallactic angle ($\theta_{PA}$) can be described with only six parameters ($q_A$, $u_A$, $q_B$, $u_B$, $IP_Q$, and $IP_U$) by expanding the rotation matrix in Equation (3) and adding instrumental polarization:

$$q_A = q_{sky,A} \cos 2\theta_{PA} + u_{sky,A} \sin 2\theta_{PA} + IP_Q$$

$$u_A = q_{sky,A} \sin 2\theta_{PA} - u_{sky,A} \cos 2\theta_{PA} + IP_U$$

$$q_B = q_{sky,B} \cos 2\theta_{PA} + u_{sky,B} \sin 2\theta_{PA} + IP_Q$$

$$u_B = q_{sky,B} \sin 2\theta_{PA} - u_{sky,B} \cos 2\theta_{PA} + IP_U.$$

In this set of equations, we applied a negative sign to the right-hand sign for both $u_A$ and $u_B$ to compensate for the sign flips between $q_{NaCo}$ and $u_{NaCo}$ seen by Witzel et al. (2011) and deBoer et al. (2014). We fit for a joint solution for all six parameters from our full set of measurements with the SciPy `least_squares` function, using a “Cauchy” loss function to account for outlier data points. We display the results of this fit in Figure 4 as the “Simple Model,” demonstrating the detection of polarization in both A and B. The fit returned $q_A = -0.007\%$, $u_A = 0.027\%$, $q_B = 0.007\%$, $u_B = -0.006\%$, $IP_Q = -1.921\%$, and $IP_U = -0.324$. The degrees of linear polarization for A and B are $p_A = 0.028\%$ and $p_B = 0.009\%$.

While this simple model serves to highlight the detection of polarization in the two components, it fails to capture several important aspects when trying to accurately estimate the degrees and angles of linear polarization. First, in assuming a perfect system, we have failed to account for polarimetric efficiencies and crosstalks, as well as instrument angular offset, all of which will affect the estimates of the on-sky values. The polarimetric efficiencies and crosstalks would typically be measured using polarized standards observed on the same night as the observations. However, because our standards were observed after the HWP mechanism was replaced, a more complicated approach is required. From Figure 4 it is also clear that the formal photometric error bars do not represent the scatter in the data, and parameter errors estimated from the
## Elias 2-25

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** The measured $q_{NaCo}$ and $u_{NaCo}$ for the Elias 2-25 polarized standard observations. Also shown are the modeled $q_{NaCo}$ and $u_{NaCo}$ values for a random selection of walker positions from the MCMC fit, for the median posterior parameters from the MCMC fit, and for a perfect $M_{inst+tel}$ Mueller matrix (although including instrumental polarization). Error bars on the measurements are shown, including the extra error term from the MCMC fit, but are smaller than symbols. For each randomly selected walker position, a system Mueller matrix is generated and the modeled polarization is calculated by propagating the known polarization of Elias 2-25 through Equation (3) (which includes all instrumental polarization effects). Deviation from the perfect model is due to nonperfect efficiencies ($\eta$) and crosstalks ($X$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** The measured (black) and modeled (blue) instrumental polarization for the HD 162973 unpolarized data set. The expected instrumental polarization is shown as a 2D histogram and is generated by picking values of $IP_Q$ and $IP_U$ from 30,000 randomly selected walker positions in the final MCMC chain and then rotating them with the associated $\delta\theta_{HWP}$ value for that walker position. The model and the data agree to within $1\sigma$.</figcaption></figure>

least-squares approach above would be inaccurate (for this reason we do not report them). Further, because we were required to apply a nonideal correction of the $u_{NaCo}$ measurements, we expect their errors to be larger than for $q_{NaCo}$ and there may be nonstandard correlations between $q_{NaCo}$ and $u_{NaCo}$. In order to fully account for these effects, we developed an MCMC-based strategy to simultaneously fit for the polarization of Luhman 16A and B and instrumental effects, and to account for increased errors and correlations in the data.

### 3.1. Polarimetric Instrument Model

We first built up a Mueller matrix model for the combination of the instrument ($M_{inst}$) and the telescope ($M_{tel}$) to translate on-sky polarization to the $q_{NaCo}$ and $u_{NaCo}$ values from each observation. Although Mueller matrices for NaCo have been measured before both in the H and K bands (Witzel et al. 2011;

deBoer et al. 2014), both were obtained while NaCo was mounted on the UT4 telescope and neither used the pupil tracking mode. The move to a new telescope (UT1), the change in observing mode, and the increased depth of our data set motivated the development of a new system model.

In pupil tracking mode, the orientations of $M_{inst}$ and $M_{tel}$ do not change relative to each other, and $M_{inst}$ and $M_{tel}$ was therefore combined for all of our observations into a single Mueller matrix, which we call $M_{inst+tel}$:

$$ M_{\text{inst+tel}} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ IP_Q & \eta_Q & X_{U \to Q} & 0 \\ IP_U & X_{Q \to U} & \eta_U & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}, \quad (12) $$

where the $IP_{Q,U}$ elements represent the instrumental polarization (i.e., the total intensity that gets polarized, even for an unpolarized source), the $\eta$ terms represent efficiency terms, and the $X_{\to}$ terms represent crosstalk between $Q$ and $U$. For completeness, we included the fourth row that corresponds to circular polarization, to which NaCo is insensitive.

Technically, all 16 elements of the $4 \times 4$ matrix need to be calculated in order to back out the true on-sky Stokes vector of a source. However, a number of reasonable assumptions were made to limit the number of free variables. First, because NaCo is only sensitive to linear polarization, the bottom row of $M_{inst+tel}$ can be effectively set to zero. Second, we assumed that our target is not significantly circularly polarized, which is reasonable for brown dwarfs and most stars, allowing us to set the fourth column to be all zeros (e.g., Clarke 2010). Third, we assumed that any polarization does not affect the total intensity, setting all but the [1, 1] element of the first row to be zero. Thus, we were left with only six free parameters in $M_{inst+tel}$ that had to be found: $IP_Q$, $IP_U$, $\eta_Q$, $\eta_U$, $X_{U \to Q}$, and $X_{Q \to U}$.

We also included in our system model two additional free parameters: (i) a rotational offset between NaCo's frame of reference and the sky, $\delta\theta_{PA}$, that was included in the rotation matrix in Equation (3), $M_{rot}(\theta_{PA}+\delta\theta_{PA})$, and (ii) an extra rotation
<figure><img src="image_11.png" /><figcaption>**Figure 8.** A comparison of the measured *q* and *u* values as a function of aperture size. An aperture size of 17 pixels corresponds to the data displayed in Figure 4. The error bars represent the expected uncertainties due to photon read noise. The orange horizontal line indicates the best-fit instrumental polarization values from Section 3.2.</figcaption></figure>

angle for the HWP, $\delta\theta_{HWP}$, that was applied in Equation (4) to the HWP Mueller matrix, $M_{HWP}(\theta_{HWP}+\delta\theta_{HWP})$. This second offset, $\delta\theta_{HWP}$, represents the unknown encoder offset between the April and May/June data and was only applied to the data obtained after the intervention (i.e., observations of Elias 2-25 and HD 162973).

In our fitting procedure, the instrumental polarization is constrained mainly by the Luhman 16 observations themselves. The instrumental polarization is essentially the mean of the *q* and *u* measurements of the source, with any astrophysical signal presented as a modulation on top of that (modulated in Parallactic angle). Under this scheme, the primary function of including the unpolarized standard is to constrain the relative offset of the encoder between the April data and the later data, $\delta\theta_{HWP}$. This offset is constrained by the difference in the instrumental polarization angle measured by the Luhman 16 data set in April and the instrumental polarization angle measured by the unpolarized standard. This offset's main purpose is to connect the constraints on the *q* and *u* efficiencies and crosstalks obtained by the polarized standard to the April Luhman 16 data.

The instrumental polarization has the potential to evolve slightly on one- to two-month timescales (e.g., due to degradation of mirror coatings); this evolution could affect the

$\delta\theta_{HWP}$ value, the crosstalks, and the efficiencies, and in turn the final derived degrees and angles of polarization for A and B. If the unpolarized standard were being used to constrain the instrumental polarization in the Luhman 16 data set, the potential for this evolution would be of significant concern in interpreting our signal. However, in our case, the unpolarized standard is not being used in such a manner. The evolution of the polarimetric efficiencies and crosstalks on this timescale is likely to be on the order of a few percent or smaller. These effects are relative to the polarization signal itself, and given the low S/N of the polarization signal of Luhman 16A and B, we anticipate that the final degrees and angles of polarizations will be dominated by statistical errors rather than any systematic offsets introduced by this temporal evolution.

Residual fast-varying polarimetric systematics may unduly increase the scatter of our data beyond the read-out and photon noise, especially for the $u_{NaCo}$ measurements of Luhman 16, where double differencing was not possible. To characterize these systematics, we included two extra error terms, one for the single-differenced data (i.e., $u_{NaCo}$ for Luhman 16 and the twilight observations), and one for the double-differenced data: $\sigma_{SD}$ and $\sigma_{DD}$, respectively. In practice, we fit for $\log(\sigma_{SD}^2)$ and $\log(\sigma_{DD}^2)$. Finally, to account for possible covariance between measurements taken within the same HWP cycle in the
<table><caption>**Table 4**
Summary of Model Parameters</caption><thead><tr><th>Parameter</th><th>Description</th><th>Constraining Data Set</th><th>Affected Data Set</th></tr></thead><tbody><tr><td>p<sub>A</sub></td><td>Degree of Linear Polarization for Luhman 16A</td><td>Luhman 16</td><td>None</td></tr><tr><td>ψ<sub>A</sub></td><td>Angle of Linear Polarization for Luhman 16A</td><td>Luhman 16</td><td>None</td></tr><tr><td>p<sub>B</sub></td><td>Degree of Linear Polarization for Luhman 16B</td><td>Luhman 16</td><td>None</td></tr><tr><td>ψ<sub>B</sub></td><td>Angle of Linear Polarization for Luhman 16B</td><td>Luhman 16</td><td>None</td></tr><tr><td>IP<sub>Q</sub></td><td>Instrumental polarization—Q</td><td>Luhman 16</td><td>Twilight, HD 162973, Elias 2-25</td></tr><tr><td>IP<sub>U</sub></td><td>Instrumental polarization—U</td><td>Luhman 16</td><td>Twilight, HD 162973, Elias 2-25</td></tr><tr><td>η<sub>Q</sub></td><td>Q efficiency</td><td>Elias 2-25</td><td>Luhman 16, Twilight</td></tr><tr><td>η<sub>U</sub></td><td>U efficiency</td><td>Elias 2-25</td><td>Luhman 16, Twilight</td></tr><tr><td>X<sub>Q→U</sub></td><td>U → Q crosstalk</td><td>Elias 2-25</td><td>Luhman 16, Twilight</td></tr><tr><td>X<sub>U→Q</sub></td><td>Q → U crosstalk</td><td>Elias 2-25</td><td>Luhman 16, Twilight</td></tr><tr><td>δθ<sub>PA</sub></td><td>Rotational offset between sky and NaCo</td><td>Twilight</td><td>Luhman 16, Elias 2-25</td></tr><tr><td>δθ<sub>HWP</sub></td><td>HWP offset after intervention</td><td>HD 162973</td><td>Elias 2-25</td></tr><tr><td>σ<sub>SD</sub></td><td>Extra Systematic Error—Single Difference</td><td>Luhman 16, HD 162973</td><td>Gives more realistic parameter errors</td></tr><tr><td>σ<sub>DD</sub></td><td>Extra Systematic Error—Double Difference</td><td>Luhman 16, HD 162973</td><td>Gives more realistic parameter errors</td></tr><tr><td>C<sup>u</sup><sub>q</sub></td><td>Covariance between q<sub>NaCo</sub> and u<sub>NaCo</sub></td><td>Luhman 16</td><td>Luhman 16</td></tr><tr><td>C<sup>q,B</sup><sub>q,A</sub></td><td>Covariance between q<sup>A</sup><sub>NaCo</sub> and q<sup>B</sup><sub>NaCo</sub></td><td>Luhman 16</td><td>Luhman 16</td></tr><tr><td>C<sup>u,B</sup><sub>q,A</sub></td><td>Covariance between q<sup>A</sup><sub>NaCo</sub> and u<sup>B</sup><sub>NaCo</sub></td><td>Luhman 16</td><td>Luhman 16</td></tr><tr><td>C<sup>u,B</sup><sub>u,A</sub></td><td>Covariance between u<sup>A</sup><sub>NaCo</sub> and u<sup>B</sup><sub>NaCo</sub></td><td>Luhman 16</td><td>Luhman 16</td></tr></tbody></table>

Luhman 16 data set, we include four covariance terms in our fit: (i) $C_q^u$, covariance between $q_{NaCo}$ and $u_{NaCo}$ (the same value for both Luhman 16A and B measurements), (ii) $C_{q,A}^{q,B}$, covariance between $q_A$ and $q_B$, (iii) $C_{q,A}^{u,B}$, covariance between $q_A$ and $u_B$, and (iv) $C_{u,A}^{u,B}$, covariance between $u_A$ and $u_B$. Although typically one would expect the $q$ and $u$ measurements to be independent, the correction of the systematic bias term $\epsilon$ for the $u_{NaCo}$ measurements in a given HWP cycle relies on the $\epsilon$ measurements of $q_{NaCo}$, motivating the inclusion of the covariance terms.

## 3.2. MCMC Model Fitting and Results

To extract accurate polarization measurements and errors, we adopted a strategy of simultaneously fitting for Luhman 16’s polarization and our system model parameters with a Bayesian MCMC approach, using all four of our data sets as input (Luhman 16, Twilight polarization, Elias 2-25, and HD 162973). Our full model includes 18 parameters, summarized in Table 4. Our choice of modeling strategy was largely driven by the many relationships between the model parameters and our different data sets, in addition to the need to compensate for the extra systematic error terms and covariances. For each model parameter, Table 4 summarizes the most constraining data set, as well as the other data sets whose interpretations are affected by that model parameter.

We constructed our log-likelihood function ($\ln p(y|\Theta)$) for the MCMC fitting as the sum of four components as follows:

$$ \ln p(y|\Theta) = -0.5\sum ((y - f(\Theta))^T C^{-1}(y - f(\Theta)) \\ + \ln(\det C))_{\text{Luhman 16}} \\ - 0.5\sum ((y - f(\Theta))^T C^{-1}(y - f(\Theta)) \\ + \ln(\det C))_{\text{Elias 2-25}} \\ - 0.5\sum ((y - f(\Theta))^T C^{-1}(y - f(\Theta)) \\ + \ln(\det C))_{\text{HD 162973}} \\ - 0.5\sum ((y - f(\Theta))^T C^{-1}(y - f(\Theta)) \\ + \ln(\det C))_{\text{Twilight}}, \quad (13) $$

where each sum is over the data (y), model ($f(\Theta)$), and covariance matrix (C) associated with each subscript for a given parameter set, $\Theta$. For the Luhman 16, Elias 2-25, and HD 162973 data sets, $y$ included both $q_{NaCo}$ and $u_{NaCo}$ measurements. The model $q_{NaCo}$ and $u_{NaCo}$ measurements ($f(\Theta)$) were calculated from Equation (3) for a given input $q_{sky}$ and $u_{sky}$ and included the $M_{inst+tel}$ Mueller matrix (Equation (12)), as well as the offset parameters $\delta\theta_{PA}$ and $\delta\theta_{HWP}$. For Luhman 16A and B, [$q_{sky}$, $u_{sky}$] were calculated from the free parameters [$p_A$, $\psi_A$] and [$p_B$, $\psi_B$] as

$$ \begin{aligned} q_{\text{sky}} &= p \cos(2\psi) \\ u_{\text{sky}} &= p \sin(2\psi). \end{aligned} \quad (14) $$

Rather than fitting for $q_{sky}$ and $u_{sky}$ and calculating $p$ and $\psi$ afterwards, we chose to fit for $p$ and $\psi$ directly and forward model into $q_{sky}$ and $u_{sky}$, which allowed us to directly obtain a posterior distribution of $p$, avoiding the positive bias associated with calculating a single value of $p = \sqrt{q^2 + u^2}$ in the presence of noise on $q$ and $u$.

For the Elias 2-25 data set, the model values for $q_{NaCo}$ and $u_{NaCo}$ (i.e., ($f(\Theta)$) were calculated from $q_{sky}$ and $u_{sky}$, using the known $p = 4.13\% \pm 0.02\%$ and $\psi_{PA} = 24^\circ \pm 1^\circ$ (Whittet et al. 1992). HD 162973 has a measured B-band polarization of $p = 0.09\% \pm 0.055\%$ (Mathewson & Ford 1970). Assuming a Serkowski polarization law (Serkowski 1973), with a maximum polarization of 0.09 and a peak wavelength of 0.55 µm, we estimated a polarization of ~0.02% at 1.6 µm. Given that the error bars on the measurement are significantly larger than this value, we considered the intrinsic source polarization to be negligible and set $q_{sky}$ and $u_{sky}$ both to zero (further justifying our choice to average together the two measurements at different parallactic angles). Our choice of 0.55 µm sits roughly in the center of the range of values found in Whittet et al. (1992). If the true value of the peak polarization wavelength is smaller, then the polarization of HD 162983 in the H band will also be smaller. Whittet et al. (1992) find peak polarization wavelengths as high as ~0.8 µm. If it were that high for HD
<table><caption>**Table 5**
Best-fit Parameters</caption><thead><tr><th>Parameter</th><th>Best-fit</th><th>Uniform Prior</th></tr></thead><tbody><tr><td>p<sub>A</sub> (%)</td><td>0.031 ± 0.004</td><td>(0.0, 1.0)</td></tr><tr><td>ψ<sub>A</sub> (°)</td><td>-32 ± 4</td><td>(-180.0, 0.0)</td></tr><tr><td>p<sub>B</sub> (%)</td><td>0.010 ± 0.004</td><td>(0.0,1.0)</td></tr><tr><td>ψ<sub>B</sub> (°)</td><td>73<sup>+13</sup><sub>-11</sub></td><td>(0.0, 180.0)</td></tr><tr><td>IP<sub>Q</sub> (%)</td><td>-1.919 ± 0.002</td><td>(-2.5, 0.0)</td></tr><tr><td>IP<sub>U</sub> (%)</td><td>-0.327 ± 0.004</td><td>(-1.0, 0.0)</td></tr><tr><td>η<sub>Q</sub></td><td>-0.91 ± 0.02</td><td>(-1.0, 1.0)</td></tr><tr><td>η<sub>U</sub></td><td>0.70 ± 0.02</td><td>(-1.0, 1.0)</td></tr><tr><td>X<sub>Q→U</sub></td><td>0.06 ± 0.07</td><td>(-1.0, 1.0)</td></tr><tr><td>X<sub>U→Q</sub></td><td>0.17<sup>+0.09</sup><sub>-0.08</sub></td><td>(-1.0, 1.0)</td></tr><tr><td>δθ<sub>PA</sub> (°)</td><td>-6<sup>+2</sup><sub>-2</sub></td><td>(-25.0, 25.0)</td></tr><tr><td>δθ<sub>HWP</sub> (°)</td><td>27.9<sup>+0.3</sup><sub>-0.3</sub></td><td>(-180.0, 180.0)</td></tr><tr><td>σ<sub>SD</sub> (%)</td><td>0.02 ± 0.02</td><td>(4.5 × 10<sup>-3</sup>, 0.7)</td></tr><tr><td>σ<sub>DD</sub> (%)</td><td>0.03 ± 0.03</td><td>(4.5 × 10<sup>-3</sup>, 0.7)</td></tr><tr><td>C<sup>u</sup><sub>q</sub></td><td>0.15 ± 0.14</td><td>(-1, 1)</td></tr><tr><td>C<sup>q,B</sup><sub>q,A</sub></td><td>0.0 ± 0.2</td><td>(-1, 1)</td></tr><tr><td>C<sup>u,B</sup><sub>q,A</sub></td><td>0.0 ± 0.1</td><td>(-1, 1)</td></tr><tr><td>C<sup>u,B</sup><sub>u,A</sub></td><td>0.4 ± 0.1</td><td>(-1, 1)</td></tr></tbody></table>

162983, the polarization in the *H* band would be ~0.05%, which is still smaller than our error bars.

In the twilight portion of Equation (13), the data y is the measured position angle and the model position angle ($f(Θ)$) was calculated from the model $q_{\text{NaCo}}$ and $u_{\text{NaCo}}$, where the input $q_{\text{sky}}$ and $u_{\text{sky}}$ were calculated with Equation (14) assuming $p = 1.0$ (this has no effect on the results) and $\psi_{\text{PA}}$ is 90° away from the solar azimuth at the time of observation.

For all data sets except the twilight position angle measurements, the covariance matrices were populated with diagonal variance terms calculated as the square of the original photon and read-out noise estimates for $q_{\text{NaCo}}$ and $u_{\text{NaCo}}$ plus either $\sigma_{\text{SD}}^2$ or $\sigma_{\text{DD}}^2$, depending on whether the data were measured as a single difference or double difference. For the Luhman 16 data set, the four covariance terms ($C_q^u$, $C_{q,A}^{q,B}$, $C_{q,A}^{u,B}$, $C_{u,A}^{u,B}$) multiplied by the square of the two associated diagonal terms were included as off-diagonal terms. With this definition, the expected range of the covariance parameters should range from -1 to 1. The covariance matrix for the twilight data was populated with diagonal terms corresponding to the square of errors on the angles of polarization.

To sample our parameter posterior distributions we employed the python-based ensemble-sampling MCMC package `emcee` (Foreman-Mackey et al. 2013). Priors for all the model parameters can be found in Table 5. The `emcee` ensemble-sampler was run for 50,000 steps with 256 walkers, after a burn-in of 1000 steps. After the run, we measured a maximum autocorrelation across all parameters of 48.4 steps, verifying that the chains should have reached equilibrium within the burn-in phase (~O(10) autocorrelation times are needed for convergence; Foreman-Mackey et al. 2013).

Posterior distributions were estimated using 1 out of every 49 chains, to ensure statistical independence. The full “corner” plot showing the marginalized and joint probability distributions can be seen in Figure A4 (Foreman-Mackey 2016). All parameters appear single-peaked and mostly Gaussian in the marginalized posterior distributions, with the exception of $\eta_U$ and $\eta_Q$, which both show a slightly skewed tail to smaller absolute values. The joint posterior distributions show

correlations between $\eta_Q$, $\eta_U$, $X_{Q→U}$, $X_{U→Q}$, $\delta\theta_{\text{PA}}$, and $\delta\theta_{\text{HWP}}$, but all parameters appear well constrained in the marginalized posteriors. Table 5 summarizes the fitting results, where the median value from each marginalized posterior distribution is taken as the best-fit value, and the upper and lower $1\sigma$ values were taken as the 16% and 84% percentile values (corresponding to a confidence interval of 68%). As demonstrated in Figures 4–7, the model appears to fit all of our input data well.

From the results of the fitting, we detect polarization in both Luhman 16A and B at the $8\sigma$ and $2.5\sigma$ levels, respectively. For Luhman 16A, we find a linear degree of polarization of $p_A = 0.031\% \pm 0.004\%$ with an angle of $\psi_A = -32^\circ \pm 4^\circ$. For Luhman 16B, we find a linear degree of polarization of $p_B = 0.010\% \pm 0.004\%$ with an angle of $\psi_B = 73^{+13}_{-11}^\circ$. As expected, we find a higher value of $\sigma_{\text{SD}}$ than $\sigma_{\text{DD}}$. In all cases, we find that the data is well fit by the model. Our fitting process recovers a weak covariance between the $q_{\text{NaCo}}$ and $u_{\text{NaCo}}$ measurements in the Luhman 16 data set, $C_q^u$, that we attribute to the correction of the $\epsilon$ systematic for the $u_{\text{NaCo}}$ measurements. We also recover to higher significance a covariance between $u_A$ and $u_B$, which we attribute to residual uncorrected systematics in both $u$ measurements as a result of the nonoptimal correction of $\epsilon$. Random walker positions were selected from the final parameter chains and have been overplotted on the measurements of all four data sets in Figures 4–7 and A1 and A2.

Our Mueller matrix system model parameters qualitatively agree well with those found by previous characterization (Witzel et al. 2011; deBoer et al. 2014), modulo several sign flips that may be attributable to different sign conventions. Exact agreement was never expected, due to the change of telescope and aging mirror coatings. Our modeling strategy is able to achieve errors on the efficiencies and crosstalks similar to those of deBoer et al. (2014), where they are able to constrain them. However, the accuracy of our constraints on the instrumental polarization is two orders of magnitude greater than either of these two works. This can be attributed to the depth and stability of our data set (see Section 2.4.2), the simplicity of our instrument setup, and our self-calibration strategy. In contrast to previous modeling efforts that had to model M3 and NaCo (upstream of the HWP) separately, operating in pupil-tracking mode has allowed us to consider only a single Mueller matrix when modeling our system, therefore simplifying the observations required to constrain it. Our data-driven modeling strategy also contrasts with that of Witzel et al. (2011; and in part deBoer et al. 2014) in that we leave the relevant elements of our Mueller matrix as completely free parameters that we fit to the data, whereas they rely on specific function forms and assumed material parameters, giving our model more freedom to fit the data.

Although the $2.5\sigma$ detection of Luhman 16B is of low S/N, we believe that it is real to within the accuracy of our system model architecture. To verify this, we explored several different model comparison metrics for different model setups. For all model setups, we considered the reduced $\chi^2$, the Aikake Information Criterion (AIC), and the Bayesian Information Criterion (BIC). For our default model, we calculated all three parameters for the fit described above, as well as the same fit, but setting $p_B = 0.0$. For this second case, we reduced the number of parameters by two, as both $p_B$ and $\psi_B$ need to be removed. We also considered a model where the instrumental polarizations ($IP_Q$ and $IP_U$) were separate free parameters for
<table><caption>**Table 6**
Model Comparison</caption><thead><tr><th>Test Model</th><th>red - χ<sup>2</sup></th><th>AIC</th><th>BIC</th><th>N<sub>Parameters</sub></th></tr></thead><tbody><tr><td>Default Model</td><td><b>1.08</b></td><td>-2033</td><td>-1984</td><td>18</td></tr><tr><td>......p<sub>B</sub> = 0</td><td>1.13</td><td>-2025</td><td>-1974</td><td>16</td></tr><tr><td>Independent IPs</td><td><b>1.08</b></td><td>-2031</td><td>-1975</td><td>20</td></tr><tr><td>......p<sub>B</sub> = 0</td><td>1.13</td><td>-2022</td><td>-1967</td><td>18</td></tr><tr><td>No B polarization</td><td>1.09</td><td>-2018</td><td>-1963</td><td>18</td></tr></tbody></table>

**Note.** The best-fit models are shown in bold.

the A and B measurements. After rerunning the MCMC fit under this new assumption, we calculated all three metrics both with and without $p_B$. Under this new assumption, we measured the same polarizations of A and B as those shown in Table 5. Finally, we consider a model where the polarization of B is forced to be zero from the start of the fitting procedure. The goodness-of-fit parameters for all model setups can be found in Table 6, with the best-fit models shown in bold. When fitting for separate IP values for A and B, we find the best-fit instrumental polarization values for each component to be within $1 - \sigma$ of the values reported in Table 5, and the instrumental polarization measured from A and B are within $2 - \sigma$ of each other. For all model comparison metrics, we find that the default model that includes a polarization of B is the best model to describe the data.

Using the best-fit system model, we inverted Equation (3) to obtain the sky-plane Q and U values as a function of time for A and B. The residuals displayed similar features in both Luhman 16A and B, which we attribute to uncorrected time-varying systematics likely related to the larger systematics in *u* (and encapsulated in $\sigma_{SD}$). We searched for variability using a Lomb–Scargle periodogram, but found no significant peaks. We conclude that we have not detected any polarimetric variability.

As a secondary check on our fitting results, we split our data in half and reran the fitting procedure. The data was split such that one of the fits included the first and third quarters of the data, and the second fit included the second and fourth quarters. Rather than splitting the data at the halfway mark, splitting the data up by quarters was done to ensure that each data set had sufficient diversity in parallactic angle. The new fits resulted in values of $p_A$, $\psi_A$, $p_B$, $\psi_B$, $IP_Q$, and $IP_U$ that agreed with the mean values presented in Table 5 to within $1 - \sigma$ (using the newly fit error bars). As expected, the new error bars were $\sim\sqrt{2}$ times worse than the original error bars.

### 3.3. Spatial Correlations

Due to our pupil tracking observing mode, the two binary components rotate around each other's center of light throughout our observation set. This naturally introduces a correlation between detector position and parallactic angle, and in turn, between detector position and the instrument-frame *q* and *u* measurements, $q_{NaCo}$ and $u_{NaCo}$. Here we revisit the question of whether or not the spatial correlation of $q_A$ measured in Section 2.4.2 can be explained by an astrophysical signal.

We began by generating fake *q* data sets for A and B, given the best-fit model of Table 5, sampling the model at the parallactic angles corresponding to our observations. We then injected noise into the fake data sets by replacing each data

point with a draw from a Gaussian distribution centered on the model value and with a width equal to the photometric errors and best-fit $\sigma_{DD}$ value added in quadrature. This procedure was repeated 1000 times, and Spearman-r and Kendall-Tau coefficients were measured for each fake data set with respect to the real detector *x* and *y* positions. Figure 9 displays the resulting distributions of the coefficients for this “noisy model.” These histograms represent the range of correlation coefficients one might expect given the final errors in our model fitting (e.g., Curran 2014).

Next, we estimated the distribution of the Spearman-r and Kendall-Tau values for the *q* measurements of Luhman 16A/B against the detector *x* and *y* positions by perturbing the original *q* measurements 1000 times. For each iteration, each *q* measurement (for both A and B) was replaced with a new value drawn from a Gaussian distribution with a mean equal to the measurement value and a width equal to the photometric error (Figure 9). As with the model data, this perturbation approach should represent the range of values expected in the two correlation coefficients given our errors. The Spearman-r and Kendall-Tau distributions for the perturbed data overlap significantly with the model distributions, suggesting that the posterior distributions recovered from the model fitting procedure are sufficient to explain the measured spatial correlation of the data.

## 4. Atmospheric Modeling

In this section, we consider what possible physical phenomena could result in the measured polarizations of Luhman 16A and B. The constant polarization that we measure is integrated over our entire 7 hr observing window, corresponding to greater than one rotation period for each component, or nearly a full period if Luhman 16A's period is 8 hr. The measured signals cannot be attributed to longitudinally varying features, such as patchy clouds or spatially varying bands, because they would cause the degrees and angles of polarization to change in time as the features rotated in and out of view and would manifest as variability in our residuals. However, we cannot rule out lower-level variable features below our detection limits that may be superimposed upon the constant signal. Polarization from a circum-brown dwarf disk can also be ruled out due to the lack of any previous evidence of extra dust in the SED. Thus, oblateness (Section 4.1) and constant cloud bands (Section 4.2) are the only remaining possible sources of polarization.

### 4.1. Polarization due to Oblateness

We first considered the polarization signal that would be caused by oblateness in the case of a homogeneous atmosphere and cloud cover. Oblateness as a function of effective temperature was calculated using updated evolutionary models for objects with masses of 27.2, 29.3, 31.4, and 34.6 $M_{Jup}$ (M. S. Marley et al. 2020, in preparation). Each model track provided the radius, moment of inertia, and effective temperature as a function of age for a given mass. The spin angular velocity was calculated for 5 and 8 hr periods, for each radius in each evolutionary track. Oblateness as a function of effective temperature was calculated using the Darwin-Radau formula that relates an object's spin angular velocity, mass, radius, and moment of inertia to its oblateness (Barnes & Fortney 2003). These curves were then interpolated to the known mass ranges
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Spearman-r and Kendall-Tau correlation values for $q_A$ and $q_B$ against the x and y detector position. The green and red distributions were calculated from our best-fit model, and the blue and orange distributions were calculated from the data. The solid and dashed black lines represent the values reported in Table 3 for $q_A$ and $q_B$, respectively.</figcaption></figure>

for Luhman 16A and B (Figure 10). The 5 hr period was considered for both Luhman 16A and B, and the 8 hr period was only considered for Luhman 16A. Using the measured effective temperatures from Faherty et al. (2014), we constrained Luhman 16B’s oblateness to $0.0131_{-0.0015}^{+0.0016}$ and Luhman 16A’s oblateness to either $0.0092 \pm 0.0009$ for a 5 hr period or $0.0036_{-0.0004}^{+0.0003}$ for an 8 hr period. We note that all of these values are smaller than Jupiter’s oblateness of 0.0649 for a rotation period of 9.9 hr (e.g., Dutta et al. 2009).

We then calculated polarization as a function of oblateness and line-of-sight inclination using a radiative transfer code previously applied to model polarized brown dwarfs (Sengupta & Marley 2009, 2010; Marley & Sengupta 2011). We define an inclination of $i = 0^\circ$ to be when a brown dwarf is viewed equator-on and $i = 90^\circ$ is when viewed pole-on. To explore a representative parameter space for Luhman 16A and B, we

considered models for 1200 and 1300 K objects, each with a sedimentation parameter $f_{sed} = 1$ (thicker clouds) and $f_{sed} = 3$ (thinner clouds), generated from the Ackerman & Marley (2001) cloud code. Figure 11 displays the expected polarization due to oblateness for the $f_{sed} = 1$ case at the inclinations measured for Luhman 16A (Karalidi et al. 2016) and Luhman 16B (Crossfield 2014). For all inclinations and masses, the $f_{sed} = 3$ models result in a polarization an order of magnitude less than what was measured. Given the oblateness ranges calculated above, the $f_{sed} = 1$ models predict a polarization in the range $p_B = 0.026\%–0.033\%$ for Luhman 16B, and either $p_A = 0.008\%–0.010\%$ or $p_A = 0.007\%–0.009\%$ for Luhman 16A for a rotation period of 8 hr and 5 hr, respectively.

It is clear that for both Luhman 16A and B, the $f_{sed} = 1$ models do not reproduce the measurements; the models
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Oblateness as a function of effective temperature for the range of masses and rotation periods relevant to Luhman 16A and B. The dashed colored lines show the oblateness curves for the upper and lower mass ranges for A and B. Luhman 16B is assumed to have a period of 5 hr and A's period is either 5 hr or 8 hr. The shaded regions highlight the effective temperatures measured by Faherty et al. (2014). Shown in gray are the original evolutionary models that were interpolated to obtain the dashed colored regions.</figcaption></figure>

underpredict the polarization of Luhman 16A and overpredict the polarization of Luhman 16B. For Luhman 16B, this discrepancy could potentially be attributed to the presence of clouds thinner than those produced by the $f_{sed} = 1$ models. Models with higher values of $f_{sed}$ (i.e., thinner clouds) result in a decrease in the expected polarization (Sengupta & Marley 2010), and previous fits to the observed spectra of Luhman 16B do indeed suggest values higher than $f_{sed} = 1$. The best-fit model to the HST/WFC3 0.8–1.6 µm time-resolved varying spectra is a linear combination of two models with $f_{sed} = 1$ and $f_{sed} = 3$ (each component with its own effective temperature), with the thinner cloud component ($f_{sed} = 3$) always having a higher coverage fraction than the thicker clouds (Buenzli et al. 2015a, 2015b). Although our polarimetric observations of Luhman 16B are consistent with thinner clouds, our data alone cannot definitively distinguish between oblateness-induced polarization with thinner clouds and polarization due to banded clouds (see Section 4.2).

On the other hand, we cannot resolve the underprediction of the oblateness model (with $f_{sed} = 1$) relative to the Luhman 16A measurement by invoking thicker clouds (i.e., $f_{sed} < 1$) without contradicting previous spectroscopic measurements. Model fitting to HST/WFC3 spectra always results in $f_{sed} \ge 1$; depending on the wavelength region analyzed, HST/WFC3 data of Luhman 16A are best fit by either homogeneous atmospheres with $f_{sed} = 2$ (Buenzli et al. 2015b), or a patchy model with a combination of $f_{sed} = 1$ and $f_{sed} = 3$ (Buenzli et al. 2015a). There is no evidence that a thicker, $f_{sed} < 1$, cloud model applies to Luhman 16A. If anything, the $f_{sed} > 1$ values found in these studies should reduce the expected polarization due to oblateness. Another mechanism other than oblateness is needed to explain the polarization of Luhman 16A.

## 4.2. Polarization due to Banded Clouds

Cloud banding in brown dwarfs has recently been inferred from modeling of quasi-periodic photometric variability (Apai et al. 2017) and has been predicted by general circulation models (Zhang & Showman 2014). Within the solar system, cloud banding can easily be seen in Jupiter, Saturn, and Neptune. de Kok et al. (2011) first showed that a single cloud band in a giant exoplanet/brown dwarf can easily cause

degrees of linear polarizations greater than 0.5%. Further modeling by Stolker et al. (2017) demonstrated that depending on the cloud band distributions and oblateness of the source, one could achieve degrees of linear polarizations up to ~1.33%.

For a given object, the net polarization at a given wavelength will depend on the line-of-sight inclination, the oblateness, and the properties of the cloud bands (i.e., the number of bands, their widths, their longitudes, and their relative polarization compared to the background atmosphere). Polarization from banded clouds can either supplement the net polarization from oblateness or work against it. For example, for a brown dwarf viewed equator-on, if polarizing clouds are found in an equatorial band, the angle of polarization from the clouds would be parallel to the direction of the pole projected on the sky and would “cancel out” some of the polarized intensity coming from oblateness, which is oriented perpendicular to the spin axis (Stolker et al. 2017). Depending on the relative strength of the polarization from the cloud bands versus that from oblateness, the polarization can (a) align itself with the spin axis, if equatorial clouds dominate the polarization, (b) be zero if the polarization from clouds and oblateness perfectly match, or (c) align itself perpendicular to the spin axis if oblateness dominates. Alternatively, if the polarizing clouds are concentrated at the poles (e.g., similar to the hazes in Jupiter poles), the clouds will only add to the oblateness-induced polarization, and the angle of polarization will always be perpendicular to the spin axis (Stolker et al. 2017). In Jupiter, these polar hazes are thought to originate from interactions with its aurora (West et al. 2004). For Luhman 16, $H\alpha$ (Faherty et al. 2014) and radio (Osten et al. 2015) emission nondetec-tions suggest negligible auroral activity, and as a result, we rule out polar hazes or clouds as possible sources of polarization.

Here we explore banding in a spherical atmosphere and demonstrate that using known characteristics of the system, we can recreate the measured polarization, without the need for significant model tuning. To estimate the number of bands to consider, we calculated the expected number of atmospheric jets given a characteristic horizontal wind speed, assumed radius, and rotation rate using a simple atmospheric scaling from Showman et al. (2010) based on 2D models. We assumed a radius of $1R_{Jup}$, and considered both 5 and 8 hr rotation periods. We calculated the number of jets for wind speeds of 10, 50, 100, and $200\ m\ s^{-1}$, consistent with the expected ranges (Table 7; Showman & Kaspi 2013). We find that for a five-hour period, we can expect at least a single jet (i.e., more than 0.5 jets) for wind speeds <200 m s⁻¹ and over two jets if the wind speeds are <10 m s⁻¹. For an eight-hour rotation period, we get one jet for wind speeds of <120 m s⁻¹ for rotation periods of 8 hr and 1.8 jets for wind speeds of 10 m s⁻¹.

Given these results, we calculated the expected degree of linear polarization for cloud configurations with one, two, and three bands at brown dwarf inclinations of 20° and 56° (roughly corresponding to the possible inclinations of Luhman 16A and B). Although the above calculations do not predict three jets for the wind speeds we explored, we include a three-band model as a proxy for the cloud bands inferred for 2MASS J13243553+6358281 (Apai et al. 2017). We assumed that the atmospheres of Luhman 16A and B are composed of two components: a “background” atmosphere at 1100 K with thicker clouds ($f_{sed} = 1$) and a “bands” atmosphere component at 1300 K with thinner clouds ($f_{sed} = 3$; Buenzli et al. 2015b).
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Polarization as a function of oblateness for $T_{\text{eff}} = 1200$ K (dashed lines) and 1300 K (solid lines) for the inclinations relevant for Luhman 16A (left) and Luhman 16B (right). All models use a sedimentation parameter $f_{\text{sed}} = 1$. Thickened solid lines highlight the oblateness ranges determined in Figure 10 for the 1300 K curves. Also shown as horizontal bands are the $\pm 1\sigma$ ranges for the polarization measured from Luhman 16A (orange) and B (blue).</figcaption></figure>

<table><caption>**Table 7**
Number of Jets</caption><thead><tr><th rowspan="2">Period (hr)</th><th colspan="4">Wind Speed</th></tr><tr><th>10 m s<sup>-1</sup></th><th>50 m s<sup>-1</sup></th><th>100 m s<sup>-1</sup></th><th>200 m s<sup>-1</sup></th></tr></thead><tbody><tr><td>5</td><td>2.2</td><td>1.0</td><td>0.7</td><td>0.5</td></tr><tr><td>8</td><td>1.8</td><td>0.78</td><td>0.6</td><td>0.4</td></tr></tbody></table>

We first used the Marley et al. (1996) radiative transfer code to calculate the temperature-pressure and composition profiles of the two atmosphere components, as well as their cloud properties. We injected these profiles into the polarization-sensitive radiative transfer code used in de Kok et al. (2011), which uses a doubling-adding method that fully includes all orders of scattering and polarization to calculate the polarization of the two atmosphere components in the H band. Clouds were modeled using Mie theory, as is standard in brown dwarf and exoplanet cloud modeling.

We modeled the atmospheres of Luhman 16A and B, assuming either one band spanning $-10^\circ$ to $+10^\circ$ of latitude, or two bands spanning $-50^\circ$ to $-30^\circ$, and $+30^\circ$ to $+50^\circ$ degrees of latitude. The models with three bands included the bands from both the one-band and two-band models. Net degrees of polarization were calculated for an inclination of $\sim 20^\circ$ for the models with one, two, or three bands and $\sim 56^\circ$ for the two-band model. In order to explain the quasi-periodic photometric variability, previous modeling relied on a sinusoidal modulation of the band brightness with longitude (Apai et al. 2017). Because our observations were averaged over a rotation period, we assumed that the bands were rotationally invariant to the extent of our time resolution and sensitivity. Our model atmosphere was split in a grid of $2^\circ \times 2^\circ$ resolution in latitude-longitude, and each pixel was assigned one model ("band" or "background atmosphere"). We assumed that our pixels were large enough to be able to ignore adjacency effects, i.e., light that is scattered within more than one pixel. For each configuration, we integrated the spatially resolved polarized fluxes across the observed disk and divided by the total intensity fluxes to get the degree of linear polarization (Table 8, Figure 12).

With these models, we can recreate the measured values to within a factor of $\sim 1.5$. By adjusting our two-banded model to have bands with widths of $32^\circ$ centered at latitudes of $-35^\circ$ and $+35^\circ$, we can reproduce the measured degree of linear polarization of 0.03% (shown in bold in Table 8). This configuration was found via a nonexhaustive, manual exploration of possible cloud configurations, and we expect that it is not a unique solution. Further, these models do not currently simultaneously include oblateness and cloud bands. Nonetheless, these models demonstrate that banded clouds are easily able to explain the level of polarization measured in Luhman 16A.

Considering the disagreement between the measurements and the homogeneous oblate polarization model, we conclude that cloud banding is the dominant source of the measured polarization in Luhman 16A. For Luhman 16B, we can neither rule out nor confirm the presence of banded clouds; the measurements could be explained as due to oblateness alone with clouds thinner than in the $f_{\text{sed}} = 1$ models, or by a combination of oblateness and banded clouds. Although we model the bands and fit the polarization as constant over longitude and time, there is the possibility of variability with amplitudes below our sensitivity limits. Indeed, the bands of Jupiter display variability, yet the bulk features remain relatively constant over time (Ge et al. 2019). Similar variability superimposed on the bands discovered here may account for the previously detected photometric variability seen in Luhman 16A.

### 4.3. Spin Axis Orientations from the Angle of Linear Polarization

The angle of polarization for Luhman 16A and B differs by $\sim 105^\circ$ (or equivalently, $\sim 75^\circ$, because of a $180^\circ$ degeneracy in the definition of the angle of polarization). For each object, the angle of polarization could be aligned either parallel or perpendicular to the position angle of the projected on-sky spin axis, as discussed in Section 4.2. The misalignment of the angle of polarization immediately suggests a misalignment of the two components' spin axes of either $\sim 15^\circ$ or $\sim 75^\circ$, depending on the source of polarization of Luhman 16B. Neither component has an angle of polarization aligned parallel or perpendicular to the longitude of the ascending node of the
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Schematic images displaying a selection of the different cloud banding and oblateness scenarios that were explored. The modeled degree of linear polarization *p* and angle of linear polarization *ψ* (and blue arrows) are shown above each model. The degrees of polarization calculated for the oblateness models (Figure 10) are displayed as upper limits, as thinner clouds could reduce the polarization. The measured values for A and B are displayed for reference, and a banded cloud model that matches the polarization of Luhman 16A has been highlighted. The angles of polarization for the measurements are given in the sky frame measured east of north, but the models are given relative to the on-sky projected rotation axis.</figcaption></figure>

<table><caption>**Table 8**
Degree and Angle of Polarization from the Banded Cloud Models Considered Here</caption><thead><tr><td>Inclination</td><td>Band Latitudes</td><td>Degree of Linear Polarization</td><td>Ψ Relative to Projected Spin Axis</td></tr></thead><tbody><tr><td>20°</td><td>[-10, 10]</td><td>0.016%</td><td>Perpendicular</td></tr><tr><td>20°</td><td>[-50, -30], [30, 50]</td><td>0.006%</td><td>Perpendicular</td></tr><tr><td>20°</td><td>[-51, -19], [19, 51]</td><td><b>0.030%</b></td><td>Perpendicular</td></tr><tr><td>20°</td><td>[-50, -30], [-10, 10], [30, 50]</td><td>0.012%</td><td>Perpendicular</td></tr><tr><td>56°</td><td>[-50, -30], [30, 50]</td><td>0.015%</td><td>Parallel</td></tr></tbody></table>

**Note.** The bold value highlights the measured degree of linear polarization of 0.03%.

binary orbit, $\Omega = 130°.12 \pm 0°.12$ (Garcia et al. 2017). Luhman 16A’s polarization is at best ~20° misaligned, and Luhman 16B is ~40°–60°, indicating that neither object has a spin axis aligned with the orbital plane, regardless of the oblateness or banded interpretation. The misalignment may imply that the two objects did not form alone together within the same disk, but instead may have instead experienced a more raucous dynamical history such as having originally been formed in a triple system (Reipurth & Mikkola 2015). To the best of the authors’ knowledge, this is the first measurement of the projected spin axis for any brown dwarf (albeit with a 90° degeneracy). Given the low-S/N detection of Luhman 16B, we consider any interpretation of its position angle to be tentative.

## 5. Total Intensity Variability Interpretation

We were unable to measure absolute photometric variability, but the flux ratio of the two objects exhibits quasi-periodic variations, with a period of 1.64 hr. We attribute this signal to Luhman 16B alone (*H*-band variability has never been detected in A; e.g., Biller et al. 2013; Buenzli et al. 2015b) and suggest

that it may represent a similar three-band scenario to that inferred for other L/T transition objects, where slight differences in the rotational periods of each band results in photometric variations that beat over time (Apai et al. 2017). In this scenario our data would represent a time when the three bands were nearly completely out of phase, but each band still rotates with a ~5 hr period. However, unambiguously distinguishing between a model composed of only cloud patches and that which includes cloud bands would require a baseline of >2 rotation periods.

## 6. Conclusions

Here we have presented new *H*-band linear polarization measurements of the Luhman 16 binary brown dwarf system obtained with VLT/NaCo. The measurements of Luhman 16B could be explained by oblateness or cloud banding, but the polarization of Luhman 16A can only be explained by bands of clouds, similar to those seen throughout the solar system. Previous photometric and spectroscopic variability studies of Luhman 16A have either suggested patchy clouds or have been
unable to constrain its cloud morphology due to nondetections of variability. In contrast, our polarimetric measurements detect bands that could not have been found using these techniques, as both methods rely on cloud morphologies or variations that rotate in and out of view. Although our data are inconclusive about the presence of bands in Luhman 16B, bands cannot be ruled out by our data or previous studies; the Doppler-imaging technique used with Luhman 16B was unable to recover longitudinal bands artificially injected into the data (Crossfield et al. 2014) and cloud spot mapping using *HST* data did not consider banded structures (Karalidi et al. 2016).

The polarization measurements of Luhman 16A represent the state of the cloud morphology at a fixed point in time. The absolute variability of previous polarization measurements of brown dwarfs suggest that cloud morphologies may vary significantly over time. Whether or not the polarization of Luhman 16A varies in time is yet to be determined. We have obtained follow-up observations of Luhman 16 over four nights in 2019 April in order to search for both short-term and long-term variability in the polarization. These data have not yet been analyzed and will be published in a follow-up study.

Cloud bands in brown dwarfs have only been previously inferred for three L/T transition brown dwarfs via photometric monitoring (Apai et al. 2017). The discovery of cloud bands in Luhman 16A provides a critical independent confirmation of banded clouds near the L/T transition and suggests that their presence is common, if not ubiquitous. Whether or not these bands persist outside of the L/T transition remains an open question that should be pursued with further studies. We have demonstrated that polarimetry provides a promising avenue to answer this question; variability occurrence rates and amplitudes decline outside of the L/T transition, and therefore many targets are unsuitable to cloud mapping via variability monitoring. High-accuracy targeted polarimetric studies sensitive to cloud bands similar to those found in this study may therefore be critical to furthering our understanding of the cloud dynamics in substellar objects. The characterization presented herein relies heavily on previous characterization of Luhman 16. Similar polarimetric studies should be pursued for other brown dwarfs with comparable levels of characterization (e.g., brown dwarfs with well-constrained masses and/or inclinations).

We note that the exquisite accuracy of our polarimetric measurements is due not only to the brightness of Luhman 16, but also to the specific instrumental setup that allowed us to

easily distinguish between instrumental polarization and astrophysical signal. Unfortunately, NaCo was decommissioned in 2019, preempting further studies with the same instrument. Other NIR polarimeters, such as SPIRou (Artigau et al. 2014), MMT-Pol (Packham & Jones 2008), and/or WIRC+Pol (Tinyanont et al. 2018) (among others) may provide other opportunities for new brown dwarf polarimetric studies. The discovery of cloud bands in Luhman 16A via polarimetry sets the stage for a wide range of new 3D atmospheric dynamics studies, not only for brown dwarfs like Luhman 16, but also for giant exoplanets.

We thank D. Saumon for updated brown dwarf evolutionary models. Support for this work was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51378.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. The research of F.S. and J.d.B. leading to these results has received funding from the European Research Council under ERC Starting Grant agreement 678194 (FALCONER). Based on observations collected at the European Southern Observatory under ESO programmes 0101.C-0561(A) and 0101.C-0561(B). (Part of) The research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. This work is supported by the National Science Foundation under grant No. AAG-1816341.

*Facility:* VLT(NaCo).

*Software:* *astropy* (Astropy Collaboration et al. 2013), *photutils* (Bradley et al. 2017), *emcee* (Foreman-Mackey et al. 2013), *NumPy*, *SciPy*, *Matplotlib*.

# Appendix
## Supplementary Figures and Tables

Figures A1 and A2 show the Luhman 16A and B *q* and *u* measurements as a function of time. Figure A3 shows the measurements of Luhman 16A and B *p* and *θ* as a function of parallactic angle, while Figure A4 displays marginalized and joint probability distributions for the 18 model parameters. Tables A1 and A2 list the data presented in Figures 1 and 3; and Figures 4, 8, A1, and A2, respectively.
<figure><img src="image_16.png" /><figcaption>**Figure A1.** Luhman 16A $q_{\text{NaCo}}$ (bottom right) and $u_{\text{NaCo}}$ (top left) measurements as a function of time. The bottom-left plot displays the measurements in the $q-u$ plane. For each data point, both the original error bars (black) and the error bars increased by $\sigma_{\text{SD}}$ (for U) and $\sigma_{\text{DD}}$ (for Q) are shown (gray). The multicolored lines are a random selection of accepted models from the MCMC fitting with the cyan lines showing only the IP components. In the $q-u$ plot, the instrumental polarization is shown as a 2D histogram and the nonperfect efficiencies and crosstalk terms cause the MCMC model tracks to deviate from a perfect model, manifesting as ovals rather than a perfect circle.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure A2. Same as Figure A1, but for Luhman 16B.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure A3.** The $p_{\text{NaCo}}$ (left) and $\theta_{\text{NaCo}}$ (right) measurements as a function of parallactic angle of Luhman 16A (top row) and B (bottom) from Figure 4, after subtracting off an estimate of the instrumental polarization. The linear trend in Luhman 16A's $\theta_{\text{NaCo}}$ measurements (overlaid with an orange offset straight line with unity slope) indicates the presence of an astrophysical signal that rotates with the sky relative to the fixed instrument frame. The black error bars represent the photometric errors.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure A4. The marginalized and joint probability distributions for the 18 model parameters.</figcaption></figure>
<table><caption>Table A1
The Flux Ratio Measurements of Luhman 16 AB as a Function of Time with an Aperture Size of 17, Corresponding to the Data Presented in Figures 1 and 3</caption><thead><tr><th>Time from Start (hr)</th><th>F<sub>A</sub>/F<sub>B</sub></th><th>Time from Start (hr)</th><th>F<sub>A</sub>/F<sub>B</sub></th><th>Time from Start (hr)</th><th>F<sub>A</sub>/F<sub>B</sub></th><th>Time from Start (hr)</th><th>F<sub>A</sub>/F<sub>B</sub></th></tr></thead><tbody><tr><td>0.00</td><td>0.95120 ± 0.00005</td><td>2.00</td><td>0.96234 ± 0.00003</td><td>4.94</td><td>0.95335 ± 0.00002</td><td>6.91</td><td>0.97187 ± 0.00002</td></tr><tr><td>0.02</td><td>0.95765 ± 0.00003</td><td>2.02</td><td>0.96288 ± 0.00003</td><td>4.96</td><td>0.95367 ± 0.00002</td><td>6.93</td><td>0.97321 ± 0.00002</td></tr><tr><td>0.04</td><td>0.95930 ± 0.00003</td><td>2.04</td><td>0.96604 ± 0.00003</td><td>4.98</td><td>0.95326 ± 0.00002</td><td>6.96</td><td>0.97260 ± 0.00002</td></tr><tr><td>0.06</td><td>0.95803 ± 0.00003</td><td>2.06</td><td>0.96748 ± 0.00003</td><td>5.00</td><td>0.95570 ± 0.00002</td><td>6.98</td><td>0.97443 ± 0.00002</td></tr><tr><td>0.09</td><td>0.95868 ± 0.00002</td><td>2.08</td><td>0.96635 ± 0.00003</td><td>5.03</td><td>0.95816 ± 0.00002</td><td>7.00</td><td>0.97523 ± 0.00002</td></tr><tr><td>0.11</td><td>0.96131 ± 0.00002</td><td>2.10</td><td>0.96764 ± 0.00003</td><td>5.05</td><td>0.95868 ± 0.00002</td><td>7.02</td><td>0.97448 ± 0.00002</td></tr><tr><td>0.13</td><td>0.96430 ± 0.00002</td><td>2.20</td><td>0.96757 ± 0.00003</td><td>5.07</td><td>0.96090 ± 0.00002</td><td>7.04</td><td>0.97303 ± 0.00002</td></tr><tr><td>0.15</td><td>0.96218 ± 0.00002</td><td>2.22</td><td>0.96569 ± 0.00003</td><td>5.09</td><td>0.96304 ± 0.00002</td><td>7.06</td><td>0.97042 ± 0.00002</td></tr><tr><td>0.20</td><td>0.96316 ± 0.00002</td><td>2.24</td><td>0.96614 ± 0.00003</td><td>5.11</td><td>0.95574 ± 0.00002</td><td>7.09</td><td>0.97583 ± 0.00002</td></tr><tr><td>0.22</td><td>0.96497 ± 0.00002</td><td>2.26</td><td>0.96961 ± 0.00003</td><td>5.13</td><td>0.95754 ± 0.00002</td><td>7.11</td><td>0.97639 ± 0.00002</td></tr><tr><td>0.24</td><td>0.96421 ± 0.00002</td><td>1.78</td><td>0.95961 ± 0.00002</td><td>4.11</td><td>0.95000 ± 0.00003</td><td>5.80</td><td>0.97172 ± 0.00002</td></tr><tr><td>0.26</td><td>0.96365 ± 0.00002</td><td>1.80</td><td>0.96045 ± 0.00003</td><td>4.13</td><td>0.94655 ± 0.00003</td><td>5.83</td><td>0.97291 ± 0.00002</td></tr><tr><td>0.28</td><td>0.96569 ± 0.00002</td><td>1.82</td><td>0.96174 ± 0.00003</td><td>4.15</td><td>0.94960 ± 0.00003</td><td>5.85</td><td>0.97107 ± 0.00002</td></tr><tr><td>0.30</td><td>0.96570 ± 0.00002</td><td>1.84</td><td>0.96080 ± 0.00003</td><td>4.17</td><td>0.95096 ± 0.00003</td><td>5.87</td><td>0.97273 ± 0.00002</td></tr><tr><td>0.32</td><td>0.96620 ± 0.00002</td><td>1.87</td><td>0.96401 ± 0.00003</td><td>4.19</td><td>0.95097 ± 0.00003</td><td>5.89</td><td>0.97510 ± 0.00002</td></tr><tr><td>0.35</td><td>0.96685 ± 0.00003</td><td>1.89</td><td>0.96557 ± 0.00003</td><td>4.21</td><td>0.94981 ± 0.00003</td><td>5.91</td><td>0.97492 ± 0.00002</td></tr><tr><td>0.40</td><td>0.96788 ± 0.00002</td><td>1.91</td><td>0.96505 ± 0.00003</td><td>4.24</td><td>0.94415 ± 0.00003</td><td>5.93</td><td>0.97614 ± 0.00002</td></tr><tr><td>0.42</td><td>0.96860 ± 0.00002</td><td>1.93</td><td>0.96629 ± 0.00003</td><td>4.26</td><td>0.94963 ± 0.00003</td><td>5.95</td><td>0.97493 ± 0.00002</td></tr><tr><td>0.44</td><td>0.96895 ± 0.00002</td><td>1.95</td><td>0.96484 ± 0.00003</td><td>4.28</td><td>0.94082 ± 0.00003</td><td>5.98</td><td>0.97256 ± 0.00002</td></tr><tr><td>0.46</td><td>0.96973 ± 0.00002</td><td>1.98</td><td>0.96165 ± 0.00003</td><td>4.30</td><td>0.94006 ± 0.00003</td><td>6.00</td><td>0.97267 ± 0.00002</td></tr><tr><td>0.48</td><td>0.97174 ± 0.00002</td><td>2.00</td><td>0.96234 ± 0.00003</td><td>4.33</td><td>0.94781 ± 0.00003</td><td>6.02</td><td>0.97214 ± 0.00002</td></tr><tr><td>0.50</td><td>0.97323 ± 0.00002</td><td>2.02</td><td>0.96288 ± 0.00003</td><td>4.35</td><td>0.94567 ± 0.00003</td><td>6.04</td><td>0.97165 ± 0.00002</td></tr><tr><td>0.52</td><td>0.97374 ± 0.00002</td><td>2.04</td><td>0.96604 ± 0.00003</td><td>4.37</td><td>0.94955 ± 0.00003</td><td>6.06</td><td>0.97342 ± 0.00002</td></tr><tr><td>0.55</td><td>0.97476 ± 0.00002</td><td>2.06</td><td>0.96748 ± 0.00003</td><td>4.39</td><td>0.95908 ± 0.00003</td><td>6.08</td><td>0.97218 ± 0.00002</td></tr><tr><td>0.57</td><td>0.97503 ± 0.00002</td><td>2.08</td><td>0.96635 ± 0.00003</td><td>4.41</td><td>0.95653 ± 0.00003</td><td>6.11</td><td>0.97382 ± 0.00002</td></tr><tr><td>0.59</td><td>0.97513 ± 0.00002</td><td>2.10</td><td>0.96764 ± 0.00003</td><td>4.43</td><td>0.95084 ± 0.00003</td><td>6.13</td><td>0.97398 ± 0.00002</td></tr><tr><td>0.61</td><td>0.97583 ± 0.00002</td><td>2.20</td><td>0.96757 ± 0.00003</td><td>4.45</td><td>0.95144 ± 0.00003</td><td>6.15</td><td>0.96999 ± 0.00002</td></tr><tr><td>0.63</td><td>0.97677 ± 0.00002</td><td>2.22</td><td>0.96569 ± 0.00003</td><td>4.47</td><td>0.96476 ± 0.00003</td><td>6.17</td><td>0.96719 ± 0.00002</td></tr><tr><td>0.66</td><td>0.97785 ± 0.00002</td><td>2.24</td><td>0.96614 ± 0.00003</td><td>4.50</td><td>0.96600 ± 0.00003</td><td>6.19</td><td>0.96839 ± 0.00002</td></tr><tr><td>0.68</td><td>0.97897 ± 0.00003</td><td>2.26</td><td>0.96961 ± 0.00003</td><td>4.52</td><td>0.96128 ± 0.00003</td><td>6.21</td><td>0.97053 ± 0.00002</td></tr><tr><td>0.70</td><td>0.98005 ± 0.00003</td><td>2.29</td><td>0.97036 ± 0.00003</td><td>4.54</td><td>0.95981 ± 0.00003</td><td>6.24</td><td>0.97219 ± 0.00002</td></tr><tr><td>0.72</td><td>0.98116 ± 0.00003</td><td>2.31</td><td>0.97242 ± 0.00003</td><td>4.56</td><td>0.95559 ± 0.00003</td><td>6.26</td><td>0.96914 ± 0.00002</td></tr><tr><td>0.74</td><td>0.98045 ± 0.00002</td><td>2.33</td><td>0.97202 ± 0.00003</td><td>4.77</td><td>0.94654 ± 0.00002</td><td>6.28</td><td>0.97049 ± 0.00002</td></tr><tr><td>0.76</td><td>0.98241 ± 0.00003</td><td>2.35</td><td>0.97660 ± 0.00003</td><td>4.79</td><td>0.94873 ± 0.00003</td><td>6.30</td><td>0.97203 ± 0.00002</td></tr><tr><td>0.78</td><td>0.98339 ± 0.00003</td><td>2.37</td><td>0.97159 ± 0.00003</td><td>4.81</td><td>0.94954 ± 0.00003</td><td>6.32</td><td>0.96940 ± 0.00002</td></tr><tr><td>0.81</td><td>0.98004 ± 0.00002</td><td>2.40</td><td>0.97228 ± 0.00003</td><td>4.83</td><td>0.94895 ± 0.00003</td><td>6.34</td><td>0.96627 ± 0.00002</td></tr><tr><td>0.83</td><td>0.98147 ± 0.00002</td><td>2.42</td><td>0.97286 ± 0.00003</td><td>4.85</td><td>0.95406 ± 0.00002</td><td>6.37</td><td>0.96543 ± 0.00002</td></tr><tr><td>0.85</td><td>0.97865 ± 0.00003</td><td>2.44</td><td>0.97542 ± 0.00003</td><td>4.87</td><td>0.95410 ± 0.00002</td><td>6.39</td><td>0.96511 ± 0.00002</td></tr><tr><td>0.87</td><td>0.98107 ± 0.00003</td><td>2.46</td><td>0.97557 ± 0.00002</td><td>4.89</td><td>0.95655 ± 0.00002</td><td>6.41</td><td>0.96831 ± 0.00002</td></tr><tr><td>0.89</td><td>0.98294 ± 0.00003</td><td>2.48</td><td>0.97510 ± 0.00002</td><td>4.92</td><td>0.95621 ± 0.00002</td><td>6.43</td><td>0.96984 ± 0.00002</td></tr><tr><td>0.92</td><td>0.97903 ± 0.00002</td><td>2.50</td><td>0.97604 ± 0.00002</td><td>4.94</td><td>0.95335 ± 0.00002</td><td>6.45</td><td>0.97152 ± 0.00002</td></tr><tr><td>0.94</td><td>0.97979 ± 0.00002</td><td>2.52</td><td>0.97788 ± 0.00003</td><td>4.96</td><td>0.95367 ± 0.00002</td><td>6.47</td><td>0.97007 ± 0.00002</td></tr><tr><td>0.96</td><td>0.98140 ± 0.00003</td><td>2.55</td><td>0.97729 ± 0.00002</td><td>4.98</td><td>0.95326 ± 0.00002</td><td>6.52</td><td>0.96569 ± 0.00002</td></tr><tr><td>0.98</td><td>0.97940 ± 0.00003</td><td>2.57</td><td>0.97790 ± 0.00002</td><td>5.00</td><td>0.95570 ± 0.00002</td><td>6.55</td><td>0.96637 ± 0.00002</td></tr><tr><td>1.00</td><td>0.98288 ± 0.00002</td><td>2.59</td><td>0.97703 ± 0.00003</td><td>5.03</td><td>0.95816 ± 0.00002</td><td>6.57</td><td>0.96439 ± 0.00002</td></tr><tr><td>1.02</td><td>0.97978 ± 0.00002</td><td>2.61</td><td>0.97642 ± 0.00002</td><td>5.05</td><td>0.95868 ± 0.00002</td><td>6.59</td><td>0.96562 ± 0.00002</td></tr><tr><td>1.04</td><td>0.97776 ± 0.00002</td><td>2.63</td><td>0.97901 ± 0.00002</td><td>5.07</td><td>0.96090 ± 0.00002</td><td>6.61</td><td>0.96711 ± 0.00002</td></tr><tr><td>1.07</td><td>0.97962 ± 0.00002</td><td>2.65</td><td>0.97800 ± 0.00003</td><td>5.09</td><td>0.96304 ± 0.00002</td><td>6.63</td><td>0.96677 ± 0.00002</td></tr><tr><td>1.09</td><td>0.97725 ± 0.00002</td><td>2.68</td><td>0.97934 ± 0.00002</td><td>5.11</td><td>0.95574 ± 0.00002</td><td>6.65</td><td>0.96693 ± 0.00002</td></tr><tr><td>1.11</td><td>0.97761 ± 0.00002</td><td>2.70</td><td>0.97933 ± 0.00002</td><td>5.13</td><td>0.95754 ± 0.00002</td><td>6.67</td><td>0.96625 ± 0.00002</td></tr><tr><td>1.13</td><td>0.97776 ± 0.00002</td><td>2.72</td><td>0.97794 ± 0.00002</td><td>5.15</td><td>0.95811 ± 0.00002</td><td>6.70</td><td>0.96538 ± 0.00002</td></tr><tr><td>1.15</td><td>0.97534 ± 0.00002</td><td>2.74</td><td>0.97695 ± 0.00003</td><td>5.17</td><td>0.95797 ± 0.00002</td><td>6.72</td><td>0.96625 ± 0.00003</td></tr><tr><td>1.18</td><td>0.97673 ± 0.00002</td><td>2.76</td><td>0.97866 ± 0.00003</td><td>5.20</td><td>0.95922 ± 0.00002</td><td>6.74</td><td>0.96764 ± 0.00002</td></tr><tr><td>1.20</td><td>0.97556 ± 0.00002</td><td>2.78</td><td>0.98052 ± 0.00002</td><td>5.22</td><td>0.96042 ± 0.00002</td><td>6.76</td><td>0.96849 ± 0.00002</td></tr><tr><td>1.22</td><td>0.97499 ± 0.00002</td><td>2.81</td><td>0.98002 ± 0.00002</td><td>5.24</td><td>0.96067 ± 0.00002</td><td>6.78</td><td>0.96972 ± 0.00002</td></tr><tr><td>1.24</td><td>0.97484 ± 0.00002</td><td>2.83</td><td>0.97816 ± 0.00002</td><td>5.26</td><td>0.96152 ± 0.00002</td><td>6.81</td><td>0.97035 ± 0.00002</td></tr><tr><td>1.26</td><td>0.96978 ± 0.00002</td><td>2.85</td><td>0.97983 ± 0.00002</td><td>5.29</td><td>0.96194 ± 0.00002</td><td>6.83</td><td>0.97026 ± 0.00002</td></tr><tr><td>1.28</td><td>0.96914 ± 0.00002</td><td>2.87</td><td>0.98074 ± 0.00002</td><td>5.31</td><td>0.96119 ± 0.00002</td><td>6.85</td><td>0.97206 ± 0.00002</td></tr><tr><td>1.30</td><td>0.96786 ± 0.00002</td><td>2.89</td><td>0.97812 ± 0.00002</td><td>5.33</td><td>0.96355 ± 0.00002</td><td>6.87</td><td>0.96975 ± 0.00002</td></tr><tr><td>1.33</td><td>0.96674 ± 0.00002</td><td>2.91</td><td>0.97761 ± 0.00002</td><td>5.35</td><td>0.96348 ± 0.00002</td><td>6.89</td><td>0.97095 ± 0.00002</td></tr><tr><td>1.35</td><td>0.96743 ± 0.00002</td><td>2.93</td><td>0.97794 ± 0.00002</td><td>5.37</td><td>0.96784 ± 0.00002</td><td>6.91</td><td>0.97187 ± 0.00002</td></tr><tr><td>1.37</td><td>0.96590 ± 0.00002</td><td>2.96</td><td>0.97795 ± 0.00002</td><td>5.39</td><td>0.96819 ± 0.00002</td><td>6.93</td><td>0.97321 ± 0.00002</td></tr><tr><td>1.39</td><td>0.96557 ± 0.00002</td><td>2.98</td><td>0.97619 ± 0.00002</td><td>5.41</td><td>0.96767 ± 0.00002</td><td>6.96</td><td>0.97260 ± 0.00002</td></tr></tbody></table>
<table><caption>**Table A1**
(Continued)</caption><thead><tr><th>Time from Start (hr)</th><th>F<sub>A</sub>/F<sub>B</sub></th><th>Time from Start (hr)</th><th>F<sub>A</sub>/F<sub>B</sub></th><th>Time from Start (hr)</th><th>F<sub>A</sub>/F<sub>B</sub></th><th>Time from Start (hr)</th><th>F<sub>A</sub>/F<sub>B</sub></th></tr></thead><tbody><tr><td>1.41</td><td>0.96432 ± 0.00002</td><td>3.00</td><td>0.97772 ± 0.00002</td><td>5.43</td><td>0.97012 ± 0.00002</td><td>6.98</td><td>0.97443 ± 0.00002</td></tr><tr><td>1.44</td><td>0.96030 ± 0.00002</td><td>3.02</td><td>0.97641 ± 0.00002</td><td>5.46</td><td>0.96944 ± 0.00002</td><td>7.00</td><td>0.97523 ± 0.00002</td></tr><tr><td>1.46</td><td>0.96349 ± 0.00002</td><td>3.04</td><td>0.97666 ± 0.00002</td><td>5.48</td><td>0.96883 ± 0.00002</td><td>7.02</td><td>0.97448 ± 0.00002</td></tr><tr><td>1.48</td><td>0.95945 ± 0.00002</td><td>3.07</td><td>0.97286 ± 0.00002</td><td>5.50</td><td>0.96943 ± 0.00002</td><td>7.04</td><td>0.97303 ± 0.00002</td></tr><tr><td>1.50</td><td>0.96110 ± 0.00003</td><td>3.09</td><td>0.97311 ± 0.00002</td><td>5.52</td><td>0.97027 ± 0.00002</td><td>7.06</td><td>0.97042 ± 0.00002</td></tr><tr><td>1.52</td><td>0.96428 ± 0.00002</td><td>3.11</td><td>0.97151 ± 0.00002</td><td>5.55</td><td>0.97161 ± 0.00002</td><td>7.09</td><td>0.97583 ± 0.00002</td></tr><tr><td>1.54</td><td>0.96320 ± 0.00003</td><td>3.13</td><td>0.97068 ± 0.00002</td><td>5.57</td><td>0.97221 ± 0.00002</td><td>7.11</td><td>0.97639 ± 0.00002</td></tr></tbody></table>

Note. The error bars represent the 1σ photometric errors, not necessarily the ultimate photometric stability of the relative measurements.

<table><caption>**Table A2**
The Polarization Measurements as a Function of Time with an Aperture Size of 17, Corresponding to the Data Presented in Figures 4, 8, 13, and 14</caption><thead><tr><th>Time from Start (hr)</th><th>q<sub>A</sub> (%)</th><th>u<sub>A</sub> (%)</th><th>q<sub>B</sub> (%)</th><th>u<sub>B</sub> (%)</th></tr></thead><tbody><tr><td>0.07</td><td>-1.927 ± 0.007</td><td>-0.335 ± 0.008</td><td>-1.946 ± 0.007</td><td>-0.385 ± 0.008</td></tr><tr><td>0.27</td><td>-1.914 ± 0.005</td><td>-0.343 ± 0.008</td><td>-1.927 ± 0.005</td><td>-0.361 ± 0.008</td></tr><tr><td>0.47</td><td>-1.927 ± 0.005</td><td>-0.327 ± 0.008</td><td>-1.917 ± 0.005</td><td>-0.327 ± 0.008</td></tr><tr><td>0.64</td><td>-1.949 ± 0.005</td><td>-0.279 ± 0.008</td><td>-1.912 ± 0.005</td><td>-0.310 ± 0.008</td></tr><tr><td>0.82</td><td>-1.963 ± 0.005</td><td>-0.276 ± 0.008</td><td>-1.925 ± 0.005</td><td>-0.317 ± 0.008</td></tr><tr><td>0.99</td><td>-1.966 ± 0.005</td><td>-0.305 ± 0.008</td><td>-1.932 ± 0.005</td><td>-0.326 ± 0.008</td></tr><tr><td>1.16</td><td>-1.944 ± 0.005</td><td>-0.324 ± 0.008</td><td>-1.918 ± 0.005</td><td>-0.311 ± 0.008</td></tr><tr><td>1.34</td><td>-1.944 ± 0.005</td><td>-0.327 ± 0.008</td><td>-1.892 ± 0.005</td><td>-0.372 ± 0.007</td></tr><tr><td>1.51</td><td>-1.954 ± 0.005</td><td>-0.332 ± 0.008</td><td>-1.887 ± 0.005</td><td>-0.326 ± 0.008</td></tr><tr><td>1.68</td><td>-1.942 ± 0.005</td><td>-0.349 ± 0.008</td><td>-1.907 ± 0.005</td><td>-0.381 ± 0.008</td></tr><tr><td>1.86</td><td>-1.954 ± 0.005</td><td>-0.303 ± 0.010</td><td>-1.917 ± 0.005</td><td>-0.309 ± 0.009</td></tr><tr><td>2.03</td><td>-1.893 ± 0.005</td><td>-0.288 ± 0.009</td><td>-1.910 ± 0.005</td><td>-0.255 ± 0.009</td></tr><tr><td>2.28</td><td>-1.963 ± 0.005</td><td>-0.363 ± 0.009</td><td>-1.910 ± 0.005</td><td>-0.358 ± 0.009</td></tr><tr><td>2.45</td><td>-1.953 ± 0.005</td><td>-0.350 ± 0.008</td><td>-1.927 ± 0.005</td><td>-0.338 ± 0.008</td></tr><tr><td>2.62</td><td>-1.938 ± 0.005</td><td>-0.373 ± 0.008</td><td>-1.927 ± 0.005</td><td>-0.357 ± 0.008</td></tr><tr><td>2.79</td><td>-1.932 ± 0.005</td><td>-0.379 ± 0.008</td><td>-1.937 ± 0.005</td><td>-0.360 ± 0.008</td></tr><tr><td>2.97</td><td>-1.918 ± 0.005</td><td>-0.353 ± 0.008</td><td>-1.894 ± 0.005</td><td>-0.317 ± 0.008</td></tr><tr><td>3.05</td><td>-1.922 ± 0.005</td><td>-0.340 ± 0.008</td><td>-1.888 ± 0.005</td><td>-0.279 ± 0.008</td></tr><tr><td>4.05</td><td>-1.918 ± 0.006</td><td>-0.343 ± 0.011</td><td>-1.910 ± 0.006</td><td>-0.305 ± 0.010</td></tr><tr><td>4.23</td><td>-1.912 ± 0.006</td><td>-0.336 ± 0.011</td><td>-1.897 ± 0.006</td><td>-0.328 ± 0.010</td></tr><tr><td>4.40</td><td>-1.867 ± 0.007</td><td>-0.318 ± 0.011</td><td>-1.926 ± 0.007</td><td>-0.356 ± 0.011</td></tr><tr><td>4.66</td><td>-1.916 ± 0.006</td><td>-0.310 ± 0.010</td><td>-1.916 ± 0.006</td><td>-0.264 ± 0.010</td></tr><tr><td>4.93</td><td>-1.904 ± 0.005</td><td>-0.315 ± 0.008</td><td>-1.925 ± 0.005</td><td>-0.317 ± 0.007</td></tr><tr><td>5.10</td><td>-1.904 ± 0.005</td><td>-0.321 ± 0.008</td><td>-1.942 ± 0.005</td><td>-0.296 ± 0.007</td></tr><tr><td>5.27</td><td>-1.914 ± 0.005</td><td>-0.350 ± 0.008</td><td>-1.910 ± 0.005</td><td>-0.311 ± 0.007</td></tr><tr><td>5.45</td><td>-1.885 ± 0.005</td><td>-0.304 ± 0.008</td><td>-1.924 ± 0.005</td><td>-0.297 ± 0.008</td></tr><tr><td>5.62</td><td>-1.883 ± 0.005</td><td>-0.322 ± 0.008</td><td>-1.928 ± 0.005</td><td>-0.335 ± 0.008</td></tr><tr><td>5.79</td><td>-1.867 ± 0.005</td><td>-0.288 ± 0.008</td><td>-1.923 ± 0.005</td><td>-0.341 ± 0.008</td></tr><tr><td>5.97</td><td>-1.901 ± 0.005</td><td>-0.307 ± 0.008</td><td>-1.926 ± 0.005</td><td>-0.381 ± 0.008</td></tr><tr><td>6.14</td><td>-1.912 ± 0.005</td><td>-0.276 ± 0.008</td><td>-1.915 ± 0.005</td><td>-0.292 ± 0.008</td></tr><tr><td>6.31</td><td>-1.878 ± 0.005</td><td>-0.321 ± 0.008</td><td>-1.932 ± 0.005</td><td>-0.325 ± 0.008</td></tr><tr><td>6.50</td><td>-1.904 ± 0.005</td><td>-0.307 ± 0.008</td><td>-1.929 ± 0.005</td><td>-0.374 ± 0.008</td></tr><tr><td>6.69</td><td>-1.901 ± 0.005</td><td>-0.291 ± 0.008</td><td>-1.936 ± 0.005</td><td>-0.362 ± 0.008</td></tr><tr><td>6.86</td><td>-1.894 ± 0.005</td><td>-0.292 ± 0.008</td><td>-1.930 ± 0.005</td><td>-0.330 ± 0.008</td></tr><tr><td>7.03</td><td>-1.908 ± 0.005</td><td>-0.319 ± 0.008</td><td>-1.945 ± 0.005</td><td>-0.390 ± 0.008</td></tr></tbody></table>

Note. The errors presented here are the 1σ photometric errors (i.e., the black error bars in Figure 4) and do not include the extra error term that was fit for (i.e., the red error bars in Figure 4).

ORCID iDs

Maxwell A. Millar-Blanchaer https://orcid.org/0000-0001-6205-9233

Julien H. Girard https://orcid.org/0000-0001-8627-0404

Theodora Karalidi https://orcid.org/0000-0001-7356-6652

Mark S. Marley https://orcid.org/0000-0002-5251-2943

Sujan Sengupta https://orcid.org/0000-0002-6176-3816

Dimitri Mawet https://orcid.org/0000-0002-8895-4735

Tiffany Kataria https://orcid.org/0000-0003-3759-9080

Jos de Boer https://orcid.org/0000-0002-1028-7522

Rebecca Jensen-Clem https://orcid.org/0000-0003-0054-2953

Sasha Hinkley https://orcid.org/0000-0001-8074-2562
References

Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872

Apai, D., Karalidi, T., Marley, M. S., et al. 2017, Sci, 357, 683

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Artigau, É 2018, Handbook of Exoplanets (New York: Springer), 94

Artigau, É, Kouach, D., Donati, J.-F., et al. 2014, Proc. SPIE, 9147, 914715

Barnes, J. W., & Fortney, J. J. 2003, ApJ, 588, 545

Biller, B. 2017, AstRv, 13, 1

Biller, B. A., Crossfield, I. J. M., Mancini, L., et al. 2013, ApJL, 778, L10

Bowler, B. P. 2016, PASP, 128, 102001

Bradley, L., Sipocz, B., & Robitaille, T. 2017, astropy/photutils, v. 0.7.2, Zenodo, doi:10.5281/zenodo.1039309

Buenzli, E., Marley, M. S., Apai, D., et al. 2015a, ApJ, 812, 163

Buenzli, E., Saumon, D., Marley, M. S., et al. 2015b, ApJ, 798, 127

Burgasser, A. J., Gillon, M., Faherty, J. K., et al. 2014, ApJ, 785, 48

Burgasser, A. J., Marley, M. S., Ackerman, A. S., et al. 2002, ApJL, 571, L151

Burgasser, A. J., Sheppard, S. S., & Luhman, K. L. 2013, ApJ, 772, 129

Clarke, D. 2010, Stellar Polarimetry (New York: Wiley)

Crossfield, I. J. M. 2014, A&A, 566, A130

Crossfield, I. J. M., Biller, B., Schlieder, J. E., et al. 2014, Natur, 505, 654

Curran, P. A. 2014, arXiv:1411.3816

deBoer, J., Girard, J. H., Mawet, D., et al. 2014, Proc. SPIE, 9147, 914787

de Juan Ovelar, M. 2013, PhD thesis, Leiden Univ.

de Kok, R. J., Stam, D. M., & Karalidi, T. 2011, ApJ, 741, 59

Dutta, G., Ajay Kumar, M. C., Vinay Kumar, P., et al. 2009, JGRD, 114, D18109

Faherty, J. K., Beletsky, Y., Burgasser, A. J., et al. 2014, ApJ, 790, 90

Foreman-Mackey, D. 2016, JOSS, 1, 24

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Garcia, E. V., Ammons, S. M., Salama, M., et al. 2017, ApJ, 846, 97

Ge, H., Zhang, X., Fletcher, L. N., et al. 2019, AJ, 157, 89

Gillon, M., Triaud, A. H. M. J., Jehin, E., et al. 2013, A&A, 555, L5

Ginski, C., Benisty, M., van Holstein, R. G., et al. 2018, A&A, 616, A79

Goldman, B., Pitann, J., Zapatero Osorio, M. R., et al. 2009, A&A, 502, 929

Goldstein, D. 2003, Polarized Light, Revised and Expanded, Optical Engineering (Boca Raton, FL: CRC Press)

Harrington, D. M., Kuhn, J. R., & Hall, S. 2011, PASP, 123, 799

Hashimoto, J., Tamura, M., Kandori, R., et al. 2009, in AIP Conf. Ser. 1158, Exoplanets and Disks: Their Formation and Diversity, ed. T. Usuda, M. Tamura, & M. Ishii (Melville, NY: AIP), 111

Karalidi, T., Apai, D., Marley, M. S., & Buenzli, E. 2016, ApJ, 825, 90

Kniazev, A. Y., Vaisanen, P., Mužić, K., et al. 2013, ApJ, 770, 124

Lazorenko, P. F., & Sahlmann, J. 2018, A&A, 618, A111

Lenzen, R., Hartung, M., Brandner, W., et al. 2003, Proc. SPIE, 4841, 944

Lomb, N. R. 1976, Ap&SS, 39, 447

Luhman, K. L. 2013, ApJL, 767, L1

Mancini, L., Giacobbe, P., Littlefair, S. P., et al. 2015, A&A, 584, A104

Marley, M. S., Saumon, D., & Goldblatt, C. 2010, ApJL, 723, L117

Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Sci, 272, 1919

Marley, M. S., & Sengupta, S. 2011, MNRAS, 417, 2874

Mathewson, D. S., & Ford, V. L. 1970, MmRAS, 74, 139

Ménard, F., Delfosse, X., & Monin, J.-L. 2002, A&A, 396, L35

Miles-Páez, P. A., Zapatero Osorio, M. R., & Pallé, E. 2015, A&A, 580, L12

Miles-Páez, P. A., Zapatero Osorio, M. R., Pallé, E., & Peña Ramírez, K. 2013, A&A, 556, A125

Miles-Páez, P. A., Zapatero Osorio, M. R., Pallé, E., & Peña Ramírez, K. 2017, MNRAS, 466, 3184

Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, ApJ, 756, 172

Osten, R. A., Melis, C., Stelzer, B., et al. 2015, ApJL, 805, L3

Packham, C., & Jones, T. J. 2008, Proc. SPIE, 7014, 70145F

Radigan, J., Jayawardhana, R., Lafrenière, D., et al. 2012, ApJ, 750, 105

Reipurth, B., & Mikkola, S. 2015, AJ, 149, 145

Rousset, G., Lacombe, F., Puget, P., et al. 2003, Proc. SPIE, 4839, 140

Scargle, J. D. 1982, ApJ, 263, 835

Sengupta, S., & Krishan, V. 2001, ApJL, 561, L123

Sengupta, S., & Marley, M. S. 2009, ApJ, 707, 716

Sengupta, S., & Marley, M. S. 2010, ApJL, 722, L142

Serkowski, K. 1973, in IAU Symp. 52, Interstellar Dust and Related Topics, ed. J. M. Greenberg & H. C. van de Hulst (Dordrecht: Reidel), 145

Showman, A. P., Cho, J. Y.-K., & Menou, K. 2010, in Exoplanets, ed. S. Seager (Tucson, AZ: Univ. Arizona Press), 471

Showman, A. P., & Kaspi, Y. 2013, ApJ, 776, 85

Stolker, T., Min, M., Stam, D. M., et al. 2017, A&A, 607, A42

Tata, R., Martín, E. L., Sengupta, S., et al. 2009, A&A, 508, 1423

Tinyanont, S., Millar-Blanchaer, M., Nilsson, R., et al. 2018, Proc. SPIE, 10702, 107023J

West, R. A., Baines, K. H., Friedson, A. J., et al. 2004, in Jupiter: the planet, satellites and magnetosphere, ed. F. Bagenal et al. (Cambridge: Cambridge Univ. Press), 79

Whittet, D. C. B., Martin, P. G., Hough, J. H., et al. 1992, ApJ, 386, 562

Witzel, G., Eckart, A., Buchholz, R. M., et al. 2011, A&A, 525, A130

Zapatero Osorio, M. R., Béjar, V. J. S., Goldman, B., et al. 2011, ApJ, 740, 4

Zapatero Osorio, M. R., Caballero, J. A., & Béjar, V. J. S. 2005, ApJ, 621, 445

Zhang, X., & Showman, A. P. 2014, ApJL, 788, L6