Characterizing 51 Eri b from 1-5 $\mu$m: a partly-cloudy exoplanet

TL;DR

This study provides detailed spectro-photometric characterization of 51 Eri b across 1-5 μm, revealing its atmospheric properties, cloud features, and formation scenario, suggesting a partly-cloudy atmosphere and a cold-start formation via core accretion.

Contribution

First comprehensive multi-wavelength spectro-photometry of 51 Eri b, combining new data with models to determine atmospheric composition, cloud structure, and formation history.

Findings

Effective temperature between 605-737 K.

Atmosphere requires patchy clouds for accurate modeling.

Consistent with cold-start formation via core accretion.

Abstract

We present spectro-photometry spanning 1-5 $\mu$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 $\mu$m) and $K2$ (2.10-2.40 $\mu$m) spectra taken with the Gemini Planet Imager as well as an updated $L_P$ (3.76 $\mu$m) and new $M_S$ (4.67 $\mu$m) photometry from the NIRC2 Narrow camera. The new data were combined with $J$ (1.13-1.35 $\mu$m) and $H$ (1.50-1.80 $\mu$m) spectra from the discovery epoch with the goal of better characterizing the planet properties. 51 Eri b photometry is redder than field brown dwarfs as well as known young T-dwarfs with similar spectral type (between T4-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. Model fits suggest that 51 Eri b has an effective temperature ranging between 605-737 K, a solar metallicity, a surface gravity of $\log$(g) = 3.5-4.0 dex, and the atmosphere requires a patchy cloud atmosphere to model the SED. From the model atmospheres, we infer a luminosity for the planet of -5.83 to -5.93 ($\log L/L_{\odot}$), leaving 51 Eri b in the unique position as being one of the only directly imaged planet consistent with having formed via cold-start scenario. Comparisons of the planet SED against warm-start models indicates that the planet luminosity is best reproduced by a planet formed via core accretion with a core mass between 15 and 127 M$_{\oplus}$.

Figures

Tables

Table 1. Table 1: Observations of 51 Eri b

Date	Instrument	Filter	Total Int.	Field	Averaged	Averaged	Averaged
			time (min)	Rot. (deg)	airmass	DIMM seeing    (as)	MASS ${\tau }_0$ (ms)
2015 Jan 30	GS/GPI	$J$aaMacintosh et al. (2015)	70	23.8	1.15	0.52	3.26
2014 Dec 18	GS/GPI	$H$aaMacintosh et al. (2015)	38	37.7	1.14	–	–
2015 Nov 06	GS/GPI	$K1$bbThis work	55	30.5	1.17	0.38	1.56
2015 Dec 18	GS/GPI	$K2$bbThis work	103	71.7	1.22	0.69	0.94
2016 Jan 28	GS/GPI	$K1$bbThis work	97	55.5	1.15	0.86	4.40
2015 Oct 27	Keck/NIRC2	$L_P$bbThis work	100	74.2	1.10	–	–
2016 Jan 02	Keck/NIRC2	$M_S$bbThis work	139	115.7	1.18	–	–
2016 Jan 21	Keck/NIRC2	$M_S$bbThis work	174	116.0	1.21	–	–
2016 Feb 04	Keck/NIRC2	$M_S$bbThis work	148	101.4	1.21	–	–
2016 Feb 05	Keck/NIRC2	$M_S$bbThis work	142	102.1	1.21	–	–

Table 2. Table 2: System properties.

Property	51 Eri A	51 Eri b
Distance (pc)	$29.43\pm 0.29$aaHipparcos catalog (van Leeuwen 2007)
Age (Myr)	$26\pm 3$bbNielsen et al. (2016)
Spectral type	F0IV	T6.5$\pm$1.5
$\log (L/L_{\odot }$)	${0.85}_{-0.07}^{+0.06}$ ccMacintosh et al. (2015) using hot-start predictions.	$-{5.83}_{-0.12}^{+0.15}$ to $-{5.93}_{-0.14}^{+0.19}$ddThis work
$T_{\text{eff}}$	7331$\pm 30$ KeeStellar photometry estimated using SED fit	605–737 KddThis work
$\log g$	3.95$\pm$0.04eeStellar photometry estimated using SED fit	3.5–4.0ddThis work
$J_{\text{MKO}}$	4.690$\pm$0.020ddThis work	19.04$\pm$0.40d,fd,ffootnotemark:
$H_{\text{MKO}}$	4.562$\pm$0.031ddThis work	18.99$\pm$0.21ddThis work
$K_{S,MKO}$	4.546$\pm$0.024ddThis work	18.49$\pm$0.19ddThis work
$K_{\text{MKO}}$	4.600$\pm$0.024eeStellar photometry estimated using SED fit	18.67$\pm$0.19ddThis work
$L_P$	4.604$\pm$0.014eeStellar photometry estimated using SED fit	16.20$\pm$0.11d,gd,gfootnotemark:
$M_S$	4.602$\pm$0.014eeStellar photometry estimated using SED fit	16.1$\pm$0.5ddThis work

Table 3. Table 3: Model grid parameters

Model	Effective	Surface Gravity	Metallicity	Cloud	Cloud Hole
Name	Temperature (K)	[$\log g$] (dex)	[M/H] (dex)	Parameter (f${}_{\text{sed}}$)	Fraction (%)
Iron/Silicate Cloud Grid	600–1000	3.25	0.0	2	0–75
Sulfide/Salt Cloud Grid	450–900	3.5–5.0	0.0, 0.5, 1.0	1, 2, 3, 5	–
Cloudless Grid	450–900	3.5–5.0	0.0, 0.5, 1.0	no cloud	–

Table 4. Table 4: Correlation model parameters

Band	$A_{\rho }$	$A_{\lambda }$	$A_{\delta }$	${\sigma }_{\rho }$	${\sigma }_{\lambda }$
$J$	0.43	0.43	0.14	0.44	0.05
$H$	0.70	0.27	0.03	0.45	0.01
$K1$	0.51	0.41	0.07	0.68	0.004
$K2$	0.30	0.62	0.08	0.43	0.004

Equations

$${\chi }_j^2=\sum _{i=1}^4{\left({\mathbf{S}}_i-{\alpha }_{i,j}{\mathbf{F}}_{i,j}\right)}^T{\mathbf{C}}_i^{-1}\left({\mathbf{S}}_i-{\alpha }_{i,j}{\mathbf{F}}_{i,j}\right),$$

$$\begin{gathered} {\chi }_j^2=\sum _{i=1}^4[{\left({\mathbf{S}}_i-{\alpha }_j{\beta }_{i,j}{\mathbf{F}}_{i,j}\right)}^T{\mathbf{C}}_i^{-1}\left({\mathbf{S}}_i-{\alpha }_j{\beta }_{i,j}{\mathbf{F}}_{i,j}\right) \\ +N_i{\left(\frac{{\beta }_{i,j}-1}{{\sigma }_i}\right)}^2] \end{gathered}$$

$${\psi }_{ij}=\frac{⟨I_i{I}_j⟩}{\sqrt{⟨I_i^2⟩⟨I_j^2⟩}}$$

$${\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}$$

$${\psi }_{ij}\equiv \frac{C_{ij}}{\sqrt{C_{ii}{C}_{jj}}}$$

Full text

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

††journal: AJ††facilities: Gemini:South (GPI), Keck:II (NIRC2)††software: GPI Data Reduction Pipeline (v1.3.0; Perrin et al. 2014), pyKLIP (Wang et al. 2015).

\watermark

DRAFT

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; Wahhaj et al. 2016; Millar-Blanchaer et al. 2016; Bonnefoy et al. 2017) showcase the advances made by these next generation AO systems. Direct imaging, unlike non-direct 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 dependant on complete understanding of the stellar properties. These new AO instruments combine excellent image stability and high throughput with 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 non-equilibrium chemistry when undertaking comprehensive modeling of the exoplanet spectral energy distribution.

In this study we focus on the planetary companions, 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 $\beta$ Pic moving group (Zuckerman et al. 2001). In this study, we adopt an age of of $26\pm 3$ Myr for the $\beta$ 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 e.g. $21\pm 4$ (Binks & Jeffries 2014), $23\pm 3$ (Mamajek & Bell 2014), $20\pm 6$ (Macintosh et al. 2015), $24\pm 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 $\sim$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 Herschel Space Observatory 70 and 100 $\mu$m bands with very low IR luminosity of L${}_{\text{IR}}$/L${}_{\star}=2\times 10^{-6}$ and an a lower limit on the inner radius of 82 au (Riviere-Marichalar et al. 2014) as well as a detection at 24 $\mu$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 $\mu$m) and Keck $L_{P}$ photometry (3.76 $\mu$m), using two different model atmosphere grids to estimate planet properties. While the models agreed on the temperature and luminosity, they were highly discrepant in terms of best fitting 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 clear atmosphere. These discrepancies indicate that more data is required to fully constrain the planet parameters.

In this paper, we present new observations and 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 $\mu$m) and $K2$ (2.10–2.40 $\mu$m) spectrum 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 $\mu$m). In Section 3, we present new near-IR photometry of the star and revise the stellar spectral energy distribution (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. 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 450K to 1000K with deep iron/silicate clouds or sulfide/salt clouds. The 1–5 $\mu$m spectral energy distribution in combination with these two model grids with 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 meridian passage. All the GPI datasets 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, 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 datacubes (Maire et al. 2014). During the process of generating the 3D (x, y, $\lambda$) cubes, the microspectra data are resampled to $\lambda/\delta\lambda$ = 65, and 75 at $K1$ and $K2$, respectively, after which they are interpolated to a common wavelength scale and corrected for geometric distortion (Konopacky et al. 2014). The datacubes 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 (Sivaramakrishnan & Oppenheimer 2006; Marois et al. 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 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 datacube 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 star point spread function (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 full width at half maximum (FWHM), $g=0.5$, and $N_{\delta}=0.5-0.75$ FWHM for the three datasets. 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_{\delta}$ 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$ broad-band 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 broad-band images were used to extract the astrometry of the planet in each dataset thanks to higher signal-to-noise ratio (SNR) 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 the contrast of the planet in each dataset. The same procedure was executed in the non-collapsed wavelength residual images but 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 SNR spectral slices. To measure uncertainties, we injected the template PSF with the measured planet contrast into each datacube 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 Nov 06 and 2016 Jan 28 $K1$ datasets agreed within the uncertainties, the latter having significantly better SNR, and were combined with weighted mean to provide the final planet contrasts.

2.1.1 Spectral covariances

Estimation of a directly imaged planets 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 A.

We use the spectral covariance when carrying out 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 towards cooler temperatures when compared to using the variance directly.

2.2 Keck $L_{P}$

We observed the 51 Eri system on 2015 Oct 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 $\sim$1 hour prior to meridian crossing to maximize the field of view rotation. The target was observed for $\sim$3 hours total, with 100 min of on-source integration. The observations were acquired using the 400mas focal plane mask and the circular undersized “incircle” 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 flatfield and masked hot and bad pixels. As these observations were taken after the April 2015 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$\pm$0.004 mas pixel*-1* and the offset angle $\beta$ ($0.262\pm 0.020$) which is required when calculating the position angle prior to rotating the images to put north up (Yelda et al. 2010). Post-processing of the data was carried out using the Python version of the Karhunen-Loève Image Projection algorithm (KLIP, Soummer et al. 2012; Amara & Quanz 2012), pyKLIP (Wang et al. 2015). As part of this study, we included a NIRC2 module in the pyKLIP codebase that is publicly available for users. 111https://bitbucket.org/pyKLIP/pyklip 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, there are two main parameters that were adjusted, the number of modes used in the Karhunen?-Loève (KL) transform and an exclusion criterion for reference PSFs, similar to $N_{\delta}$ mentioned above, that determines 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 signal to noise 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$\pm$0.15 mag which agrees very well with the photometry in the original epoch, 11.62$\pm$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 Jan 02, 21 and 2016 Feb 04, 05 with Keck/NIRC2 Narrow camera. The details of the observations are presented in Table 1. Each night the target was observed for a period of $\sim$6 hours, as part of two separate NASA and UC Keck observing programs (Program ID - 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.3s with 200 co-adds, 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 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 datasets is the background subtraction. Since the thermal background at 5$\mu$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, a) 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 pre-determined pixel value to align all the images and b) 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 and 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 the $L_{P}$ data. The planet is not detected in each of the individual half-night datasets, requiring a combination of all four half-nights to increase the signal to noise ratio 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 dataset, the planet rotated $\sim$0.48 degrees or $\sim$0.4 pixel, which is a sufficiently large correction that it must be included in the data reduction. Each night’s data was reduced individually to generate 603 PSF subtracted images. These images were then combined by dividing each image into 13 annuli which were combined 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 $\sim$2–3 sigma. 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 planet-star contrast in the $M_{S}$ is 11.5$\pm$0.5 mag.

3 Results

To estimate 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 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_{\text{eff}}$=3.35 $\mu$m, $\Delta\lambda$=1.11 $\mu$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 re-observe the star in the $JHK_{S}$ filters and fit all the available photometry to estimate the photometry in filters where no calibrated stellar data exists.

3.1 Revised Stellar Photometry at $J$,$H$,$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 a $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 Feb 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 airmass 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 pixels correction. Aperture photometry was performed on both targets, with the curve of growth used to select an aperture which 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 relations222http://www.astro.caltech.edu/~jmc/2mass/v3/transformations/ yields $K_{S,\text{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.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 $\sim$10% in the $J$-band and $\sim$15% in the H-band, which is significant given the high SNR of the $H$-band data.

3.2 Fitting the Spectral Energy Distribution of 51 Eri A

To mitigate the effects of incorrect photometry, rather than scale the spectrum in pieces using the relevant broadband photometry, we decided to fit the full SED of 51 Eri A using literature photometry and colors, including Geneva $U,B_{1},B,B_{2},V_{1},V,G$ (Rufener & Nicolet 1988), Tycho2 $B_{T},V_{T}$ / Hipparcos $H_{P}$ (Høg et al. 2000; ESA 1997), MKO $JHK_{S}$ (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 $\sim$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 grid333https://phoenix.ens-lyon.fr/Grids/BT-NextGen/SPECTRA/ (Allard et al. 2012), we estimated the stellar spectrum using a five parameter MCMC grid search. The best fit atmosphere was found with $T_{\text{eff}}=7331\pm 30$ K, $\log g=3.95\pm 0.04$, [$M/H$] = $-0.12\pm 0.06$, and a stellar radius, $R=1.45\pm 0.02$ $R_{\odot}$ (assuming a parallax of $33.98\pm 0.34$ mas; van Leeuwen 2007). No correction for extinction is performed as the extinction in the direction of 51 Eri is negligible ($A_{V}=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$, NIRC2 $L_{P}$ and $M_{S}$ 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_{P}$ photometry

51 Eri b emits a substantial amount of flux in the mid-IR and $L_{P}$ photometry in Macintosh et al. (2015) was used to constrain the effective temperature of the planet. There exists no $L_{P}$ flux measurement for the star and thus they used the $W1$ magnitude reported in the AllWISE catalog ($W1=4.543\pm 0.210$; Cutri et al. 2013), and assumed a color of $W1-L_{P}=0$ based on the F0IV spectral type of 51 Eri (Abt & Morrell 1995). The $L_{P}$ photometry we estimated via the SED fits for 51 Eri is $L_{P}$ = $4.604\pm 0.014$ mag, which is consistent with the value reported in Macintosh et al. (2015) (4.52$\pm$0.21 mag) but with significantly smaller uncertainties.

As a final check for consistency, the 2MASS $K_{S}$ magnitude of 51 Eri ($K_{\text{S, 2MASS}}=4.537\pm 0.024$) was used instead as a starting point. The $K_{S}-L_{P}$ color for early F-type dwarfs and subgiants was estimated by folding model stellar spectra ($7200\leq T_{\text{eff}}/\text{K}\leq 7400$, $4.0\leq\log g\leq 4.5$, $[M/H]=0$) from the BT-Settl model grid through the relative spectral response of the 2MASS $K_{S}$ (Cohen et al. 2003) and NIRC2 $L_{P}$ filters. Over this range of temperatures and surface gravities, the color was calculated as $K_{S}-L_{P}=-0.001\pm 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_{S}-L_{P}=0.014\pm 0.055$. A color of $K_{S}-L_{P}=-0.001\pm 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_{S,\text{2MASS}}$ magnitude of 51 Eri, gives an $L_{P}$ apparent magnitude of $4.538\pm 0.060$. Each estimate for the stellar $L_{P}$ magnitude are within 1-$\sigma$ of each other, and thus we adopt the value derived from the SED fit i.e. $L_{P}$ = $4.604\pm 0.014$ mag.

3.3 51 Eri b Spectral Energy Distribution

We present the final spectral energy distribution 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_{P}M_{S}$ 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 color magnitude diagrams (CMD) for ultracool objects in Figure 5, and compare the photometry of field M, L, and T dwarfs and young brown dwarfs and imaged companions to that of 51 Eri b (red star). The colors of 51 Eri b seems to match the phase space of the late-T dwarfs. To classify the spectral type of 51 Eri b, we do a chi-square 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 sub-sample of the brown dwarfs have corresponding mid-IR photometry and thus 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-square 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 unrestricted fit, we compute the $\chi^{2}$ statistic for the $j^{\text{th}}$ object within the library as

[TABLE]

where ${\bm{S}}_{i}$ is the spectrum of the planet, ${\bm{C}}_{i}$ is the covariance matrix calculated in Section A, and ${\bm{F}}_{i,j}$ is the spectrum of the $j^{\text{th}}$ comparison brown dwarf, all for the $i^{\text{th}}$ 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 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}$,

[TABLE]

where $N_{i}$ is the number of spectral channels in the 51 Eri b spectrum for the $i^{\text{th}}$ filter, and $\sigma_{i}$ is the uncertainty on 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-square 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 (Kirkpatrick et al. 2010; Burgasser et al. 2006; 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, i.e. 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 i.e. T$6.3\pm 1.3$ and $6.1\pm 1.4$ for unrestricted and restricted fits, see Figure 6. We adopt a spectral type for 51 Eri b of T$6.5\pm 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}_{\nu}=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\approx 4.5$; Knapp et al. 2004) and spectroscopic measurements ($\log g=4.7$–$4.9$; Burgasser et al. 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 that of G 204-39 B.

Additional good matches to the 51 Eri b spectrum include 2MASS J22282889–4310262 (2M 2228–43, $\chi^{2}_{\nu}=1.07$ and $1.26$ for the two fits) and 2MASS J10073369–4555147 (2M 1007–45, $\chi^{2}_{\nu}=1.07$ and $1.33$). 2M 2228–43 a well-studied T6 brown dwarf that exhibits spectrophotometric variability in multiple wavelengths 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\pm 2$ pc (Smart et al. 2013). It was identified by Looper et al. (2007) as a low surface gravity object based on its H2O$-J$ vs $K/H$ spectral ratios defined in Burgasser et al. (2006); comparisons against solar-metallicity models imply an age of between 200 and 400 Myrs (Looper et al. 2007).

The best fit object 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, and/or metallicity could be the cause (e.g., Knapp et al. 2004). The poor fit to the color of 51 Eri b is especially 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 hundred of million year old L-type brown dwarf and planetary mass companions as well as the identification of L-dwarf sub-classes 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 similar manner to what was done above for field brown dwarfs. The chi-square 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 similar spectral type 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 mass (10–12 M${}_{\text{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 of the AB Doradus moving group, however Delorme et al. (2017) find that the parallax and kinematics of the free-floating object rule out its membership to 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$\mu$m suggests that the object has low surface gravity, and is most likely a young planetary mass object (2–13 MJup). An alternative solution is that it is a higher mass, 2–40 MJup, 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 under-luminous in the $K$-band. A likely reason for this is that 51 Eri b is much younger than both the comparison companions and thus has the lowest surface gravity amongst 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_{P}$ CMD shows little reddening, which is natural if clouds are causing the effect. The effect of clouds is negligible in the $K$ and $L_{P}$ 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_{P}$ vs $H-L_{P}$ color for the planet and noted that it was clearly redder than the field. Rather than simply being redder than the field T-dwarfs due to the presence of clouds, we present a second possible 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 dependant 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 re-plot 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 over-plot two low mass, 5 and 14 M${}_{\text{Jup}}$, evolutionary model tracks (assuming hot-start conditions) from Saumon & Marley (2008); 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\text{M}_{\text{Jup}}$ planet the L to T transition begins and ends about 1 magnitude fainter in $J$ band than 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.

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 space (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 silicates 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-square 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-square statistic. Note that we use the restricted fit equation in the computation of the best fitting model. This equation permits each of 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 earlier 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 Silicates Cloud Models

In sec. 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 silicates 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 cloud hole fraction and the unique aspect of this grid is the presence of iron/silicate clouds in an atmosphere with clear indications of methane absorption. The clouds are modelled using the prescription presented in Ackerman & Marley (2001), where cloud thickness is parameterized via an efficiency factor (f${}_{\text{sed}}$). Where small values of f${}_{\text{sed}}$ indicate atmospheres with thick clouds while large values of f${}_{\text{sed}}$ are for atmospheres with large particles that rain out of the atmosphere leaving optically thinner clouds. As mentioned early 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, thus for each T${}_{\text{eff}}$, the models went from fully cloudy i.e. f${}_{\text{sed}}$ = 2 and 0% holes to an atmosphere with f${}_{\text{sed}}$ = 2 and 75% holes (patchy clouds). The methodology used to calculated the flux emitted from the patchy cloud atmosphere include both cloud and cloud-free regions simultaneously in the atmosphere using a single, global temperature-pressure profile and are not created via a linear combination of two models as is sometimes done in the literature Marley et al. (2010). The iron and silicates 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 600K to 1000K 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 10pc from the object. One of the free parameters in most model fitting codes is the term R2/d2 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) and thus 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, e.g. 1.33–1.14 R${}_{\text{Jup}}$ for a 2–10 M${}_{\text{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${}_{\text{Jup}}$ and T${}_{\text{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${}_{\text{Jup}}$, and T${}_{\text{eff}}$ $\sim$ 690K, biasing the luminosity of the planet 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 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–750K). 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 CH4 in the atmosphere and could explain higher flux at 1.6$\mu$m and in the $L_{P}$ band. It would also introduce CO, accounting for 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, Na2S, 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 which includes salt and sulfide clouds to test additional parameters such as the surface gravity and metallicity (which were varied, unlike the iron/silicates grid) and the properties of clouds typically associated with T-dwarfs. The grid was designed specifically for lower temperature objects ($450\sim 900$ K Morley et al. 2012, 2014) and has been 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) which 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 range of parameters varied are presented in Table 3, including temperatures, surface gravities, metallicities, and sedimentation factor (f${}_{\text{sed}}$) ranging from cloudy (f${}_{\text{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). Between effective temperatures of 450–775 K, the grid is complete with models available for every step of the varied parameters. For effective temperatures between 800–900K, the temperature steps switch from increments of 25K to 50K and there are no models with f${}_{\text{sed}}$ values of 1 and 2. This grid does not include opacity effect due to iron and silicates condensates. A future series of paper describing an extended atmosphere model grid will describe the updates, however the present grid extends the models to greater than solar metallicites.

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-square in one of four cases, namely, solar and cloudless (top-left), solar and cloudy (top-right), non-solar and cloudless (bottom-left), non-solar 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 the 51 Eri b and the best fitting atmosphere for both solar and non-solar metallicity have very similar reduced chi-square 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${}_{\text{eff}}$ (${605}_{-66}^{+61}$ K) is much cooler in comparison to the iron/silicates grid, but the values are within 2-$\sigma$ of each other. We also note out that the median might not be the best estimate for the effective temperature PDF in the iron/silicates 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 does 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 the iron/silicates 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/silicates 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/silicates 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 be favoring 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 require cloud opacity. This can be achieved in one of two ways: either make the deep iron/silicates 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 under luminous in the $J$ and over luminous $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 $\mu$m SED, a further test of current models. In Figure 15, we show ten of the best fitting models assuming cloudy (sulfide/salt clouds) or cloudless atmospheres extended out to 20$\mu$m. It is clear from these models that observations with the coronagraph on Near InfraRed Camera (NIRCam), spanning the 3–5 $\mu$m wavelength will add significant constraints on the atmosphere of the planet. If the planet can be studied with the Mid Infrared Instrument (MIRI), it could be used to apply constraints on chemical disequilibrium in the atmosphere through observations NH3 in the 10–11$\mu$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/silicates 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, 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 post-formation 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–2 M${}_{\text{Jup}}$. However, for the cold-start case the planet mass could lie anywhere between 2–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 models444http://www.astro.princeton.edu/~burrows/ 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${}_{\text{B}}$/baryon, where k${}_{\text{B}}$ is Boltzmann’s constant, with steps of 0.25 k${}_{\text{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$\times$ solar metallicity (Burrows et al. 2011). The bolometric luminosity of each point in the grid for each of the four atmosphere scenario was computed by integrating the SED over the wavelength range. Because of the sparse sampling of the grid, we linearly interpolate the evolutionary tracks with steps of 0.06 k${}_{\text{B}}$/baryon and 0.2 M${}_{\text{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 super-critical 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*⊕*, 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 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-start 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–127 M*⊕*.

Conversely 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 and/or higher signal to noise data will help to reduce the width of the two branches and independent mass constraints, from dynamical measurements, will enable to infer the initial entropy and possible formation route. Atmospheric retrievals and/or higher resolution spectra aimed at exploring and characterizing the planets chemical composition might also help understand whether the planet has higher C/O ratios compared to 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 Gemini Planet Imager (K1 and K2 bands) as well as the first photometric measurement of the planet at $M_{S}$ obtained with the NIRC2 Narrow camera. We also obtained an additional $L_{P}$ photometric point that agrees very well with the $L_{P}$ 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_{P}$ photometry to present the spectral energy distribution spanning 1–5 $\mu$m for the planet.

As part of the data analysis, we calculated the covariance for each of the spectral datasets i.e. $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 T$6.5\pm 1.5$. Due 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 amongst 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 color magnitude diagrams 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 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/silicates grid, we estimate that the planet has a patchy atmosphere with 10–25 % hole fraction in the surface cloud cover, which is consistent with the f${}_{\text{sed}}$ values of 2–3 resulting from the sulfide/salt grid. The median effective temperature from the two grids is ${737}_{-46}^{+39}$ K and ${605}_{-66}^{+61}$ K for iron/silicates and sulfide/salt respectively. This value is slightly cooler, compared to Macintosh et al. (2015), where the best fit models had temperatures of 700K and 750K 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 planet 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_{P}$ 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 are 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 datasets using a common pipeline for data processing and analysis to understand whether the differences arise from the algorithms or due to other causes.

With future space missions such as the James Webb Space Telescope, the 3–24 $\mu$m SED of this planet could be observed at higher SNR, 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.

Acknowledgements: The authors thank Gabriel Marleau for 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). 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. 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 benefitted from collaborations and/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 (and/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é.

Appendix A Derivation of Spectral Covariance

We follow the method described in Greco & Brandt (2016) to measure the inter-pixel correlation within the PSF-subtracted images, and 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

[TABLE]

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 parametrized correlation model of Greco & Brandt (2016),

[TABLE]

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 do not contribute to the off-diagonal terms of the correlation matrix. The amplitude 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, $\lambda_{c}$ is the central wavelength of the spectrum (1.65 $\mu$m for $H$). The colored lines in the plot are the best fits to Equation A2.

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 is 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

[TABLE]

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ρ) is significantly higher than 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 at $J$ band at wavelengths separated by up to five spectral channels.