# Individual, Model-independent Masses of the Closest Known Brown Dwarf Binary to the Sun

E. Victor Garcia, S. Mark Ammons, Maissa Salama, Ian Crossfield, Eduardo

Bendek, Jeffrey Chilcote, Vincent Garrel, James R. Graham, Paul Kalas,

Quinn Konopacky, et al.

► To cite this version:

E. Victor Garcia, S. Mark Ammons, Maissa Salama, Ian Crossfield, Eduardo Bendek, et al. Individual, Model-independent Masses of the Closest Known Brown Dwarf Binary to the Sun. The Astrophysical Journal, 2017, 846 (2), 10.3847/1538-4357/aa844f. [DOI: https://doi.org/10.3847/1538-4357/aa844f](https://doi.org/10.3847/1538-4357/aa844f). [arXiv:1707.0849](https://arxiv.org/abs/1707.0849)

HAL Id: hal-01777519

https://hal.science/hal-01777519v1

Submitted on 28 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Individual, Model-independent Masses of the Closest Known Brown Dwarf Binary to the Sun

E. Victor Garcia¹, S. Mark Ammons¹, Maissa Salama², Ian Crossfield³, Eduardo Bendek⁴, Jeffrey Chilcote⁵,
Vincent Garrel⁶, James R. Graham⁷, Paul Kalas⁷, Quinn Konopacky⁸, Jessica R. Lu²,
Bruce Macintosh⁹, Eduardo Marin⁶, Christian Marois¹⁰,¹¹, Eric Nielsen⁹,¹², Benoît Neichel¹³, Don Pham¹,
Robert J. De Rosa⁷, Dominic M. Ryan⁷, Maxwell Service², and Gaetano Sivo⁶

¹Lawrence Livermore National Laboratory, L-210, 7000 East Avenue, Livermore, CA 94550, USA; eugenio.v.garcia@gmail.com

²Institute for Astronomy, University of Hawaii at Manoa, Honolulu, HI 96822, USA

³Department of Physics, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA

⁴NASA Ames Research Center, Moffett Field, CA 94035, USA

⁵Dunlap Institute for Astronomy and Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada

⁶Gemini Telescope, Colina el Pino S/N, La Serena, Chile

⁷Astronomy Department, University of California, Berkeley, CA 94720, USA

⁸Center for Astrophysics and Space Sciences, University of California, San Diego, La Jolla, CA 92093, USA

⁹Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Stanford, CA 94305, USA

¹⁰National Research Council of Canada Herzberg, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada

¹¹Department of Physics and Astronomy, University of Victoria, 3800 Finnerty Road, Victoria, BC V8P 5C2, Canada

¹²SETI Institute, Carl Sagan Center, 189 Bernardo Avenue, Mountain View, CA 94043, USA

¹³Aix Marseille Universit, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France

Received 2016 August 26; revised 2017 August 1; accepted 2017 August 2; published 2017 September 5

## Abstract

At a distance of ~2 pc, our nearest brown dwarf neighbor, Luhman 16 AB, has been extensively studied since its discovery 3 years ago, yet its most fundamental parameter—the masses of the individual dwarfs—has not been constrained with precision. In this work, we present the full astrometric orbit and barycentric motion of Luhman 16 AB and the first precision measurements of the individual component masses. We draw upon archival observations spanning 31 years from the European Southern Observatory (ESO) Schmidt Telescope, the Deep Near-Infrared Survey of the Southern Sky (DENIS), public FORS2 data on the Very Large Telescope (VLT), and new astrometry from the Gemini South Multiconjugate Adaptive Optics System (GeMS). Finally, we include three radial velocity measurements of the two components from VLT/CRIRES, spanning one year. With this new data sampling a full period of the orbit, we use a Markov chain Monte Carlo algorithm to fit a 16-parameter model incorporating mutual orbit and barycentric motion parameters and constrain the individual masses to be $27.9_{-1.0}^{+1.1} M_\odot$ for the T dwarf and $34.2_{-1.1}^{+1.3} M_\odot$ for the L dwarf. Our measurements of Luhman 16 AB's mass ratio and barycentric motion parameters are consistent with previous estimates in the literature utilizing recent astrometry only. The GeMS-derived measurements of the Luhman 16 AB separation in 2014–2015 agree closely with *Hubble Space Telescope* (HST) measurements made during the same epoch, and the derived mutual orbit agrees with those measurements to within the HST uncertainties of 0.3–0.4 mas.

**Key words:** astrometry – binaries: visual – brown dwarfs – stars: fundamental parameters – stars: individual (WISE J104915.57–531906.1)

## 1. Introduction

Individual mass measurements of brown dwarfs are critical to calibrating the mass–luminosity relationship for substellar objects. To date, there is only a single field brown dwarf binary with direct individual mass measurements of its components (SDSS 105213.51+442255.7AB, Dupuy et al. 2015). In the substellar mass regime of ≈0.01–0.08 $M_\odot$ (Dieterich et al. 2014; Auddy et al. 2016), there are only four other systems with model-independent, directly measured masses with uncertainties <10%: the young, eclipsing brown dwarf M6+M6 binary with masses of $0.054 \pm 0.005 M_\odot$ and $0.034 \pm 0.003 M_\odot$ 2MASS J05352184-0546085AB from Stassun et al. (2006); a substellar ($0.063 \pm 0.005 M_\odot$) companion GL 802B with a directly observed mass via barycentric astrometry from Ireland et al. (2008); a substellar companion to HR 7672B with directly observed masses of $0.0655 \pm 0.0029 M_\odot$ and $0.069_{-0.012}^{+0.003} M_\odot$ from combined orbital astrometry and radial velocity measurements from Crepp et al. (2012); and GJ 569 Bab with directly observed masses of $0.073 \pm 0.008 M_\odot$ and $0.053 \pm 0.006 M_\odot$

with combined radial velocity and orbital astrometry from Konopacky et al. (2010). There are also a number of brown dwarfs with mass determinations that depend on a stellar model–derived mass of the primary (host) star (Deleuil et al. 2008; Anderson et al. 2011; Bouchy et al. 2011a, 2011b; Johnson et al. 2011; Siverd et al. 2012; Díaz et al. 2013; Montet et al. 2015, 2016). Given that there are only seven completely model-independent mass determinations of individual brown dwarfs (as compared to more than 30 planets and about 100 or more stars), more measurements are sorely needed.

At a distance of ~2 pc, Luhman 16 AB (WISE J104915.57–531906.1) is the closest known binary brown dwarf system. Given its proximity, Luhman 16 represents a unique opportunity to directly test brown dwarf evolution models. There has been a flurry of literature on Luhman 16 AB since it was first identified as a high proper motion L/T brown dwarf binary by Luhman (2013). The binary is composed of a L7.5 primary and a T0.5 secondary (Burgasser et al. 2013; Kniazev et al. 2013). Both components show Li I (6707 Å) absorption, confirming
Luhman 16 AB as the closest brown dwarf binary discovered (Faherty et al. 2014), and exhibit significant photometric and spectroscopic variability, suggesting patchy clouds (Biller 2013; Gillon et al. 2013; Burgasser et al. 2014; Buenzli et al. 2015; Mancini et al. 2015; Street et al. 2015). Lodieu et al. (2015) obtained VLT X-shooter spectroscopy of Luhman 16 AB, finding evidence for lithium, placing an upper limit on the masses of $\lesssim 0.06 M_\odot$. Crossfield et al. (2014) and Karalidi et al. (2016) find evidence for spots, using Doppler imaging and light-curve inversion techniques, respectively. Sahlmann & Lazorenko (2015) derived a distance of $1.998 \pm 0.0004$ pc and further confirmed the high proper motion of $\approx 2''$ 8 yr$^{-1}$ of Luhman (2013). Sahlmann & Lazorenko (2015) used the FORS2/VLT I-band imaging to obtain barycentric motion, finding a mass ratio of $q \equiv M_B/M_A = 0.78 \pm 0.10$ for the system. Bedin et al. (2017) use astrometry recently obtained from the *Hubble Space Telescope (HST)* Wide Field Camera 3 (WFC3) to constrain the system's barycentric motion and mass ratio, obtaining $q \equiv M_B/M_A = 0.848 \pm 0.023$. Bedin et al. (2017) also obtain constraints on the orbital parameters and total mass of Luhman 16 AB, although these are relatively weak due to the limited coverage of the system's full period.

Despite the characterization of Luhman 16 AB's surfaces, cloud structures, temperature, atmospheric composition, photometric variability, parallax, proper motion, and mass ratio over the last 3 years, we have yet to directly measure Luhman 16 AB's most fundamental parameter, the individual component masses, with precision. Previous authors have noted a degeneracy between barycentric proper motion and the mass ratio of the binary (Sahlmann & Lazorenko 2015), frustrating mass measurements that rely on modern epoch data alone. In this work, we constrain the total mass $M_{\text{tot}}$ by combining high-precision astrometry with archival imaging data from 1984 and 1999 to sample the full orbital period. We measure positions of the individual components of Luhman 16 at the 1984 epoch, place these data on the ICRS system using publicly available GAIA DR1 data of field reference stars (Gaia Collaboration et al. 2016a, 2016b), and include the 1984 positions as explicit constraints in the model fit. We then fit a global, 16-parameter model that includes seven parameters describing the mutual orbit of the binary (total mass $M_{\text{tot}} = M_A + M_B$, semimajor axis log a, eccentricity e, longitude of ascending node Ω, argument of periastron ω, time of periastron passage τ, and inclination cos i), six parameters describing the space motion of the system barycenter (parallax $\bar{\omega}_{\text{rel}}$, proper motion $\mu_\alpha$, $\delta$, position $\alpha_0$, $\delta_0$, and mass ratio q), and three parameters modeling chromatic differential atmospheric refraction effects ($\rho_{\text{fors2}}$, $d_{\text{fors2}}$, $\rho_{\text{gems}}$). By jointly fitting for the mutual orbit and barycentric motion parameters, we account for any correlated uncertainty between these. We present comparisons between our measurements and those of Bedin et al. (2017) in the conclusion.

## 2. Data

The astrometry and relative radial velocity data used in this work span 31 years, roughly the orbital period previously estimated for this binary (Luhman 2013; Mamajek 2013). We obtained astrometry of Luhman 16 AB using the following facilities: Gemini Multiconjugate Adaptive Optics System (Neichel et al. 2014; Rigaut et al. 2014) in 2014–2015; archival data from VLT/FORS2 data from 2013–2015; a photographic plate taken with the 1 m Schmidt telescope at the European Southern Observatory (ESO) La Silla Observatory in

1984; and an image from the Deep Near-Infrared Survey of the Southern Sky (DENIS, Epchtein et al. 1999) in 1999. The elongated 1984 ESO R-band image of Luhman 16 AB, the unresolved I-band DENIS 1999 image, a sample reduced GeMS 2014 image, and a VLT/FORS2 2013 epoch are shown in Figure 1. We also include archival VLT/CRIRES radial velocity measurements from 2013–2014, described at the end of this section.

### 2.1. GeMS Multi-conjugate Adaptive Optics (MCAO) Imaging

Observations of Luhman 16 AB and NGC 1851 were obtained with the Gemini South telescope using the Gemini South Adaptive Optics Imager (GSAOI; McGregor et al. 2004; Carrasco et al. 2012) and the Gemini Multiconjugate Adaptive Optics System (GeMS; Rigaut et al. 2014). Together, the two instruments deliver near-diffraction-limited images from 0.9 to 2.4 µm, with a field of view (FOV) of 85 × 85". This AO system uses five sodium Laser Guide Stars (LGSs) to correct for atmospheric turbulence and up to three Natural Guide Stars (NGSs) brighter than $R = 15.5$ mag to compensate for tip-tilt and plate mode variations<fnref n="14" /> over a 2' FOV of the AO bench unit (CANOPUS; Rigaut et al. 2014). With a resolution typically exceeding 0".1 in H and $K_s$, GeMS is capable of resolving Luhman 16 AB throughout its entire orbit.

We measure the GeMS static distortion map with observations of NGC 1851 obtained as part of the GeMS Science Verification program. Stellar positions from GeMS are compared with known stellar positions measured with the *HST* Advanced Camera for Surveys (ACS; Ford et al. 1998). The distortion measurement uses two epochs of undithered frames of NGC 1851 taken one month apart, on 2012 December 30 and 2013 January 28. Images are flat-fielded with combined twilight flats and masked for bad pixels using IRAF's GAREDUCE reduction routines. Images are not mosaicked into 4196 × 4196 frames.

Images of Luhman 16 AB were obtained with GeMS on the nights of 2014 February 12, April 12, April 13, May 27, December 3, and 2015 March 5 in the J-continuum narrowband filter.<fnref n="15" /> A sample image is shown in the lower left panel of Figure 1. Exposure times were 15 s, and a position angle of 200° was used. AO-corrected Strehl ratios as measured on Luhman 16 AB were 0.5%–5%. The data were reduced following the standard procedures for near-infrared imaging provided by the Gemini/GSAOI package inside IRAF (Tody 1986). Images were flat-fielded with combined twilight flats and masked for bad pixels using IRAF's GAREDUCE reduction routines.

### 2.2. FORS2 Imaging

We analyzed FORS2 (Appenzeller et al. 1998) images of Luhman 16 AB taken between 2013 April 4 and 2015 February 3 that were retrieved from the ESO archive (Boffin et al. 2014; Programs 291.C-5004 and 593.C-0314; PI: Boffin). A sample image is shown in the lower right panel of Figure 1. The images were obtained in the I-band Bessel filter with exposures of 15 s. The images were dithered in a 3 × 3 box pattern with 1"0 steps. We analyzed the dark images, stacked the twilight

<footnote n="14">¹⁴ For details, see “Tip-Tilt and Plate Scales modes” here: http://www.gemini.edu/sciops/instruments/gems/introduction-gems/loops-and-offloads.</footnote>

<footnote n="15">¹⁵ The J-continuum central wavelength is 1.207 µm. See here: https://www.gemini.edu/sciops/instruments/gsaoi/instrument-description/filters.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sample images of Luhman 16 AB used to derive barycenter and orbital astrometry listed in Tables 2 and 1. All images are the same 9" × 9" cutout and rotated so north is up. Luhman 16 B is brighter than Luhman 16 A for the GeMS images (lower left panel) due to a J-band flux reversal across the L-T transition, as noted by Faherty et al. (2014).</figcaption></figure>

flat-field image, and flat-fielded the science images with the ESO Reflex Data Reduction Pipeline (Freudling et al. 2013). Image FWHMs varied between 0.′′51 and 0.′′85, and the observations span 661 days. Each of 27 epochs consists of 16–42 frames taken over 1800 s on average, resulting in 825 exposures. We exclude 107 frames with point-spread function (PSF) FWHM > 0.′′75, resulting in 718 frames. Although more archival FORS2 images are available after 2015 February 3, we choose not to analyze these because the PSFs of the individual components of Luhman 16 AB begin to overlap after this point in time in the seeing-limited FORS2 images, making precision astrometry difficult to recover.

## 2.3. ESO R-band Photometric Plate Image

We make use of a scanned plate image of Luhman 16 AB taken by the ESO 1 m Schmidt Telescope on 1984 March 5, first identified by Mamajek (2013). From Figure 1 (upper left), Luhman 16 AB appears well resolved in 1984, with the brighter component in the southeast. SUPERCOSMOS (Hambly et al. 2001) obtained a scan of ESO red plate R5562 with a plate scale of 0.′′675 pixel⁻¹. The filter used for the photographic plate was the RG630 filter, and the plate used IIIaF emulsion. The exposure time was 120 s.

## 2.4. DENIS I-band Image

Luhman 16 AB was identified by Luhman (2013) in an I-band (0.83 µm) image from DENIS on 1999 February 11. From Figure 1 (upper right), Luhman 16 AB appears unresolved in very good seeing, implying the binary was at small projected separation at that time. The plate scale was 1.′′0 pixel⁻¹ for the I-band filter and $t_{exp} = 9$ s.

## 2.5. VLT/CRIRES Relative Radial Velocities

The Luhman 16 AB system was observed several times by the high-resolution infrared spectrograph CRIRES (Kaeufl et al. 2004) fed by Multi-Applications Curvature Adaptive optics, as part of programs 291.C-5006 and 093.C-0335 (PI Crossfield). The former program observed Luhman 16 for five hours on UT 2013 May 05, and the latter observed for several hours each on UT 2014 May 02 and 20. The spectra span wavelengths from 2.288 to 2.345 µm; the VLT adaptive optics system spatially resolved the two brown dwarfs, so the CRIRES slit was aligned along the binary position angle to observe both brown dwarfs simultaneously.

We use standard ESO tools to reduce these data and extract 1D spectra, following exactly the same approach laid out by Crossfield et al. (2014). Again following that work, we
<table><caption>**Table 1**
Position Measurements for Luhman 16 AB</caption><thead><tr><th>Instrument</th><th>MJD<br>(days)</th><th>α<sub>1</sub><br>(°, mas)</th><th>δ<sub>1</sub><br>(°, mas)</th><th>α<sub>2</sub><br>(°, mas)</th><th>δ<sub>2</sub><br>(°, mas)</th></tr></thead><tbody><tr><td>ESO Schmidt (Red)</td><td>45764.95</td><td>162.3487979 ± 179.51</td><td>-53.3212708 ± 169.87</td><td>162.3482743 ± 233.46</td><td>-53.3208797 ± 226.13</td></tr><tr><td rowspan="14">FORS2/VLT</td><td>56397.00</td><td>162.3111711 ± 1.65</td><td>-53.3182805 ± 1.42</td><td>162.3106864 ± 1.65</td><td>-53.3180021 ± 1.46</td></tr><tr><td>56402.06</td><td>162.3111447 ± 1.52</td><td>-53.3182687 ± 1.36</td><td>162.3106612 ± 1.53</td><td>-53.3179911 ± 1.36</td></tr><tr><td>56408.04</td><td>162.3111090 ± 1.53</td><td>-53.3182555 ± 1.34</td><td>162.3106269 ± 1.52</td><td>-53.3179791 ± 1.34</td></tr><tr><td>56418.03</td><td>162.3110542 ± 1.72</td><td>-53.3182314 ± 1.53</td><td>162.3105750 ± 1.70</td><td>-53.3179576 ± 1.65</td></tr><tr><td>56424.09</td><td>162.3110300 ± 2.33</td><td>-53.3182171 ± 1.99</td><td>162.3105528 ± 2.25</td><td>-53.3179446 ± 1.93</td></tr><tr><td>56432.97</td><td>162.3109810 ± 2.24</td><td>-53.3181949 ± 1.45</td><td>162.3105053 ± 1.69</td><td>-53.3179242 ± 1.48</td></tr><tr><td>56438.03</td><td>162.3109655 ± 1.53</td><td>-53.3181826 ± 1.36</td><td>162.3104909 ± 1.51</td><td>-53.3179129 ± 1.37</td></tr><tr><td>56443.99</td><td>162.3109410 ± 1.89</td><td>-53.3181675 ± 1.45</td><td>162.3104679 ± 1.81</td><td>-53.3178995 ± 1.45</td></tr><tr><td>56448.00</td><td>162.3109277 ± 1.54</td><td>-53.3181581 ± 1.35</td><td>162.3104557 ± 1.56</td><td>-53.3178908 ± 1.34</td></tr><tr><td>56452.97</td><td>162.3109107 ± 1.50</td><td>-53.3181461 ± 1.32</td><td>162.3104397 ± 1.50</td><td>-53.3178798 ± 1.32</td></tr><tr><td>56457.00</td><td>162.3109010 ± 1.58</td><td>-53.3181373 ± 1.44</td><td>162.3104316 ± 1.58</td><td>-53.3178717 ± 1.44</td></tr><tr><td>56459.97</td><td>162.3108918 ± 1.59</td><td>-53.3181301 ± 1.43</td><td>162.3104225 ± 1.55</td><td>-53.3178656 ± 1.38</td></tr><tr><td>56465.96</td><td>162.3108774 ± 1.49</td><td>-53.3181169 ± 1.78</td><td>162.3104097 ± 1.55</td><td>-53.3178537 ± 2.28</td></tr><tr><td>56693.25</td><td>162.3103384 ± 1.51</td><td>-53.3182177 ± 1.35</td><td>162.3099316 ± 1.53</td><td>-53.3180095 ± 1.37</td></tr><tr><td>56702.19</td><td>162.3102721 ± 1.52</td><td>-53.3182173 ± 1.32</td><td>162.3098669 ± 1.51</td><td>-53.3180112 ± 1.32</td></tr><tr><td>56726.14</td><td>162.3100969 ± 1.55</td><td>-53.3182010 ± 1.35</td><td>162.3097001 ± 1.56</td><td>-53.3180015 ± 1.38</td></tr><tr><td>56735.19</td><td>162.3100353 ± 1.52</td><td>-53.3181895 ± 1.37</td><td>162.3096414 ± 1.62</td><td>-53.3179926 ± 1.34</td></tr><tr><td>56746.12</td><td>162.3099561 ± 1.56</td><td>-53.3181720 ± 1.35</td><td>162.3095644 ± 1.64</td><td>-53.3179776 ± 1.33</td></tr><tr><td>56757.11</td><td>162.3098827 ± 1.52</td><td>-53.3181511 ± 1.33</td><td>162.3094951 ± 1.53</td><td>-53.3179596 ± 1.35</td></tr><tr><td>56773.00</td><td>162.3097794 ± 1.51</td><td>-53.3181165 ± 1.34</td><td>162.3093965 ± 1.51</td><td>-53.3179293 ± 1.35</td></tr><tr><td>56782.99</td><td>162.3097251 ± 1.82</td><td>-53.3180926 ± 1.31</td><td>162.3093454 ± 1.88</td><td>-53.3179084 ± 1.78</td></tr><tr><td>56795.98</td><td>162.3096628 ± 1.54</td><td>-53.3180604 ± 1.31</td><td>162.3092875 ± 1.60</td><td>-53.3178793 ± 1.32</td></tr><tr><td>57013.31</td><td>162.3092730 ± 1.68</td><td>-53.3180357 ± 1.48</td><td>162.3089733 ± 2.39</td><td>-53.3179157 ± 1.70</td></tr><tr><td>57018.33</td><td>162.3092495 ± 1.63</td><td>-53.3180422 ± 1.35</td><td>162.3089518 ± 1.70</td><td>-53.3179239 ± 1.38</td></tr><tr><td>57025.34</td><td>162.3092125 ± 1.50</td><td>-53.3180504 ± 1.34</td><td>162.3089167 ± 1.51</td><td>-53.3179339 ± 1.34</td></tr><tr><td>57034.35</td><td>162.3091610 ± 1.53</td><td>-53.3180597 ± 1.33</td><td>162.3088688 ± 1.56</td><td>-53.3179459 ± 1.35</td></tr><tr><td>57044.36</td><td>162.3090986 ± 1.58</td><td>-53.3180676 ± 1.42</td><td>162.3088113 ± 1.67</td><td>-53.3179570 ± 1.53</td></tr><tr><td rowspan="7">GeMS</td><td>56701.22</td><td>162.3102827 ± 1.72</td><td>-53.3182209 ± 1.72</td><td>162.3098775 ± 1.72</td><td>-53.3180149 ± 1.72</td></tr><tr><td>56759.18</td><td>162.3098694 ± 1.72</td><td>-53.3181510 ± 1.72</td><td>162.3094820 ± 1.72</td><td>-53.3179602 ± 1.72</td></tr><tr><td>56760.17</td><td>162.3098617 ± 1.72</td><td>-53.3181480 ± 1.72</td><td>162.3094745 ± 1.72</td><td>-53.3179576 ± 1.72</td></tr><tr><td>56804.07</td><td>162.3096275 ± 1.72</td><td>-53.3180436 ± 1.72</td><td>162.3092542 ± 1.72</td><td>-53.3178649 ± 1.72</td></tr><tr><td>57000.34</td><td>162.3093359 ± 1.72</td><td>-53.3180168 ± 1.72</td><td>162.3090309 ± 1.72</td><td>-53.3178932 ± 1.72</td></tr><tr><td>57086.32</td><td>162.3087965 ± 1.72</td><td>-53.3180638 ± 1.72</td><td>162.3085243 ± 1.72</td><td>-53.3179656 ± 1.72</td></tr></tbody></table>

Note. Observed astrometry for Luhman 16 AB measured as described in Section 3. The CDAR is not subtracted from these data.

calibrate the extracted spectra using a forward-modeling approach (Blake et al. 2007; Bean et al. 2010).

In chronological order, we measure a relative velocity ($v_A - v_B$) for this system of 2740, 1940, and 1850 m s⁻¹, as plotted in Figure 11 and shown in Table 3. Though the relative radial velocity uncertainties are quite small (tens of m s⁻¹), the uncertainties are in fact dominated by systematics resulting from the inhomogeneous brown dwarf surface (and the induced nonuniformities in the brown dwarf line profiles). We conservatively estimate the uncertainty at each of the three epochs to be 200 m s⁻¹, consistent with the internal scatter in these measurements reported by Crossfield et al. (2014).

### 3. Astrometry and Uncertainty Analysis

The processed Luhman 16 astrometry from GeMS, VLT/FORS2, and the ESO photographic plate are shown in Table 1. This astrometry still contains chromatic differential atmospheric refraction (CDAR), which is dealt with as detailed in Section 4.1. The CDAR-subtracted separations are shown in Table 2.

All astrometry in this work is ultimately tied to the ICRS system via GAIA DR1 measurements of background stars (Gaia

Collaboration et al. 2016a, 2016b). The linking of reference frames is described in Section 3.5. The first linking of reference frames is between a distortion-corrected FORS2 reference frame and GAIA. The other data sets, including the ESO photographic plate and the distortion-corrected GeMS data, are linked to the FORS2 reference frame via matched background stars. In each data set, individual frames are never stacked before measuring stellar positions; positions are only averaged over frames following correction for static distortion. Within a particular data set, reference frames are defined by averaging positions of background stars over the entire data set following a rough alignment, as described below in Section 3.5. These unique reference frames are used when transforming from camera to camera. Throughout this work, tangent plane coordinates are converted to spherical coordinates using the reference point 10:49:14.765–53:18:15.998, which is on the ICRS system.

#### 3.1. GeMS Distortion Map Measurement with NGC 1851 Observations

In order to obtain Luhman 16 AB astrometry from GeMS data, we first obtain the distortion map over the GeMS field of view. We use an image of globular cluster NGC 1851 to derive this distortion map. We run *Starfinder* (Diolaiti et al. 2000)
<table><caption>**Table 2**
Orbital Astrometry</caption><thead><tr><th>Instrument</th><th>MJD<br>(days)</th><th>ΔX<br>(mas)</th><th>ΔY<br>(mas)</th></tr></thead><tbody><tr><td>ESO<br>Schmidt (Red)</td><td>45764.95</td><td>-1125.97 ± 160.00</td><td>1408.10 ± 135.00</td></tr><tr><td rowspan="16">FORS2/VLT</td><td>56397.00</td><td>-1042.56 ± 1.41</td><td>1002.32 ± 2.03</td></tr><tr><td>56402.06</td><td>-1039.95 ± 0.47</td><td>999.61 ± 0.81</td></tr><tr><td>56408.04</td><td>-1036.84 ± 0.47</td><td>995.03 ± 1.02</td></tr><tr><td>56418.03</td><td>-1030.59 ± 1.10</td><td>986.03 ± 2.30</td></tr><tr><td>56424.09</td><td>-1026.21 ± 0.65</td><td>981.09 ± 1.63</td></tr><tr><td>56432.97</td><td>-1023.17 ± 0.81</td><td>974.80 ± 1.60</td></tr><tr><td>56438.03</td><td>-1020.58 ± 0.66</td><td>970.98 ± 1.12</td></tr><tr><td>56443.99</td><td>-1017.40 ± 0.86</td><td>965.15 ± 1.37</td></tr><tr><td>56448.00</td><td>-1014.90 ± 0.65</td><td>962.26 ± 1.07</td></tr><tr><td>56452.97</td><td>-1012.86 ± 0.53</td><td>958.55 ± 0.99</td></tr><tr><td>56457.00</td><td>-1009.39 ± 0.64</td><td>956.42 ± 1.39</td></tr><tr><td>56459.97</td><td>-1009.20 ± 1.06</td><td>952.25 ± 0.69</td></tr><tr><td>56465.96</td><td>-1005.63 ± 0.78</td><td>947.62 ± 1.98</td></tr><tr><td>56693.25</td><td>-874.99 ± 1.44</td><td>749.51 ± 2.28</td></tr><tr><td>56702.19</td><td>-871.52 ± 0.89</td><td>741.83 ± 1.68</td></tr><tr><td>56726.14</td><td>-853.54 ± 2.36</td><td>718.19 ± 3.34</td></tr><tr><td>56735.19</td><td>-847.04 ± 1.08</td><td>709.08 ± 1.35</td></tr><tr><td>56746.12</td><td>-842.25 ± 1.00</td><td>699.73 ± 1.29</td></tr><tr><td>56757.11</td><td>-833.66 ± 1.24</td><td>689.32 ± 1.95</td></tr><tr><td>56773.00</td><td>-823.64 ± 0.81</td><td>673.97 ± 1.25</td></tr><tr><td>56782.99</td><td>-816.79 ± 0.97</td><td>663.35 ± 2.68</td></tr><tr><td>56795.98</td><td>-807.20 ± 1.24</td><td>652.02 ± 1.32</td></tr><tr><td>57013.31</td><td>-644.66 ± 1.85</td><td>432.15 ± 1.25</td></tr><tr><td>57018.33</td><td>-640.39 ± 2.46</td><td>426.14 ± 1.95</td></tr><tr><td>57025.34</td><td>-636.18 ± 1.33</td><td>419.43 ± 1.68</td></tr><tr><td>57034.35</td><td>-628.30 ± 1.11</td><td>409.74 ± 1.78</td></tr><tr><td>57044.36</td><td>-617.60 ± 1.34</td><td>398.22 ± 2.12</td></tr><tr><td rowspan="6">GeMS</td><td>56701.22</td><td>-871.33 ± 0.26</td><td>741.46 ± 0.28</td></tr><tr><td>56759.18</td><td>-833.18 ± 0.26</td><td>686.77 ± 0.28</td></tr><tr><td>56760.17</td><td>-832.71 ± 0.26</td><td>685.66 ± 0.28</td></tr><tr><td>56804.07</td><td>-802.82 ± 0.26</td><td>643.57 ± 0.28</td></tr><tr><td>57000.34</td><td>-655.87 ± 0.26</td><td>444.92 ± 0.28</td></tr><tr><td>57086.32</td><td>-585.38 ± 0.26</td><td>353.41 ± 0.28</td></tr><tr><td>DENIS</td><td>51220.50</td><td>&lt;675.00</td><td>&lt;675.00</td></tr></tbody></table>

Note. Observed separations for Luhman 16 AB measured as described in Section 3. We present an upper limit on the separation at the epoch of the DENIS observations as Luhman 16 AB is unresolved. The maximum-likelihood CDAR is subtracted from these orbital astrometry observations, as detailed in Sections 4.1 and 4.2.

on individual frames, using 20 bright, unsaturated stars distributed throughout the field to create a master PSF for each frame. A first extraction is performed down to a $20\sigma$ limit, secondary stars in the PSF images are subtracted, and a new master PSF free from secondary sources is created. Stars are then re-extracted down to $7\sigma$ with a PSF correlation threshold of 0.7. The *Starfinder* function *matchCoord.pro* is used to approximately register each frame to the first in each set. Frames with fewer than 1500 stars matching the reference frame's list to within three pixels are rejected. This step trims Set 1 from 18 to 13 frames and Set 2 from 26 to 24 frames.

To establish a set of common stars in all frames, stars falling more than three pixels from the reference frame's position in any of the frames are rejected from the list, leaving about 1000 stars in each set in the full GSAOI field. Stars located within 10 pixels of any other common stars are also eliminated. Stars within two pixels of a bad pixel or within 10 pixels of the edge of the $85'' \times 85''$ field of view are rejected. If *Starfinder*

splits a single star into multiple detections in a single frame, those measurements are excluded, but other frames with a single detection are retained for that star. Following these exclusions, ≈1200 stars remain for each set.

For each frame, a four-parameter linear transformation is performed to transform each GeMS chip to the appropriate area on the ACS frame to enable matching between GeMS stars and ACS stars. The deviations between the ACS and original, untransformed GeMS positions are measured and recorded for each star in each frame (a total of $2 \times 10^5$ stellar measurements). A global 120-parameter fourth-order static polynomial distortion pattern (30 parameters per chip) in addition to a six-parameter linear transformation for each frame is fitted to all stellar deviation measurements using the IDL *MPFITFUN* minimization routine. The nonlinear components of the model remain constant to within 2 mas across the full field when the position angle of the instrument is rotated and the dither position is changed by up to 30''.

The resulting nonlinear components of the model distortion map are shown in Figure 2. These vectors have been obtained by applying the distortion correction model to a regularly spaced grid of points and fitting and subtracting a six-parameter linear transformation to each chip.

## 3.2. GeMS Luhman 16 AB Position Measurements

The GSAOI images of Luhman 16 AB were analyzed by first running *Starfinder* on individual frames. A $16 \times 16$ pixel ($0''32 \times 0''32$) image of Luhman A is used as the reference PSF for cross-correlation with other stars in the field. The correlation between the Luhman 16 A and Luhman B PSFs is consistently above 0.99, and the median value over all frames is 0.994. We establish a set of common stars in all frames using a procedure similar to that described above for the NGC 1851 analysis. If *Starfinder* splits a single star into multiple detections in a single frame, those measurements are excluded (but other frames with a single detection are retained for that star).

The frames are distortion corrected using the distortion map measured from NGC 1851 images, as described in the section above. Following distortion correction, the stellar positions for all epochs are transformed to the same astrometric frame with six-parameter linear transformations measured from ≈12 background stars.

## 3.3. FORS2 Luhman 16 AB Position Measurements

We use *Starfinder* to measure stellar positions in the FORS2 chip 1 images. We use 20 bright, unsaturated stars distributed throughout the field to create a master PSF for each frame. A first extraction is performed down to a $20\sigma$ limit, secondary stars in the PSF images are subtracted, and a new master PSF free from secondary sources is created. Stars are then re-extracted down to $7\sigma$ with a PSF correlation threshold of 0.7. The *Starfinder* function *matchCoord.pro* is used to approximately register each frame to the first in each set. Frames with less than 300 stars matching the reference frame's list to within three pixels are rejected.

To establish a set of common stars in all frames, stars falling more than three pixels from the reference frame's position in any of the frames are rejected from the list, leaving about 200 stars in each set in the full FORS2 field. Stars located within 10 pixels of any other common stars are also eliminated. Stars
<figure><img src="image_4.png" /><figcaption>Figure 2. Nonlinear components of GSAOI distortion map. Plotted vectors and distortion map are derived as described in Section 3.1.</figcaption></figure>

within two pixels of a bad pixel or within 10 pixels of the edge of the field of view are rejected. If Starfinder splits a single star into multiple detections in a single frame, those measurements are excluded (but other frames with a single detection are retained for that star).

### 3.3.1. Generation of FORS2 Static Distortion Map

For each dithered frame, a six-parameter linear transformation is performed to match the FORS2 chip to the GAIA DR1 positions using ≈200 background stars. This transformation removes the linear effects of achromatic differential atmospheric refraction (ADAR; Chapter 4, Lu 2008), but nonlinear chip distortions have not yet been removed. The deviations between the GAIA positions and the original, untransformed FORS2 positions are measured and recorded for each star in each frame (a total of 1.3 × 10⁵ stellar measurements). A global 20-parameter third-order static polynomial distortion pattern in addition to a six-parameter linear transformation for each frame is fitted to all stellar deviation measurements. The six additional parameters for each frame fit the 9 × 9 dithering pattern used in the original observations. The output of this procedure is a set of polynomial coefficients that can be used to correct the nonlinear components of FORS2 camera distortion.

### 3.4. ESO R-band Plate Position Measurements

We run Starfinder on the digitized ESO 1984 R-band plate to detect background stars. A 5 × 5 pixel (3.4 × 3.4) image of a bright, unsaturated star is used as the reference PSF for cross-correlation with other stars in the field. As expected for a digitized photographic plate with partially saturated background stars, correlations between stars and the reference PSF are not high. The median correlation is 0.85, and the threshold for inclusion is 0.6.

### 3.5. Transforming GeMS, FORS2, and ESO R-band Plate Astrometry to a Common Reference Frame

Reference frames for each camera (GeMS and FORS2) are established according to the following procedure. First, following distortion correction on the position measurements, each frame is roughly aligned to the frame with the highest signal-to-noise ratio (S/N) with two shift parameters and a

rotation angle. Positions of reference stars are then averaged to define that camera's unique reference frame. For the ESO R-band plate, there is only one image in the data set, so the raw positions of the background stars define that reference frame.

The GeMS and ESO R-band data sets are then transformed to the FORS2 reference frame, and then GeMS, FORS2, and ESO R-band data are transformed to the GAIA DR1 ICRS system (Gaia Collaboration et al. 2016a, 2016b). The GeMS system is transformed to the FORS2 system with a six-parameter linear transformation of the general form of Equation (1) via matching of 14 isolated background stars (15 < R < 18). The per-coordinate standard deviation of the residuals of that fit is 2.6 mas. These low residuals suggest that the high-order camera distortion corrections for both the GeMS and FORS2 instruments are accurate to within 2–3 mas over the 1.4 field of view of the GeMS camera.

Next, we remove ADAR from our FORS2 reference frame. Rather than calculate a theoretical ADAR as a function of zenith angle, temperature, pressure, relative humidity, and wavelength of observation (Gubler & Tytler 1998; Lu 2008), we transform the FORS2 system to the ADAR-free GAIA DR1 ICRS system (Gaia Collaboration et al. 2016a, 2016b), using a six-parameter linear fit:

$$ 
\begin{aligned}
x' &= Ax + By + C \\
y' &= Dx + Ey + F,
\end{aligned}
\qquad (1) $$

where the coefficients A, B, C, and D account for the linear ADAR and a measurement of the FORS2 plate scale, and parameters E and F account for the absolute shift and rotation in coordinates between the FORS2 system and the GAIA DR1 ICRS system. Thus, we empirically remove the linear ADAR from our FORS2 and GeMS observations of Luhman 16 AB. The quadratic component of ADAR is negligible in our case. We transformed the FORS2 system to the GAIA DR1 ICRS with a six-parameter linear transformation above fitted to 194 matching isolated stars (15 < R < 18). The standard deviation of the residuals of that fit is 1.2 mas. Following the transformations derived from background star positions, the Luhman 16 AB separations as measured by FORS2 and GeMS agree to within the FORS2 error bars (1–2 mas) for comparable epochs.

The ESO R plate is linked to the FORS2 system via a transformation consisting of a pair of second-order polynomials. These polynomials are fitted to a set of 27 common background stars within a 1.0 radius of Luhman 16 AB. This single step accounts for DAR and camera distortion in the ESO R plate. The standard deviation of the residuals about that fit is 84 mas.

### 3.6. ESO R-band Plate Orbital Astrometry

Given the long time baseline between the 1984 ESO Schmidt photometric plate image of Luhman 16 AB and our 2013–2015 VLT/FORS2 and GeMS data, the position of either star in 1984 plays a critical role in breaking the mass-period–eccentricity degeneracies in our orbit model. Therefore, we paid special attention to deriving a separation and position angle for Luhman 16 AB from our ESO 1984 plate.

First, we determined that the best PSF model using background stars is a Lorentzian PSF. We fit Moffat, Gaussian, and Lorentzian PSF functions to isolated, unsaturated stars in
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top left: ESO 1984 photometric plate image of Luhman 16 AB. Top right: minimum-variance binary PSF Lorentzian model (see Section 3.6). Bottom left: minimum-variance “data-model” residuals, with the upper right few pixels masked by the presence of another star. Bottom right: The histogram of the residuals is Gaussian distributed, with median absolute residuals of 7.4%.</figcaption></figure>

the plate image $\lesssim 2'$ of Luhman 16 AB, finding Lorentzian PSFs to yield the lowest “model-data” residuals for fits to each nearby background star. The PSF models are constructed at $10 \times$ the 0.67 pixel$^{-1}$ plate scale.

Next, to model the binary Luhman 16 AB, we fit a binary model of two co-added, shifted, scaled Lorentzian PSFs to the 1984 ESO image of Luhman 16 AB using a variant of IDL's MPFITFUN, which minimizes the $\chi^2$ using a Levenberg-Marquardt implementation in IDL (MPFIT, Markwardt 2009). We fit for nine free parameters: the sky background, the fluxes of Luhman 16 A and B, the positions X and Y of Luhman 16 A and B on the photometric plate, and the X axis and Y axis Lorentzian widths of the Lorentzian PSFs. We fit our binary model to an $11 \times 11$ pixel (7.37 $\times$ 7.37) region centered on Luhman 16 AB. We masked a few pixels on the upper right-hand corner of the image that contained a nearby star (Figure 3).

The above procedure produces a $\chi^2$ map of parameter space for each binary model as a function of the parameters of interest. Our final contour map of $\Delta\chi^2 = \chi^2 - \chi_{\min}^2$ is shown in Figure 4 (left panel) as dark green, medium green, and light green solid lines for $1\sigma$, $2\sigma$, and $3\sigma$ contours corresponding to $\Delta\chi^2 < 10.4, 17.2, \text{and } 25.2$, respectively (Press et al. 2002).

We compare our minimum-variance binary Lorentzian PSF model to our 1984 ESO image in Figure 3. We find the median “data-model” residuals are 7.5% (bottom right panel, Figure 3).

Our minimum-variance model gives a reduced chi-squared of $\chi_{\text{red}}^2 = 0.45$ for 112 degrees of freedom. Our minimum-variance PA 319 $\pm$ 6$^\circ$ matches within uncertainty the 138$^\circ$ elongation reported by Mamajek (2013), albeit with a 180$^\circ$ degeneracy. Our minimum-variance separation is $1801 \pm 200$ mas.

We took several steps to ensure our minimum-variance binary PSF model was not a local optimum in $\chi^2$ space. We ran MPFITFUN twice with a maximum of 3000 iterations, with the second run starting at the minimum-variance parameters of the first run. We tested the sensitivity of our minimum-variance binary PSF model to the initial guess, performing $10^6$ MPFITFUN runs as detailed above, with random initializations for each initial guess. MPFITFUN converged to our best-fit parameters the majority of the time, meaning our minimum-variance solution is robust to the initial guess.

The southeastern component of the binary was identified as Luhman 16 A with a spectral type of L7.5 (Kniazev et al. 2013; Luhman 2013; Burgasser et al. 2014). We identify Luhman 16 B as the northwestern component of the binary, which has a spectral type of T0.5 (Kniazev et al. 2013; Luhman 2013; Burgasser et al. 2014). Similarly, the brighter, southeastern component of the binary in 2013–2014 was identified by Street et al. (2015) in a 700–950 nm filter (see Figure 3). We find the southeast component of the binary in the ESO Schmidt plate to have a flux ratio $F_1/F_2 \approx 1.2$ relative to the northwest component of the binary at the R band. Thus we conclude
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left: $\Delta\chi^2$ 1$\sigma$, 2$\sigma$, and 3$\sigma$ contours for the binary separation from our binary Lorentzian PSF model fits to the 1984 ESO resolved photometric plate of Luhman 16 AB (see Section 3.6). We find a minimum-variance separation of 1801 ± 200 mas and PA of 319 ± 6°, in good agreement with the 138° elongation measured from this photometric plate by Mamajek (2013). Right: contours for the binary separation from the DENIS I-band unresolved image of Luhman 16 AB (see Section 3.7), which places a constraint on a maximum $\Delta$R.A. $\lesssim$ 675 mas and $\Delta$decl. $\lesssim$ 675 mas, due to excellent seeing conditions.</figcaption></figure>

that the orientation of the binary in the 1984 epoch is the same as in 2013; that is, Luhman 16 A is in the southeastern position, and Luhman 16 B is in the northwestern position.

### 3.7. I-band DENIS 1999 Orbital Astrometry Upper Limit

Luhman 16 AB is unresolved in 1999, as shown in Figure 1 (upper right panel). We quantify this upper limit on the binary separation at all position angles as follows. First, we use the same procedure detailed in Section 3.6, to find a minimum-variance circular Lorentzian PSF model for background stellar PSFs in the image. We use Lorentzian models rather than Gaussian models as the mean reduced chi-squared $\chi^2_{red}$ was significantly lower (1.17 versus 1.96). We fit our model Lorentzian PSF to a $7 \times 7$ pixel (1"/pixel$^{-1}$ plate scale) region centered on each stellar PSF, finding an average Lorentzian FWHM of 0"79 ± 0"05 for reference stars in the field near Luhman 16 AB. We used these average Lorentzian parameters as input for the next step.

Similar to Section 3.6, we use a binary PSF of two co-added, scaled, and shifted Lorentzians to place a constraint on the maximum possible separation of the Luhman 16 AB binary while remaining consistent with the unresolved DENIS image. We minimized the $\chi^2$ using MPFIT as above, with five parameters of interest: the position $X_1$, $Y_1$ of the primary, the position $X_2$, $Y_2$ of the secondary, and the total flux of the system. We started our binary fitting routine with initial guesses from a grid of 40 × 40 different positions, with step sizes of 0"1, for a total of 1600 runs. A contour map of $\Delta\chi^2 = \chi^2 - \chi^2_{min}$ is shown as the solid dark blue, blue, and light blue lines in Figure 4 (right panel) for 1$\sigma$, 2$\sigma$, and 3$\sigma$ confidence. We find the best-fit position angle can vary by ±180° (dotted lines in Figure 4), and we use both solutions in fitting the orbit of Luhman 16 AB (see Section 4.2).

### 3.8. Uncertainties on Stellar Positions

There is significant systematic epoch-to-epoch variation in the observations of Luhman 16 AB with GeMS and FORS2. We cannot characterize systematic uncertainty for Luhman 16 using Luhman 16's position measurements directly, given that the binary has orbital motion across multiple epochs. However,

we can at least characterize the random terms on a per-epoch basis by computing the standard deviation of the Luhman 16 AB position measurements in a set of frames within an epoch $t$ (we denote this parameter $\sigma_{\alpha_A}(t)$). To estimate the systematic component of the position uncertainty, we measure the epoch-to-epoch variation of the mean position of background stars as described below.

Section 3.8.1 describes the measurement of position uncertainties in the FORS2 reference frame, Section 3.8.2 details the measurement of position uncertainties in the GeMS reference frame, and Section 3.9 describes the methods used to estimate the uncertainties on the Luhman 16 AB separation. Later in the paper, a subsection gives the methodology for estimating the Luhman 16 AB barycenter uncertainties from the position uncertainties.

#### 3.8.1. Stellar Position Uncertainties for FORS2

We assign uncertainties on the mean position measurements of the individual components of Luhman 16 AB in the FORS2 reference frame at an epoch $t$ as follows:

$$ 
\begin{aligned}
\sigma^2_{F,\alpha_A,\text{tot}}(t) &= \sigma^2_{F,\alpha_A}(t)/N(t) + \sigma^2_{F,\text{sys,L16}} \\
\sigma^2_{F,\delta_A,\text{tot}}(t) &= \sigma^2_{F,\delta_A}(t)/N(t) + \sigma^2_{F,\text{sys,L16}} \\
\sigma^2_{F,\alpha_B,\text{tot}}(t) &= \sigma^2_{F,\alpha_B}(t)/N(t) + \sigma^2_{F,\text{sys,L16}} \\
\sigma^2_{F,\delta_B,\text{tot}}(t) &= \sigma^2_{F,\delta_B}(t)/N(t) + \sigma^2_{F,\text{sys,L16}},
\end{aligned}
\quad (2) $$

where $\sigma_{F,\alpha_A/B}(t)$ and $\sigma_{F,\delta_A/B}(t)$ are the standard deviations of the FORS2 position measurements of Luhman 16 A/B over all frames within an epoch $t$, $N(t)$ is the number of frames in epoch $t$, and $\sigma_{\text{sys,L16}}$ is an estimate of the systematic component of the uncertainty of the positions of the individual components of Luhman 16.

Stars in the background have predictable celestial motion and can serve as ideal references whose position residuals, following subtraction of a fitted proper motion and parallax model, reflect a more accurate measure of the total uncertainty. We assume that the total uncertainties on the Luhman 16 positions are comparable to the uncertainties of background
<figure><img src="image_7.png" /><figcaption>Figure 5. Position uncertainties over all epochs for VLT/FORS2 observations of background stars (black circles) and Luhman 16 AB (red triangles). The plotted uncertainties for Luhman 16 AB are the median values calculated over all epochs.</figcaption></figure>

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

stars of similar brightness:

$$ \text{median}(\sigma_{F,\text{tot},\text{L16}}(t)) \sim \sigma_{F,\text{std},\text{bkgd}}. \quad (3) $$

To estimate $\sigma_{F,sys,bkgd}$, we first obtain the mean positions of isolated background stars averaged over frames in each epoch t. Next, we fit and subtract a minimum-variance parallax/proper-motion model for each background star. For the set of background stars near Luhman 16 AB in the image and of comparable brightness, we compute the standard deviation of the residual positions over all epochs (black circles in Figure 5).

Substituting into Equation (2), the scalar $\sigma_{F,sys,L16}$ can then be approximated as

$$ \sigma_{F,sys,L16}^2 \sim \sigma_{F,\text{std},\text{bkgd}}^2 - \text{median}(\sigma_{F,\alpha/\delta_{A/B}}^2(t)/N(t)), \quad (4) $$

which is a function of quantities measurable from the data itself. The scalar $\sigma_{F,sys,L16}$ can then be used in Equation (2) to calculate the vectors of uncertainties $\sigma_{F,\alpha_A/B,\text{tot}}(t)$ and $\sigma_{F,\delta_A/B,\text{tot}}(t)$.

To estimate the errors in the ICRS reference frame, we add in quadrature a term $\sigma_{F-1}$ accounting for uncertainties introduced by converting the FORS2 astrometric system to the ICRS system. We estimate $\sigma_{F-1}$ as the standard deviation of the reference star positions in 2015 as measured by GAIA minus the transformed FORS2 positions at the 2015 epoch. This value, averaged over both axes, is 1.14 mas. Our final position uncertainties in the ICRS system are

$$ \begin{aligned} \sigma_{L,\alpha_A/B,\text{tot}}^2(t) &= \sigma_{F,\alpha_A/B,\text{tot}}^2(t) + \sigma_{F-1}^2 \\ \sigma_{L,\delta_A/B,\text{tot}}^2(t) &= \sigma_{F,\delta_A/B,\text{tot}}^2(t) + \sigma_{F-1}^2 \end{aligned} \quad (5) $$

and are given in Table 1.

### 3.8.2. Stellar Position Uncertainties for GeMS

Given that we have only one to four frames per epoch, we are unable to estimate a per-epoch uncertainty using the methods above (Section 3.8.1). The PSF FWHMs delivered by the GeMS AO system do not vary significantly across epochs (~100 mas FWHM, ~20% variation), and the signal-to-noise ratios of the detections are typically quite high (500–1000), so the contribution from random noise sources like shot noise and read noise is expected to be small (~0.1 mas). The uncertainties are dominated by terms associated with

conversion between reference frames rather than poor S/N. Thus, we assign the same Luhman 16 AB position uncertainty to all epochs.

We estimate the position uncertainty for Luhman 16 AB in the ICRS reference frame in both R.A. and decl. as

$$ \sigma_{L,\text{tot}}^2 = \sigma_{F-1}^2 + \sigma_{\text{Ge-F}}^2 + \sigma_{\text{GeMS-internal}}^2, \quad (6) $$

where the additional term $\sigma_{\text{Ge-F}} = 1.22$ mas, as above (Section 3.8.1), models uncertainties introduced by transforming the GeMS astrometric system to the FORS2 system. We estimate the uncertainty $\sigma_{\text{GeMS-internal}} = 0.41$ mas by taking the standard deviation of the mean of background star positions over all epochs. We estimate $\sigma_{\text{Ge-F}}$ by computing the standard deviation of the reference star positions as measured by FORS2 minus the transformed GeMS positions. The sum of all uncertainties in quadrature is $\sigma_{L,\text{tot}} = 1.72$ mas.

### 3.9. Luhman 16 AB Separation Uncertainties for FORS2 and GeMS

The uncertainty on the mean separation between components of Luhman 16 AB is expected to be smaller than the uncertainties on the mean positions of the individual components. Most sources of uncertainty are correlated for closely spaced stars, including the differential tip/tilt jitter (DTTJ), PSF anisotropy, and other systematic terms (Ammons et al. 2011, 2012). Thus, we make estimates of the uncertainties on the separations $\sigma_{\alpha_A-\alpha_B}(t)$, $\sigma_{\delta_A-\delta_B}(t)$ at epoch t independently from the uncertainties on the mean positions of the individual components, as done in Section 3.8.

In the FORS2 reference frame, the uncertainties on the mean separation measurements of the individual components of Luhman 16 AB at an epoch t are

$$ \begin{aligned} \sigma_{F,\alpha_A-\alpha_B,\text{tot}}^2(t) &= \sigma_{F,\alpha_A-\alpha_B}^2(t)/N(t) + \sigma_{F,\text{sys,L16}}^2 \\ \sigma_{F,\delta_A-\delta_B,\text{tot}}^2(t) &= \sigma_{F,\delta_A-\delta_B}^2(t)/N(t) + \sigma_{F,\text{sys,L16}}^2, \end{aligned} \quad (7) $$

where $\sigma_{F,\alpha_A-\alpha_B}(t)$ and $\sigma_{F,\delta_A-\delta_B}(t)$ are the standard deviations of the separation measurements over frames within an epoch t. We estimate the systematic component of the uncertainty on the separation of Luhman 16 AB $\sigma_{F,sys,L16}$ using a calibration background double star. The background double star has a separation of 1"8 and a position of 10:49:11.81–53:18:57. This
<table><caption>**Table 3**
CRIRES Relative Radial Velocity</caption><thead><tr><th>MJD<br>(days)</th><th>ΔV<br>(m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>56417.5</td><td>2740 ± 200</td></tr><tr><td>56779.5</td><td>1940 ± 200</td></tr><tr><td>56797.5</td><td>1850 ± 200</td></tr></tbody></table>

**Note.** VLT/CRIRES observed the relative radial velocities of Luhman 16 AB used in fitting the orbital parameters, as detailed in Section 2.5.

calibration star is ≈0.4 away from Luhman 16 AB. We chose this double star because it has a separation and flux ratio similar to the components of Luhman 16 AB. However, the double star is fainter by ≈2 mag in the I band and thus has more Poisson noise. We assume that the true space motion of this double star is defined by pure parallax and proper motion for the two components. With a component magnitude of V ≈ 18.5, the double star is likely composed of unbound stars at distances beyond 1 kpc.

For the GeMS data, the double star used above does not yield a usable calibration of the uncertainty as the S/N is lower. There are only between one and three frames per epoch in the Luhman 16 separation data from GeMS, so there are not enough data to estimate the uncertainties from the data itself for each epoch. Instead, for epochs with more than one frame, we conservatively compute the standard deviation of the separation measurements and average these over all epochs to estimate a single uncertainty on the separation measurement, 0.26 mas for R.A. and 0.28 mas for decl.

## 4. MCMC Analysis

We constrain a joint 16-parameter astrometry model using a Markov chain Monte Carlo (MCMC) searching algorithm. This model includes the orbital astrometry parameters (Keplerian orbital parameters and total dynamical mass) as well as the barycentric astrometry parameters (mass ratio, parallax, barycentric proper motion, position at a reference epoch, and CDAR parameters for both FORS2 and GeMS data sets). The starting priors for each parameter are given in Table 4. With our simultaneous fit of all parameters, we can inspect joint posterior probability distributions for correlated uncertainty between model parameters derived from orbital and barycentric astrometry. We present the results of this analysis as well as the individual masses $M_A$ and $M_B$ for Luhman 16 A and B in Section 5.

We perform the following procedure in every step of our MCMC chains, detailed in the next few sections:

1. Subtract the chromatic differential atmospheric refraction from the coordinate positions for both brown dwarfs (see Section 4.1).

2. Compute the Luhman 16 separation for all data from the CDAR-corrected coordinate positions, and then compute the logarithm of the likelihood (“log likelihood”) for the relative orbit terms $L_{rel}$ (see Section 4.2).

3. Compute relative radial velocity log likelihood $L_{RV}$ (see Section 4.3).

4. Compute the barycenter positions and the orbit log likelihood $L_{abs}$ from the CDAR-corrected coordinate positions for ESO 1984, FORS2, and GeMS data (see Section 4.4).

<table><caption>**Table 4**
Uniform Priors for Our MCMC Analysis</caption><thead><tr><th>Parameter</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td colspan="3"><strong>Orbital Parameters</strong></td></tr><tr><td>Semimajor Axis, a (au)</td><td>2.0</td><td>50.0</td></tr><tr><td>Total Mass, M<sub>tot</sub> (M<sub>J</sub>)</td><td>10.5</td><td>125.7</td></tr><tr><td>Eccentricity, e</td><td>0.0</td><td>0.9</td></tr><tr><td>Inclination, cos i</td><td>-0.9</td><td>1.0</td></tr><tr><td>Longitude of Ascending Node, Ω (°)</td><td>0.0</td><td>360.0</td></tr><tr><td>Argument of Periastron, ω (°)</td><td>0.0</td><td>360.0</td></tr><tr><td>Time of Periastron Passage, τ (Reduced Julian Day)</td><td>10,000</td><td>80,000</td></tr><tr><td colspan="3"><strong>Barycenter Parameters</strong></td></tr><tr><td>Mass Ratio, q</td><td>0.2</td><td>1.5</td></tr><tr><td>Ref R.A., α<sub>0</sub> (°)</td><td>162.310622</td><td>162.311178</td></tr><tr><td>Ref Decl., δ<sub>0</sub> (°)</td><td>-53.318413</td><td>-53.317857</td></tr><tr><td>Proper Motion R.A., μ<sub>α</sub> (")</td><td>-2.8</td><td>-2.7</td></tr><tr><td>Proper Motion Decl., μ<sub>δ</sub> (")</td><td>0.3</td><td>0.4</td></tr><tr><td>Parallax, ω<sub>rel</sub> (")</td><td>0.49</td><td>0.51</td></tr><tr><td colspan="3"><strong>Instrument Parameters</strong></td></tr><tr><td>CDAR FORS2, ρ<sub>FORS2</sub> (mas)</td><td>25.0</td><td>50.0</td></tr><tr><td>CDAR FORS2, d<sub>FORS2</sub> (mas)</td><td>-58.9</td><td>-38.5</td></tr><tr><td>CDAR GeMS, ρ<sub>GeMS</sub> (mas)</td><td>0.0</td><td>5.0</td></tr></tbody></table>

**Note.** Uniformly sampled prior probability density distributions used in our MCMC analysis (see Section 4).

5. Compute the sum of the negative log likelihood for the orbital, barycentric motion, and radial velocity parameters $L_{rel} + L_{abs} + L_{RV}$.

We use uniformly distributed priors for all 16 parameters, with ranges shown in Table 4. The maximum-likelihood parameters and 68.3% and 95.4% confidence intervals from our MCMC analysis are shown in Table 5.

### 4.1. Modeling Chromatic Differential Atmospheric Refraction

Chromatic differential atmospheric refraction produces a change in the measured separation vector between any two stars that depends on the zenith angle, the parallactic angle, atmospheric pressure, atmospheric temperature, and the two stellar spectra. Astrometric measurements in VLT/FORS2 have the additional complication that a majority of the CDAR offset is corrected with the Longitudinal Atmospheric Dispersion Corrector (LADC), a maneuvering prism with a single degree of freedom (Avila et al. 1997).

At each step in our MCMC chains, we model the CDAR for both the FORS2 and GeMS astrometry following Lazorenko et al. (2009), Sahlmann et al. (2013), and Sahlmann & Lazorenko (2015). In measuring proper motion and parallaxes of single stars in FORS2 data, these authors find that the offsets produced by CDAR from the atmosphere and the LADC are sufficiently modeled by two parameters ρ and d, in units of mas. The first of these, ρ, accounts for the atmospheric CDAR offset, and the second, d, accounts for the LADC correction of CDAR. The correction of CDAR by the LADC is imperfect because the LADC corrects for CDAR at the zenith angle of the beginning of a set of observations, but does not update during the observation. There is a small offset between the actual CDAR ρ and the correction of CDAR by the LADC d.

<footnote>ª CDAR is chromatic differential atmospheric refraction.</footnote>
<table><caption>**Table 5**
MCMC Posterior Probability Distributions</caption><thead><tr><th>Parameter</th><th>Sample Fit<sup>a</sup></th><th>Median</th><th>68.3% c.i.</th><th>95.4% c.i.</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-weight: bold;">Orbital Parameters</td></tr><tr><td>Orbital Period, P (years)</td><td>27.4</td><td>27.4</td><td>27.6, 27.1</td><td>27.6, 26.7</td></tr><tr><td>Semimajor Axis, a (au)</td><td>3.55</td><td>3.54</td><td>3.59, 3.49</td><td>3.64, 3.43</td></tr><tr><td>Total Mass, M<sub>tot</sub> (M<sub>J</sub>)</td><td>62.5</td><td>62.0</td><td>64.1, 60.3</td><td>66.2, 58.9</td></tr><tr><td>Eccentricity, e</td><td>0.35</td><td>0.35</td><td>0.37, 0.33</td><td>0.39, 0.31</td></tr><tr><td>Inclination, i (°)</td><td>79.6</td><td>79.5</td><td>79.3, 79.8</td><td>79.0, 80.1</td></tr><tr><td>Longitude of Ascending Node, Ω (°)</td><td>130.17</td><td>130.12</td><td>130.24, 130.00</td><td>130.34, 129.87</td></tr><tr><td>Argument of Periastron, ω (°)</td><td>129.6</td><td>130.4</td><td>133.9, 126.9</td><td>137.6, 123.6</td></tr><tr><td>Time of Periastron Passage, τ (day)</td><td>48050</td><td>48030</td><td>48340, 47780</td><td>48720, 47580</td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Barycenter Parameters</td></tr><tr><td>Mass Ratio, q</td><td>0.82</td><td>0.82</td><td>0.85, 0.78</td><td>0.88, 0.75</td></tr><tr><td>Ref R.A., α<sub>0</sub> (°)</td><td>162.31093</td><td>162.31093</td><td>162.31093, 162.31092</td><td>162.31094, 162.31092</td></tr><tr><td>Ref Decl., δ<sub>0</sub> (°)</td><td>-53.31815</td><td>-53.31815</td><td>-53.31814, -53.31815</td><td>-53.31814, -53.31815</td></tr><tr><td>Proper Motion R.A., μ<sub>α</sub> (")</td><td>-2.761</td><td>-2.762</td><td>-2.760, -2.765</td><td>-2.757, -2.767</td></tr><tr><td>Proper Motion Decl., μ<sub>δ</sub> (")</td><td>0.358</td><td>0.358</td><td>0.361, 0.354</td><td>0.365, 0.351</td></tr><tr><td>Parallax, ϖ<sub>rel</sub> (")</td><td>0.50150</td><td>0.50101</td><td>0.50153, 0.50048</td><td>0.50204, 0.49995</td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Instrument Parameters</td></tr><tr><td>CDAR<sup>b</sup> FORS2, ρ<sub>FORS2</sub> (mas)</td><td>29.92</td><td>38.06</td><td>46.09, 29.74</td><td>49.34, 26.13</td></tr><tr><td>CDAR FORS2, d<sub>FORS2</sub> (mas)</td><td>-42.70</td><td>-48.73</td><td>-41.93, -55.32</td><td>-39.06, -58.04</td></tr><tr><td>CDAR GeMS, ρ<sub>GeMS</sub> (mas)</td><td>0.14</td><td>0.49</td><td>1.12, 0.13</td><td>1.94, 0.02</td></tr></tbody></table>

Notes. Luhman 16 AB parameters in our MCMC analysis (see Section 4).

<sup>a</sup> CDAR is chromatic differential atmospheric refraction.

<sup>b</sup> The sample fit is taken directly from the accepted zero-temperature walker MCMC chain. The sample 16-parameter fit is chosen to be within the 68.3% confidence interval. The median parameters are not representative of a single actual 16-parameter fit (i.e., the set of parameters accepted by the MCMC sampler).

Thus, we include the ρ for FORS2, ρ<sub>FORS2</sub>; the d for FORS2, d<sub>FORS2</sub>; and the ρ for GeMS, ρ<sub>GeMS</sub> in our MCMC analysis. The prior constraints on these parameters are shown in Table 4. We place a prior on ρ<sub>FORS2</sub> and d<sub>FORS2</sub> to be opposite signs, choosing (p > 0, d < 0), given that we define the atmospheric CDAR offset to always be positive and the correction to always be negative. Similarly, we place a positive prior constraint on ρ<sub>GeMS</sub> > 0. Finally, we constrain both ρ<sub>FORS2</sub> and d<sub>FORS2</sub> to be within reasonable physical ranges for total possible atmospheric CDAR and its correction.

Both ρ and d include dependencies on the spectrum of the target star. For the case of FORS2 imaging of Luhman 16 AB, the individual stellar components have different spectra, and a separate ρ and d would have to be fitted for each star individually. Rather than add these new terms to the global model fit, we explicitly measure the color index CI<sub>A,B</sub> of each star using MagE medium-resolution (λ/Δλ ≈ 4000) spectra and FIRE near-infrared (λ/Δλ ≈ 8000) spectra (Faherty et al. 2014, their Figure 1).

We write the CDAR parameters for Luhman B as a function of the ρ and d for Luhman A and the ratio of color indices CI<sub>A</sub>/CI<sub>B</sub>. The CDAR parameters ρ and d are expressed in units of mas (tan z)<sup>-1</sup>, where z is the zenith angle of the observation. We compute ρ and d for both Luhman 16 A and B as

$$ 
\begin{aligned}
\rho_A &= \rho \\
\rho_B &= \rho(CI_B/CI_A) \\
d_A &= d \\
d_B &= \rho(CI_B/CI_A).
\end{aligned} 
\qquad (8) 
$$

The color indices CI<sub>AB</sub> are defined as the shift in the position of the stellar photocenter due to the effect of CDAR. We model this

using formula (6) from Gubler & Tytler (1998) as follows:

$$ CI = \frac{\int_{\lambda_{\min}}^{\lambda_{\max}} A(\lambda)T(\lambda)I(\lambda)d\lambda}{\int_{\lambda_{\min}}^{\lambda_{\max}} T(\lambda)I(\lambda)d\lambda}, \qquad (9) $$

where T(λ) is the wavelength-dependent throughput of the atmosphere, telescope, and filters; I(λ) is the spectrum of the star; and A(λ) is the chromatic component of the refraction equation that is linear with the tangent of the zenith angle. As described in Gubler & Tytler (1998), the nonlinear component of the refraction B(λ) is three orders of magnitude smaller than the linear component and negligible for the zenith angles of concern here. We use the function refco within the Python implementation of SLALIB to calculate A(λ) for observatory parameters corresponding to the VLT, for an altitude of 2635 m, pressure of 800 mbar, temperature of 278 K, and relative humidity of 10%. Stellar spectra I(λ) in the I band are obtained from Faherty et al. (2014). T(λ) is calculated as the product of the I-band filter function, the QE curve of the detectors in FORS2, and the atmospheric transmission.<fnref n="16" /> The ratio of the color indices between Luhman 16 A and B in the FORS2 passband is (CI<sub>B</sub>/CI<sub>A</sub>) = 0.985763 mas (tan z)<sup>-1</sup>.

Similarly, we calculated the color indices for GeMS observations using the J-continuum narrowband filter using Equation (9), detector characteristics corresponding to GSAOI,

<footnote n="16">¹⁶ We use the Cerro Pachon model with 1.5 air masses and a water vapor column of 4.3 mm. See here: http://www.gemini.edu/sciops/telescopes-and-sites/observing-condition-constraints/ir-transmission-spectra.</footnote>
<table><caption>**Table 6**
Chromatic Differential Atmospheric Refraction Parameters</caption><thead><tr><th>Instrument</th><th>MJD</th><th>f<sub>3m</sub></th><th>tan z<sub>m</sub></th><th>tan z<sub>L,m</sub></th><th>cos γ</th><th>sin γ</th></tr></thead><tbody><tr><td rowspan="20">FORS2/VLT</td><td>56396.9951520</td><td>0.9791000</td><td>0.7546000</td><td>0.9286000</td><td>-0.6232000</td><td>-0.7816000</td></tr><tr><td>56402.0627130</td><td>0.9915000</td><td>0.5493000</td><td>0.6744000</td><td>-0.9973000</td><td>-0.0604000</td></tr><tr><td>56408.0405880</td><td>0.9809000</td><td>0.5521000</td><td>0.6636000</td><td>-0.9910000</td><td>-0.1282000</td></tr><tr><td>56418.0314980</td><td>0.9882000</td><td>0.5505000</td><td>0.6652000</td><td>-0.9941000</td><td>0.0876000</td></tr><tr><td>56424.0852670</td><td>0.9971000</td><td>0.7135000</td><td>0.8958000</td><td>-0.6789000</td><td>0.7336000</td></tr><tr><td>56432.9703430</td><td>1.0054000</td><td>0.5540000</td><td>0.6837000</td><td>-0.9855000</td><td>-0.1509000</td></tr><tr><td>56438.0316910</td><td>1.0052000</td><td>0.6561000</td><td>0.8121000</td><td>-0.7741000</td><td>0.6323000</td></tr><tr><td>56443.9929650</td><td>1.0101000</td><td>0.5929000</td><td>0.7314000</td><td>-0.8962000</td><td>0.4426000</td></tr><tr><td>56447.9978380</td><td>0.9892000</td><td>0.6350000</td><td>0.7693000</td><td>-0.8124000</td><td>0.5820000</td></tr><tr><td>56452.9719010</td><td>0.9842000</td><td>0.6028000</td><td>0.7154000</td><td>-0.8769000</td><td>0.4730000</td></tr><tr><td>56456.9999390</td><td>1.0011000</td><td>0.7344000</td><td>0.9070000</td><td>-0.6485000</td><td>0.7599000</td></tr><tr><td>56459.9735320</td><td>0.9914000</td><td>0.6628000</td><td>0.8087000</td><td>-0.7625000</td><td>0.6462000</td></tr><tr><td>56465.9621110</td><td>1.0012000</td><td>0.6804000</td><td>0.8611000</td><td>-0.7322000</td><td>0.6810000</td></tr><tr><td>56693.2499970</td><td>0.9981000</td><td>0.5626000</td><td>0.6850000</td><td>-0.9664000</td><td>-0.2378000</td></tr><tr><td>56702.1862400</td><td>0.9887000</td><td>0.6508000</td><td>0.7781000</td><td>-0.7880000</td><td>-0.6114000</td></tr><tr><td>56726.1444170</td><td>0.9856000</td><td>0.5884000</td><td>0.7062000</td><td>-0.9101000</td><td>-0.4007000</td></tr><tr><td>56735.1855440</td><td>0.9945000</td><td>0.5759000</td><td>0.6892000</td><td>-0.9348000</td><td>0.3368000</td></tr><tr><td>56746.1216870</td><td>0.9966000</td><td>0.5497000</td><td>0.6716000</td><td>-0.9968000</td><td>-0.0488000</td></tr><tr><td>56757.1135790</td><td>0.9882000</td><td>0.5595000</td><td>0.6750000</td><td>-0.9730000</td><td>0.2073000</td></tr><tr><td>56772.9970680</td><td>0.9914000</td><td>0.6361000</td><td>0.7689000</td><td>-0.8148000</td><td>-0.5743000</td></tr><tr><td>56782.9852030</td><td>0.9883000</td><td>0.5937000</td><td>0.7066000</td><td>-0.8974000</td><td>-0.4399000</td></tr><tr><td>56795.9827440</td><td>0.9915000</td><td>0.5522000</td><td>0.6678000</td><td>-0.9914000</td><td>-0.0771000</td></tr><tr><td>57013.3143040</td><td>0.9976000</td><td>0.7249000</td><td>0.8769000</td><td>-0.6650000</td><td>-0.7466000</td></tr><tr><td>57018.3266580</td><td>0.9896000</td><td>0.6318000</td><td>0.7469000</td><td>-0.8224000</td><td>-0.5647000</td></tr><tr><td>57025.3373360</td><td>0.9899000</td><td>0.5669000</td><td>0.6854000</td><td>-0.9567000</td><td>-0.2760000</td></tr><tr><td>57034.3479320</td><td>0.9853000</td><td>0.5530000</td><td>0.6668000</td><td>-0.9876000</td><td>0.1321000</td></tr><tr><td>57044.3577360</td><td>0.9873000</td><td>0.6153000</td><td>0.7421000</td><td>-0.8512000</td><td>0.5177000</td></tr><tr><td>57061.2201260</td><td>0.9642000</td><td>0.4498000</td><td>0.0000000</td><td>-0.9218000</td><td>-0.3867000</td></tr><tr><td>56759.1760450</td><td>0.9890000</td><td>0.6395000</td><td>0.0000000</td><td>-0.5429000</td><td>0.8396000</td></tr><tr><td>56760.1718400</td><td>0.9747000</td><td>0.6243000</td><td>0.0000000</td><td>-0.5675000</td><td>0.8224000</td></tr><tr><td>56804.0736350</td><td>0.9871000</td><td>0.7281000</td><td>0.0000000</td><td>-0.4182000</td><td>0.9083000</td></tr><tr><td>57000.3432270</td><td>0.9695000</td><td>0.6273000</td><td>0.0000000</td><td>-0.5403000</td><td>-0.8413000</td></tr><tr><td>57086.3244570</td><td>0.9653000</td><td>0.8635000</td><td>0.0000000</td><td>-0.2703000</td><td>0.9625000</td></tr></tbody></table>

Note. Chromatic differential atmospheric refraction parameters used to calculate parameters $f_{1,x,m}$, $f_{1,y,m}$, $f_{2,x,m}$, and $f_{2,y,m}$ as detailed in Section 4.1. The parameter $\tan z_{L,m}$ is zero for the GeMS data given that GeMS does not have an atmospheric differential corrector like VLT/FORS2.

and atmospheric parameters appropriate for the Gemini South Observatory. Due to the similarity of the J-continuum spectra of Luhman 16 A and B as measured in Faherty et al. (2014) and the narrow width of the filter band, the ratio in color indices for the GeMS observations evaluates to $0.9999$ mas ($\tan z)^{-1}$. We neglect this correction as it is significantly smaller than the error bars of the GeMS observations.

Following Lazorenko et al. (2009), Sahlmann et al. (2013), and Sahlmann & Lazorenko (2015), we express the CDAR correction in R.A. and decl. as a function of parameters $f_{1,x,m}$, $f_{1,y,m}$, $f_{2,x,m}$, and $f_{2,y,m}$. We use Equations (2) and (3) of Sahlmann et al. (2013) to compute $f_{1,x,m}$, $f_{1,y,m}$, $f_{2,x,m}$, and $f_{2,y,m}$ using parameters related to zenith angle, atmospheric pressure, atmospheric temperature, and parallactic angle from Table 6. We recompute and then subtract off CDAR from the coordinate positions at each step of our MCMC chains by fitting for $\rho$ and $d$ for both the FORS2 and GeMS data:

$$ \begin{aligned} \alpha'_{A,FORS2} &= \alpha_{A,FORS2} + (\rho_{A,FORS2} f_{1,x,m} + d_{A,FORS2} f_{2,x,m}) \\ \delta'_{A,FORS2} &= \delta_{A,FORS2} - (\rho_{A,FORS2} f_{1,y,m} + d_{A,FORS2} f_{2,y,m}) \\ \alpha'_{B,FORS2} &= \alpha_{B,FORS2} + (\rho_{B,FORS2} f_{1,x,m} + d_{B,FORS2} f_{2,x,m}) \\ \delta'_{B,FORS2} &= \delta_{B,FORS2} - (\rho_{B,FORS2} f_{1,y,m} + d_{B,FORS2} f_{2,y,m}), \end{aligned} \quad (10) $$

where $(\alpha'_{AB,FORS2}, \delta'_{AB,FORS2})$ are the CDAR-corrected coordinate positions for a given MCMC step and given $\rho_{AB,FORS2}d_{AB,FORS2}$ parameters, and $(\alpha_{AB,FORS2}, \delta_{AB,FORS2})$ is the uncorrected FORS2 coordinate positions. For the GeMS data, as mentioned above, we neglect the difference in color indices between Luhman 16 AB and use a single $\rho_{GeMS}$ parameter to model CDAR:

$$ \begin{aligned} \alpha'_{A,GeMS} &= \alpha_{A,GeMS} + (\rho_{GeMS} f_{1,x,m}) \\ \delta'_{A,GeMS} &= \delta_{A,GeMS} - (\rho_{GeMS} f_{1,y,m}) \\ \alpha'_{B,GeMS} &= \alpha_{B,GeMS} + (\rho_{GeMS} f_{1,x,m}) \\ \delta'_{B,GeMS} &= \delta_{B,GeMS} - (\rho_{GeMS} f_{1,y,m}), \end{aligned} \quad (11) $$

where $(\alpha'_{AB,GeMS}, \delta'_{AB,GeMS})$ are the CDAR-corrected positions for a given MCMC step and a given $\rho_{GeMS}$, and $(\alpha_{AB,GeMS}, \delta_{AB,GeMS})$ are the uncorrected GeMS positions. We do not subtract the effects of CDAR at the 1984 and 1999 epochs, which are negligible relative to the substantial astrometric uncertainty for those data.

The photocenter shift due to CDAR expected for the 1984 ESO Schmitt observations is <50 mas, which is 4× smaller than the uncertainties on the barycenter position due to other sources of error. Therefore, we neglect the CDAR correction in
the 1984 barycenter position as well as the correction on the relative separation. For the unresolved 1999 DENIS data, we do not attempt to constrain the barycenter position, and the change in the relative separation due to CDAR is far smaller than the error bars, so we again neglect CDAR.

## 4.2. Orbital Astrometry Model

At each step in our MCMC chains, we compute the log likelihood of our orbital astrometry model parameters $L_{rel}$. As the CDAR-subtracted Luhman 16 separations are themselves dependent on the three CDAR parameters, we recompute the CDAR magnitudes and thus the corrected separations at each step in our MCMC chains:

$$ 
\begin{aligned}
\Delta\alpha(t) &= (\alpha'_{B}(t) - \alpha'_{A}(t))\cos\frac{\delta'_{B}(t) + \delta'_{A}(t)}{2} \\
\Delta\delta(t) &= \delta'_{B}(t) - \delta'_{A}(t),
\end{aligned}
\quad (12) $$

where $\alpha'_{A,B}(t)$, $\delta'_{A,B}(t)$ are the CDAR-subtracted coordinates of Luhman 16 AB from Equations (10) and (11) at epoch $t$ and $\Delta\alpha(t)$, and $\Delta\delta(t)$ are the relative separations. Uncertainties on the separations $\sigma_{\Delta\alpha}$, $\sigma_{\Delta\delta}$ for GeMS and FORS2 wide-field imaging are computed as detailed in Sections 3.2 and 3.3, respectively.

We model orbital astrometry from Equation (12) with six Keplerian orbit parameters: semimajor axis $a$, eccentricity $e$, inclination $i$, longitude of ascending node $\Omega$, argument of periastron $\omega$, and time of periapsis $\tau$ (Green 1985). For our MCMC analysis, we adopt uniform priors in total mass $M_{tot} = M_A + M_B$, semimajor axis log $a$, eccentricity $e$, longitude of ascending node $\Omega$, argument of periastron $\omega$, time of periastron passage $\tau$, and inclination $\cos i$, as shown in Table 4.

For the 1999 DENIS data, we explicitly compute probability density functions (PDF) as a function of binary separation and position angle, as described in Section 3.7. The confidence intervals derived from these PDFs are shown in Figure 4. At each iteration of the MCMC run, the orbital model predicts the separation and P.A. at the 1999 epoch. We add the log PDF (log 1 - $p$) at the model separation and P.A. to the overall likelihood function $L$.

## 4.3. Relative Radial Velocity Model

The VLT/CRIRES relative radial velocities for Luhman 16 AB (see Section 2.5) provide information about the mutual orbit complementary to our orbital astrometry above. At each step in our MCMC chains, we compute the relative radial velocity model $V_A - V_B$ at each of three epochs using Equation (2.44), page 42 from Hilditch (2001). We add the likelihood function $L_{RV}$ to the overall likelihood $L$. Based on the evolution in the model posteriors before and after the addition of the radial velocity data, we assess that these data are primarily constraining the longitude of ascending node $\Omega$, the argument of periastron $\omega$, and the eccentricity $e$. These constraints disfavored sections of the mutual orbit parameter space that were not ruled out with the orbital astrometry alone.

## 4.4. Barycentric Motion Model

Luhman 16 AB was unresolved by DENIS in 1999, and the barycenter would be difficult to derive from the photocenter. Thus we exclude the DENIS I-band 1999 epoch from our barycentric motion model. At each step in our MCMC chains,

we compute the log likelihood of the barycenter model parameters $L_{abs}$. We compute the position of the barycenter (Sahlmann & Lazorenko 2015) at each epoch $t$ from the CDAR-subtracted coordinate positions (i.e., relative to background stars):

$$ 
\begin{aligned}
\alpha_{\gamma}(t) &= \left(\frac{1}{1+q}\right)(\alpha'_{A} + q\alpha'_{B}) \\
\delta_{\gamma}(t) &= \left(\frac{1}{1+q}\right)(\delta'_{A} + q\delta'_{B}),
\end{aligned}
\quad (13) $$

where $q = \frac{M_B}{M_A}$ is the mass ratio; $\alpha_{\gamma}(t)$, $\delta_{\gamma}(t)$ are the R.A. and decl. of the barycenter at epoch $t$; and $\alpha'_{A,B}$, $\delta'_{A,B}$ are the Luhman 16 AB positions from Equations (10) and (11). Our model for the position of the barycenter at epoch $t$ is given as

$$ 
\begin{aligned}
\alpha_{\gamma}(t) &= \alpha_0 + \mu_{\alpha}t + \bar{\omega}\Pi_{\alpha} \\
\delta_{\gamma}(t) &= \delta_0 + \mu_{\delta}t + \bar{\omega}\Pi_{\delta},
\end{aligned}
\quad (14) $$

where $\alpha_{\gamma}(t)$, $\delta_{\gamma}(t)$ are the coordinates of the barycenter; $\bar{\omega}_{rel}$ is the parallax factor; and $\mu_{\alpha}$, $\mu_{\delta}$ are the proper motions in R.A. and decl. To maximize $L_{abs}$, we minimize the residuals between Equations (13) and (14). Thus, from Equations (13) and (14), we model six parameters associated with the location of the barycenter—and importantly, the barycenter position changes significantly as a function of mass ratio. For our MCMC analysis, we adopt uniform priors for the reference epoch positions $\alpha_0$, $\delta_0$, proper motion $\mu_{\alpha}$ and $\mu_{\delta}$, parallax $\bar{\omega}_{rel}$, and mass ratio $q$, as shown in Table 4.

### 4.4.1. Derivation of Barycentric Position Errors

Unlike the relative separation uncertainties, the uncertainties on the barycenter position have a dependence on the mass ratio $q$. Propagated uncertainties on the barycenter $\sigma_{\alpha_{\gamma}}(t)$, $\sigma_{\delta_{\gamma}}(t)$ are computed for the FORS2, GeMS, and ESO 1984 plate as

$$ 
\begin{aligned}
\sigma_{\alpha_{\gamma}}^{2}(t) &= \left(\frac{1}{1+q}\right)^{2}\sigma_{\alpha_{A}}(t)^{2} + \left(\frac{q}{1+q}\right)^{2}\sigma_{\alpha_{B}}(t)^{2} \\
&\quad + 2\frac{q}{1+q}\frac{1}{1+q}\sigma_{\alpha_{A},\alpha_{B}}(t) \\
\sigma_{\delta_{\gamma}}^{2}(t) &= \left(\frac{1}{1+q}\right)^{2}\sigma_{\delta_{A}}(t)^{2} + \left(\frac{q}{1+q}\right)^{2}\sigma_{\delta_{B}}(t)^{2} \\
&\quad + 2\frac{q}{1+q}\frac{1}{1+q}\sigma_{\delta_{A},\delta_{B}}(t),
\end{aligned}
\quad (15) $$

where $\sigma_{\alpha_{A}}(t)$, $\sigma_{\alpha_{B}}(t)$ are the position uncertainties for Luhman 16 AB in the ICRS system at epoch $t$ as derived in Section 3.8, $q = \frac{M_B}{M_A}$ is the mass ratio, and $\sigma_{\alpha_{A},\alpha_{B}}(t)$, $\sigma_{\delta_{A},\delta_{B}}(t)$ are the correlated error between the positions of Luhman 16 AB. Given that we have only one to three frames per epoch for our GeMS astrometry, we do not include correlated error in the GeMS barycenter uncertainty calculation.

We include the correlated error in our computation of the FORS2 barycenter uncertainty (Equation (15)). The correlated errors are calculated for each epoch from the frame-by-frame
<figure><img src="image_9.png" /><figcaption>Figure 6. Posterior probability distributions for the orbital astrometry parameters (see Section 4 and Table 5). There are clear degeneracies between several of the parameters (ω, e, a, Ω, and M<sub>tot</sub>). However, the confidence intervals on the total mass M<sub>tot</sub> given the degeneracies remain small at ≈3.25%.</figcaption></figure>

positions of both stars:

$$
\begin{align}
\sigma_{X_A, X_B}(t) &= \frac{1}{1-N} \sum_{i=1}^{i=N} [(X_{A,i} - \bar{X}_A)(X_{B,i} - \bar{X}_B)] \nonumber \\
\sigma_{\alpha_A, \alpha_B}(t) &= 125\sigma_{X_A, X_B}(t) \nonumber \\
\sigma_{Y_A, Y_B}(t) &= \frac{1}{1-N} \sum_{i=1}^{i=N} [(Y_{A,i} - \bar{Y}_A)(Y_{B,i} - \bar{Y}_B)] \nonumber \\
\sigma_{\delta_A, \delta_B}(t) &= 125\sigma_{Y_A, Y_B}(t), \tag{16}
\end{align}
$$

where σ<sub>X<sub>A</sub>,x<sub>B</sub></sub>(t), σ<sub>Y<sub>A</sub>,y<sub>B</sub></sub>(t) are the correlated errors in pixels for
epoch t using N frames; X<sub>A</sub>, i, Y<sub>A,i</sub> are the positions of Luhman
16 A at frame i in units of pixels; and X̄<sub>A</sub> is the mean position
over N frames. We convert the correlated errors σ<sub>X<sub>A</sub>,x<sub>B</sub></sub>(t),

σ<sub>Y<sub>A</sub>, y<sub>B</sub></sub>(t) pixels to σ<sub>α<sub>A</sub>, α<sub>B</sub></sub>(t), σ<sub>δ<sub>A</sub>, δ<sub>B</sub></sub>(t) arcseconds using the
FORS2 pixel scale of 125 mas pixel<sup>-1</sup>.

4.5. Running the MCMC

To determine the posterior probability distributions for the
orbit and barycenter astrometry parameters, we performed a
parallel-tempered MCMC analysis (Earl & Deem 2005) to
maximize the likelihood function L = L<sub>rel</sub> + L<sub>abs</sub> + L<sub>RV</sub> of the
model given the data (Gregory 2011; Kalas et al. 2013). We
compute the CDAR-subtracted coordinate positions as detailed
in Section 4.1. We compute the likelihood for the orbit
parameters L<sub>rel</sub> as described in Section 4.2, the likelihood for
barycentric orbit parameters L<sub>abs</sub> as described in Section 4.4,
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Posterior probability distributions for the barycentric astrometry and instrument parameters (see Section 4 and Table 5). There are obvious degeneracies between several of the parameters, in particular between proper motion $\mu_\alpha$, $\mu_\delta$ and the mass ratio $q$. Nevertheless, our mass ratio $0.82^{+0.03}_{-0.03}$ is in good agreement with $0.78 \pm 0.10$ from Sahlmann & Lazorenko (2015).</figcaption></figure>

and the likelihood for the relative radial velocity data as described in Section 4.3.

We minimize the negative log likelihood function using the Python implementation of the affine-invariant parallel-tempered ensemble sampler `PTsampler` from the Python module `emcee v2.1.0`<fnref n="17" /> (Foreman-Mackey et al. 2013). We use 300 walkers and $2 \times 10^6$ steps, discarding the first 20% of steps as the burn-in for each walker. We set temperature $T = 14$ to maximize the search of parameter space for secondary minima. We visually inspect the highest temperature chain, finding that the walkers cover the entirety of the prior uniformly distributed parameter ranges in Table 4. This means the

convergence of our MCMC analysis to the best-fit parameters and confidence intervals in Table 5 are not the result of a local minimum. We also visually inspect the lowest temperature chain, finding the walkers to stabilize over the $2 \times 10^6$ steps, which confirms the MCMC run convergence to our maximum-likelihood parameters. The Gelman–Rubin statistics, which check for convergence of MCMC chains run in parallel, were $\approx 1.0004$ for all 16 parameters, indicating successful convergence (Gelman & Rubin 1992).

## 5. Results

The posterior probability distributions for the individual masses $M_A$, $M_B$, mass ratio $q \equiv \frac{M_B}{M_A}$, and sum of masses

<footnote n="17">¹⁷ http://dan.iel.fm/emcee/current/</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Posterior probability distributions from our MCMC analysis for the cross section between our orbital astrometry parameters (X axis) or barycentric astrometry parameters (Y axis). There are no obvious degeneracies between the parameters derived from Luhman 16 separation data and the barycentric motion parameters.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Posterior probability distributions for the mass ratio, total mass, and individual masses (see Section 4 and Table 5). The individual masses $M_A$ and $M_B$ are largely decorrelated. We marginalize over the posterior probability distributions to obtain a well-constrained total mass $M_{\text{tot}} = 62.0_{-1.7}^{+2.0} M_{\odot}$ with a 3.25% (1σ) uncertainty.</figcaption></figure>

$M_A + M_B$ are shown in Figure 9. The posterior probability distributions are smooth, indicating no multiple local minima. As shown in Figures 6–8, we compute posterior probability distributions of each of the 16 model parameters by marginalizing over the other 15 parameters from the zero-temperature walker distributions. Our 1σ (68.3% confidence) and 2σ (95.4%) credible intervals, derived from our MCMC posterior distributions, are shown in Table 5. Our individual masses derived for Luhman 16 AB are $34.2_{-1.1}^{+1.3} M_J$ and $27.9_{-1.0}^{+1.1} M_J$, respectively.

As shown in Figure 6, we find slightly elliptical 2D posterior PDFs when plotting orbit parameters $\Omega$ versus $a$, $e$ versus $i$, $M_{\text{tot}}$ versus $e$, and $M_{\text{tot}}$ versus $i$—those parameters have small degeneracies. Nevertheless, we have a strong constraint on the total mass $M_{\text{tot}} = 62.0_{-1.7}^{+2.0}$. We rule out alternative orbit models with larger total mass $M_{\text{tot}} > 65 M_J$ by starting multiple MCMC runs, as detailed in Section 4, with initial guesses of $M_{\text{tot}} = 70–120 M_J$. We find in all of these runs that the zero-temperature walkers for $M_{\text{tot}}$ steadily move toward lower masses (~62 $M_J$) over the duration of the run.

As shown in Figure 7, proper motions $\mu_\alpha$, $\mu_\delta$ and reference epoch R.A. $\alpha_0$ and decl. $\delta_0$ are degenerate with mass ratio $q$. Nevertheless, we constrain the mass ratio to $q = 0.82_{-0.03}^{+0.03}$, in good agreement with the mass ratio of $q = 0.78 \pm 0.10$ from Sahlmann & Lazorenko (2015), which uses only the VLT/FORS2 2013–2014 astrometry included in this analysis. Our relative parallax of $\omega_{\text{rel}} = 0.50101_{-0.00053}^{+0.00052}$ arcsec is $\sim 1.5$σ off from the relative parallax of $\omega_{\text{rel}} = 0.50023 \pm 0.00011$ arcsec measured by Sahlmann & Lazorenko (2015). This

discrepancy likely results from the addition of new 2015 FORS2 data, 2014–2015 GeMS data, and explicit astrometric analysis of the 1984 ESO plate in our analysis. Our median proper motions from the posterior probability density functions of $\mu_\alpha = -2.76216_{-0.00243}^{+0.00243}$ arcsec per year and $\mu_\delta = 0.35779_{-0.00344}^{+0.00344}$ arcsec per year are in good agreement with Sahlmann & Lazorenko (2015) proper motions of $\mu_\alpha = 2.754 \pm 0.00625$ arcsec per year and $\mu_\delta = 0.3587 \pm 0.0095$ arcsec per year.

For earlier runs of the MCMC analysis, we placed much wider constraints of $\pm 200$ mas on CDAR parameters $\rho_{\text{FORS2}}$ and $d_{\text{FORS2}}$, finding the two parameters to be tightly correlated (see Figure 3 from Lazorenko et al. 2007, which uses different units). This linear relationship between $\rho_{\text{FORS2}}$ and $d_{\text{FORS2}}$ is shown in our Figure 7. This trend is expected as the LADC parameter $d_{\text{FORS2}}$ models a near cancellation of atmospheric DAR, as modeled by $\rho_{\text{FORS2}}$. We narrowed the priors for these parameters to physically meaningful values as calculated from Equation (3). These limits are given in Table 4 and reflected in Figure 7. Sahlmann & Lazorenko (2015) find CDAR parameters $\rho_{\text{FORS2}} = 35.23$ mas and $d_{\text{FORS2}} = -48.71$ mas, which agree with the linear relationship between those parameters shown in Figure 7. Furthermore, we find $\rho_{\text{GeMS}}$ to be nearly zero, which agrees with our expectation that the usage of a narrowband $J_{\text{cont}}$ filter minimizes the effect (Cameron et al. 2009).

In this study, we explored the possibility that the orbit model and the barycentric motion model could have degenerate parameters, given that both are derived from the same data sets.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Panel (A): $10^4$ randomly sampled orbits from our posterior probability distributions (gray lines) using an MCMC analysis technique (see Section 4.2). The uncertainties for the FORS2 (red circles) and GeMS (teal squares) orbital astrometry are smaller than the symbol sizes. The $1\sigma$, $2\sigma$, and $3\sigma$ contours for our DENIS 1999 are shown in dark blue, blue, and light blue, respectively. The dotted lines correspond to the 180-degree degeneracy in the binary fitting routine (see Section 3.7). The ESO 1984 epoch is shown as a green downward-facing triangle. Panels (B) and (C): mutual orbit of Luhman 16 AB in delta R.A. and decl. with the same models shown in panel (A). Panels (D)–(F): The “data-model” residuals for our best-fit orbit, excluding the 1984 and 1999 epochs for clarity. We find a small ≲0.8 mas systematic offset between our FORS2 and GeMS astrometry and our maximum-likelihood orbit, likely due to small systematic errors in the CDAR parameters (see Section 4.1).</figcaption></figure>

From Figure 8, however, we find no degeneracies between the
orbit parameters and the barycenter astrometry parameters—all
of the 2D posterior PDFs are circularly symmetric. Figure 9
shows that the posterior PDFs for *q* and *M*<sub>tot</sub> are circularly
symmetric, indicating little degeneracy between the mass of the
primary *q* and *M*<sub>tot</sub>. We do find a small degeneracy between the
mass of the primary *M*<sub>A</sub> and secondary *M*<sub>B</sub>, which is
incorporated into our confidence intervals. These masses are
within range of the model-based masses of 21–62 *M*<sub>J</sub> estimated
from the observed luminosities and effective temperatures of

Luhman 16 AB and an upper limit of <62 *M*<sub>J</sub> (for both stars)
from lithium detection by Faherty et al. (2014) and Lodieu
et al. (2015).

In panels (B) and (C) of Figure 10, we compare $10^4$
randomly sampled orbits from our zero-temperature walker
posterior probability distributions to our CDAR-subtracted
separation data. The standard deviation of the absolute value of
the FORS2 and GeMS residuals (“data-model”) is ≈0.8 mas
for the separation fits, as shown in panels (D)–(F) of Figure 10.
From panel (D), we see a small systematic offset of ≈0.5 mas
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Shown are $10^4$ randomly sampled relative radial velocity orbits from our posterior probability MCMC chains (same as Figure 10), compared to the VLT/CRIRES relative radial velocities $V_A - V_B$, as detailed in Sections 2.5 and 4.3. The relative radial velocities are in good agreement with our maximum-likelihood orbit.</figcaption></figure>

between the mean residuals for FORS2 orbital astrometry (red points) and the GeMS astrometry (teal points). This small offset may be due to small systematic errors in computation of the FORS2 CDAR parameters; however, the overall uncertainty on total mass (3.25%) remains small. In multiple MCMC runs with and without the 1984 ESO and 1999 DENIS astrometry, we find the combination of both of these epochs plays a crucial role in constraining a mass–eccentricity–period degeneracy observed when using only the GeMS and FORS2 data to fit the orbital parameters of Luhman 16 AB. When using only the GeMS and FORS2 astrometry, we find that dynamical mass could only be constrained to 60–90 $M_J$. Observations in 2017–2018 correspond to the time of closest separation for Luhman 16, which could render the DENIS 1999 data unnecessary. Similarly, observations in 2020–2022 would constrain the most uncertain parts of the orbit, possibly rendering the ESO 1984 data unnecessary.

In Figure 11, we compare the same $10^4$ randomly sampled orbits to the CRIRES relative radial velocity points, finding good agreement. The relative radial velocities had little effect on further constraining the total mass beyond the constraints provided by the orbital astrometry. Using the RV data, we rule out orbits with argument of periastron $\omega > 150^\circ$ and longitude of ascending node $\Omega < 129.5$. Those ranges of parameter space could not be ruled out by the astrometry alone.

In Figure 12, we compare randomly sampled barycenter models from our posterior probability distributions to the observed barycenter motion for Luhman 16 AB. We find residuals of 1.41 mas in R.A. and 0.99 in decl. for our barycentric astrometry with no statistically significant systematic trends.

## 5.1. Comparison to Bedin et al. (2017) Results

Bedin et al. (2017) present HST Wide Field Camera 3 (WFC3) measurements of the separation of the individual

<table><caption>**Table 7**

Masses and Luminosities of Luhman 16 AB</caption><thead><tr><th></th><th>Mass (M<sub>J</sub>)</th><th>log L<sub>☉</sub><sup>a</sup></th><th>log L<sub>☉</sub><sup>b</sup></th></tr></thead><tbody><tr><td>Luhman 16 A</td><td>34.2 ± 1.2</td><td>-4.67 ± 0.04</td><td>-4.66 ± 0.08</td></tr><tr><td>Luhman 16 B</td><td>27.9 ± 1.0</td><td>-4.71 ± 0.10</td><td>-4.68 ± 0.13</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Faherty et al. (2014).

<sup>b</sup> Lodieu et al. (2015).

components of Luhman 16 AB. The authors also present derivations of the orbit and barycentric motion parameters of the system. Due to the limited temporal coverage of the data relative to the full orbital period, the constraints on the orbit, total system mass, and individual masses are relatively weak. Here we compare our measured separation data for the 2014.5–2015.5 epoch, in which both of our observational programs have coverage. We also compare our derived model parameters.

Our best-fitting model of the Luhman 16 AB orbit compares favorably to the raw separation measurements from Bedin et al. (2017). The standard deviations of the residuals (i.e., our minimum $\chi^2$ model minus separation data from Bedin et al. 2017 over the longest common epoch of 2014.5–2015.5) are 0.318 mas in X and 0.431 mas in Y, which compare favorably with the 0.3–0.4 measurement errors given in Table 8 of that work.

All six barycentric motion parameters (parallax $\pi_{rel}$, proper motion $\mu_{\alpha,\delta}$, position $\alpha_0$, $\delta_0$, and mass ratio $q$) measured by Bedin et al. (2017) fall within our 68.3% confidence intervals. The widths of our 68.3% confidence intervals are generally similar to the $1\sigma$ error bars quoted by Bedin et al. (2017) for these parameters.

For the orbital parameters, given our more complete sampling of the orbital period, the widths of our confidence intervals are generally improved over the Bedin et al. (2017) $1\sigma$ error bars. Our measurements of the total mass $M_{tot} = M_A + M_B$, semimajor axis $a$, eccentricity $e$, longitude of ascending node $\Omega$, argument of periastron $\omega$, time of periastron passage $\tau$, inclination $i$, and period $P$ are more precise by factors of 2–30. Although our measurements of the longitude of ascending node $\Omega$ and argument of periastron $\omega$ are different from those of Bedin et al. (2017) by a $\pi$ phase shift and the inclination $i$ has an additional negative sign, our measurements of all orbital parameters are consistent to within the $1\sigma$ error bars quoted in Bedin et al. (2017), with the exception of the eccentricity $e$. Our measurement of this parameter has a 68.3% confidence interval of [0.33, 0.37], and the value from Bedin et al. (2017) has confidence limits of [0.399, 0.527], so the difference is $1.76\sigma$ discrepant from our most likely value of 0.35. Given that this is the only parameter with a discrepancy above $1\sigma$, and the probability of such a difference occurring by chance alone is not small, we conclude that all of our orbital parameters are consistent to within the errors.

## 5.2. Evolutionary Models

Our masses and the bolometric luminosities of Luhman 16 AB due to Lodieu et al. (2015) and Faherty et al. (2014) are shown in Table 7 and plotted in Figure 13 (we adopt the Lodieu et al. luminosities). We compare our observations to cloudy ($f_{sed} = 2$), cloud-free, and hybrid ($f_{sed} = 2$ cloudy-to-cloud-free) evolutionary models from Saumon & Marley (2008; SM08).
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Panel (A): Our FORS2 and GeMS barycentric astrometry (red and teal circles and squares) are accurate enough to distinguish between different mass ratios. The black and gray lines are the barycenter locations corresponding to MCMC chains within $1\sigma$ and $2\sigma$ confidence intervals of our maximum-likelihood mass ratio $q = 0.82_{-0.03}^{+0.03}$. While the true location of the barycenter relative to either star does not change, the barycenter will vary as a function of the mass ratio (Equation (13), Section 4.4). Panels (B) and (C): “data-model” residuals for our maximum-likelihood barycentric astrometry, excluding the 1984 epoch for clarity. We find there is no systematic error between our GeMS (teal squares) and FORS2 (red circles) astrometry and our astrometric model.</figcaption></figure>

The evolutionary models use solar metallicity ($\frac{M}{H} = 0$). The uncertainties on the bolometric luminosities include the known few percent photometric variability of both brown dwarfs (see Faherty et al. 2014; Lodieu et al. 2015, for details). We find the hybrid models give a common age of $\approx 600–800$ Myr for both Luhman 16 AB to within $1\sigma$ uncertainty. Similar to the findings of Dupuy et al. (2015) for L–T transition binary SDSSJ 1052 AB, we find a mass–luminosity relationship across the L–T transition consistent with zero. We also find that SM08 cloudless models cannot reproduce our masses and the observed luminosities of Luhman 16 AB to within $1\sigma$. Stronger conclusions to differentiate between the

SM08 hybrid and SM08 cloudy models require greater precision in the bolometric luminosities for Luhman 16 AB.

## 6. Conclusion

In this work, we derive the individual masses for the closest known brown dwarf binary, Luhman 16 AB. We draw upon archival observations from the European Southern Observatory Schmidt Telescope, the Deep Near-Infrared Survey of the Southern Sky (DENIS), public FORS2 data on the Very Large Telescope (VLT), and new astrometry from the Gemini South Multiconjugate Adaptive Optics System (GeMS). We use an
<figure><img src="image_17.png" /><figcaption>**Figure 13.** The hybrid (cloudy-to-cloud-free) or cloudy models of Saumon & Marley (2008) best reproduce the flat mass–luminosity relationship across the L–T transition observed from directly measured masses and bolometric luminosities of brown dwarf binaries Luhman 16 AB (Lodieu et al. 2015 for luminosities; this work for masses) and SDSS J105213.51+442255.7AB (Dupuy et al. 2015) to within uncertainty. The uncertainties on the bolometric luminosity due to the few percent inherent photometric variability of Luhman 16 AB are included in the error bars shown here.</figcaption></figure>

MCMC analysis to simultaneously fit the barycentric space motion of the binary as well as the mutual Keplerian orbit. We see little evidence for degeneracies between parameters derived from the separation data and those describing the barycentric motion. Our individual masses derived for Luhman 16 AB of $34.2_{-1.1}^{+1.3} M_J$ and $27.9_{-1.0}^{+1.1} M_J$ respectively represent the second derivation of masses for the L and T components of a dwarf binary with uncertainties <10%.

Luhman 16 AB mass measurements will also play an important role in constraining the mass–luminosity relationship across the L/T transition. Future work may find that our individual mass measurements of Luhman 16 AB complement the observations by Dupuy et al. (2015) of L–T transition binary SDSS J105213.51+442255.7, who found that hybrid cloudy-to-cloud-free brown dwarf evolutionary models of Saumon & Marley (2008) pass the coevality test of binarity. On the other hand, cloudless models of Saumon & Marley (2008) are unable to reproduce the observed masses and luminosities of SDSS J1052+44 while maintaining coevality, that is, derive ages for both stars from the masses and luminosities that are within $1\sigma$ uncertainty. Future work includes comparing the masses and luminosities of Luhman 16 AB to every set of the most recent cloudy and cloud-free

evolutionary models, to test whether the models can reproduce the observed masses and luminosities.

This work is based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), Ministerio da Ciencia, Tecnologia e Inovacao (Brazil), and Ministerio de Ciencia, Tecnologia e Innovacion Productiva (Argentina). Data are acquired through the Gemini Science Archive and processed using the Gemini IRAF package. This work has made use of data from the European Space Agency (ESA) mission GAIA (http://www.cosmos.esa.int/gaia), processed by the GAIA Data Processing and Analysis Consortium (DPAC; http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the GAIA Multilateral Agreement. This research is based on the data obtained from the ESO Science Archive Facility under program IDs 291.C-5004 and 593.C-0314. This research has made use of data obtained from the SuperCOSMOS Science
Archive, prepared and hosted by the Wide Field Astronomy Unit, Institute for Astronomy, University of Edinburgh, which is funded by the UK Science and Technology Facilities Council. This work was partially funded by NASA/NEXSS NNX15AD95G. This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France. This work is performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 with document release number LLNL-JRNL-701044. The posterior probability distributions in Figures 6-9 are made with the python function "corner" (Foreman-Mackey et al. 2016).

## ORCID iDs

E. Victor Garcia https://orcid.org/0000-0001-7779-8573

S. Mark Ammons https://orcid.org/0000-0001-5172-7902

Maissa Salama https://orcid.org/0000-0002-5082-6332

Jeffrey Chilcote https://orcid.org/0000-0001-6305-7272

Quinn Konopacky https://orcid.org/0000-0002-9936-6285

Jessica R. Lu https://orcid.org/0000-0001-9611-0009

Bruce Macintosh https://orcid.org/0000-0003-1212-7538

Christian Marois https://orcid.org/0000-0002-4164-4182

Eric Nielsen https://orcid.org/0000-0001-6975-9056

Robert J. De Rosa https://orcid.org/0000-0002-4918-0247

Gaetano Sivo https://orcid.org/0000-0002-1296-5868

## References

Ammons, S. M., Bendek, E. A., & Guyon, O. 2011, Proc. SPIE, 8151, 81510T
Ammons, S. M., Bendek, E. A., Guyon, O., Macintosh, B., & Savransky, D. 2012, Proc. SPIE, 8447, 84470P

Anderson, D. R., Collier Cameron, A., Hellier, C., et al. 2011, ApJL, 726, L19

Appenzeller, I., Fricke, K., Fürtig, W., et al. 1998, Msngr, 94, 1

Auddy, S., Basu, S., & Valluri, S. R. 2016, AdAst, arXiv:1607.04338

Avila, G., Rupprecht, G., & Beckers, J. M. 1997, Proc. SPIE, 2871, 1135

Bean, J. L., Miller-Ricci Kempton, E., & Homeier, D. 2010, Natur, 468, 669

Bedin, L. R., Pourbaix, D., Apai, D., et al. 2017, MNRAS, 470, 1140

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

Blake, C. H., Charbonneau, D., White, R. J., Marley, M. S., & Saumon, D. 2007, ApJ, 666, 1198

Boffin, H. M. J., Pourbaix, D., Mužić, K., et al. 2014, A&A, 561, L4

Bouchy, F., Bonomo, A. S., Santerne, A., et al. 2011a, A&A, 533, A83

Bouchy, F., Deleuil, M., Guillot, T., et al. 2011b, A&A, 525, A68

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

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

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

Cameron, P. B., Britton, M. C., & Kulkarni, S. R. 2009, AJ, 137, 83

Carrasco, E. R., Edwards, M. L., McGregor, P. J., et al. 2012, Proc. SPIE, 8447, 84470N

Crepp, J. R., Johnson, J. A., Fischer, D. A., et al. 2012, ApJ, 751, 97

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

Deleuil, M., Deeg, H. J., Alonso, R., et al. 2008, A&A, 491, 889

Díaz, R. F., Damiani, C., Deleuil, M., et al. 2013, A&A, 551, L9

Dieterich, S. B., Henry, T. J., Jao, W.-C., et al. 2014, AJ, 147, 94

Diolaiti, E., Bendinelli, O., Bonaccini, D., et al. 2000, Proc. SPIE, 4007, 879

Dupuy, T. J., Liu, M. C., Leggett, S. K., et al. 2015, ApJ, 805, 56

Earl, D. J., & Deem, M. W. 2005, PCCP, 7, 3910

Epchtein, N., Deul, E., Derriere, S., et al. 1999, A&A, 349, 236

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

Ford, H. C., Bartko, F., Bely, P. Y., et al. 1998, Proc. SPIE, 3356, 234

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

Foreman-Mackey, D., Vousden, W., Price-Whelan, A., et al. 2016, corner.py: corner.py, v2.0.0, doi:10.5281/zenodo.53155

Freudling, W., Romaniello, M., Bramich, D. M., et al. 2013, A&A, 559, A96

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016a, A&A, 595, A2

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016b, A&A, 595, A1

Gelman, A., & Rubin, D. B. 1992, StaSc, 7, 457

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

Green, R. M. 1985, Spherical Astronomy (Cambridge: Cambridge Univ. Press)

Gregory, P. C. 2011, MNRAS, 415, 2523

Gubler, J., & Tytler, D. 1998, PASP, 110, 738

Hambly, N. C., MacGillivray, H. T., Read, M. A., et al. 2001, MNRAS, 326, 1279

Hilditch, R. W. 2001, An Introduction to Close Binary Stars (Cambridge: Cambridge Univ. Press)

Ireland, M. J., Kraus, A., Martinache, F., Lloyd, J. P., & Tuthill, P. G. 2008, ApJ, 678, 463

Johnson, J. A., Apps, K., Gazak, J. Z., et al. 2011, ApJ, 730, 79

Kaeufl, H.-U., Ballester, P., Biereichel, P., et al. 2004, Proc. SPIE, 5492, 1218

Kalas, P., Graham, J. R., Fitzgerald, M. P., & Clampin, M. 2013, ApJ, 775, 56

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

Konopacky, Q. M., Ghez, A. M., Barman, T. S., et al. 2010, ApJ, 711, 1087

Lazorenko, P. F., Mayor, M., Dominik, M., et al. 2007, A&A, 471, 1057

Lazorenko, P. F., Mayor, M., Dominik, M., et al. 2009, A&A, 505, 903

Lodieu, N., Zapatero Osorio, M. R., Rebolo, R., et al. 2015, A&A, 581, A73

Lu, J. R. 2008, PhD thesis, Univ. California, Los Angeles

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

Mamajek, E. E. 2013, On the Nearby Binary Brown Dwarf WISE J104915.57-531906.1 (Luhman 16), arXiv:1303.5345

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

Markwardt, C. B. 2009, in ASP Conf. Ser. 411, Astronomical Data Analysis Software and Systems XVIII, ed. D. A. Bohlender, D. Durand, & P. Dowler (San Francisco, CA: ASP), 251

McGregor, P., Hart, J., Stevanovic, D., et al. 2004, Proc. SPIE, 5492, 1033

Montet, B. T., Johnson, J. A., Fortney, J. J., & Desert, J.-M. 2016, ApJL, 822, L6

Montet, B. T., Johnson, J. A., Muirhead, P. S., et al. 2015, ApJ, 800, 134

Neichel, B., Rigaut, F., Vidal, F., et al. 2014, MNRAS, 440, 1002

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical Recipes in C++: The Art of Scientific Computing (Cambridge: Cambridge Univ. Press)

Rigaut, F., Neichel, B., Boccas, M., et al. 2014, MNRAS, 437, 2361

Sahlmann, J., & Lazorenko, P. F. 2015, MNRAS, 453, L103

Sahlmann, J., Lazorenko, P. F., Ségransan, D., et al. 2013, A&A, 556 A133

Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327

Siverd, R. J., Beatty, T. G., Pepper, J., et al. 2012, ApJ, 761, 123

Stassun, K. G., Mathieu, R. D., & Valenti, J. A. 2006, Natur, 440, 311

Street, R. A., Fulton, B. J., Scholz, A., et al. 2015, ApJ, 812, 161

Tody, D. 1986, Proc. SPIE, 627, 733