# Characterizing 51 Eri b from 1 to 5 µm: A Partly Cloudy Exoplanet

Abhijith Rajan, Julien Rameau, Robert J. de Rosa, Mark S. Marley, James R.

Graham, Bruce Macintosh, Christian Marois, Caroline Morley, Jennifer

Patience, Laurent Pueyo, et al.

► To cite this version:

Abhijith Rajan, Julien Rameau, Robert J. de Rosa, Mark S. Marley, James R. Graham, et al. Characterizing 51 Eri b from 1 to 5 µm: A Partly Cloudy Exoplanet. The Astronomical Journal, 2017, 154, 10.3847/1538-3881/aa74db. insu-03692496

HAL Id: insu-03692496

https://insu.hal.science/insu-03692496v1

Submitted on 3 Sep 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>

Characterizing 51 Eri b from 1 to 5 µm: A Partly Cloudy Exoplanet

Abhijith Rajan¹, Julien Rameau², Robert J. De Rosa³, Mark S. Marley⁴, James R. Graham³, Bruce Macintosh⁵, Christian Marois⁶,⁷,
Caroline Morley⁸, Jennifer Patience¹, Laurent Pueyo⁹, Didier Saumon¹⁰, Kimberly Ward-Duong¹, S. Mark Ammons¹¹,
Pauline Arriaga¹², Vanessa P. Bailey⁵, Travis Barman¹³, Joanna Bulger¹⁴, Adam S. Burrows¹⁵, Jeffrey Chilcote¹⁶, Tara Cotten¹⁷,
Ian Czekala⁵, Rene Doyon², Gaspard Duchêne³,¹⁸, Thomas M. Esposito³, Michael P. Fitzgerald¹², Katherine B. Follette⁵,
Jonathan J. Fortney¹⁹, Stephen J. Goodsell²⁰, Alexandra Z. Greenbaum²¹, Pascale Hibon²², Li-Wei Hung¹², Patrick Ingraham²³,
Mara Johnson-Groh⁷, Paul Kalas³,²⁴, Quinn Konopacky²⁵, David Lafrenière², James E. Larkin¹², Jérôme Maire²⁵, Franck Marchis²⁴,
Stanimir Metchev²⁶,²⁷, Maxwell A. Millar-Blanchaer²⁸,<fnref n="35" />, Katie M. Morzinski²⁹, Eric L. Nielsen⁵,²⁴, Rebecca Oppenheimer³⁰,
David Palmer¹¹, Rahul I. Patel³¹, Marshall Perrin⁹, Lisa Poyneer¹¹, Fredrik T. Rantakyrö²², Jean-Baptiste Ruffio⁵,
Dmitry Savransky³², Adam C. Schneider¹, Anand Sivaramakrishnan⁹, Inseok Song¹⁷, Rémi Soummer⁹, Sandrine Thomas²³,
Gautam Vasisht²⁸, J. Kent Wallace²⁸, Jason J. Wang³, Sloane Wiktorowicz³³, and Schuyler Wolff³⁴

¹School of Earth and Space Exploration, Arizona State University, P.O. Box 871404, Tempe, AZ 85287, USA; arajan6@asu.edu

²Institut de Recherche sur les Exoplanètes, Département de Physique, Université de Montréal, Montréal QC H3C 3J7, Canada

³Astronomy Department, University of California, Berkeley; Berkeley CA 94720, USA

⁴NASA Ames Research Center, Mountain View, CA 94035, USA

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

⁶National Research Council of Canada Herzberg, 5071 West Saanich Rd, Victoria, BC V9E 2E7, Canada
⁷University of Victoria, 3800 Finnerty Rd, Victoria, BC V8P 5C2, Canada

⁸Department of Astronomy, Harvard University, Cambridge, MA 02138, USA

⁹Space Telescope Science Institute, Baltimore, MD 21218, USA

¹⁰Los Alamos National Laboratory, P.O. Box 1663, MS F663, Los Alamos, NM 87545, USA

¹¹Lawrence Livermore National Laboratory, Livermore, CA 94551, USA

¹²Department of Physics & Astronomy, University of California, Los Angeles, CA 90095, USA

¹³Lunar and Planetary Laboratory, University of Arizona, Tucson AZ 85721, USA

¹⁴Subaru Telescope, NAOJ, 650 North A'ohoku Place, Hilo, HI 96720, USA

¹⁵Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

¹⁶Dunlap Institute for Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada

¹⁷Department of Physics and Astronomy, University of Georgia, Athens, GA 30602, USA

¹⁸Univ. Grenoble Alpes/CNRS, IPAG, F-38000 Grenoble, France

¹⁹Department of Astronomy, UC Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

²⁰Gemini Observatory, 670 N. A'ohoku Place, Hilo, HI 96720, USA

²¹Department of Astronomy, University of Michigan, Ann Arbor MI 48109, USA

²²Gemini Observatory, Casilla 603, La Serena, Chile

²³Large Synoptic Survey Telescope, 950N Cherry Av, Tucson, AZ 85719, USA

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

²⁵Center for Astrophysics and Space Science, University of California San Diego, La Jolla, CA 92093, USA

²⁶Department of Physics and Astronomy, Centre for Planetary Science and Exploration, The University of Western Ontario, London, ON N6A 3K7, Canada

²⁷Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794-3800, USA

²⁸Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91125, USA

²⁹Steward Observatory, University of Arizona, Tucson AZ 85721, USA

³⁰Department of Astrophysics, American Museum of Natural History, New York, NY 10024, USA

³¹Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

³²Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY 14853, USA
³³The Aerospace Corporation, El Segundo, CA 90245, USA

³⁴Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA

Received 2017 March 24; revised 2017 May 8; accepted 2017 May 8; published 2017 June 16

Abstract

We present spectrophotometry spanning 1–5 µm of 51 Eridani b, a 2–10 $M_{\text{Jup}}$ planet discovered by the Gemini Planet Imager Exoplanet Survey. In this study, we present new K1 (1.90–2.19 µm) and K2 (2.10–2.40 µm) spectra taken with the Gemini Planet Imager as well as an updated $L_P$ (3.76 µm) and new $M_S$ (4.67 µm) photometry from the NIRC2 Narrow camera. The new data were combined with J (1.13–1.35 µm) and H (1.50–1.80 µm) spectra from the discovery epoch with the goal of better characterizing the planet properties. The 51 Eri b photometry is redder than field brown dwarfs as well as known young T-dwarfs with similar spectral type (between T4 and T8), and we propose that 51 Eri b might be in the process of undergoing the transition from L-type to T-type. We used two complementary atmosphere model grids including either deep iron/silicate clouds or sulfide/salt clouds in the photosphere, spanning a range of cloud properties, including fully cloudy, cloud-free, and patchy/intermediate- opacity clouds. The model fits suggest that 51 Eri b has an effective temperature ranging between 605 and 737 K, a solar metallicity, and a surface gravity of $\log(g) = 3.5–4.0$ dex, and the atmosphere requires a patchy cloud atmosphere to model the spectral energy distribution (SED). From the model atmospheres, we infer a luminosity

<footnote n="35">³⁵ NASA Hubble Fellow.</footnote>
for the planet of −5.83 to −5.93 (log L/L$_{⊙}$), leaving 51 Eri b in the unique position of being one of the only directly imaged planets consistent with having formed via a cold-start scenario. Comparisons of the planet SED against warm-start models indicate that the planet luminosity is best reproduced by a planet formed via core accretion with a core mass between 15 and 127 M$_{⊙}$.

*Key words:* instrumentation: adaptive optics – planets and satellites: atmospheres – planets and satellites: composition – planets and satellites: gaseous planets – stars: individual (51 Eridani)

*Supporting material:* data behind figures

# 1. Introduction

Until recently, most of the imaged planetary-mass companions detected were typically orbiting their parent star at large orbital separations, >30 au. However, new instrumentation with second-generation adaptive optics such as the Gemini Planet Imager (GPI; Macintosh et al. 2014) and Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE; Beuzit et al. 2008) are now routinely obtaining deep contrasts (>10$^{5}$–10$^{6}$) in the inner arcsecond (5–30 au). The recent detection of new companions (Macintosh et al. 2015; Konopacky et al. 2016; Wagner et al. 2016; Milli et al. 2017) and debris disks (Currie et al. 2015; Millar-Blanchaer et al. 2016; Wahhaj et al. 2016; Bonnefoy et al. 2017) showcase the advances made by these next-generation Adaptive Optics (AO) systems. Direct imaging, unlike nondirect methods such as radial velocity and transits, measures light from companions directly, which permits measuring the atmospheric spectrum, with the caveat that the final calibration is dependent on a complete understanding of the stellar properties. These new AO instruments combine excellent image stability and high throughput with Integral Field Unit (IFU) spectrographs, enabling the measurement of a spectrum of the planet in the near-infrared (IR) wavelength range. Combining the near-IR spectra with mid-IR photometry from instruments such as Keck/NIRC2, MagAO/Clio, or LBT/LMIRCam provides valuable constraints on the effective temperature and nonequilibrium chemistry when undertaking comprehensive modeling of the exoplanet spectral energy distribution (SED).

In this study we focus on the planetary companion 51 Eridani b (51 Eri b; Macintosh et al. 2015). 51 Eri b is the first planet discovered by the Gemini Planet Imager Exoplanet Survey (GPIES), a survey targeting 600 young and nearby stars using GPI to search for exoplanets. The planet orbits 51 Eri A, a young F0IV star that is part of the β Pic moving group (Zuckerman et al. 2001). In this study, we adopt an age of 26 ± 3 Myr for the β Pic moving group (Nielsen et al. 2016). However, the age of the group is a topic of considerable debate and has been revised several times, for example, 21 ± 4 (Binks & Jeffries 2014), 23 ± 3 (Mamajek & Bell 2014), 20 ± 6 (Macintosh et al. 2015), and 24 ± 3 (Bell et al. 2015). The primary is part of a hierarchical triple with two M-star companions, GJ 3305AB, separated from the primary by ~2000 au (Feigelson et al. 2006; Kasper et al. 2007; Montet et al. 2015). 51 Eri A is known to have an IR excess, and a debris disk was detected in the Herschel Space Observatory 70 and 100 µm bands with a very low IR luminosity of L$_{IR}$/L$_{*}$ = 2 × 10$^{-6}$ and a lower limit on the inner radius of 82 au (Riviere-Marichalar et al. 2014), as well as a detection at 24 µm with the Spitzer Space Telescope (Rebull et al. 2008). The debris disk was not detected in Macintosh et al. (2015), which, given the low fractional luminosity, would be extremely challenging. The analysis of the atmosphere of 51 Eri b by Macintosh et al. (2015) was based on GPI JH spectra

(1.1–1.8 µm) and Keck L$_{P}$ photometry (3.76 µm), using two different model atmosphere grids to estimate the planet properties. While the models agreed on the temperature and luminosity, they were highly discrepant in terms of best fitting the surface gravity with one grid, suggesting low surface gravity and youth, while the other required a high surface gravity and an old planet. Similarly, one grid best fit the atmosphere when using a linear combination of cloudy and clear models, while the other best fit the data with a clear atmosphere. These discrepancies indicate that more data are required to fully constrain the planet parameters.

In this paper, we present new observations and a revised data analysis that can be used to discriminate between some of the disagreements. In Section 2, we present the first K1 (1.90–2.19 µm) and K2 (2.10–2.40 µm) spectra of the planet taken with GPI. We also present updated L$_{P}$ photometry and new observations of the planet in the M$_{S}$ band (4.67 µm). In Section 3, we present new near-IR photometry of the star and revise the stellar SED used in the rest of the analysis. In Section 4, we examine the near- and mid-IR photometry of 51 Eri b in relation to that of other field and young brown dwarfs through the brown dwarf color–magnitude diagram (CMD). We also compare the near-IR spectrum of 51 Eri b to field brown dwarfs and planetary-mass companions to estimate the best-fitting spectral type of the planet. Finally, in Section 5 we model the planet SED using two different grids spanning effective temperatures from 450 to 1000 K with deep iron/silicate clouds or sulfide/salt clouds. The 1–5 µm SED in combination with these two model grids will help refine the planet properties and clarify whether the atmosphere is best fit by clouds or not.

# 2. Observations and Data Reduction

## 2.1. GPI K1 and K2

51 Eri b was observed with the Integral Field Spectrograph (IFS) of GPI through the K1 filter on 2015 November 06 UT and 2016 January 28 and through the K2 filter on 2015 December 18 UT (see Table 1). Standard procedures, namely using an argon arc lamp, were used to correct the data for instrumental flexure. To maximize the parallactic rotation for angular differential imaging (ADI; Marois et al. 2006), the observations were centered on the meridian passage. All of the GPI data sets underwent the same initial data processing steps using the GPI Data Reduction Pipeline v1.3.0 (DRP; Perrin et al. 2014). The processing steps included dark-current subtraction and bad-pixel identification and interpolation; this is followed by compensating for instrument flexure using the argon arc spectrum (Wolff et al. 2014). Following this step, the microspectra are extracted to generate the IFS data cubes (Maire et al. 2014). During the process of generating the 3D (x, y, λ) cubes, the microspectra data are resampled to λ/δλ = 65 and 75 at K1 and K2, respectively, after which they are interpolated to a common wavelength scale and
<table><caption>**Table 1**
Observations of 51 Eri b</caption><thead><tr><th>Date</th><th>Instrument</th><th>Filter</th><th>Total Int. time (minutes)</th><th>Field Rot. (°)</th><th>Averaged air mass</th><th>Averaged DIMM seeing (as)</th><th>Averaged MASS τ<sub>0</sub> (ms)</th></tr></thead><tbody><tr><td>2015 Jan 30</td><td>GS/GPI</td><td>J<sup>a</sup></td><td>70</td><td>23.8</td><td>1.15</td><td>0.52</td><td>3.26</td></tr><tr><td>2014 Dec 18</td><td>GS/GPI</td><td>H<sup>a</sup></td><td>38</td><td>37.7</td><td>1.14</td><td>...</td><td>...</td></tr><tr><td>2015 Nov 06</td><td>GS/GPI</td><td>K1<sup>b</sup></td><td>55</td><td>30.5</td><td>1.17</td><td>0.38</td><td>1.56</td></tr><tr><td>2015 Dec 18</td><td>GS/GPI</td><td>K2<sup>b</sup></td><td>103</td><td>71.7</td><td>1.22</td><td>0.69</td><td>0.94</td></tr><tr><td>2016 Jan 28</td><td>GS/GPI</td><td>K1<sup>b</sup></td><td>97</td><td>55.5</td><td>1.15</td><td>0.86</td><td>4.40</td></tr><tr><td>2015 Oct 27</td><td>Keck/NIRC2</td><td>L<sub>P</sub><sup>b</sup></td><td>100</td><td>74.2</td><td>1.10</td><td>...</td><td>...</td></tr><tr><td>2016 Jan 02</td><td>Keck/NIRC2</td><td>M<sub>S</sub><sup>b</sup></td><td>139</td><td>115.7</td><td>1.18</td><td>...</td><td>...</td></tr><tr><td>2016 Jan 21</td><td>Keck/NIRC2</td><td>M<sub>S</sub><sup>b</sup></td><td>174</td><td>116.0</td><td>1.21</td><td>...</td><td>...</td></tr><tr><td>2016 Feb 04</td><td>Keck/NIRC2</td><td>M<sub>S</sub><sup>b</sup></td><td>148</td><td>101.4</td><td>1.21</td><td>...</td><td>...</td></tr><tr><td>2016 Feb 05</td><td>Keck/NIRC2</td><td>M<sub>S</sub><sup>b</sup></td><td>142</td><td>102.1</td><td>1.21</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.**

corrected for geometric distortion (Konopacky et al. 2014). The data cubes are then aligned to a common center calculated using the four satellite spots (Wang et al. 2014). The satellite spots are copies of the occulted central star, generated by the use of a regular square grid printed on the apodizer in the pupil plane (Marois et al. 2006; Sivaramakrishnan & Oppenheimer 2006; Macintosh et al. 2014). The satellite spots also help convert the photometry from contrast units to flux units. No background subtraction was performed since the following steps of high-pass filtering and point-spread function (PSF) subtraction efficiently remove this low-frequency component.

Further steps to remove quasi-static speckles and large-scale structures were executed outside the DRP. Each data cube was filtered using an unsharp mask with a box width of 11 pixels. The four satellite spots were then extracted from each wavelength slice and averaged over time to obtain templates of the star PSF. The Linear Optimized Combination of Images algorithm (LOCI; Lafrenière et al. 2007) was used to suppress the speckle field in each frame using a combination of aggressive parameters: $dr = 5$ px, $N_A = 200$ PSF FWHM, $g = 0.5$, and $N_δ = 0.5–0.75$ FWHM for the three data sets, where $dr$ is the radial width of the optimization zone, $N_A$ is the number of PSF FWHM that can be included in the zone, $g$ is the ratio of the azimuthal and radial widths of the optimization zone, and $N_δ$ defines the maximum separation of a potential astrophysical source in FWHM between the target and the reference PSF. The residual image of each wavelength slice was built from a trimmed (10%) temporal average of the sequence.

Final K1 and K2 broadband images were created using a weighted mean of the residual wavelength frames according to the spectrum of the planet, examples of which can be found in Figure 1. These broadband images were used to extract the astrometry of the planet in each data set thanks to a higher signal-to-noise ratio (S/N) than in individual frames. To do so, a negative template PSF was injected into the raw data at the estimated position and flux of the planet before applying LOCI and reduced using the same matrix coefficients as the original reduction (Marois et al. 2010). The process was iterated over these three parameters (x position, y position, flux) with the amoeba simplex optimization (Nelder & Mead 1965) until the integration squared pixel noise in a wedge of $2 \times 2$ FWHM was minimized. The best-fit position was then used to extract

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Final PSF-subtracted images of 51 Eri b. Top: LOCI-reduced GPIES images at the K1 (2016 January 28, left) and K2 bands (2015 December 18, right). Bottom: pyKLIP-reduced NIRC2 images, smoothed with a box of width 2 pixels, at $L_P$ (2015 October 27, left), and a combined image of all four $M_S$ data sets (right). The images are scaled linearly, but are different in each panel in order to saturate the core of the planet PSF.</figcaption></figure>

the contrast of the planet in each data set. The same procedure was executed in the noncollapsed wavelength residual images, varying only the flux of the negative template PSF and keeping the position fixed to prevent the algorithm from catching nearby brighter residual speckles in the lower S/N spectral slices. To measure uncertainties, we injected the template PSF with the measured planet contrast into each data cube at the same separation and 20 different position angles. We measured the fake signal with the same extraction procedure. The contrasts measured in the 2015 November 06 and 2016 January 28 K1 data sets agreed within the uncertainties, the latter having significantly better S/N, and were combined with the weighted mean to provide the final planet contrasts.

<footnote>a Macintosh et al. (2015).</footnote>

<footnote>b This work.</footnote>
### 2.1.1. Spectral Covariances

The estimation of a directly imaged planet's properties from its measured spectrum is complicated by the fact that spectral covariances are present within the extracted spectra. In the GPI data, these are caused by the residual speckle noise in the final PSF-subtracted image and the oversampling of the individual microspectra during the initial data reduction process. Atmosphere modeling without properly accounting for these covariances can lead to biased results. We present the derivation of the correlation using the parameterization of Greco & Brandt (2016) in the Appendix.

We use the spectral covariance when carrying out a comparison of the planet spectrophotometry against other field and young dwarfs as well as during model fitting. The covariance helps correctly account for the correlation in the spectra while also increasing the importance of the photometry, and thus the use of the covariance tends to move the best fits toward cooler temperatures when compared to using the variance directly.

### 2.2. Keck $L_p$

We observed the 51 Eri system on 2015 October 27 in the $L_p$ filter with the NIRC2 camera (McLean & Sprayberry 2003) at the Keck II observatory (Program ID U055N2). The observations were taken in ADI mode, starting ~1 hr prior to meridian crossing to maximize the field-of-view rotation. The target was observed for ~3 hr total, with 100 minutes of on-source integration. The observations were acquired using the 400 mas focal plane mask and the circular undersized "in-circle" cold stop. To calibrate the planet brightness, unsaturated observations of the star were taken at the end of the observing sequence. The images were dark and flat-field corrected. We used the $K_s$-band lamp flats to build the flat-field and masked hot and bad pixels. As these observations were taken after the 2015 April servicing of NIRC2, the geometric distortion was corrected using the solution presented in Service et al. (2016, updating the original Yelda et al. 2010 solution), with an updated plate scale of 9.971 ± 0.004 mas pixel⁻¹ and the offset angle $\beta$ (0.262 ± 0.020) that is required when calculating the position angle prior to rotating the images to put north up (Yelda et al. 2010). Postprocessing of the data was carried out using the Python version of the Karhunen–Loève (KL) Image Projection algorithm (KLIP; Amara & Quanz 2012; Soummer et al. 2012), *pyKLIP* (Wang et al. 2015). As part of this study, we included an NIRC2 module in the *pyKLIP* code base that is publicly available for users.<fnref n="36" /> The algorithm accepts aligned images and performs PSF subtraction using KLIP where the image can be divided into sections both radially and azimuthally. Aside from the choice of zones, two main parameters were adjusted, the number of modes used in the KL transform and an exclusion criterion for reference PSFs, similar to $N_\delta$ mentioned above, that determine the number of pixels an astrophysical source would move due to the rotation of the reference stack. We carried out a parameter search where the four parameters mentioned were varied to optimize the S/N in the planet signal. The planet photometry was estimated using the method described above for the K1 and K2 filters, using a negative template PSF. The $L_p$ magnitude contrast for the star/ planet is 11.58 ± 0.15 mag, which agrees very well with the

photometry in the original epoch, 11.62 ± 0.17 mag. The weighted mean of both measurements is used in the rest of the analysis.

### 2.3. Keck $M_S$

Observations of 51 Eri b were taken in the $M_S$-band filter over four separate half nights on 2016 January 02 and 21 and 2016 February 04 and 05 with the Keck/NIRC2 Narrow camera. The details of the observations are presented in Table 1. Each night the target was observed for a period of ~6 hr, as part of two separate NASA and UC Keck observing programs (Program IDs N179N2, U117N2). The data were obtained in ADI mode, with the field of view rotating at the sidereal rate. To reduce the effects of persistence and enable accurate thermal background correction, the star was nodded across the detector in four large dithers centered in each quadrant of the detector. Furthermore, to prevent saturation of the detector by the thermal background, the exposures were limited to 0.3 s with 200 coadds, without using an occulting spot. The images were dark and flat-field corrected with twilight sky flats, followed by hot and bad pixel correction. As with the $L_p$ data, the solution provided by Service et al. (2016) was used to correct the NIRC2 Narrow camera geometric distortion. Finally, all of the images were rotated to put north up.

An additional step required for the $M_S$-band data that is not as critical for the other data sets is the background subtraction. Since the thermal background at 5 µm is large and highly time variable, rather than median combine or high-pass filter to remove the background, we adopted the least-squares sky-subtraction algorithm proposed in Galicher et al. (2011). For each point in the dither pattern, the algorithm uses the images where the star is in one of the other three positions to construct a reference library. We used a ring centered on the star to estimate the thermal background in each image, with an inner annulus of 24 pixels and an outer annulus of 240 pixels. The final calibration step involved aligning the background-corrected PSFs. Since the core of the PSF is saturated in the data, we aligned the data using two different methods: (1) fitting a 2D Gaussian to the wings of the stellar PSF to estimate the center of the star and then shifting the PSF to a predetermined pixel value to align all the images and (2) using the rotation symmetry of the PSF using the method described in Morzinski et al. (2015). To compare the two methods, we calculated the residuals between images aligned using the methods and compared the noise in the residuals and found them to be similar, so we chose to go with the 2D Gaussian, which is computationally faster.

The procedure used for the PSF subtraction for the $M_S$ data was similar to that for the $L_p$ data. The planet is not detected in each of the individual half-night data sets, requiring a combination of all four half nights to increase the S/N to detect the planet flux. To correctly combine the planet flux across the multiple epochs, we adjusted the PA to account for the astrophysical motion of the planet around the star, for which we used the best-fitting orbit presented in De Rosa et al. (2015). In the month between the first and last data sets, the planet rotated ~0.48 or ~0.4 pixel, which is a sufficiently large correction that it must be included in the data reduction. Each night's data were reduced individually to generate 603 PSF-subtracted images. These images were then combined by dividing each image into 13 annuli, which were combined

<footnote n="36">³⁶ https://bitbucket.org/pyKLIP/pyklip</footnote>
using a weighted mean, where the weights are the inverse variance in each annulus. As seen in Figure 1, we detect the planet signal at ~2–3σ. To confirm that we are detecting the planet, we rotated the data to match the PA value of the $L_P$ epoch to find that the flux peak in the $M_S$ band matches the location of the planet in $L_P$. We measured a star-to-planet contrast of 11.5 mag using the same procedure as described for the $L_P$ data. We injected 25 fake PSFs that were scaled to match the contrast measured for the planet and detected the fakes at the same contrast as the planet. The final magnitude of the star-to-planet contrast in the $M_S$ is 11.5 ± 0.5 mag.

## 3. Results

To estimate the stellar parameters of 51 Eri A, Macintosh et al. (2015) made use of Two Micron All-Sky Survey photometry (2MASS; Cutri et al. 2003; Skrutskie et al. 2006). However, the J- and H-band photometry for the star are flagged as “E,” indicating that the photometry is of the poorest quality and potentially unreliable (as compared to an “A” flag for the K-band photometry). Further, the study used photometry taken with the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) in the W1 filter ($\lambda_{eff} = 3.35$ µm, $\Delta\lambda = 1.11$ µm) as an approximation for the $L_P$-band magnitude of the primary star. The photometry for 51 Eri A in W1, from the AllWISE catalog (Cutri et al. 2013), has large errors and contributes to more than half the error budget of the final planet photometry. In this study, we thus chose to reobserve the star in the JHK$_S$ filters and fit all the available photometry to estimate the photometry in filters where no calibrated stellar data exist.

### 3.1. Revised Stellar Photometry at J, H, and K$_S$

The 2MASS near-IR colors of 51 Eri A were compared to empirical colors for young F0 stars taken from Kenyon & Hartmann (1995), where an F0IV star should have $J - H = 0.13$ mag and $H - K = 0.03$ mag. The colors of 51 Eri A estimated using the 2MASS photometry are however discrepant, with $J - H = -0.03 \pm 0.08$ and $H - K = 0.23 \pm 0.08$ mag. The discrepant near-IR colors combined with poor-quality flags suggest that the published photometry is potentially incorrect.

We observed the star 51 Eri A using the 6.5 m MMT on Mt. Hopkins with the ARIES instrument (McCarthy et al. 1998) on 2016 February 28 UT under photometric conditions. We obtained data in the MKO JHK$_S$ broadband filters (Tokunaga et al. 2002), for a total of 3.4 minutes in each filter. To flux calibrate these observations, we observed a photometric standard star at a similar air mass as 51 Eri A, HR 1552 (Carter 1990). The raw images for both targets were processed through a standard near-IR reduction pipeline, performing dark-current subtraction, flat-field calibration, and bad-pixel correction. Aperture photometry was performed on both targets, with the curve of growth used to select an aperture that minimized the error on the measured flux. The measured brightness of 51 Eri A is presented in Table 2.

Converting the MKO $K_S$-band measurement into the 2MASS system using empirical relations<fnref n="37" /> yields $K_{S,2MASS} = 4.551 \pm 0.032$ mag, which is within $1\sigma$ of the published 2MASS photometry. Furthermore, the $J - H$ and $H - K$ colors estimated from the revised photometry are $0.128 \pm 0.001$ mag.

<table><caption>**Table 2**
System Properties</caption><thead><tr><th>Property</th><th>51 Eri A</th><th>51 Eri b</th></tr></thead><tbody><tr><td>Distance (pc)</td><td>...</td><td>29.43 ± 0.29<sup>a</sup></td></tr><tr><td>Age (Myr)</td><td>...</td><td>26 ± 3<sup>b</sup></td></tr><tr><td>Spectral Type</td><td>F0IV</td><td>T6.5 ± 1.5</td></tr><tr><td>log(L/L<sub>☉</sub>)</td><td>0.85<sup>+0.06</sup><sub>-0.07</sub><sup>c</sup></td><td>-5.83<sup>+0.15</sup><sub>-0.12</sub> to -5.93<sup>+0.19</sup><sub>-0.14</sub><sup>d</sup></td></tr><tr><td>T<sub>eff</sub></td><td>7331 ± 30 K<sup>e</sup></td><td>605–737 K<sup>d</sup></td></tr><tr><td>log g</td><td>3.95 ± 0.04<sup>e</sup></td><td>3.5–4.0<sup>d</sup></td></tr><tr><td>J<sub>MKO</sub></td><td>4.690 ± 0.020<sup>d</sup></td><td>19.04 ± 0.40<sup>d,f</sup></td></tr><tr><td>H<sub>MKO</sub></td><td>4.562 ± 0.031<sup>d</sup></td><td>18.99 ± 0.21<sup>d</sup></td></tr><tr><td>K<sub>S,MKO</sub></td><td>4.546 ± 0.024<sup>d</sup></td><td>18.49 ± 0.19<sup>d</sup></td></tr><tr><td>K<sub>MKO</sub></td><td>4.600 ± 0.024<sup>e</sup></td><td>18.67 ± 0.19<sup>d</sup></td></tr><tr><td>L<sub>P</sub></td><td>4.604 ± 0.014<sup>e</sup></td><td>16.20 ± 0.11<sup>d,g</sup></td></tr><tr><td>M<sub>S</sub></td><td>4.602 ± 0.014<sup>e</sup></td><td>16.1 ± 0.5<sup>d</sup></td></tr></tbody></table>

**Notes.**

a *Hipparcos* catalog (van Leeuwen 2007).

b Nielsen et al. (2016).

c Macintosh et al. (2015) using hot-start predictions.

d This work.

e Stellar photometry estimated using SED fit.

f Distance modulus = 2.34 ± 0.02 mag.

g Weighted mean of the two $L_P$ observations.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison of the JH spectra of 51 Eri b using the literature 2MASS values against the new photometry measured in this study. The updated photometry increases the planet flux by ~10% in J and ~15% in the H band. The updated stellar photometry is used in the remainder of this study. However, the final stellar spectrum used to correct the planet spectrum does not depend on individual filter photometry, as in Macintosh et al. (2015) and shown in this plot, but is generated by modeling the full stellar SED prior to converting the planet spectra from contrast to flux units.</figcaption></figure>

0.037 mag and $0.016 \pm 0.039$ mag, which are consistent with the empirical expectations.

The published 51 Eri b spectrum in Macintosh et al. (2015) was calibrated using the Pickles stellar models (Pickles 1998) to estimate the spectrum of the primary, where each band was scaled using the published 2MASS photometry. In Figure 2 we present a comparison between the published spectrum and one scaled using the new MKO photometry, using the same stellar models. The revised photometry scales the planet spectrum higher by ~10% in the J band and ~15% in the H band, which is significant given the high S/N of the H-band data.

### 3.2. Fitting the SED of 51 Eri A

To mitigate the effects of incorrect photometry, rather than scale the spectrum in pieces using the relevant broadband

<footnote n="37">³⁷ http://www.astro.caltech.edu/~jmc/2mass/v3/transformations/</footnote>
<figure><img src="image_5.png" /></figure>

**Figure 3.** Top panel: photometry of 51 Eri A from the literature, and from the results presented in this study (solid symbols). One hundred models were randomly selected from the MCMC search and are plotted (translucent black curves). For each model, the synthetic magnitude was calculated for each filter. The median value for each filter is shown as an open square. The 2MASS photometry points are plotted to illustrate the offset relative to the new MKO measurements and are not included in the fit. For the plotted Geneva photometry, we computed the Geneva V-band photometry using the best-fit spectrum and then used the color relations to calculate the photometry in the remaining filters. Bottom panel: the fractional residuals relative to the median model. The data used to create this figure are available.

photometry, we decided to fit the full SED of 51 Eri A using literature photometry and colors, including Geneva U, B<sub>I</sub>, B<sub>2</sub>, V<sub>I</sub>, V, G (Rufener & Nicolet 1988), *Tycho2* B<sub>T</sub>, V<sub>T</sub>/ Hipparcos H<sub>P</sub> (ESA 1997; Høg et al. 2000), MKO JHK<sub>S</sub> (this work), and WISE W1, W2 (Cutri et al. 2013) measurements. We made use of the Geneva color relations as constraints to the full SED fit since the published Geneva V magnitude, which anchors the colors to estimate the remaining photometry, appears to be offset by ~5% when compared to the *Tycho2* photometry. The WISE W2 photometry was corrected using the Cotten & Song (2016) relation for bright stars. We combine the photometry with model stellar atmospheres from the BT-NextGen grid<sup>38</sup> (Allard et al. 2012), and we estimated the stellar spectrum using a five-parameter Markov chain Monte Carlo (MCMC) grid search. The best-fit atmosphere was found with T<sub>eff</sub> = 7331 ± 30 K, log g = 3.95 ± 0.04, [M/H] = −0.12 ± 0.06, and a stellar radius R = 1.45 ± 0.02 R<sub>☉</sub> (assuming a parallax of 33.98 ± 0.34 mas; van Leeuwen 2007). No correction for extinction is performed as the extinction in the direction of 51 Eri is negligible (A<sub>V</sub> = 0.00; Guarinos 1992). These values are consistent with previous literature estimates (e.g., Koleva & Vazdekis 2012). The final SED of 51 Eri A is shown in Figure 3, which highlights the significantly discrepant 2MASS JH-band photometry that was used previously to calibrate the spectrum of 51 Eri b. We extracted MKO K and NIRC2 L<sub>P</sub> and M<sub>S</sub> photometry from the SED fit using the filter response functions presented in Tokunaga et al. (2002); see Table 2.

### 3.2.1. Confirming the Stellar L<sub>P</sub> Photometry

51 Eri b emits a substantial amount of flux in the mid-IR, and L<sub>P</sub> photometry in Macintosh et al. (2015) was used to constrain the effective temperature of the planet. There exists no L<sub>P</sub> flux measurement for the star, so they used the W1 magnitude reported in the AllWISE catalog (W1 = 4.543 ± 0.210; Cutri et al. 2013) and assumed a color of W1 − L<sub>P</sub> = 0

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

**Figure 4.** Final spectral energy distribution of the directly imaged exoplanet 51 Eri b. The new K1 and K2 GPI spectra along with the updated L<sub>P</sub> and new M<sub>S</sub> photometry are shown with red squares. The GPI *J* and *H* spectra, updated to account for the revised stellar flux, from the discovery paper (Macintosh et al. 2015), are plotted with blue circles. The filter extent is shown with the horizontal line over each band. To reduce crowding in the spectra, the errors for one out of every two data points are plotted. The data used to create this figure are available.

based on the F0IV spectral type of 51 Eri (Abt & Morrell 1995). The L<sub>P</sub> photometry we estimated via the SED fits for 51 Eri is L<sub>P</sub> = 4.604 ± 0.014 mag, which is consistent with the value reported in Macintosh et al. (2015, 4.52 ± 0.21 mag) but with significantly smaller uncertainties.

As a final check for consistency, the 2MASS K<sub>S</sub> magnitude of 51 Eri (K<sub>S,2MASS</sub> = 4.537 ± 0.024) was used instead as a starting point. The K<sub>S</sub> − L<sub>P</sub> color for early F-type dwarfs and subgiants was estimated by folding model stellar spectra (7200 ≤ T<sub>eff</sub>/K ≤ 7400, 4.0 ≤ log g ≤ 4.5, [M/H] = 0) from the BT-SETTL model grid through the relative spectral response of the 2MASS K<sub>S</sub> (Cohen et al. 2003) and NIRC2 L<sub>P</sub> filters. Over this range of temperatures and surface gravities, the color was calculated as K<sub>S</sub> − L<sub>P</sub> = −0.001 ± 0.001. In order to realistically assess the uncertainties on this color, the near- to thermal-IR spectra of F-type dwarfs and subgiants within the IRTF library (Rayner et al. 2009) were processed in the same fashion, resulting in a K<sub>S</sub> − L<sub>P</sub> = 0.014 ± 0.055. A color of K<sub>S</sub> − L<sub>P</sub> = −0.001 ± 0.055 was adopted based on the color calculated from the model grid and the uncertainty calculated from the empirical IRTF spectra. This color, combined with the K<sub>S,2MASS</sub> magnitude of 51 Eri, gives an L<sub>P</sub> apparent magnitude of 4.538 ± 0.060. The estimates for the stellar L<sub>P</sub> magnitude are within 1σ of each other, so we adopt the value derived from the SED fit: L<sub>P</sub> = 4.604 ± 0.014 mag.

### 3.3. 51 Eri b SED

We present the final SED of the planet 51 Eri b in Figure 4 and use it to analyze the system properties in the following sections. Using the stellar SED estimated earlier, we have updated the *J* and *H* spectra that were published in Macintosh et al. (2015). In Table 2, we present the properties of the system, including updated MKO JHK and NIRC2 L<sub>P</sub>M<sub>S</sub> photometry for both the star and the planet. A future study will refine the orbital solution presented in De Rosa et al. (2015).

## 4. Analysis

### 4.1. Comparison against Field Brown Dwarfs

We plot a series of CMDs for ultracool objects in Figure 5 and compare the photometry of field M, L, and T dwarfs and young

<footnote n="38">38 https://phoenix.ens-lyon.fr/Grids/BT-NextGen/SPECTRA/</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The brown dwarf and imaged exoplanet color–magnitude diagram. 51 Eri b is shown with the red star. The colors of 51 Eri b place it among late-T dwarfs, where it is redder than most comparable-temperature brown dwarfs, likely indicative of greater cloud opacity in the atmosphere. The photometry for the field M-dwarfs (black circles), young M-dwarfs (blue triangles), field L-dwarfs (dark gray circles), young L-dwarfs (light blue triangles), and T-dwarfs (light gray circles) is taken from the compilation of Dupuy & Liu (2012) and Liu et al. (2016). We used a linear fit to convert WISE W1 photometry to L<sub>p</sub>, similar to what was done in Macintosh et al. (2015). The photometry for the directly imaged planets and young brown dwarfs were taken from Males et al. (2014), Bonnefoy et al. (2014), Bowler et al. (2017), Marois et al. (2010), Chauvin et al. (2005), Rameau et al. (2013), Naud et al. (2014), Leggett et al. (2007), Delorme et al. (2017), Goldman et al. (2010), Janson et al. (2011), and Kuzuhara et al. (2013).</figcaption></figure>

brown dwarfs and imaged companions to that of 51 Eri b (red
star). The colors of 51 Eri b seem to match the phase space of the
late-T dwarfs. To classify the spectral type of 51 Eri b, we do a
chi-squared comparison of the GPI JHK1K2 spectrum of
51 Eri b to a library of brown dwarf spectra compiled from the
IRTF (Cushing et al. 2005), SpeX (Burgasser 2014), and Montreal
(e.g., Gagné et al. 2015; Robert et al. 2016) spectral libraries. Only
a small subsample of the brown dwarfs have corresponding

mid-IR photometry, so we choose to restrict our comparison to the
near-IR. The spectra within the library were convolved with a
Gaussian kernel to match the spectral resolution of GPI.

To compute the chi-squared between the spectrum of
51 Eri b and the objects within the library, we use two
different equations. The first method permits each individual
filter spectrum to vary freely (unrestricted fit). In the unrest-
RICTED fit, we compute the χ² statistic for the jth object within
the library as

$$ \chi_j^2 = \sum_{i=1}^{4} (S_i - \alpha_{i,j} F_{i,j})^T C_i^{-1} (S_i - \alpha_{i,j} F_{i,j}), \quad (1) $$

where $S_i$ is the spectrum of the planet, $C_i$ is the covariance matrix calculated in Section 6, and $F_{i,j}$ is the spectrum of the jth comparison brown dwarf, all for the ith filter. For each object, the scale factor $\alpha_{i,j}$ that minimizes $\chi^2$ is found using a downhill simplex minimization algorithm. In this method, the scale factor for each object, $\alpha_{i,j}$, is allowed to vary between the four filters (JHK1K2). This is equivalent to allowing the near-IR colors to vary freely up and down in order to better fit the object (e.g., Burningham et al. 2011).

In the second method, the individual filter spectra are still allowed to vary, only within the satellite spot brightness ratio uncertainty (restricted fit), thereby restricting the scale factor for each filter. For the restricted fit, the scale factor is split into two components. The first, $\alpha_j$, is independent of the filter and accounts for the bulk of the difference in flux between 51 Eri b and the comparison object due to differing distances and radii. The second, $\beta_{i,j}$, is a filter-dependent factor that accounts for uncertainties in the satellite spot ratios given in Maire et al. (2014). Equation (1) is modified to include an additional cost term restricting the possible values of $\beta_{i,j}$:

$$ \chi_j^2 = \sum_{i=1}^{4} \left[ (S_i - \alpha_j \beta_{i,j} F_{i,j})^T C_i^{-1} (S_i - \alpha_j \beta_{i,j} F_{i,j}) + N_i \left( \frac{\beta_{i,j} - 1}{\sigma_i} \right)^2 \right] \quad (2) $$

where $N_i$ is the number of spectral channels in the 51 Eri b spectrum for the ith filter, and $\sigma_i$ is the uncertainty in the satellite spot flux ratio given in Maire et al. (2014) for the same filter. The second term in Equation (2) penalizes values of the scale factor, $\beta_{i,j}$, that are very different from the satellite spot uncertainty and thus increases the chi-squared for objects significantly different from 51 Eri b.

The spectral type of 51 Eri b was estimated for both fits from the $\chi^2$ of the L5–T9 near-IR spectral standards (Burgasser et al. 2006; Kirkpatrick et al. 2010; Cushing et al. 2011). To compute the weighted mean and standard deviation of 51 Eri b, we converted the spectral type to a numerical value for the standard brown dwarfs, that is, L5 = 75, T5 = 85. Each numerical spectral type when compared to 51 Eri b is weighted according to the ratio of its $\chi^2$ to the minimum $\chi^2$ for all standards (e.g., Burgasser et al. 2010), and the lowest value was adopted as the spectral type of 51 Eri b. A systematic uncertainty of one-half subtype was assumed for the standards. We find that the two estimates are consistent with one another, that is, T6.3±1.3 and 6.1±1.4 for unrestricted and restricted fits, respectively; see Figure 6. We adopt a spectral type for 51 Eri b of T6.5±1.5 from the unrestricted fit, rounded to the nearest half subtype.

The best-fit object for both the unrestricted and restricted fits was G 204-39 B (SDSS J175805.46+463311.9; $\chi^2 = 1.033$ and 1.209), a T6.5 brown dwarf common proper motion companion to the nearby M3 star G 204-39 A (Faherty et al. 2010). G 204-39 B has marginally low surface gravity based on photometric (log g ≈ 4.5; Knapp et al. 2004) and spectroscopic measurements (log g = 4.7–4.9; Burgasser et al.

<figure><img src="image_8.png" /><figcaption>Figure 6. Comparison of L5 to T9 field (gray circles) and young (yellow stars) brown dwarf JHK spectra to 51 Eri b using the reduced $\chi^2$. The standard brown dwarf for each spectral bin is plotted with a red cross (Burgasser et al. 2006; Kirkpatrick et al. 2010; Cushing et al. 2011). The dashed and dotted vertical lines give the best-fitting spectral type and corresponding uncertainty. Top: Each spectral band of the comparison was allowed to float to find the lowest chi-squared while fitting the planet spectrum. Bottom: The spectrum was allowed to float up or down in flux, but was penalized by the spot ratio uncertainty in each respective band.</figcaption></figure>

2006), indicative of it being younger than the field population. While the binary system is not thought to be a member of any known young moving group (Gagné et al. 2014), the stellar primary can be used to provide a constraint on the age of the system. Combining the X-ray and chromospheric activity indicators for the M-dwarf primary and a comparison of the luminosity of the secondary with evolutionary models, Faherty et al. (2010) adopt an age of 0.5–1.5 Gyr for the system. 51 Eri b is redder than the spectrum of G 204-39 B (Figure 7), especially in terms of the H-K color, which is a photometric diagnostic of low surface gravity among T-dwarfs (e.g., Knapp et al. 2004). This is consistent with the younger age of 51 Eri b, and the most likely cause for this is that it has lower surface gravity than does G 204-39 B.

Additional good matches to the 51 Eri b spectrum include 2MASS J22282889–4310262 (2M 2228–43, $\chi^2 = 1.07$ and 1.26 for the two fits) and 2MASS J10073369–4555147 (2M 1007–45, $\chi^2 = 1.07$ and 1.33). 2M 2228–43 is a well-studied T6 brown dwarf that exhibits a spectrophotometric variability in multiple wavelengths that is indicative of patchy clouds in the photosphere (Buenzli et al. 2012; Yang et al. 2016). 2M 1007–45 is a T5 brown dwarf at a distance of 17 ± 2 pc (Smart et al. 2013). It was identified by Looper et al. (2007) as a low surface gravity object based on its H$_2$O –J versus K/H spectral ratios defined in Burgasser et al. (2006); comparisons against solar-metallicity models imply an age of between 200 and 400 Myr (Looper et al. 2007).

The best-fit objects for each spectral type between spectral types T4.5 and T7.5 using the restricted fit are plotted in Figure 7. While the quality of the fits were generally good, none of the objects were able to provide a good match across all of the bands simultaneously, being too luminous in either the J or K bands. Differences in surface gravity, effective temperature, or metallicity could be the cause (e.g., Knapp et al. 2004). The poor fit to the color of 51 Eri b is especially
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of the spectra, using the restricted fit, of the best-fitting T4.5 to T7.5 field brown dwarfs to 51 Eri b. The spectra shown in this figure are a subset of the data plotted in Figure 6. The brown dwarf spectral fits plotted here use the restricted chi-squared equations presented in Equation (2). The T4.5 and T5.0 spectra are from Looper et al. (2007), the T5.5 is from Burgasser et al. (2008), the T6.0 is from Burgasser et al. (2004), the T6.5 and T7.5 are from Burgasser et al. (2006), and the T7.0 is from Dupuy & Liu (2012).</figcaption></figure>

apparent in the CMDs plotted in Figure 5, with 51 Eri b having unusually red near-IR colors relative to similar spectral type objects.

## 4.2. Comparison against Young Brown Dwarfs

Searches for young companions and moving group objects have resulted in detections of several tens to hundreds of million year old L-type brown dwarf and planetary-mass companions as well as the identification of L-dwarf subclasses based on youth (e.g., Allers & Liu 2013; Filippazzo et al. 2015; Faherty et al. 2016; Liu et al. 2016). In comparison, there exist relatively few known (or suspected) young T-dwarf brown dwarfs. In Figure 8, we plot the known young T-dwarfs and compare them in a manner similar to what was done above for field brown dwarfs. The chi-squared for the fits is not much better than what is seen for the field dwarfs, which is likely due to the absence of young T-dwarfs of spectral type similar to 51 Eri b.

The brown dwarf SDSS J1110+0116 with a spectral type of T5/T5.5 is the best-fitting young comparison object. It has been identified as a bona fide member of the AB Doradus moving group and is thus young (110–130 Myr) and low in mass (10–12 $M_{Jup}$; Gagné et al. 2015). The other young field object that closely matches the near-IR spectrum of 51 Eri b is the T7 peculiar brown dwarf, CFBDSIR J2149–0403 (Delorme et al. 2012). CFBDSIR J2149–0403 was originally suggested to be a member of the AB Doradus moving group, but Delorme et al. (2017) find that the parallax and kinematics of the free-floating object rule out its membership in any known young moving group. However, despite the lack of proof of youth, medium-resolution spectroscopy examining the equivalent width of the KI doublet at 1.25 µm suggests that the object has low surface gravity and is most likely a young planetary-mass object (2–13 $M_{Jup}$). An alternative solution is that it is a higher mass, 2–40 $M_{Jup}$, brown dwarf with high metallicity. CFBDSIR J2149–0403 shows stronger methane absorption features in the red end of the H-band spectrum as compared to 51 Eri b. However, it is worth pointing out that while both young objects, SDSS J1110+0116 and CFBDSIR J2149–0403, are reasonable matches across the *J* and *H* spectra of 51 Eri b, they appear to be underluminous in the K band. A likely reason for this is that 51 Eri b is much younger than both of the comparison companions and thus has the lowest surface gravity among the three objects (Burgasser et al. 2006).

## 4.3. A Very Red T6 or an L-T Transition Planet?

Based on the position of 51 Eri b in Figure 5, it appears that the trend of planetary-mass objects having redder colors compared to the field, seen in young L-type brown dwarfs and planetary-mass companions (Faherty et al. 2016; Liu et al. 2016), possibly continues for the T-type companions. Note that the *K* − *L*<sub>P</sub> CMD shows little reddening, which is natural if clouds are causing the effect. The effect of clouds is negligible in the *K* and *L*<sub>P</sub> bands. Across both the near- and mid-IR CMDs, 51 Eri b is one of the reddest T-type objects, and within its spectral classification it has the reddest colors. This trend in the 51 Eri b colors was originally noted in Macintosh et al. (2015), where they compared the *L*<sub>P</sub> versus *H* − *L*<sub>P</sub> color for the planet and noted that it was clearly redder than the field. Rather than simply being redder than the field T-dwarfs because of the presence of clouds, we present a second possible
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison of the spectra of known young T-dwarfs to that of 51 Eri b. Similar to the field sequence, the fits presented here were computed using the restricted chi-squared. From top to bottom, the four spectra were sourced from Luhman et al. (2007), Naud et al. (2014), Burgasser et al. (2006), and Delorme et al. (2012).</figcaption></figure>

interpretation for the red colors of 51 Eri b, of the planet still undergoing the process of transitioning from L-type to T-type. This hypothesis assumes that the evolutionary track followed is gravity dependent, with examples for higher mass objects shown in Figure 9. In this scenario, 51 Eri b transitions at fainter magnitudes than that seen for field L-T transition brown dwarfs, and it has not yet completed its evolutionary transition to reach the blue colors typical of field mid-T dwarfs.

In Figure 9, we replot the *J* vs. *J* − *H* panel from the series of CMDs shown in Figure 5. In addition to the photometry of 51 Eri b and the field and young brown dwarfs, we also

<figure><img src="image_11.png" /><figcaption>**Figure 9.** The *J* vs. *J* − *H* brown dwarf and imaged exoplanet color-magnitude diagram reproduced from Figure 5. The photometry for 51 Eri b is shown with the red star. Also plotted on the CMD are the evolutionary tracks for 5 and 14 *M*<sub>Jup</sub> objects (Marley et al. 2012), with the solid red line and dashed black line, respectively. The models assume a simple gravity dependence for the initiation of the transition. A few ages for the 5 *M*<sub>Jup</sub> track have been overplotted. The L-T transition for the 5 *M*<sub>Jup</sub> planet starts at approximately 900 K and 20 Myr, but for a lower-mass planet such as 51 Eri b, it will occur at younger ages.</figcaption></figure>

overplot two low-mass, 5 and 14 *M*<sub>Jup</sub>, evolutionary model tracks (assuming hot-start conditions) from Saumon & Marley (2008) and Marley et al. (2012). If the L-T transition is gravity dependent, as multiple lines of evidence now suggest (Leggett et al. 2008; Dupuy et al. 2009; Stephens et al. 2009), then lower mass objects may turn blue at fainter absolute magnitudes than field objects. In Figure 9, we show a simple model in which the L-to-T transition begins at 900 K at log g = 4 (solid red line) instead of 1200 K at log g = 5.3 (dashed black line). In the case of a 5*M*<sub>Jup</sub> planet, the L-to-T transition begins and ends about one magnitude fainter in the *J* band than is observed for the field population. Furthermore, the congruence of the spectrum of SDSS J1110+0116 with 51 Eri b (Figure 8) is interesting as SDSS J1110+0116 lies just short of the blue end of the field L-to-T transition, although it does so at an absolute magnitude just slightly fainter than the field transition magnitude. While these simple models explain the fainter absolute magnitudes of the transition, their colors are too blue and appear to miss the younger brown dwarf and free-floating planets. Similarly, the models are too blue to match the T-dwarf sequence. Clearly, more sophisticated modeling of evolution through the L-to-T transition, accounting for inhomogeneous cloud cover and a gravity-dependent transition mechanism as well as a range of initial conditions, is required. Testing this hypothesis is difficult and would require knowledge of the true mass of the companion as well as the formation mechanism. If this hypothesis is true, then the only objects that are brighter on the CMD should be higher mass objects. There should not be any lower mass objects above and to the left of 51 Eri b on the *J* vs. *J* − *H* CMD shown in Figure 9.
<table><caption>**Table 3**
Model Grid Parameters</caption><thead><tr><th>Model Name</th><th>Effective Temperature (K)</th><th>Surface Gravity [log g] (dex)</th><th>Metallicity [M/H] (dex)</th><th>Cloud Parameter (f<sub>sed</sub>)</th><th>Cloud Hole Fraction (%)</th></tr></thead><tbody><tr><td>Iron/Silicate Cloud Grid</td><td>600–1000</td><td>3.25</td><td>0.0</td><td>2</td><td>0–75</td></tr><tr><td>Sulfide/Salt Cloud Grid</td><td>450–900</td><td>3.5–5.0</td><td>0.0, 0.5, 1.0</td><td>1, 2, 3, 5</td><td>...</td></tr><tr><td>Cloudless Grid</td><td>450–900</td><td>3.5–5.0</td><td>0.0, 0.5, 1.0</td><td>no cloud</td><td>...</td></tr></tbody></table>

## 5. Modeling the Atmosphere of 51 Eri b

For the purpose of modeling the complete SED of 51 Eri b, we made use of two updated atmospheric model grids from the same group, focusing on different parameter spaces (see Table 3). The first grid, described in Marley et al. (1996, 2002, 2010), focused on the higher effective temperature atmospheres (L-dwarfs) and includes iron and silicate clouds in the atmosphere. The second grid, described in Morley et al. (2012, 2014) and Skemer et al. (2016), is designed for lower effective temperatures (T- and Y-dwarfs) and include salt and sulfide clouds in the atmosphere, which are expected to condense in the atmospheres of mid-to-late-T dwarfs.

The methodology used to fit the models to the data is the same for both model grids. To fit the models to the data, we bin the model spectra to match the spectral resolution of the GPIES spectra across each of the JHK1K2 filters. For the photometry, we integrated the model flux through the Keck/NIRC2 $L_p$ and $M_s$ filter profiles, respectively. The estimation of the best-fitting model is done by computing the chi-squared value for each model in the grid compared to the data using Equation (2). We made use of the covariance matrices estimated for the four spectral channels described in the Appendix and also included the variance for each of the two photometric data points to compute the chi-squared statistic. Note that we use the restricted fit equation in the computation of the best-fitting model. This equation permits each of the individual filters to scale within the $1\sigma$ error of the satellite spot ratios. We also did the fitting without the scaling factor and found that the results are similar.

As stated in Section 2.1.1, the use of the covariance affects the model fitting where the peak of the posterior distribution occurs at slightly cooler effective temperatures, consistent within the errors. Due to the high spectral correlation in the J band (see Figure 19), when using the covariance, the best-fitting models are not models that pass through the data but rather models that have lower flux in the J band than the data. We present the specific modeling details in the following text.

### 5.1. Iron and Silicate Cloud Models

In Section 4.3, we suggested that 51 Eri b, rather than having completely evolved to T-type, could be transitioning from L to T. In this scenario, the cloud composition of the planetary atmosphere might still be influenced by the deep iron and silicate condensate grains and patchy cloud atmosphere. Therefore, we compared the planet SED to a grid of models with a fixed low surface gravity and solar metallicity, where the key variable is the cloud hole fraction, and the unique aspect of this grid is the presence of iron or silicate clouds in an atmosphere with clear indications of methane absorption. The clouds are modeled using the prescription presented in Ackerman & Marley (2001), where cloud thickness is parameterized via an efficiency factor ($f_{sed}$). Small values of

$f_{sed}$ indicate atmospheres with thick clouds, whereas large values of $f_{sed}$ are for atmospheres with large particles that rain out of the atmosphere, leaving optically thinner clouds. As mentioned earlier, the primary condensate species in this grid are iron, silicate, and corundum clouds, molecules that are expected to dominate clouds in L-dwarfs (Saumon & Marley 2008; Stephens et al. 2009). At the L-T transition, clouds are expected to be patchy, so for each $T_{eff}$, the models went from fully cloudy, with $f_{sed} = 2$ and 0% holes, to an atmosphere with $f_{sed} = 2$ and 75% holes (patchy clouds). The methodology used to calculate the flux emitted from the patchy cloud atmosphere includes both cloud and cloud-free regions simultaneously in the atmosphere using a single global temperature-pressure profile and is not created via a linear combination of two models, as is sometimes done in the literature (Marley et al. 2010). The iron and silicate cloud grid models use solar metallicity (Lodders 2003). The opacity database used for the absorbers are described in Freedman et al. (2008), including updated molecular line lists for ammonia and methane (Yurchenko et al. 2011; Yurchenko & Tennyson 2014). The models span effective temperatures from 600 to 1000 K for solar metallicity ($[M/H] = 0.0$) and low surface gravity ($log g = 3.25, 3.50$; see Table 3).

Presented in Figure 10 is the best-fitting model to the SED of 51 Eri b. Stated in the figure are the model parameters along with the radius of the planet required to scale the model spectrum to match the planet SED. This scaling factor is required since the model spectra are typically computed to be the emission at the photosphere or at 10 pc from the object. One of the free parameters in most model fitting codes is the term $R^2/d^2$ to scale the model flux to match the SED, where R is the radius of the planet and d is the distance to the object. For 51 Eri, the distance is known to better than 2% (see Table 2), so we only fit the radius term. Shown in Figure 11 is the posterior distribution for the radius, where we find that the best-fitting radii are significantly smaller than that predicted by evolutionary models, for example, 1.33–1.14 $R_{Jup}$ for a 2–10 $M_{Jup}$ hot/cold-start planet at the age of 51 Eri (Marley et al. 2007; Fortney et al. 2008). This discrepancy has been noted previously as well for the HR8799 planets (Marois et al. 2008; Bowler et al. 2010; Barman et al. 2011; Currie et al. 2011; Marley et al. 2012), $\beta$ Pic b (Morzinski et al. 2015), and for 51 Eri b itself in the discovery paper (Macintosh et al. 2015). In an attempt to circumvent this issue, while modeling the SED we adopted a Bayesian prior probability density function for the radius in the form of a Gaussian centered on the expected radius from evolutionary models (green line in Figure 12), with the width chosen to include the radius of Jupiter. Without the prior (i.e., using a uniform prior), the median radius is $0.68 R_{Jup}$ and $T_{eff} \sim 740$ K; with the prior, the median radius value is forced closer to the predictions of evolutionary models (red line in Figure 12) at $0.98 R_{Jup}$ and $T_{eff} \sim 690$ K, biasing the luminosity of the planet
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Spectral energy distribution of 51 Eri b with the best-fitting iron and silicate cloudy model.</figcaption></figure>

to larger values. When fitting the SED, the term that is conserved is the luminosity rather than the effective temperature or the radius. Adopting the evolutionary radius and marginalizing over the uncertainty in radius raises the luminosity ($\log L/L_\odot$) from −5.83 to −5.65. Since observational constraints on the radius for young planets are unavailable, we chose to use an uninformative prior.

Plotted in Figure 11 are the normalized posterior distributions for each of the model parameters varied in the model fit, along with the covariances to show how each of the parameters are affected. Since the grid only had a few models with $\log g = 3.5$, with the majority being 3.25, we marginalized over the surface gravity. The irregular shape of the effective temperature posterior is caused by the missing models in the grid. The median effective temperature, 737 K, estimated from the grid falls right in between the range of best-fitting temperatures from the models in the Macintosh et al. (2015) paper (700–750 K). However, based on the shape of the posterior and the covariances, the peak of the effective temperature distribution extends to cooler temperatures. Since the L-to-T transition has been suggested to arise from holes or low-opacity patches appearing in an initially more uniform cloud deck (Ackerman & Marley 2001; Burgasser et al. 2002; Marley et al. 2010), our finding here that partly cloudy models best fit the 51 Eri b spectrum is consistent with this interpretation. In general, however, the models struggled to fit the entire planet SED, typically being able to fit either the near- or mid-IR portions of the SED. The inability to fit mid-IR photometry suggests that chemical equilibrium models are not appropriate. Disequilibrium chemistry predicts less CH$_4$ in the atmosphere and could explain the higher flux at 1.6 µm and in the $L_P$ band. It would also introduce CO, accounting for the lower flux in the $M_S$ band.

## 5.2. Sulfide and Salt Cloud Models

In Section 4.1, we showed that the best-fitting spectral type of 51 Eri b is a mid-to-late-T dwarf. At the effective temperatures of mid-to-late-T dwarfs, Cr, MnS, Na$_2$S, ZnS, and KCl are expected to condense and form clouds high in the photosphere. The second grid we tested the planet SED against made use of a model grid that includes salt and sulfide clouds to test additional parameters such as the surface gravity and metallicity (which were varied, unlike the iron/silicate grid) and the properties of clouds typically associated with T-dwarfs. The grid was designed specifically for lower temperature objects (450 ~ 900 K; Morley et al. 2012, 2014) and has been

used successfully to reproduce the SED of GJ 504 b (Skemer et al. 2016), a cool low-mass companion with a similar spectral type (late-T) that is comparable to 51 Eri b (Kuzuhara et al. 2013). Note that the use of this cloud grid does not preclude the possibility of the planet transitioning from L to T.

Also included as part of this grid are the clear atmosphere models from Saumon & Marley (2008), the ranges for which are presented in Table 3. The ranges of parameters varied are presented in Table 3, including temperatures, surface gravities, metallicities, and sedimentation factor ($f_{sed}$) ranging from cloudy ($f_{sed} = 1$) to cloud-free. The cloud model used in the sulfide/salt grid is the same as the one described above. In addition to the opacity updates mentioned above, opacity effects due to alkali metals (Li, Na, K) have been included using the results from Allard et al. (2005). For effective temperatures between 450 and 775 K, the grid is complete with models available for every step of the varied parameters. For effective temperatures between 800 and 900 K, the temperature steps switch from increments of 25–50 K, and there are no models with $f_{sed}$ values of 1 and 2. This grid does not include the opacity effect that is due to iron and silicate condensates. A future series of papers describing an extended atmosphere model grid will describe the updates, but the present grid extends the models to greater than solar metallicities.

In Figure 13, we present the four best-fitting model atmospheres for 51 Eri b. Presented in each panel are the atmosphere with the lowest reduced chi-squared in one of four cases, namely, solar and cloudless (top left), solar and cloudy (top right), nonsolar and cloudless (bottom left), and nonsolar and cloudy (bottom right). Both cloudless model atmospheres are warmer and thus fit the near-IR spectrum of the planet while completely missing the $L_P$ photometry. The cloudy atmosphere model fits are cooler and do a much better job of fitting the overall SED of 51 Eri b, and the best-fitting atmospheres for both solar and nonsolar metallicity have very similar reduced chi-squared values.

The normalized posterior distributions for the different parameters varied as part of the model fitting are shown in Figure 14. The best-fitting $T_{eff}$ (605$^{+61}_{-66}$ K) is much cooler in comparison to the iron/silicate grid, but the values are within $2\sigma$ of each other. We also note that the median might not be the best estimate for the effective temperature PDF in the iron/silicate grid where the peak extended to cooler temperatures. For the surface gravity and metallicity posterior distributions, we present the median values and error bar assuming a Gaussian distribution, though they may not be Gaussian. The surface gravity PDF suggests that the planet has high surface gravity. However, 51 Eri b is clearly a low-mass companion, indicating that the data do not constrain the gravity. A prior might help constrain the distribution, but there are currently no physically motivated priors available for the surface gravity of young planets. Similarly, the PDF for the metallicity is also unconstrained, and higher resolution spectra in the K band might help provide greater constraints on the metallicity of 51 Eri b (Konopacky et al. 2013).

A difference between the iron/silicate and salt/sulfide atmosphere grids is in the planet radius, where the best-fit radii for the cloudy models and the median radius of the PDF for the salt/sulfide models are much closer to evolutionary model predictions. A possible explanation for this discrepancy is that fitting the lower effective temperatures while still matching the bolometric luminosity requires a larger radius. If
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Normalized posterior distributions for the iron and silicate model grid. The PDFs are for the parameters varied in our fit along with the inferred distribution of the luminosity of 51 Eri b. The lines on the 1D histogram indicate the 16th, 50th, and 84th percentile values, while those on the 2D histogram are the 1σ (solid red), 2σ (dashed green), and 3σ (dotted blue) values of the distribution. The values printed above each histogram are the median value along with the 1σ error on it.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** The figure shows the effect of applying a Gaussian radius prior when modeling with the iron/silicate grid. The prior shown by the green line is centered on the radius given by evolutionary models: 1.29 $R_{\text{Jup}}$ (Marley et al. 2007; Fortney et al. 2008). Also plotted are the likelihood (black) and posterior distribution (red).</figcaption></figure>

the iron/silicate models extended to lower temperatures, assuming the continued presence of these clouds at these colder temperatures, it is likely that the radius discrepancy would not be as apparent. The sedimentation factor was fixed (at $f_{\text{sed}} = 2$) in the iron/silicate grids but had varying hole

fractions ($h_{\text{frac}}$). In the sulfide/salt grid, $f_{\text{sed}}$ was varied, and the median value for the distribution is $f_{\text{sed}} = 2.48$. If we equate the $h_{\text{frac}}$ from the iron/silicate model with the $f_{\text{sed}}$ as the physics controlling the emission of flux from the photosphere, then for both model grids the best-fitting models tend to favor the presence of clouds over cloud-free atmospheres. Furthermore, in both cases, the best-fitting models were not the fully cloudy atmospheres, with the smallest $h_{\text{frac}}/f_{\text{sed}}$. While the cloud compositions in both models are different, fitting either grid requires the cloud opacity. This can be achieved in one of two ways: either make the deep iron/silicate clouds be very vertically extended (small $f_{\text{sed}}$) or introduce a new cloud layer in the form of the sulfide/salt clouds.

The cloudy model atmosphere fits presented in Figure 13 match the *H* through *K* spectrum while being slightly underluminous in the *J* and overluminous $M_S$ bands. Given the large photometric errors in the $M_S$ data, the model photometry lies within $2\sigma$ of the data. JWST and other future low-background mid-IR instruments will better constrain the 3–24 µm SED, a further test of current models. In Figure 15, we show 10 of the best-fitting models assuming cloudy (sulfide/salt clouds) or cloudless atmospheres extended out to 20 µm. It is clear from these models that observations with the
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Spectral energy distribution of 51 Eri b with the best-fitting salt and sulfide cloud models. Each panel shows the best-fitting model under the specific conditions: the top two panels show the best-fitting *solar*-metallicity models with a cloudless atmosphere on the left and cloudy atmospheres on the right. The bottom two panels show the best-fitting *nonsolar* models with a cloudless atmosphere on the left and cloudy atmospheres on the right.</figcaption></figure>

coronagraph on the Near InfraRed Camera, spanning the 3–5 µm wavelengths, will add significant constraints on the atmosphere of the planet. If the planet can be studied with the Mid Infrared Instrument, it could be used to apply constraints on the chemical disequilibrium in the atmosphere through observations of NH₃ in the 10–11 µm range.

## 5.3. Luminosity of the Planet

The two different grids used in this study have produced similar luminosity predictions for the planet despite the different cloud compositions. From the iron/silicate grid we infer a bolometric luminosity of $log L/L_\odot = -5.83^{+0.15}_{-0.12}$, and $log L/L_\odot = -5.93^{+0.19}_{-0.14}$ from the sulfide/salt model atmospheres. We compare these luminosity estimates to predictions of evolutionary models to infer the planet mass and discuss its initial formation conditions.

### 5.3.1. Standard Cold- and Hot-start Models

In Figure 16 we compare the bolometric luminosity to evolutionary models for planets formed via the two extreme scenarios, namely, the hot-start and cold-start models (Burrows et al. 1997; Marley et al. 2007). In the hot-start scenario, planets are formed with high initial entropy and are very luminous at birth. This scenario is usually associated with rapid formation in the circumstellar disk through disk instabilities. Alternatively, in the cold-start scenario, which is often associated with current 1D models of the core-accretion mechanism, planets start with a solid core that accretes gas from the stellar disk. The accreting gas loses energy via a radiatively efficient accretion shock and form with low initial entropy and thereby lower postformation luminosity.

The other directly imaged companions plotted in Figure 16 can all be considered as having formed via the hot-start

scenario. Despite the older age assessment for the companion in this study, $26 \pm 3$ Myr (Nielsen et al. 2016) compared to $20 \pm 6$ Myr (Macintosh et al. 2015), the revised luminosity when compared to the system age places 51 Eri b in a location where either cold or hot initial conditions are possible. Based on the hot-start tracks, it would have an inferred mass between 1 and $2 M_{\text{Jup}}$. However, for the cold-start case, the planet mass could lie anywhere between 2 and $12 M_{\text{Jup}}$, since the model luminosity is largely independent of mass at the age of 51 Eri b. Dynamical mass estimates for the planet could help clarify the formation mechanism, especially if the planet mass > $2 M_{\text{Jup}}$.

### 5.3.2. Warm-start Models

Spiegel & Burrows (2012) proposed a complete family of solutions existing between the hot- and cold-start extreme cases. Warm-start models<fnref n="39" /> explore a wide range of initial entropies aimed at covering the possible range of initial parameters that govern the formation of planets. In Figure 17, we compare the inferred bolometric luminosity and the planet SED to models from Spiegel & Burrows (2012). The Spiegel & Burrows (2012) models are evolutionary tracks calculated assuming different initial entropies for the planet, between 8 and 13 $k_B$/baryon, where $k_B$ is Boltzmann's constant, with steps of 0.25 $k_B$/baryon and masses between 1 and 15 $M_{\text{Jup}}$ with steps of 1 $M_{\text{Jup}}$. Four different model atmospheres are considered in combination with the evolutionary model: cloud-free and solar metallicity to fully cloudy with 3× solar metallicity (Burrows et al. 2011). The bolometric luminosity of each point in the grid for each of the four atmosphere scenarios was computed by integrating the SED over the wavelength range. Because of the sparse sampling of

<footnote n="39">³⁹ http://www.astro.princeton.edu/~burrows/</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 14. Normalized posterior distributions for the sulfide and salt model grid, same as Figure 11.</figcaption></figure>

the grid, we linearly interpolate the evolutionary tracks with steps of 0.06 $k_B$/baryon and 0.2 $M_{Jup}$.

In the top row of Figure 17, we plot the probabilities for each grid point measured by comparing the average of the inferred bolometric luminosities from the SED fit ($\log L/L_\odot = -5.87 \pm 0.15$) to the predictions of the Spiegel & Burrows (2012) models with the four atmosphere conditions. For the bottom row in Figure 17, the surface is calculated by fitting the planet SED to the Spiegel & Burrows (2012) model atmosphere grid, using Equation (2). For both comparisons, luminosity and SED, we chose the age of the evolutionary grid best matching the age of 51 Eri (25 Myr), to minimize the number of interpolations, and only varied the mass of the planet and initial entropy for the models.

Mordasini (2013) find that the luminosity of a planet that underwent accretion through a supercritical shock (the standard cold-start core accretion hypothesis) is highly dependent on the mass of the core, $M_{\text{core}}^{2-3}$. Therefore, the continuum of warm-start models can also be explained by similar bulk-mass planets with increasing core mass. These models suggest that the entropy of 51 Eri b can be explained via core accretion, with a core mass ranging between 15 and 127 $M_\odot$, which can reproduce the planet luminosity with various initial entropies.

The four panels generated by fitting the inferred luminosity (upper four panels) appear highly consistent and in agreement with the results from Figure 16. The $1\sigma$ contour encompasses the entire available entropy space, where for intermediate and
<figure><img src="image_17.png" /><figcaption>**Figure 15.** The ten best-fitting cloudy (red) and cloudless (blue) atmospheres over the wavelength range of the James Webb Space Telescope. The median of the models is plotted with a thicker line. The models indicate the divergence between the model fits over the wavelength covered by JWST.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Luminosity of imaged planetary-mass companions as a function of age. For 51 Eri b, the red star is the inferred luminosity from the iron/silicate model grid, while the green square is the inferred luminosity from the sulfide/salt model grid. Also plotted are evolutionary tracks assuming different starting conditions: the “hot-start” (dotted lines) and “cold-start” (solid lines) models of Marley et al. (2007) and Fortney et al. (2008). 51 Eri b is consistent with both hot- and cold-start formation models. A subset of known, directly imaged companions are plotted in the figure, illustrating the difference between 51 Eri b and other imaged planets. Data for the companions, 2M1207 b, HR8799 bcde, HD 95086 b, and GJ 504 b, were taken from Patience et al. (2010), Rajan et al. (2015), Zurlo et al. (2016), De Rosa et al. (2016), and Skemer et al. (2016), respectively.</figcaption></figure>

high entropies the most likely mass for the planet is between 2 and 3 $M_{\text{Jup}}$ and for low initial entropy the most likely mass for the planet increases, making distinguishing between cold, warm, and hot starts difficult.

When we compare the model spectra directly to the planet SED, the surface is qualitatively similar to that made with the luminosity but shifted to higher mass and with the $1\sigma$ contours and best-fit models favoring lower entropy. According to the Mordasini (2013) models, the fits presented here would be consistent with a planet having core masses ranging from 15 to 127 $M_{\oplus}$.

Compared to other directly imaged companions (see figures in Marleau & Cumming 2014), 51 Eri b is the only planet compatible with very low initial entropy and the cold-start case. Tighter constraints on the bolometric luminosity or higher S/N data will help to reduce the width of the two branches, and independent mass constraints from dynamical measurements will enable us to infer the initial entropy and possible formation route. Atmospheric retrievals or higher resolution spectra aimed at exploring and characterizing the planet's chemical composition might also help us understand whether the planet has higher C/O ratios than the star, since planetary C/O can be used to understand planet formation (Öberg et al. 2011; Konopacky et al. 2013).

## 6. Conclusion

In this paper, we have presented the first spectrum of 51 Eridani b in the *K* band obtained with the GPI (*K1* and *K2* bands) as well as the first photometric measurement of the planet at *M*<sub>S</sub> obtained with the NIRC2 Narrow camera. We also obtained an additional *L*<sub>P</sub> photometric point that agrees very well with the *L*<sub>P</sub> measurement taken in the discovery paper (Macintosh et al. 2015). In addition, we revised the stellar photometry by observing the star in the near-IR and estimating its photometry in the mid-IR through an SED fit. The new data are combined with the published *J* and *H* spectra and the *L*<sub>P</sub> photometry to present the SED spanning 1–5 µm for the planet.

As part of the data analysis, we calculated the covariance for each of the spectral data sets, that is, *J*, *H*, *K1*, and *K2*, using the formalism presented in Greco & Brandt (2016). The spectral covariance was used in all the chi-squared minimization performed as part of this study, in combination with the photometric variance. Using the covariance ensured that the photometric points were weighted in a suitable manner and resulted in cooler effective temperatures for the best fits.

We compared the planet photometry to field and young brown dwarfs by fitting their near-IR spectra to 51 Eri b to estimate a spectral type of T6.5 ± 1.5. Due to the relative paucity of known young T-dwarfs, our comparison of the planet spectrum to young T-dwarfs only included a handful of objects, and among the sample, 51 Eri b appears to have the lowest surface gravity based on a comparison of their spectral shape and amplitude.

In a comparison of the near- and mid-IR photometry for the planet to the field and young brown dwarf population via a range of CMDs, we note that 51 Eri b is redder than brown dwarfs of similar spectral types. This was also noted in the discovery paper, and it was proposed that this might be due to the presence of clouds, similar to young L-type planetary-mass companions. In this study, we extended this idea to suggest that a possible reason for the presence of clouds (compared to the field) is that the planet is still transitioning from the L-type to the T-type. This would occur at a lower *J* magnitude than field brown dwarfs due to its lower mass when including a gravity-dependent transition in the evolution (Saumon & Marley 2008).

We also fit the planet SED with two different model atmosphere grids that varied in the composition of molecules that could condense in the atmosphere. The best-fitting models in both cases were those that contained large amount of condensates in the atmosphere as compared to cloud-free atmospheres. Through the iron/silicate grid, we estimate that the planet has a patchy atmosphere with 10%–25% hole
<figure><img src="image_19.png" /></figure>

**Figure 17.** Comparison of the planet spectrum and luminosity to a combination of initial entropy ($k_B$/baryon) and planet mass ($M_{Jup}$) from the warm-start evolutionary models of Spiegel & Burrows (2012). The four different atmospheres tested include cloud-free and hybrid cloud models, with both solar and supersolar metallicity. Also plotted are the $1\sigma$ (solid red), $2\sigma$ (dotted cyan), and $3\sigma$ (dashed yellow) contours. The entropy plotted in the figure and used in the modeling is not the entropy for the evolved object but rather the entropy at formation. The best-fitting model fit is indicated by the large circle (white and red circle). The orange filled circles show the hot-start model limits, while the blue filled circles show the cold-start limits, which are presented as the boundary cases in Spiegel & Burrows (2012). The top row compares the model luminosity to the inferred luminosity for 51 Eri b, and the bottom row compares directly the SED to the evolutionary model spectra.

fraction in the surface cloud cover, which is consistent with the $f_{sed}$ values of 2–3 resulting from the sulfide/salt grid. The median effective temperature from the two grids is 737$^{+39}_{-46}$ K and 605$^{+61}_{-66}$ K for iron/silicate and sulfide/salt, respectively. This value is slightly cooler than that in Macintosh et al. (2015), where the best-fit models had temperatures of 700 K and 750 K, respectively. The surface gravity and metallicity both appear to be unconstrained by the data, but empirical fits to young T-dwarfs suggest that the planet has lower surface gravity.

The two atmosphere grids provide similar luminosity estimates, which were compared to hot-, warm-, and cold-start models. 51 Eri b appears to be one of the only directly imaged planets that is consistent with the cold-start scenario, and a comparison of the planet SED to a range of initial entropy models indicates that cloudy atmospheres with low initial entropies provide the best fit to the planet SED.

Following the submission of this study for publication, a paper on 51 Eri b using spectrophotometry taken with the VLT/SPHERE was published by Samland et al. (2017). Their study includes new *YJH* spectra as well as *K1K2* photometry, in addition to the *H* spectrum and *L*<sub>P</sub> photometry from Macintosh et al. (2015). Their results are consistent in parts with ours, although we note that the SPHERE *J*-band spectrum is fainter than the GPI *J* spectrum, while their *K1K2* photometry is brighter than the GPI spectrum (and corresponding integrated GPI photometry). These differences could very well be caused by the application of different algorithms, where Samland et al. (2017) demonstrate that different algorithms can result in spectra with a range of flux values, including ones that agree with the GPI *J* spectrum. Future

studies will need to analyze all the available data sets using a common pipeline for data processing and analysis to understand whether the differences arise from the algorithms or are due to other causes.

With future space missions such as the *James Webb Space Telescope*, the 3–24 µm SED of this planet could be observed at higher S/N, providing tests of current atmospheric models. The best-fitting atmosphere models further indicate that the planet might have a cloudy atmosphere with patchy clouds, making 51 Eri b a prime candidate for atmospheric variability studies that might be possible with future instrumentation. Further analysis of this data using methods such as atmosphere retrievals could permit an exploration of other planet parameters that were not considered in this study, such as chemical composition of the atmosphere and the thermal structure.

The authors thank Gabriel Marleau for the discussion on warm-start models. The Gemini Observatory is operated by the AURA 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), Ministério da Ciência, Tecnologia e Inovação (Brazil), and Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina). This work was supported by NSF grants AST-1411868 (A.R., J.L.P., B.M.), AST-1518332 (R.J.D.R.), and DGE-1311230 (K.W.D.). F.M. and E.N. are supported by NASA Grant NNX14AJ80G. This work was supported by Fonds de Recherche du Québec (J.R., R.D., D.L.). K.M.M. and T.S.B. are supported by the NASA Exoplanets Research
Program (XRP) by cooperative agreement NNX16AD44G. G.V. and J.K.W. acknowledge JPL's ESI program for GPI-related funding. The results reported herein benefited from collaborations or information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate and the NExSS grant NNX15AD95G. Portions of this work were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. This research has benefited from the SpeX Prism Library (or SpeX Prism Library Analysis Toolkit) maintained by Adam Burgasser at http://www.browndwarfs.org/spexprism, the IRTF Spectral Library maintained by Michael Cushing, and the Montreal Brown Dwarf and Exoplanet Spectral Library maintained by Jonathan Gagné.

*Facilities:* Gemini:South (GPI), Keck:II (NIRC2).

*Software:* GPI Data Reduction Pipeline (v1.3.0; Perrin et al. 2014), pyKLIP (Wang et al. 2015).

## Appendix
### Derivation of Spectral Covariance

We follow the method described in Greco & Brandt (2016) to measure the interpixel correlation within the PSF-subtracted images, and we convert these into a covariance matrix. For each image (J, H, K1, and K2), the correlation $\psi_{ij}$ between pixel values at wavelengths $\lambda_i$ and $\lambda_j$ within a 1.5 $\lambda/D$ annulus was estimated as

$$ \psi_{ij} = \frac{\langle I_i I_j \rangle}{\sqrt{\langle I_i^2 \rangle \langle I_j^2 \rangle}}, \quad (3) $$

where $\langle I_i \rangle$ is the average intensity within the annulus at wavelength $\lambda_i$. This was repeated for all wavelength pairs and at five different separations: 350, 454 (the separation of 51 Eri b), 550, 650, and 750 mas. To avoid biasing the measurement, 51 Eri b was masked in the 454 mas annulus.

The measurements of the correlation $\psi_{ij}$ at the eight different separations within the final image were used to fit the parameterized correlation model of Greco & Brandt (2016):

$$ \psi_{ij} \approx A_{\rho} \exp \left[ -\frac{1}{2} \left( \frac{\rho}{\sigma_{\rho}} \frac{\lambda_{i} - \lambda_{j}}{\lambda_{c}} \right)^{2} \right] + A_{\lambda} \exp \left[ -\frac{1}{2} \left( \frac{1}{\sigma_{\lambda}} \frac{\lambda_{i} - \lambda_{j}}{\lambda_{c}} \right)^{2} \right] + A_{\delta} \delta_{ij}, \quad (4) $$

where the symbols are as in Greco & Brandt (2016). This model is based on the assumption that the correlation consists of three components. The first two terms model the contribution of the speckle noise and the correlation induced by the interpolation within the reduction process. The third models uncorrelated noise, such as read noise, which does not contribute to the off-diagonal terms of the correlation matrix. The amplitudes of the first two terms ($A_\rho$, $A_\lambda$) were allowed to vary with separation, while the two correlation lengths ($\sigma_\rho$, $\sigma_\lambda$) were fixed. As the sum of the amplitudes must equal unity, $A_\delta$ was derived from the other amplitudes. Figure 18 shows an example of the spectral correlation as a function of the angular separation for the H-band spectral cube, and $\lambda_c$ is the central

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Example of the correlation function at the various angular separations included in the fit for the H-band spectral cube. The different colors correspond to the angular separations, with the circles being the value of the correlation for all of the wavelength pairs, and the lines of the same color indicate the best fit to the distribution calculated using Equation (4).</figcaption></figure>

<table><caption>**Table 4**
Correlation Model Parameters</caption><thead><tr><th>Band</th><th>A<sub>ρ</sub></th><th>A<sub>λ</sub></th><th>A<sub>δ</sub></th><th>σ<sub>ρ</sub></th><th>σ<sub>λ</sub></th></tr></thead><tbody><tr><td>J</td><td>0.43</td><td>0.43</td><td>0.14</td><td>0.44</td><td>0.05</td></tr><tr><td>H</td><td>0.70</td><td>0.27</td><td>0.03</td><td>0.45</td><td>0.01</td></tr><tr><td>K1</td><td>0.51</td><td>0.41</td><td>0.07</td><td>0.68</td><td>0.004</td></tr><tr><td>K2</td><td>0.30</td><td>0.62</td><td>0.08</td><td>0.43</td><td>0.004</td></tr></tbody></table>

wavelength of the spectrum (1.65 µm for H). The colored lines in the plot are the best fits to Equation (4).

Due to the high dimensionality of the problem, we use a parallel-tempered Markov chain Monte Carlo algorithm (Foreman-Mackey et al. 2013) to find the global minimum. The best-fit parameters at the separation of 51 Eri b within the PSF-subtracted image at each band are given in Table 4. Using these parameters, the covariance matrix, C, was constructed for each band. The diagonal elements contained the square of the uncertainties of the spectrum of the planet, and the off-diagonal elements were calculated using

$$ \psi_{ij} \equiv \frac{C_{ij}}{\sqrt{C_{ii} C_{jj}}}. \quad (5) $$

The fitted parameters in Table 4 demonstrate that the primary cause of correlation at the shorter wavelengths is speckle noise, with the correlation induced by interpolation becoming more significant in the K1 and K2 images. In each case, the amplitude of the speckle noise term ($A_\rho$) is significantly higher than is seen for HD 95086 b (De Rosa et al. 2016). This can be attributed to the fact that 51 Eri A is approximately two magnitudes brighter at K1 (than HD 95086 A), leading to a significantly brighter speckle field. The typical correlation lengths in the PSF-subtracted image for each band are visualized in Figure 19, with the data being highly correlated
<figure><img src="image_21.png" /><figcaption>**Figure 19.** The correlation matrices calculated for each of the four *JHK1K2* spectra. Going from the *J* band through *K2*, the correlation length can be seen to change as a function of speckle vs. background noise. The spectra are highly correlated at *J* with up to five channels showing high correlation values, down to ~3 at *K2*, which is a consequence of the spectral resampling.</figcaption></figure>

at the *J* band at wavelengths separated by up to five spectral channels.

References

Abt, H. A., & Morrell, N. I. 1995, ApJS, 99, 135

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

Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

Allard, N. I., Allard, F., & Kielkopf, J. F. 2005, A&A, 440, 1195

Allers, K. N., & Liu, M. C. 2013, ApJ, 772, 79

Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948

Barman, T. S., Macintosh, B., Konopacky, Q. M., & Marois, C. 2011, ApJ, 733, 65

Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593

Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, Proc. SPIE, 7014, 701418

Binks, A. S., & Jeffries, R. D. 2014, MNRAS, 438, L11

Bonnefoy, M., Marleau, G.-D., Galicher, R., et al. 2014, A&A, 567, L9

Bonnefoy, M., Milli, J., Ménard, F., et al. 2017, A&A, 597, L7

Bowler, B. P., Liu, M. C., Dupuy, T. J., & Cushing, M. C. 2010, ApJ, 723, 850

Bowler, B. P., Liu, M. C., Mawet, D., et al. 2017, AJ, 153, 18

Buenzli, E., Apai, D., Morley, C. V., et al. 2012, ApJL, 760, L31

Buenzli, E., Apai, D., Radigan, J., Reid, I. N., & Flateau, D. 2014, ApJ, 782, 77

Burgasser, A. J. 2014, in International Workshop on Stellar Spectral Libraries, Astronomical Society of India Conf. Ser. 11, ed. H. P. Singh, P. Prugniel, & I. Vauglin, 7

Burgasser, A. J., Burrows, A., & Kirkpatrick, J. D. 2006, ApJ, 639, 1095

Burgasser, A. J., Cruz, K. L., Cushing, M., et al. 2010, ApJ, 710, 1142

Burgasser, A. J., Geballe, T. R., Leggett, S. K., Kirkpatrick, J. D., & Golimowski, D. A. 2006, ApJ, 637, 1067

Burgasser, A. J., Liu, M. C., Ireland, M. J., Cruz, K. L., & Dupuy, T. J. 2008, ApJ, 681, 579-93

Burgasser, A. J., Lopez, M. A., Mamajek, E. E., et al. 2016, ApJ, 820, 32

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

Burgasser, A. J., McElwain, M. W., Kirkpatrick, J. D., et al. 2004, AJ, 127, 2856

Burningham, B., Leggett, S. K., Homeier, D., et al. 2011, MNRAS, 414, 3590

Burrows, A., Heng, K., & Nampaisarn, T. 2011, ApJ, 736, 47

Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856

Carter, B. S. 1990, MNRAS, 242, 1

Chauvin, G., Lagrange, A.-M., Dumas, C., et al. 2005, A&A, 438, L25

Cohen, M., Wheaton, W. A., & Megeath, S. T. 2003, AJ, 126, 1090

Cotten, T. H., & Song, I. 2016, ApJS, 225, 15

Currie, T., Burrows, A., Itoh, Y., et al. 2011, ApJ, 729, 128

Currie, T., Lisse, C. M., Kuchner, M., et al. 2015, ApJL, 807, L7

Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al. 2011, ApJ, 743, 50

Cushing, M. C., Rayner, J. T., & Vacca, W. D. 2005, ApJ, 623, 1115

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, yCat, 2246, 0

Cutri, R. M., et al. 2013, yCat, 2328, 0

De Rosa, R. J., Nielsen, E. L., Blunt, S. C., et al. 2015, ApJL, 814, L3

De Rosa, R. J., Rameau, J., Patience, J., et al. 2016, ApJ, 824, 121

Delorme, P., Dupuy, T., Gagné, J., et al. 2017, arXiv:1703.00843

Delorme, P., Gagné, J., Malo, L., et al. 2012, A&A, 548, A26

Dupuy, T. J., & Liu, M. C. 2012, ApJS, 201, 19

Dupuy, T. J., Liu, M. C., & Ireland, M. J. 2009, ApJ, 699, 168

ESA 1997, in Hipparcos and Tycho Catalogues, ESA Special Publication 1200 (Noordwijk: ESA), http://adsabs.harvard.edu/abs/1997ESASP1200.....E

Faherty, J. K., Burgasser, A. J., West, A. A., et al. 2010, AJ, 139, 176

Faherty, J. K., Riedel, A. R., Cruz, K. L., et al. 2016, ApJS, 225, 10

Feigelson, E. D., Lawson, W. A., Stark, M., Townsley, L., & Garmire, G. P. 2006, AJ, 131, 1730

Filippazzo, J. C., Rice, E. L., Faherty, J., et al. 2015, ApJ, 810, 158

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

Fortney, J. J., Marley, M. S., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, AN, 326, 925

Fortney, J. J., Marley, M. S., Saumon, D., & Lodders, K. 2008, ApJ, 683, 1104

Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504

Gagné, J., Burgasser, A. J., Faherty, J. K., et al. 2015, ApJL, 808, L20

Gagné, J., Faherty, J. K., Cruz, K. L., et al. 2015, ApJS, 219, 33

Gagné, J., Lafrenière, D., Doyon, R., Malo, L., & Artigau, É. 2014, ApJ, 783, 12

Galicher, R., Marois, C., Macintosh, B., Barman, T., & Konopacky, Q. 2011, ApJL, 739, L41

Goldman, B., Marsat, S., Henning, T., Clemens, C., & Greiner, J. 2010, MNRAS, 405, 1140

Greco, J. P., & Brandt, T. D. 2016, ApJ, 833, 134

Guarinos, J. 1992, in Astronomy for Large Databases, II, European Southern Observatory Conf. Workshop No. 43, ed. A. Heck & F. Murtagh (Garching: ESO), 301

Hög, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27

Janson, M., Carson, J., Thalmann, C., et al. 2011, ApJ, 728, 85

Kalas, P. G., Rajan, A., Wang, J. J., et al. 2015, ApJ, 814, 32

Kasper, M., Apai, D., Janson, M., & Brandner, W. 2007, A&A, 472, 321

Kasper, M., Apai, D., Wagner, K., & Robberto, M. 2015, ApJL, 812, L33

Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117

Kirkpatrick, J. D., Looper, D. L., Burgasser, A. J., et al. 2010, ApJS, 190, 100

Knapp, G. R., Leggett, S. K., Fan, X., et al. 2004, AJ, 127, 3553

Koleva, M., & Vazdekis, A. 2012, A&A, 538, A143

Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Sci, 339, 1398

Konopacky, Q. M., Rameau, J., Duchêne, G., et al. 2016, ApJL, 829, L4

Konopacky, Q. M., Thomas, S. J., Macintosh, B. A., et al. 2014, Proc. SPIE, 9147, 84

Kraus, A. L., Ireland, M. J., Cieza, L. A., et al. 2014, ApJ, 781, 20

Kuzuhara, M., Tamura, M., Kudo, T., et al. 2013, ApJ, 774, L11

Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660, 770

Leggett, S. K., Saumon, D., Albert, L., et al. 2008, ApJ, 682, 1256

Leggett, S. K., Saumon, D., Marley, M. S., et al. 2007, ApJ, 655, 1079

Liu, M. C., Dupuy, T. J., & Allers, K. N. 2016, ApJ, 833, 96

Liu, M. C., Magnier, E. A., Deacon, N. R., et al. 2013, ApJL, 777, L20

Lodders, K. 2003, ApJ, 591, 1220

Looper, D. L., Kirkpatrick, J. D., & Burgasser, A. J. 2007, AJ, 134, 1162

Luhman, K. L., Patten, B. M., Marengo, M., et al. 2007, ApJ, 654, 570

Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Sci, 350, 64

Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, PNAS, 111, 12661

Maire, J., Ingraham, P. J., De Rosa, R. J., et al. 2014, Proc. SPIE, 9147, 914785

Males, J. R., Close, L. M., Morzinski, K. M., et al. 2014, ApJ, 786, 32

Mamajek, E. E., & Bell, C. P. M. 2014, MNRAS, 445, 2169

Marleau, G.-D., & Cumming, A. 2014, MNRAS, 437, 1378

Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541

Marley, M. S., Saumon, D., Cushing, M., et al. 2012, ApJ, 754, 135

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

Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Sci, 272, 1919
Marley, M. S., Seager, S., Saumon, D., et al. 2002, ApJ, 568, 335

Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556

Marois, C., Lafrenière, D., Macintosh, B., & Doyon, R. 2006, ApJ, 647, 612

Marois, C., Macintosh, B., Barman, T., et al. 2008, Sci, 322, 1348

Marois, C., Macintosh, B., & Véran, J.-P. 2010, Proc. SPIE, 7736, 77361J

Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Natur, 468, 1080

McCarthy, D. W., Burge, J. H., Angel, J. R. P., et al. 1998, Proc. SPIE, 3354, 750

McLean, I. S., & Sprayberry, D. 2003, Proc. SPIE, 4841, 1

Millar-Blanchaer, M. A., Wang, J. J., Kalas, P., et al. 2016, AJ, 152, 128

Milli, J., Hibon, P., Christiaens, V., et al. 2017, A&A, 597, L2

Montet, B. T., Bowler, B. P., Shkolnik, E. L., et al. 2015, ApJL, 813, L11

Mordasini, C. 2013, A&A, 558, A113

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

Morley, C. V., Marley, M. S., Fortney, J. J., et al. 2014, ApJ, 787, 78

Morzinski, K. M., Males, J. R., Skemer, A. J., et al. 2015, ApJ, 815, 108

Nakajima, T., Oppenheimer, B. R., Kulkarni, S. R., et al. 1995, Natur, 378, 463

Naud, M.-E., Artigau, É., Malo, L., et al. 2014, ApJ, 787, 5

Nelder, J. A., & Mead, R. 1965, CompJ, 7, 308

Nielsen, E. L., De Rosa, R. J., Wang, J., et al. 2016, AJ, 152, 175

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJL, 743, L16

Oppenheimer, B. R., Kulkarni, S. R., Matthews, K., & Nakajima, T. 1995, Sci, 270, 1478

Patience, J., King, R. R., de Rosa, R. J., & Marois, C. 2010, A&A, 517, A76

Perrin, M. D., Maire, J., Ingraham, P., et al. 2014, Proc. SPIE, 9147, 91473J

Pickles, A. J. 1998, PASP, 110, 863

Rajan, A., Barman, T., Soummer, R., et al. 2015, ApJL, 809, L33

Rameau, J., Chauvin, G., Lagrange, A.-M., et al. 2013, ApJL, 772, L15

Rayner, J. T., Cushing, M. C., & Vacca, W. D. 2009, ApJS, 185, 289

Rebull, L. M., Stapelfeldt, K. R., Werner, M. W., et al. 2008, ApJ, 681, 1484

Riviere-Marichalar, P., Barrado, D., Montesinos, B., et al. 2014, A&A, 565, A68

Robert, J., Gagné, J., Artigau, É., et al. 2016, ApJ, 830, 144

Rufener, F., & Nicolet, B. 1988, A&A, 206, 357

Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, arXiv:1704.02987

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

Service, M., Lu, J. R., Campbell, R., et al. 2016, PASP, 128, 095004

Sivaramakrishnan, A., & Oppenheimer, B. R. 2006, ApJ, 647, 620

Skemer, A. J., Morley, C. V., Zimmerman, N. T., et al. 2016, ApJ, 817, 166

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Smart, R. L., Tinney, C. G., Bucciarelli, B., et al. 2013, MNRAS, 433, 2054

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJL, 755, L28

Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174

Stephens, D. C., Leggett, S. K., Cushing, M. C., et al. 2009, ApJ, 702, 154

Tokunaga, A. T., Simons, D. A., & Vacca, W. D. 2002, PASP, 114, 180

van Leeuwen, F. 2007, A&A, 474, 653

Wagner, K., Apai, D., Kasper, M., et al. 2016, Sci, 353, 673

Wahhaj, Z., Milli, J., Kennedy, G., et al. 2016, A&A, 596, L4

Wang, J. J., Rajan, A., Graham, J. R., et al. 2014, Proc. SPIE, 9147, 55

Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015, pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophysics Source Code Library, ascl:1506.001

Wolff, S. G., Perrin, M. D., Maire, J., et al. 2014, Proc. SPIE, 9147, 91477H

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Yang, H., Apai, D., Marley, M. S., et al. 2016, ApJ, 826, 8

Yelda, S., Lu, J. R., Ghez, A. M., et al. 2010, ApJ, 725, 331

Yurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, MNRAS, 413, 1828

Yurchenko, S. N., & Tennyson, J. 2014, MNRAS, 440, 1649

Zuckerman, B., Song, I., Bessell, M. S., & Webb, R. A. 2001, ApJL, 562, L87

Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587, A57