Brown Dwarf Retrieval

# Uniform Atmospheric Retrieval Analysis of Ultracool Dwarfs I: Characterizing Benchmarks, Gl570D and HD3651B

## Abstract

Interpreting the spectra of brown dwarfs is key to determining the fundamental physical and chemical processes occurring in their atmospheres. Powerful Bayesian atmospheric retrieval tools have recently been applied to both exoplanet and brown dwarf spectra to tease out the thermal structures and molecular abundances to understand those processes. In this manuscript we develop a significantly upgraded retrieval method and apply it to the SpeX spectral library data of two benchmark late T-dwarfs, Gl570D and HD3651B, to establish the validity of our upgraded forward model parameterization and Bayesian estimator. Our retrieved metallicities, gravities, and effective temperatures are consistent with the metallicity and presumed ages of the systems. We add the carbon-to-oxygen ratio as a new dimension to benchmark systems and find good agreement between carbon-to-oxygen ratios derived in the brown dwarfs and the host stars. Furthermore, we have for the first time unambiguously determined the presence of ammonia in the low-resolution spectra of these two late T-dwarfs. We also show that the retrieved results are not significantly impacted by the possible presence of clouds, though some quantities are significantly impacted by uncertainties in photometry. This investigation represents a watershed study in establishing the utility of atmospheric retrieval approaches on brown dwarf spectra.

12

## 1. Introduction

The spectrum of a brown dwarf opens a series of windows into the depths of its atmosphere, revealing its composition and thermal structure. Differing wavelengths peer to a diversity of depths and are influenced by varying atmospheric species. Teasing out the chemical composition and atmospheric structure of brown dwarfs by synthesizing the information emerging from each window is key to understanding the processes acting in their atmospheres and possibly their evolutionary histories. Historically, empirical analysis techniques, such as color magnitude diagrams, the appearance of spectroscopic features, or comparisons to forward models have been used to understand broad trends in the brown dwarf population. Such analyses have led to the L-T-Y classification schemes as well as widely accepted spectral and photometric indicators of gravity and metallicity (Kirkpatrick 2005; Cushing et al. 2011). However such approaches have important limitations. The physics and chemistry of these dense molecular atmospheres is complex and subtle processes can significantly influcence thermal spectra.

The traditional method for interpreting brown dwarf spectra is to construct sets of sophisticated coupled radiative-convective-chemical equilibrium models to which the data can be compared (Marley et al. 1996; Allard et al. 1996; Tsuji et al. 1996; Burrows et al. 2001). Systematic comparisons of models to data constrain atmospheric composition, temperature, and cloud properties (e.g., Saumon et al. 2006; Cushing et al. 2008; Rice et al. 2010). While these grid modeling studies have unquestionably advanced our understanding of brown dwarf atmospheres, there is still much that we don’t understand. For instance, can dynamical processes in brown dwarf atmospheres cause deviations from radiative convective equilibrium? How do their molecular compositions vary with altitude? Can we measure the molecular abundances rather than assume them based upon chemical equilibrium or disequilibrium? Can their elemental abundances deviate from the often assumed solar composition? It is difficult to answer such questions from comparison solely with forward models, as not all processes can be modeled with fidelity and, in any case, it is not always obvious which processes are responsible for given deviations of models from data.

Line et al. (2014b) presented a novel approach for answering the above questions. Rather then rely on grid models that only constrain a few basic parameters, they used atmospheric retrieval methods common in Earth and planetary atmosphere studies to invert brown dwarf spectra for the temperature structures and abundances. Using such approaches allows the maximum amount of information to be extracted from brown dwarf spectra. Line et al. (2014b) found for, Gl570D, with few assumptions about the nature of the chemistry or the temperature-pressure profile, that the retrieved quantities were consistent with previous grid modeling studies (Saumon et al. 2006), though with some minor deviations such as a more isothermal upper atmosphere and a depleted ammonia abundance. While individual objects themselves are interesting, the key to understanding the underlying physical processes in brown dwarf atmospheres (or any class of atmospheres) requires an investigation of a large population of objects. With large populations one can identify trends and correlations of various physical parameters that can lead to insight into new phenomena.

We first aim to improve upon the techniques presented in Line et al. (2014a) and to validate our improved approach against benchmark brown dwarfs. Then in a followup paper we will apply our new approach to a small sample of T-dwarfs. In that sample we will identify trends within our retrieved results such as metallicity vs. gravity, molecular abundances vs. effective temperature etc., and also how empirical spectral properties such as color or other indices correlate with the retrieved physical parameters. Preliminarily, we choose mid- to late-T dwarfs as they, based upon previous investigations, appear to be largely free of thick silicate clouds that can complicate the interpretation of L and early-T spectra (Kirpatrick 2005).

In this manuscript we present the upgraded retrieval approach and apply it on two benchmark brown dwarfs, HD3651B and Gl570D, using the the SpeX Prism Library (Burgasser et al. 2006a). In 2 we present a modified Bayesian retrieval approach and a novel approach for inverting for temperature structures. In 3 we present our retrieved results and how they are impacted by various assumptions. Finally, in 4 we compare our retrieved values to benchmark properties such as the age and metallicity. Additionally, we present a detailed stellar abundance analysis of Gl570A and HD3651A and derive their carbon-to-oxygen ratios in order to add another dimension to the benchmark comparison.

## 2. Methods

Here we give a brief review of current state of knowledge of Gl570D and HD651B, describe the data we use, the forward radiative transfer model, and the inverse methods to retrieve temperature and abundance information. Building upon Line et al. (2014b), we have made significant upgrades to our methodology in terms of the forward model, Bayesian estimator, and treatment of the data. We highlight those differences where applicable.

### 2.1. Current State of Knowledge on Gl570D & HD3651B

The targets of this study have been selected for their suitability as robust test cases for validating our approach. Both objects are wide-orbit common proper motion companions to stars, allowing us to check our derived properties for consistency with those found for their stellar hosts by other routes. For this reason, such systems are often referred to as benchmarks, although the quality and context of the available constraints varies widely depending the mass and evolutionary phase of the stellar primary (e.g., Pinfield et al. 2006).

Gl570D was the first T-dwarf companion to a star identified following the prototypical T-dwarf Gl229B (Burgasser et al. 2006a). It is a wide component in a hierarchical quadruple system, whose inner components are an M1V+M3V spectroscopic binary (Gl570B and C) and a K4V primary (Gl570A) (Gleise 1969; Duquennory & Mayor 1988; Mariotti et al. 1990; Foreville et al. 1999), from which Gl570D lies at a projected separation of  AU (Burgasser et al. 2000). Gl570D has been subject to a number of grid-model fitting studies, which have to varying extents used the primary star to restrict the parameter space available for the models (e.g., Geballe et al. 2001; Saumon et al. 2006; Legett et al. 2007; Saumon et al. 2012) and has been used as an anchor point for applying trends seen in self-consistent grid models to estimate parameters of the wider T dwarf population (Burgasser et al. 2006b).

Liu, Leggett & Chiu (2007) used a variety of stellar age indicators for Gl570A to constrain the age of the system to the range 1–5 Gyr, whilst Saumon et a. (2006) collated literature values to estimate its metallicity as . In the Appendix we present a new measurement of .

HD3651B was identified as a wide common proper motion companion to the planet-hosting K0V star HD3651A with a projected separation of 480 AU (Murgrauer et al. 2006; Luhman et al. 2007; Liu, Leggett & Chiu 2007). Like Gl570D, it has been the subject of a number of spectroscopic studies that have been constrained by the properties of the primary star (e.g., Leggett et al. 2007; Burgasser 2007). Liu, Leggett & Chiu (2007) reviewed X-ray luminosity, chromospheric and rotation-based age indicators for HD3651A, and found the target lies in the unreliable “older” tail of each of these diagnostics. Like them, we adopt the isochronal age range of 3–12Gyrs from Velenti & Fischer (2005). As an exoplanet host star, HD3651A has been the subject of several recent composition studies. The determinations of are consistently super-Solar, ranging from (Santos et al. 2004) to (Ghezzi et al. 2010). In the Appendix we have further developed these targets potential to contribute to our understanding its substellar companions by compiling detailed abundance measurements from the literature, and measuring new abundances for both. Most significantly for this study we present new determinations of the C to O ratios for both primary stars.

### 2.2. Data

We use the data within the SpeX Prism Library (Burgasser et al. 2006) to perform our analysis on our target objects. Since a given SpeX spectrum is continuous, we avoid having to consider various instrumental systematics (e.g., as had to be done in Line et al. 2014b) and subsequent impact on their interpretation that come along with having to stitch spectra from multiple instruments together. We use the SpeX data taken in the SXD mode which cover wavelengths between 0.8 - 2.5 m with a wavelength dependent resolving power ranging from 87 to 300. The spectra within the SpeX SXD library are oversampled relative to a spectral resolution element, and are therefore each pixel is not an independent sample. Using the full oversampled data would result in over constrained results. We therefore sample every few pixels. We choose the sampling length based upon the auto-corrletion length scale of the residuals of a typical model fit to the data. This is 2.7 pixels. We thus take every other 3rd pixel to be statistically independent. An alternative approach would be to model the covariant error structure of the oversampled data through a Guassian process (e.g., Czlecka et al. 2014).

The native format of the spectral fluxes and error bars within the library are normalized spectra. In order to perform the subsequent analysis, these spectra and error bars must be converted into physical units via photometric calibration. The SpeX database provides the 2MASS photometric J, H, and K-band magnitudes for each object. We convert the 2MASS magnitudes into MKS flux units (W m m) using the Spitzer Science Center Magnitude/Flux density converter 3 which uses the zero point fluxes described in Cohen, Wheaten & Megeath (2003) . The H-band flux is used to derive the final flux-calibrated spectrum just as in Saumon et al. (2006) and Liu, Leggett & Chiu (2007). Figure 1 shows the flux calibrated spectra. Unfortunately, the calibrated spectra are different by some scale factor that is larger than the quoted photometric uncertainty, depending upon which photometric point is used to calibrate. Our retrieval model (discussed below) includes a scaling factor as a free parameter so these differences are not critical to our analysis (with the exception of the derived spectroscopic radius, see the Results section for more on this). We do note, however, that better photometry, namely (MKO) exists for these objects. Furthermore, Stephens & Leggett (2004) suggest that the 2MASS photometry is not the most accurate due to the shape of the 2MASS filter profiles with respect to telluric transmittance. They also provide correction factors for the 2MASS photometry based upon the more accurate MKO photometry. However, we choose to use the un-corrected 2MASS photometry as they shall provide the most conservative impact of inconsistencies in photometry on the derived quantities.

### 2.3. Forward Radiative Transfer Model

The forward radiative model is a derivative of the CHIMERA forward model (Line et al. 2013a; 2014a,b) which computes the upwelling 1-dimensional disk integrated thermal emission spectrum given the molecular abundances, temperature-pressure (TP) profile, and gravity . Near infrared spectra of late-T’s are typically dominated by strong absorption features from water, methane, and alkali metals and little if any obvious absorption exists due to other gases. We include constant-with-altitude volume (molar) mixing ratios for HO, CH, CO, CO, NH, HS, and alkali opacities. These are the species known to be found in cool dwarf atmospheres that have spectral signatures in the near infrared. The alkali opacities include only sodium and potassium and are treated as only one free parameter with their ratio assumed to be solar. Hydrogen/helium in solar ratio is assumed to make up the remainder of the gas.

The TP profile is also included as a free parameter (see §2.4.2). The TP profile is partitioned into 15 evenly spaced slabs (or knots) in log pressure between 315 bar and 1 mbar. Because a 15 layer atmosphere does not have enough vertical resolution to accurately compute fluxes, the 15 level profile is spline interpolated to a finer 70 level grid before before being passed to the radiative transfer. Using fewer TP points permits swifter convergence of the Bayesian estimator.

In addition to the TP profile and gas mixing ratios, we include gravity as a free parameter. Gravity controls the column optical depth. Most of the opacity database is drawn from Freedman et al. (2014) and references therein as in Line et al. (2014b), but we have also incorporated the most up-to-date methane line list (Yurchenko et al. 2014) using the line broadening coefficients from Margolis (1996). Instead of using line-by-line or correlated-K we simply sample the hi-resolution cross sections at 1 cm resolution (see Sharp & Burrows 2007 section 2) which is more than sufficient for moderate resolution spectra.

Finally, the high resolution spectra are convolved with a wavelength dependent Gaussian instrumental profile that reflects the wavelength dependent resolving power, and then binned to the data wavelength grid in order for direct data-model comparison.

Figure 2 shows the sensitivity of a model spectrum at SpeX resolutions to the various parameters. Many of these parameters are sensitive to similar wavelengths. This will result in correlations/degeneracies amongst the gases and the temperature at different levels in the atmosphere.

All of the aforementioned parameters (23 of them) control the flux at the top of a brown dwarf atmosphere. We also include additional “systematic” parameters that facilitate the direct comparison of the model to the data. These parameters account for the radius to distance ratio (and implicitly the flux calibration), uncertainties in the wavelength calibration, underestimation of the spectral error bars, and a smoothing parameter for the temperature profile-for a total of 27 free parameters. These “systematic” parameters will be discussed in more detail in the following section. A list of all of the parameters and a brief description of each is presented in Table 1. The model has many parameters, but this allows us to make very few implicit assumptions about the nature of the atmosphere. Unconstrained parameters will simply appear unconstrained. The beauty of modern retrieval approaches (below) is that they can accommodate numerous parameters and fully account for all of the correlations amongst them. A larger number of parameters will of course result in more conservative estimate of the uncertainties through marginalization.

### 2.4. Retrieval Model

#### Bayesian Implementation

The retrieval model is the Bayesian engine that optimizes the forward model to fit the data. We use the Markov chain Monte Carlo approach implemented with affine-invariant ensemble sampler, EMCEE (Foreman-Mackey et al. 2013). This is a significant advancement over the optimal estimation and bootstrap Monte Carlo approaches used in Line et al. (2014b) as we are now able to make fewer a priori assumptions about the smoothness of the temperature profile or the Gaussian shape of the parameter uncertainties. EMCEE requires only a functional form for the log of the posterior probability to perform the optimization. The posterior probability is a combination of the likelihood and the prior described as follows. Starting from Bayes theorem

 p(x|y)=L(y|x)p(x)E (1)

where x is the parameter vector described in 2.3 and y is the data vector–in our case the spectrum, is the posterior probability distribution, is the likelihood distribution which penalizes poor fits to the data, is the prior which represents any external constraints, and is a normalization factor known as the evidence, or marginal likelihood, which is required for Bayesian model comparison but not for parameter estimation. We use the following log-likelihood function:

 lnL(y|x)=−12n∑i=1(yi−Fi(x))2s2i−12ln(2πs2i) (2)

Here, the index denotes the data point, in our case some property at a single wavelength bin, is the measured flux, is the modeled flux that comes out of the forward model (2.3), and is the data error given by

 s2i=σ2i+10b (3)

where is the measured error for the data point and is a free parameter. Differing from Line et al. (2014b), we modify the standard error on the data point by the factor to account for underestimated uncertainties and/or unknown missing forward model physics (Foreman-Mackey et al. 2013, Hogg et al. 2010, Tremain et al. 2002), e.g., imperfect fits. This results in a more generous estimate of the parameter uncertainties. Note that this is similar to inflating the error bars post-facto in order to achieve reduced chi-squares of unity, except that this approach is more formal because uncertainties in this parameter are properly marginalized into the other relevant parameters. Generally, the factor takes on values that fall between the minimum and maximum of the square of the data uncertainties. The first term inside the summation in equation 2 is the familiar “chi-square”. This term penalizes large residuals. The second term in the summation is the Gaussian normalization factor that is normally excluded from standard fitting routines due to the unchanging data errors. Because the data errors include the free parameter this normalization can change, and hence has to be taken into account. Really, the purpose of this term is to provide a balance for the error bar inflation parameter to prevent it from approaching infinity.

The prior, , can be broken up into several pieces as

 p(x|y)=p(T)p(x′)p(γ) (4)

where is the prior on the log of the gas mixing ratios, the instrumental parameters, gravity, and the radius-to-distance scaling while and are the temperature profile priors. The parameter is the relative weighting of the temperature prior (see equation 5) . The prior details are shown in Table 2.

Because we have both measured fluxes and parallaxes for each object, we are able to calculate the “photometric” radius. If we can measure both radius and a gravity we can then constrain the mass. We know for brown dwarfs the mass cannot exceed (e.g., Burrows et al. 2001). Therefore, rather than place individual priors on the radius and gravity we enforce a prior on the derived mass to fall between the physical plausible values of 1 and 80 . This constraint prevents the retrieved radii and gravities from entering an unphysical region of parameter space.

#### A Novel Temperature-Pressure Profile Retrieval Approach

We present our novel method for retrieving temperature profiles in atmospheres. A common issue in planetary atmospheric retrievals is how to parameterize the temperature profile in an atmosphere. There are two philosophies. One philosophy, mostly used in the exoplanet atmosphere community when the data is sparse, is to parameterize the atmosphere with some analytic function that can be described by a small set of parameters (e.g., Madhusudhan & Seager 2009; Line et al. 2012; 2013 Benneke & Seager 2012). This is advantageous as the entire TP structure can be controlled by just a few simple parameters. It is disadvantageous because it is relying on a to infer the temperature structure. This could potentially result in biases in the retrieved profiles. For instance, if one were to use the simple Eddington approximation for the temperature profile, there would be one free parameter–the mean opacity. While just one free parameter is ideal, we know that the Eddington approximation is a poor approximation for brown dwarf atmospheres because the true opacity is not constant with pressure, nor is it gray. While parameterizations are appealing in their simplicity, they can often times be too much of an oversimplification of the physics, and are thus not appropriate.

The classic planetary science approach, the other extreme, is to retrieve the temperature at each model layer in the atmosphere (e.g., Rodgers 2000; Irwin et al. 2008; Lee et al. 2012) within an optimal estimation framework. For typical model atmosphere grids this can be anywhere between 50-100 independent temperature-pressure points. This was the approach used in Line et al. (2014b). Because many of the atmospheric levels are degenerate, the retrieved profiles often result in unphysical oscillations, or ringing (Rogers 2000). The standard remedy to this problem is to implement some covariance matrix, or a smoothing kernel (or Tikonov Regularization), given some set smoothing length scale and width (Irwin et al. 2008). While this reduces wild oscillations, this smoothing length scale and width must be chosen and cannot change during the course of a retrieval. Furthermore, these values are often case specific and must be tuned by hand. Our novel approach remedies all of the aforementioned issues, and can be readily implemented within a Bayesian framework.

We borrow much of this work from the non-parametric regression literature (Lang & Brezger 2004; Rahman 2005; Jullion & Lambert 2007). The goal is to allow flexibility to fit for each of the independent temperature-pressure points while preserving smoothness, without having to a priori set the degree of smoothness. The best way to do that is by penalizing the second derivative of the temperature structure. The second derivative is the “roughness” of a function. Temperature vectors with wild, unphysical oscillations, will have large summed second derivatives. These types of roughness-penalized non-parametric polynomial fits are known as P-splines. The degree to which this roughness is penalized is included as a free parameter (). This variable smoothing is implemented as,

 lnp(T)=−12γN∑i=1(Ti+1−2Ti+Ti−1)2−12ln(2πγ) (5)

Inside the sum is the discrete second derivative of the temperature profile at each level, weighted by . Based on experimentation (Lang & Brezger 2004; Rahman 2005; Jullion & Lambert 2007) the hyperprior on should take the form of an inverse gamma distribution with the properties shown in Table 2. A variable allows the data to dictate the degree of smoothing. If the data warrants little smoothing and there truly are oscillations in the TP profile, then will be large, lending little weight to the smoothing prior. When the data does not justify rough TP profiles, will be small resulting in a larger penalty to rough profiles. As mentioned in §2.3, 15-knots, or anchor points are used in our TP profile. This number is somewhat arbitrary as the number and location of the knots should not matter within this framework just so long as there are enough evenly spaced knots to sufficiently resolve potential structure (Eilers & Marx 1996). We have tested higher numbers of knots (30 vs. 15) and have indeed found little sensitivity to the choice.

The combined log-likelihood and log-prior are readily implemented within the EMCEE sampler. The MCMC is initialized with a tight Gaussian ball with 8 chains or walkers (Foreman-Mackey et al. 2013) per parameter about an educated starting guess in order to minimize burn in time (the time it takes for the MCMC sampler to locate the posterior distribution). We note that the final results are insensitive to the initial starting point; poor starting points result in longer burn in times. Convergence is monitored with the Gelman-Rubin statistic. The statistics are summarized by drawing thousands of random points from the last (which we take to be the posterior) of the ensemble of chains which generally entail 500-1000 independent samples (as dictated by the autocorrelation length scale). In the following sections we describe the physical properties of the two benchmark late T-dwarfs evaluated within this framework.

## 3. Retrieval Results

The fiducial retrieval results are summarized in Figures 3 and 4. The top row of Figure 3 shows an ensemble of thousands of fits derived from the posterior and their residuals. We also summarize the retrieved and . The residuals appear to be minimal and random, with exception of the peaks near Y and J band. The slight misfit at these wavelengths could potentially be due to uncertainties in the alkali (sodium, potassium) cross sections, which have the least certain cross-sections of the opacities that absorb at these wavelengths. We note that these residuals are much smaller than can be obtained by typical grid model fits suggesting that the additional parameters we include in our model are required to produce these better fits.

The retrieved temperature profiles (bottom row, Figure 3) for each object are summarized with a median, 1, and 2 credibility region from the ensemble of thousands of randomly drawn TP profiles from the posterior. The retrieved TP profiles appear to be consistent with physically based expectations. In each panel we show for comparison a self-consistent grid model profile (black, from the grid models of Saumon & Marley 2008) corresponding to the median retrieved effective temperature and gravity. The self-consistent grid models assume cloud-free 1-dimensional radiative convective equilibrium with a solar composition atmosphere in thermochemical equilibrium. The agreement is astounding and this is a point we would like to stress. The self-consistent grid model profiles generally fall well within the 2 credibility region. We have made no assumptions about the nature of the TP profiles other than smoothness, the degree of which was allowed to vary. We also tested different starting guesses (e.g., isothermal) and the results are no different. This suggests that these are the actual, true, temperature-pressure profiles in these atmospheres and that the assumption of one-dimensional radiative-convection is sufficient. However, there is some divergence at pressure levels less than about 1 bar and greater than a few 10s of bars. The retrieved profiles tend to become more isothermal near the top of the atmosphere than the radiative convective models, though this divergence is much smaller than the width of the confidence intervals and is therefore not significant (as constrained by SpeX data alone). We also note that a recent paper by Tremblin et al. (2015) predicts that condensation-induced fingering convection can result in a cooler deep atmosphere than expected from standard dry convection assumptions, consistent with what we are finding but not conclusive. Higher resolution data with more vertical resolution and altitude range or longer wavelength data will (and have, as in Line et al. 2014b) provide better constraints to the TP profile that will allow us to further test deviations from the standard radiative-convective assumptions. In this investigation, we purposefully avoid combining different datasets in this investigation due to the introduction of additional systematics that come with multiple data sets.

Figure 4 summarizes the posterior for the physical atmospheric parameters. The effective temperature, radius, mass, metallicity (), and C/O ratio are all derived quantities and were not directly retrieved. The effective temperature distribution is obtained by computing the bolometric flux (1-20 m) over thousands of spectra generated from the posterior. The radius is derived from the retrieved spectral scaling factor, , given the distance. The distance and photometric uncertainties are formally propagated into the radius uncertainty via Monte Carlo error propagation. The mass is derived from the retrieved gravity and radius. We also reiterate that the radius, and hence the derived mass, are very sensitive to the accuracy (e.g., missing systematics not accounted for in the quoted photometric errors) of the photometry. We discuss the impact of the photometric calibration on the retrieved quantities further in 3.3. The metallicity is derived by summing up the molecular mixing ratios for each species weighted by the number of metal atoms divided by the abundance of hydrogen and then comparing that to the sum of solar metals relative to hydrogen. The C to O ratio is computed by dividing the sum of the carbon bearing species by the oxygen bearing species appropriately weighted by the number of carbon/oxygen atoms in each species.

The marginalized probabilities for each parameter (the histograms) are shown along the diagonals for each object. For both objects the HO, CH, NH and Na+K are constrained (68% confidence) to better than 0.3 dex (or a factor of 2). For comparison the best constraints we have on gases in exoplanet atmospheres are to within a factor of 10 (e.g., Kreidberg et al. 2014).

Perhaps the most surprising finding is the tight constraint on ammonia for each object. The robustness of the ammonia constraints in the near-IR are discussed in 3.1. Additionally there are strong correlations of log with spectrally prominent absorbers, HO, CH, and NH . As expected for any atmosphere in hydrostatic equilibrium, there is a positive correlation between metallicity and gravity (this can also be seen in the metallicity vs. gravity panels in Figure 4). As gravity increases, the optical depth at a given layer in the atmosphere decreases (, where is the opacity), so the opacity must increase to maintain that same optical depth. The strong absorbers, HO, CH, and NH, are also correlated with each other and with temperature (not shown). These correlations are positive, which seems backwards for overlapping absorbers, but it is because as one absorber increases, the temperature first increases to maintain that same flux at the wavelength of that absorber, resulting in a higher flux at a different wavelength where another absorber is present. This absorber must increase in abundance to suppress the flux at this different wavelength. The other gases, CO, CO, and HS, are largely unconstrained; only upper limits can be obtained. This is mainly because of their relatively low thermochemical abundances (despite potential vertical mixing) and relatively weak bands in the near infrared. We would expect CO and CH to flip roles for objects hotter than 1100 K.

In order to check whether or not the retrieved abundances are chemically realistic, we compare them to thermochemical equilibrium models. To do this, we use the NASA Chemical Equilibrium with Applications model (CEA, Gordon & McBride 1996) recently used for exoplanet atmosphere studies (Line et al. 2010; Moses et al. 2011; Line et al. 2011). CEA only requires the local temperature, pressure, and elemental abundances at given model layer in order to compute the equilibrium abundances. Equilibrium condensation chemistry (no rainout) is included but the code has a difficult time with enstatite (MgSiO) condensation. For this reason we do not include magnesium or silicon species, but account for the depletion of oxygen due to presumed cloud formation in the deep atmosphere by removing 3.28 oxygen atoms for every silicon atom (Burrows & Sharp 1999). We have not accounted for perturbations to equilibrium chemistry such as horizontal or vertical mixing.

Figure 5 compares the retrieved results for the well constrained species to thermochemical equilibrium abundances along the median temperature profile. In principle there will be a spread, albeit minor, due to uncertainties in the temperature profile, but since the goal is to just check for realism in the abundances we ignore such a spread. For each object we show two cases of equilibrium abundances. The first is for solar elemental composition (solid lines) and the second is a by-hand fit of the metallicity and C to O ratio to the retrieved quantities. We stress that the metallicity and C/O are different from the metallicity and C/O. Condensate processes can deplete oxygen or other species in the atmosphere resulting in metallicity and C/O’s that differ from the or bulk values.

We find for both Gl570D and HD3651B that the assumption of solar elemental abundances overestimates the retrieved water and methane abundances but appears to do a good job for the alkali and ammonia mixing ratios. However, by hand-tuning the bulk metallicity and C/O in the thermochemical model we can better match the retrieved water and methane abundances. This is not a rigorous “fit” to the chemistry by any means, but simply an attempt to show that we have retrieved chemically plausible molecular abundances. A perhaps more rigorous approach for obtaining the allowed ranges of C to O ratios and metallicities would be to perform a “retrieval on the retrieval” where by the chemical model would fit the retrieved molecular abundances within a Bayesian framework. This is currently beyond the scope if this work. We also note that the ammonia thermochemical profiles agree well with our column averaged uniform-with-altitude retrieved values. Saumon et al. (2006) suggest that quenching of ammonia due to vertical mixing occurs in the deep atmosphere near a temperature of 2200 K. Such temperatures occur at the deepest pressures (several hundred bars) on our retrieved profiles. If ammonia indeed quenches at these deep levels then the ammonia profile would be nearly constant with altitude well within our retrieved range for both objects (Figure 5).

Another surprising find is that the retrieved water/methane abundance is lower (water/methane 0.8-0.9) than what one may expect for a solar composition atmosphere. Typically water is more abundant than methane by a factor of 1.5 at solar composition (see Saumon et al. 2006). This suggests a super-solar (greater than 0.5) atmospheric carbon-to-oxygen ratio. Accounting for the draw down of O due to silicate condensation, the methane/water abundance shown in Saumon et al. (2006) for Gl570D suggests an C/O of 0.63. Our retrieved carbon-to-oxygen ratios for both objects are higher than one. The inferred, by-hand C to O ratios are less than unity but are still slightly higher than solar (see §4 for a comparison to the host star values). We note that, in our previous study, Line et al. (2014b), we found a fairly low (0.2) carbon-to-oxygen ratio for Gl570D. These differences are likely due to the inclusion of the Akari and Spitzer Infrared Spectrometer data and the treatment of the systematics between them. Our current investigation is more straightforward as we only focus on the SpeX data set, and thus do not have to worry about potential biases due to differing unaccounted for systematics amongst different datasets. In 4 we show that our thermochemically self-consistent carbon-to-oxygen ratios for the brown dwarf’s are consistent with the those derived from the stellar primaries.

### 3.1. Ammonia in the Near-IR

One of the more remarkable findings in this investigation is the strong evidence for the presence of ammonia in low resolution near infrared spectra. Ammonia at longer wavelengths in T-dwarfs is not new. Cushing et al. (2008) and Saumon et al. (2006) convincingly demonstrated the presence of ammonia in Gl570D using the strong 9.6 m band in the Spitzer IRS data. However, spectroscopic features of ammonia are not expected to present themselves below 2.5 microns until the Y-dwarfs (Kirpatrick et al. 2005) unless observed at high resolution (Canty et al. 2015). Therefore, we were surprised to find how strongly ammonia can be constrained (e.g., an actual bounded limit as opposed to an upper limit only) with -resolution near infrared data in both objects despite the lack of obvious spectral features in this wavelength range. How are we to believe that our constraint is real? We show three lines of evidence supporting our strong ammonia constraint.

One line of evidence comes from a standard Bayesian hypothesis testing procedure. Such procedures determine whether or not a parameter within nested models is justified given the data (e.g., Trotta et al. 2008). A commonly used approach is to compare the Bayesian information criterion (BIC) between a model with and without a particular parameter. Parameters that provide better fits to the data and produce a delta-chi-square that is greater than the penalty of adding that additional parameter are justified. However, the BIC is a truncated Laplace approximation to the full Bayesian evidence, or marginal likelihood (Kass & Raferty 1995). We do not use the BIC here, rather we compute the full Bayesian evidence. This is done by numerically integrating over the entire posterior using the approach described in Weinberg et al. (2012) (see also Swain, Line & Deroo 2014 for an application to exoplanet spectra). By computing the evidence of the full model that contains all parameters to the one that removes ammonia, we can obtain a Bayes factor. Bayes factors greater than one suggest that the model containing the parameter in question is favored, while Bayes factors less than one suggest otherwise. A Bayes factor can then be converted into a confidence of detection (Trotta 2008). Table 3 shows the Bayes factors and the corresponding detection significances for three different nested models each removing only one gas (NH, HS, or HO) from the full model in Table 1.

We show the detection significances for HO as an example of a gas that is visibly obvious in the near infrared and well constrained (Figure 4), and hence we would expect an extremely high detection significance. We detect water at at an extremely high degree of confidence, 17 , in both objects. HS is an example of a poorly constrained species (Figure 4), and hence would expect, and indeed do find, a low detection significance below 2 . In fact for, HD3651B the Bayes factor is less than one ((Bayes factor) 0) or evidence against HS. Since we don’t visibly see any obvious spectroscopic features due to NH, but do indeed obtain a strong constraint (the marginalized posterior is bounded on both sides as opposed to an upper limit like HS), we may expect the detection significance to fall in between the two aforementioned extremes. This is what we do indeed find, a 6 detection of NH in both objects, which is considered strong. We should note that Bayes factors can be sensitive to the prior ranges. We found that in our case, the Bayes factor calculation is insensitive to our prior ranges.

We also show two additional, more straight forward, lines of evidence in Figure 6. For the first test we re-ran the retrieval but initialized the MCMC with a a non-detectable ammonia abundance far from the retrieved value. If the retrieved value we obtain is true, then we would expect the ensemble of Markov chains to converge towards the true value regardless of the starting point–given a long enough run time. That is indeed what we find. The left two panels in Figure 6 show the evolution of the Markov chains. They readily rebound from the poor initial starting point (essentially no ammonia) and converge to a nice tightly packed bundle within the target distribution. Less than 10 of the chains remain outliers. Had we run for even longer (weeks perhaps) these chains would likely have fallen in line with the rest and too converge to the true answer.

For the final test for ammonia we create two synthetic brown dwarf spectra for which we know the true TP profile and composition. We choose parameter values and data properties similar to Gl570D. We create two synthetic spectra: one with ammonia and one without. We then apply the retrieval to these two spectra. What is shown in the right panel Figure 6 are the retrieved ammonia distributions for these two scenarios. The blue histogram is the retrieved ammonia distribution for the synthetic spectra generated ammonia. In this case only an upper limit of ammonia can be obtained. This means that there are no spectral features present in the synthetic spectrum to suggest the existence of ammonia. When the ammonia abundance creeps up to a high enough values (in this case 1ppm) it begins to present itself in the spectrum in an undesirable fashion. The red histogram is the retrieved ammonia distribution for the synthetic spectra generated ammonia. This histogram is nicely bounded on both sides suggesting that the spectral features due to ammonia are enough to place both a lower and upper bound on the retrieved abundance. These three lines of evidence strongly suggest the presence of ammonia in the the low-resolution near infrared spectra of these two T-dwarfs.

### 3.2. Verifying Cloud Free

T-dwarfs are typically assumed to be cloud free given their blue colors (Burrows et al. 1997; Allard et al. 2001) and the success of cloud-free grid models to reasonably explain their spectra (e.g., Stephens et al. 2009), though some of the redder objects are better matched with models that include sulfide-like clouds (Morley et al. 2012). Since we assumed cloud free atmospheres for our retrieval conclusions, we need to make sure that this is indeed the case. We are not interested in the cloud properties themselves, rather the impact that some unaccounted for gray absorber may have on the spectra. Therefore we model clouds rather simplistically as a gray absorber with opacity between two defined pressure levels, and . The cloud optical depth is

 τc=κcPc,bottom−Pc,topg (6)

We assume spans one scale height so that we only need to retrieve the cloud base location () and . The gray approximation can be readily justified. From grid model investigations with clouds (e.g., Morley et al. 2012) high sedimentation values (Ackerman & Marley 2001) are required to best match the spectra. High sedimentation values generally result in a wider range of particle sizes thus washing out Mie scattering features resulting in gray, or at least, nearly gray absorption.

In order to determine whether or not clouds are present, we undergo the same Bayesian hypothesis testing procedure described in 3.1. This time the full model includes all of the parameters in Table 1, plus the two cloud parameters. We compare the evidence of original cloud free model to this new full model. Table 4 shows the Bayes factor and detection significance for the cloud. We find that for Gl570D the detection significance is below 2 suggesting a weak detection of clouds. HD3651B presents a slightly higher, or weak to moderate detection of a cloud. Figure 7 shows the parameter distributions for both the cloudy and clear atmospheres; including the clouds has an insignificant impact on all of the retrieved parameter values (e.g., the change in the median value is less than the typical width of the distributions). Therefore, we are justified in assuming cloud free atmospheres for these two objects.

### 3.3. Impact of Photometry

The flux calibrated spectra depend on the choice of photometry used as demonstrated in Figure 1. Here, we explore the impact on the retrieved quantities of the choice in photometry. For each object we calibrate the spectrum with either the J-band, H-band, or K-band 2MASS photometry, as shown in Figure 1. As mentioned earlier, Stephens & Leggett (2004) provide correction factors for the 2MASS photometry. We do not apply those correction factors here, rather the goal is to determine what effects “bad” photometry may have on the retrieved quantities. We then execute the retrieval on each of those three calibrated spectra. Figure 8 shows the resulting retrieved quantities for each of the photometric calibration scenarios. The impact is minimal for most quantities (e.g., the shift in the median is well within the 1 uncertainties) with the exception of the photometric radius. This is unsurprising as the overall scaling to the spectrum depends on . Shifts in this scaling due to photometry will result in changes in the derived radius. We also find small ( 1) shifts in the retrieved gravity due to the prior upper limit on the mass (masses cannot exceed 80 M). This is because the mass depends on both the radius and gravity therefore the small shifts in the radius propagate through the mass upper limit to the gravity. If the radius increases due to a change in photometry, then the gravity has to decrease. These shifts in derived radius and gravity stress the importance of precision photometry on these objects.

## 4. Validating Retrieval with Benchmarks

Both Gl570D and HD3651B happen to orbit stars with known properties. This makes them powerful benchmark systems for which we can test the validity of our retrieval approach.

The basic properties we use to evaluate the benchmark systems are their evolutionary derived ages, metallicities, and the carbon-to-oxygen ratios. A fundamental assumption with benchmark systems is that the primary and companion both formed out of the same nebular material at the same time and should each individually indicate the same elemental abundances and age. This is unlike planetary systems in which we believe that planet formation processes within the protoplanetary disks can alter the planetary atmosphere abundances relative to their host star (Oberg et al. 2011; Fortney et al. 2013).

A summary of the relevant stellar primary and retrieved brown dwarf properties are shown in Table 5. Since stellar photospheres are generally assumed to be well mixed and free from condensates, their photospheric abundances are representative of their values, however for the brown dwarfs we retrieve the quantities rather than the quantities, which can be different due to the aforementioned reasons. Therefore, in Table 5 we show both the elemental quantities and the thermochemically self-consistent elemental quantities derived from the hand tuned fits (Figure 5 ) of the chemical model to the retrieved abundances. Typically the atmospheric oxygen is depleted by 20-30 relative to the due to the sequestration of oxygen in condensates (depending on the metallicity and C/O). This results in a higher carbon to oxygen ratio and an overall lower atmospheric metallicity, since oxygen is the dominant metal atom. In Table 5 we list the metallicity and C/O based on the atmospheric measurements and the estimated depletion of O and silicates due to condensation. These are the nominal values against which we compare to the stellar abundances. We find that the metallicities in both objects are somewhat lower than what is measured in the primaries but are still consistent. From the C/O measurements in both the brown dwarf (A) and host star, we have an additional benchmark constraint. We find that for the Gl570 system, the 1 inferred intrinsic C/O range falls entirely within the stellar primary C/O values. We consider this an excellent agreement. For the HD3651 system, the inferred intrinsic C/O is higher than the stellar primary by 30%, however their medians are consistent at the 2 level. This represents, for the first time, a comparison of a stellar and companion brown dwarf carbon-to-oxygen ratios. This suggests that we are in fact retrieving the proper molecular abundances in the brown dwarf atmospheres.

Finally, we compare the evolution-derived age to the estimated system age. Figure 9 shows the Saumon & Marley (2008) isochrones in log -effective temperature space. The shaded regions represented the estimated system ages as described extensively in Liu, Leggett, & Chiu (2007). Our retrieved gravity and effective temperature are consistent with the presumed system ages. We also obtain photometric masses from the retrieved gravity and photometric radii. We find that the 1 range in photometric masses for Gl570D (15-58 M) is consistent with the evolution model masses (Figure 9); however for HD3651B we find somewhat higher photometric masses (45-78 M) than anticipated from the evolution models. It is unclear why this may be.

In summary, we find ages derived from our retrieved gravity and effective temperatures are consistent with the measured system age and that the retrieved metallicities and C to O ratios, after taking into account the loss of oxygen due to condensates, are in good agreement with the primaries.

## 5. Summary & Conclusions

We have established a new minimal-assumption retrieval approach for brown dwarf atmospheres that significantly advances the work of our previous retrieval study (Line et al. 2014b). The new retrieval approach relies upon a more robust Bayesian estimator, forward model, temperature profile parameterization, a single continuous spectrum, and treatment of unknown systematic uncertainties permitting generous uncertainty estimates. Differences in our results compared with Line et al. (2014b) are likely due to these major changes. From our new approach applied to two benchmark late T-dwarfs, Gl570D and HD3651B, we determined the allowed range of the thermal structures, molecular abundances, gravities, and radius-to-distance scalings directly from the data. We found that this parameter set provides very good fits to the data in the form of minimal, nearly random, residuals. We validated the chemical plausibility of the retrieved molecular abundances using a a well vetted thermochemical equilibrium model.

Perhaps the most significant highlight of our work is the robust detection of ammonia in the low resolution near-infrared spectra in these late T-dwarfs. We presented three lines of evidence to support this claim. Furthermore, we showed that clouds play a minimal role in sculpting the spectra of these two objects and their inclusion had little to no influence on the other parameters. We also suggested that large systematic uncertainties in photometry can result in biased estimates of the photometric radii.

An additional highlight, is for the first time, using the carbon-to-oxygen ratio of the host-companion as an extra dimension in establishing benchmark systems. We found a remarkable agreement of the carbon-to-oxygen ratios derived from a stellar abundance analysis for Gl570A and our retrieval analysis on Gl570D, and a consistent agreement within the HD3651 system (considering the spread in literature C/O values). This is quite the accomplishment as two completely separate techniques on two different objects for which we would expect similar abundances, are in good agreement. This further bolsters the suggestion of Fortney (2012) to explore the role of C/O in T-dwarf atmospheres and that the C to O ratio should be considered as an additional dimension when interpreting brown dwarf spectra. Finally, we found that the ages derived from the evolution models and our retrieved gravity and effective temperatures are consistent with the estimated system ages. It would be interesting in future investigations to identify systems for which the companion-primary C/O and metallicities differ by a significant amount. This could point to new physics and/or chemistry operating in substellar atmospheres.

This investigation further establishes the power of our novel retrieval approach in understanding the atmospheric and bulk properties of brown dwarfs. In a future study we plan to apply this technique to a wider range of objects with the goal of identifying trends in the thermal structures and molecular abundances and how they correlate with empirical metrics. This will undoubtably verify hypothesized physical and chemical mechanisms operating in brown dwarf atmospheres, and likely identify unknown ones as well.

## 6. Acknowledgements

The authors would like to thank Adam Burgasser, Brendan Bowler, Kelle Cruz, Mike Cushing, Michael Liu, and Emily Rice for useful discussions on benchmark systems, data treatment, and various data-model comparison approaches. The authors thank Richard Freedman and Roxana Lupu for providing gas opacities and Caroline Morley for radiative transfer code comparisons and helpful discussions. We thank Jacob Lustig-Yeager and Kyle Luther for re-writing portions of the code in python and C for significant speed improvements and also Dan Foreman-Mackey for making EMCEE available to the community. Finally, we thank the anonymous referee and statistics consultant for useful and insightful comments. JT acknowledges financial support from the Carnegie Origins Postdoctoral Fellowship Program. BB acknowledges financial support from the European Commission in the form of a Marie Curie International Outgoing Fellowship (PIOF-GA-2013-629435). JF acknowledges funding support from NSF award AST-1312545. MM acknowledges support from the NASA Astrophysics Theory and Planetary Atmospheres programs

## Appendix A Appendix

Comparing the carbon-to-oxygen ratio in stellar-brown dwarf companion systems provides an additional benchmark dimension. Both brown dwarfs and the stars that they orbit are presumed to form out of the same molecular cloud, and thus would be expected to have the same elemental abundances. Metallicity and age are usually the benchmark dimensions. Additional dimensions provide more constraints on the system properties. Determining stellar abundances is no easy task as different groups using different techniques with different data on the same objects often times report significantly different results (Hinkle et al. 2014). In this section we provide our own analysis to determine the stellar carbon-to-oxygen ratios in Gl570A and HD3651A.

### a.1. Observations and Stellar Parameter Analysis

#### Gl570A

The observations of the K4 dwarf Gl570A were conducted on 13 July 2014 (UT) with the Magellan Inamori Kyocera Echelle (MIKE) spectrograph (Bernstein et al. 2003) on the 6.5m Landon Clay (Magellan II) Telescope at Las Campanas Observatory. Three frames of 100s each were taken of the target with the 0.5”x5” slit and 1x1 binning. On 12 July 2014 (UT), three frames of 500s each with the 0.35”x5” slit and 1x1 binning were taken of Vesta, as a solar standard. On both nights calibrations (biases, quartz and milky flats, and ThAr lamp spectra) were taken at the beginning of the night. MIKE is a double echelle spectrograph, meaning a dichroic splits the light into blue (3350-5000 Å ) and red (5000-9400 Å) arms. The data were reduced, extracted, combined, and wavelength calibrated with the Carnegie Python Distribution (CarPy) MIKE pipeline, written by D. Kelson (see also Kelson 2003). The resulting S/N in the Gl570A spectrum was 160, and 200 at 6300 Å . Continuum normalization, order stitching, and Doppler-shift correction were performed with standard packages in IRAF12.

Though the stellar parameters of Gl570A have been measured previously (e.g., Feltzing & Gustafsson 1998; Thorén & Feltzing 2000; Ghezzi et al. 2010; Lee et al. 2011), for consistency we re-derived the T, log , microturbulence (), and [Fe/H] values from the MIKE spectra, based on the methods from our previous work (e.g., Teske et al. 2014). Briefly, Gl570A’s stellar parameters were derived from equivalent width (EW) measurements of Fe I and Fe II. We used the iron line list of Tsantaki et al. (2013), optimized for cool stars (T 5000 K) by matching spectroscopic and infrared flux method (IRFM) temperatures. We forced zero correlation between [Fe I/H] and lower excitation potential () to set T, zero correlation between [Fe I/H] and reduced equivalent width [log(EW/)] to set , and zero difference (within two decimal places) between [Fe I/H] and [Fe II/H] to set log . The abundances of Fe were determined using the local thermodynamic equilibrium (LTE) spectral analysis code MOOG (Sneden 1973), with model atmospheres interpolated from the Kurucz ATLAS9 NOVER grids13. Equivalent widths were measured in IRAF with the ‘splot’ task, and abundances were normalized to the solar values as measured in our Vesta spectrum on a line-by-line basis. The log(Fe) values for the Sun were determined with our Vesta spectrum and a solar Kurucz model with =5777 K, log =4.44 dex, [Fe/H]=0.00 dex, and =1.38 km s. In Gl570A, 89 Fe I and 10 Fe II lines were measured; the line properties and EWs are provided in Table 6.

Our final parameters and errors for Gl570A are listed below (Table 7, along with those of several other studies for comparison. The errors are calculated as in our previous work (Teske et al. 2013ab, 2014, 2015) – the change in () required to cause a correlation coefficient between [Fe I/H] and ([Fe I/H] and reduced EW) significant at the 1 level was adopted as the uncertainty in these parameters. The uncertainty in log was calculated differently, through an iterative process described in detail in Baubar & King (2010). Uncertainties in [Fe I/H] and [Fe II/H] are calculated from the quadratic sum of the individual uncertainties in these abundances due to the derived uncertainties in , log , and , as well as the uncertainty in the mean (17) of each abundance. The uncertainty due to the stellar parameters is measured from the sensitivity of the abundance to each parameter for changes of K in , dex in log , and km s in . The uncertainty due to each parameter is then the product of this sensitivity and the corresponding parameter uncertainty. For the abundances determined through spectral synthesis (e.g., from [O I], see below), models with this range of stellar parameters were compared to the data and the elemental abundance adjusted to determine the best fit.

While our log value is moderately lower than other some other studies, it agrees within errors. As described below, these stellar parameter errors are propogated through the other abundance measurements.

#### Hd3651a

The bright (=5.88) K0 dwarf HD3651A hosts a 0.2 M planet (Fischer et al. 2003), and has thus been the target of many spectroscopic observations and stellar parameter analyses (15, according to SIMBAD). Given the many previous stellar parameter analyses based on high-resolution, high-S/N data, we do not derive yet another set of stellar parameters here. Instead we use an archive HIRES spectrum (J. Johnson 2015, private communication), with S/N200 at the [O I] 6300 Å line, along with several reported sets of stellar parameters (Table 8) to verify the previously-measured carbon and oxygen abundances. We also assess realistic uncertainties on these measurements, as no formal uncertainties were previously published; our analysis is reported in the next section. The solar standard in this case is an archive HIRES spectrum of reflected light from the asteroid Vesta (A. Howard 2014, private communication), taken in the same configuration as the HD3651A spectrum.

### a.2. C/O Ratio Analysis

#### Gl570A

The determination of [C/H] in Gl570A required a different technique than our previous work, due to the low temperature of the star. The carbon abundances for Gl570A derived from EW measurements of widely-used high-excitation ( 7.67 eV) C I lines (5380, 7711, 7113Å ) and a line-by-line comparison to solar resulted in [C/H]=0.99, an unrealistically high value. These C I lines suffer NLTE effects, such that an LTE analysis overestimates the abundances, and corrections based on NLTE atomic models are predicted to increase in magnitude (larger negative corrections) with higher T and lower log (e.g., Asplund 2005; Takeda& Honda 2005; Fabbian et al. 2006). However, similar to the O I triplet lines at 7771/7774/7775Å (see discussion below), the NLTE corrections for cool stars are predicted to be minimal, less than the solar-type star corrections of -0.05. As discussed in Teske et al. (2013a) for the cool star (5350 K) 55 Cnc, applying the predicted NLTE corrections to the O I triplet abundances actually increases the [O/H] values, rather than decreases, because the solar corrections exceed the cooler-star corrections and thus their differences are increased. We find the same problem with the C I lines in Gl570A.

Thus, the [C/H] for Gl570A is instead derived from the low-excitation (=1.26 eV) forbidden [C I] line at 8727.13Å (Figure 10). This line is weak, and possibly blended with an Fe I line at 8727.10Å (Lambert & Swing 1967), but has a well-determined transition probability and is not susceptible to depatures from LTE (Gustafsson et al. 1999; Asplund et al. 2005). Using a spectral synthesis analysis with a line list between 8724-8730Å , gathered from the Vienna Atomic Line Database (VALD; Kupka et al. 1999), we derive A(C)=8.57. This also agrees with our derivation using the IRAF-measured EW (5.30 mÅ ) and the ‘abfind’ driver in MOOG with the Gl570A model derived from the same MIKE spectra. Via the same procedure, we measure A(C)=8.43 from the 8727 Å [C I] line. The resulting [C/H] is listed in Table 9. Note that in Table 9, we include the formal errors on each abundance measurement ([Ni/H], [C/H], [O/H]), which is the quadratic sum of the three individual parameter uncertainties (, log , ) and as described above for the case of iron.

#### Hd3651a

The [C I] line used to measure the carbon abundance of Gl570A was not available in our HD3651A HIRES spectrum, so we instead rely on the lowest excitation (=7.685 eV) available C I lines at 5052.2 Å and 5380.3 Å. These lines, as mentioned above, are predicted to have negligible NLTE corrections at the low temperature of HD3651A (5100-5200 K), on par with those predicted for the Sun (0.05 dex; Takeda & Honda 2005). Thus, any deviations from LTE should cancel with the calculation of the solar-normalized [C/H] value derived from these lines. In this case, the C I lines resulted in reasonable [C/H] values (see Tables 6 and 10). Measurements of other available C I lines used in previous work (6588, 7111, 7113Å, e.g. Teske et al. 2014) result in 0.1-0.25 dex larger [C/H] values – as expected since these higher excitation lines are more susceptible to NLTE effects – and thus are excluded for the final [C/H] for HD3651A. We use our equivalent width measurements (Table 6) with three different stellar models (Allende Prieto et al.’s, Delgado Mena et al.’s, and Petigura & Marcy’s, see Table 8) to derive carbon abundances, and use the resulting spread in [C/H] as a measure of its uncertainty. Our measured [C/H] values (0.21-0.37) agree with those of Allende Prieto et al. (2004) and Delgado Mena et al. (2010), as expected (Table 10). Note that in Table 10, we include only the standard deviation () errors, since the stellar parameter errors are not derived in this work, and are calculated differently for each source; the error calculation method used for Gl570A does not apply.

#### Oxygen Measurements

As discussed in Teske et al. (2013) and (2014), and by many others in the past (e.g., Nissen & Edvardsson 1992; King & Boesgaard 1995; Asplund et al. 2004; Schuler et al. 2006a,b; Caffau et al. 2008), oxygen abundances are notoriously difficult to derive, particularly in stars that are cooler/hotter and/or more/less metal-rich than the Sun. The oxygen abundance indicators we explored here include the forbidden [O I] line at 6300Å, which is well-described by LTE but blended with a Ni I line, and the O I triplet at 7771-7775Å , made up of three unblended and usually prominent lines amenable to direct EW measurement.

We derived the [O I] 6300.304Å line using two methods. First, we measured the EW of the feature directly from the spectrum as input for the ‘blends’ driver in MOOG, accounting for Ni I Ni I feature at 6300.335Å with log(Ni)=-2.695 and log(Ni)=-2.275 as derived by Bensby et al. (2004). This resulted in [O/H] -0.08 for Gl570A and 0.19-0.23 for HD3651A, depending on the set of stellar parameters. Second, we performed a spectral synthesis analysis on the line, using as a “known” our measured [Ni/H] abundance based on a line-by-line analysis with the Sun (as with Fe; see Tables 9 and 10). The synthesized spectra (with the stellar parameters derived or noted above) were convolved with a Gaussian profile, based on near-by unblended lines, to represent the instrument PSF, stellar macroturbulence, and rotational broadening; we also fixed the nickel abundance to our measured value. Unlike in the case of 55 Cnc (Teske et al. 2013a), the oxygen line strength did not change drastically (0.02 dex) by changing the Ni abundance within our derived [Ni/H]. The remaining free parameters were continuum normalization, wavelength shift, and oxygen abundance. The best fit to the synthesized spectra for the [O I] was determined by minimizing the deviations between the observed and synthetic spectra (see Figure 10). This resulted in [O/H] -0.13 for Gl570A and 0.17-0.27 for HD3651A, depending on the set of stellar parameters. For our final [O/H] we took the mean of these measurements, -0.110.19 dex for Gl570A and 0.18-0.25 dex for HD3651A.

The O I triplet suffers NLTE effects due to the dilution of the each line’s source function with respect to the Planck function (e.g., Kiselman 1993; Gratton et al. 1999; Kiselman 2001), so abundances derived assuming LTE are overestimated. As with C I, the predicted NLTE corrections increase with decreasing gas pressures and/or increasing temperatures, but decrease for cool stars (e.g., Takeda 2003; Ramírez et al. 2007, Fabbian et al. 2009) like both Gl570A and HD3651A. We measured the EWs of the three O I triplet lines and derived an A(O) for each line in Gl570A, HD3651A, and their respective solar standards (see Table 6); the average [O/H] of the line-by-line differences with the Sun is 0.0513 for Gl570A (line 3 of Table 9) and 0.17-0.38 for HD3651 (lines 3, 6, and 9 of Table 10). We then applied NLTE corrections from three different sources – Takeda (2003), Ramírez et al. (2007), and Fabbian et al. (2009) – to Gl570A, HD3651, and the respective solar standard measurements, and recalculated the line-by-line abundance differences (see discussion in Teske et al. 2013a for details of each of these correction schemes). Corrections in absolute abundance (not relative to solar) for the Sun range from 0.13-0.21, for Gl570A range from 0.04-0.07, and for HD3651 range from 0.52-0.85. The resulting [O/H], averaged over all three lines and all three sources of NLTE corrections, is 0.14 dex for Gl570A (line 4 of Table 9) and 0.25-0.46 dex, depending on the stellar parameters. In the case of Gl570A, based on previous NLTE abundance uncertainty calculations, we assume the same uncertainty as the LTE abundance; for HD3651, again the (standard deviations) across the three lines and three sources are listed in Table 10.

### a.3. What are the C/O Ratios of Gl570A and HD3651A?

We calculate the C/O ratio18 of Gl570A and HD3651A with the Asplund et al. (2009) solar A(C)=8.43 and A(O)=8.69 values and our [C/H] and [O/H] values as follows:

C/O 10

with the error on C/O represented by the errors of [C/H] and [O/H] added in quadrature. From our different [O/H] indicators, the C/O ratio for Gl570A ranges from 0.550.17 to 0.97 0.22, and for HD3651A ranges from 0.47-0.79, depending on the set of stellar parameters and oxygen abundance indicators.

The [O I] line has been designated as a consistently reliable oxygen abundance indicator (e.g., Lambert 1978; Allende Prieto et al. 2001; Asplund et al. 2004; Schuler et al. 2006a) due to its formation in the ground state making an LTE approximation exceedingly good (Caffau et al. 2008). As noted above, the O I triplet suffers significant NLTE effects, which are predicted to decrease with temperature. However, studies by Schuler et al. (2004,2006b) and King & Schuler (2005) of dwarf stars in the Pleiades, M34, and Hyades open clusters, and the Ursa Major moving group found that [O/H] values derived from the O I triplet significantly increased with decreasing (5400). If the assumption holds that stars within a single cluster or moving group should be chemically homogenous, the increasing [O/H] with is in direct contrast with all the available NLTE calculations. These studies do not point to a definitive cause of the NLTE correction discrepancy in cool dwarfs, though they suggest age likely plays a role. However, the studies do indicate for cool stars like Gl570A and HD3651A, the O I triplet-tempearture trend appears to contradict the predicted NLTE oxygen abundance corrections. Due to larger corrections for the Sun versus Gl570A and HD3651A, applying the NLTE-corrections of Takeda (2003), Ramírez et al. (2007), and Fabbian et al. (2009) all result in [O/H] values for Gl570A and HD3651 that are in general larger than in the LTE case, and C/O values that are in general smaller.

Given the caveats of the O I triplet NLTE values, we chose here to combine the [O/H] values from the three O I lines, and [O/H] for our final best estimate for the oxygen abundances of Gl570A and HD3651A, resulting in [O/H] -0.030.12 for Gl570A (where the uncertainty here is the errors of the O I and [O I] abundances added in quadrature) and [O/H]0.21-0.29 for HD3651A, depending on the set of stellar parameters. With [C/H]0.140.11, our final C/O for Gl570A 0.810.16. For HD3651, taking [O/H]0.23, 0.07 and [C/H]0.28, 0.08 from the three different stellar parameter analyses, the final C/O for HD3651A=0.62, 0.11.

### a.4. Comparison with the Literature

#### Gl570A

Our [O/H] differs from the [O/H] values for Gl570A reported by Feltzing & Gustafsson (1998) and Petigura & Marcy (2011), who both find [O/H]0.15. This higher oxygen abundance, combined with our measured [C/H] (none of the other studies measured carbon in Gl570A), would lower the C/O ratio to 0.54, which matches the solar value (C/O=0.550.10; Asplund et al. 2009; Caffau et al. 2011). Considering our errors, and those reported by Petigura & Marcy (2011) for [O/H], the high and low C/O values would just barely overlap within errors.

Both Feltzing & Gustafsson (1998) and Petigura & Marcy (2011) also measure higher [Ni/H] values (0.18 dex), although the Feltzing & Gustafsson abundance overlaps with ours within errors. Using a stellar model with Petigura & Marcy (2011)’s derived parameters for Gl570A (4744 K , 4.76 dex log , 0.10 dex [Fe/H], and assuming 1.00 km s) and our measured EWs for Ni I, C I, [O I], and the O I triplet in LTE results in [Ni/H]0.14 and C/O100.89. Performing the same exercise with Feltzing & Gustafsson (1998)’s derived stellar parameters (4585 K T, 4.70 dex log , 0.04 dex [Fe/H], 1.00 km s) results in [Ni/H]0.14 and C/O100.63, where in both formulae the average [O/H]([O/H][O/H])/2. Given our uncertainties and the typical uncertainties in [C/H] and [O/H], particularly in cool stars, these alternative C/O values and our reported C/O agree within errors. This exercise also demonstrates that the different [O/H] and [Ni/H] abundances derived in this work versus in Feltzing & Gustafsson (1998) and Petigura & Marcy (2011) are likely due to differences in stellar parameters, and not the quality of the spectra or the empirical measurements – using our measurements and their models results in [O/H] and [Ni/H] values very similar to what they report.

#### Hd3651a

The [C/H] values we derive for HD3651A based on our measurements of a high-S/N archive HIRES spectrum are in decent agreement with those reported by Allende Prieto et al. (2004) and Delgado Mena et al. (2010). Since our [C/H] is derived from stronger C I lines, versus Allende Prieto’s [C I] 8727 Å measured [C/H], we might expect our abundance to be higher, which it is. We have perfect agreement with Delgado Mena et al. (2010) when using their stellar parameters and our C I equivalent width measurements ([C/H]=0.25).

The most challenging aspect of the C/O measurement in this (and many) stars is pinning down the [O/H]. In comparison to Allende Prieto et al. (2004)’s [O/H] derived from the [O I] 6300 Å line, our [O I] 6300 Å measurement combined with their stellar parameters produces almost an almost identical [O/H] (0.21 versus 0.23). This is not the case when comparing [O/H] values of Delgado Mena et al. (2010) and Petigura & Marcy (2011) to those measured here, where we find oxygen abundance higher by 0.12 and 0.18, respectively. Other authors (Fortney 2012; Nissen 2013; Teske et al. 2013a; Nissen et al. 2014) have called into question the high C/O ratios (often caused by low oxygen abundances) reported in the previous papers, and in some cases have reported different [O/H] results for the same stars. Differences in [O/H] measured from the same line in the same star could arise due to several challenging aspects of the 6300 Å line, such as continuum placement, telluric contamination, weakness of the line at low metallicity, and blending with an Ni I line that can make up to 30% of the line strength in the Sun (Caffau et al. 2008). Delgado Mena et al. (2010) specifically removed spectra with obvious telluric contamination, and estimated the EW of the Ni line blended with [O I] to form the 6300 Å line using the “ewfind” driver of MOOG and the Ni abundances measured from 50 Ni I lines. The [Ni/H] value we derived (0.22) from 27 Ni I lines and the same stellar parameters at Delgado Mena et al. is slightly lower than their value (0.15). Presumably a smaller Ni abundance comes from a smaller Ni EW, meaning a smaller contribution to the 6300 Å line and thus a larger contribution from O, and yet Delgado Mena et al. also find a smaller [O/H]. The Ni I and O I line parameters that we use (Bensby et al. 2004 for Ni I; Story & Zeippen 2000 for [O I]) differ from those used by Delgado Mena et al. (Allende Prieto et al. 2001, for Ni I; Lambert et al. 1978 for [O I]), but this is not enough to account for the 0.12 dex difference. Without a published EW measurements or synthesis fits, it is unclear why our [O/H] for HD3651A differs from Delgado Mena et al.’s.

In their [O I] 6300 Å measurements, Petigura & Marcy (2011) also discard any stars with telluric line contamination, but many of their spectra are contaminated by iodine lines, which are 5% deep in the region of the line. To account for this, they shift the stellar spectrum (with iodine) to match in wavelength the most recent iodine reference spectrum, and then divide the stellar spectrum by the iodine reference spectrum. The resulting spectrum can contain artifacts at the 1% level, likely within any EW measurement. Petigura & Marcy treat the Ni blend differently: Once they derive an [O/H] from the 6300 Å line by comparing synthetic spectra to their observed spectra, they refit the oxygen line to spectra with 0.03 dex Ni, and add the resulting errors in [O/H] in quadrature to their statistical errors. Petigura & Marcy authors also use different line parameters for the Ni I line blended with oxygen at 6300 Å), determined by fitting their solar spectrum to match their adopted solar abundance distribution. The [Ni/H] abundance derived by Petigura & Marcy for HD3651A (0.24) agrees with what we derive using our EWs measured from 27 Ni I lines and their stellar parameters, but they do not explicitly use the measured nickel abundance for each star in their measurement of [O/H]. This may be the reason behind our differing [O/H] values.

Ramírez et al. (2013) do not measure [O/H] from the 6300 Å forbidden line, but instead from the O I 7771-5 Å triplet. As noted above, these lines suffer NLTE effects that are not well understood or calibrated for cool stars. Combining our triplet line EWs with the stellar parameters of HD3651A from Ramírez et al. (who do not list their measured EWs) results in [O/H]=0.11, =0.05, which is slightly higher than, but overlaps within errors of, Ramírez et al.’s [O/H]=0.050.04 (where here we have added Ramírez et al.’s line-to-line scatter, 0.01, and their uncertainty in the stellar parameters, 0.04, in quadrature for a total error). Similarly, we find [O/H]=0.19, =0.04, using our triplet EWs, Ramírez et al.’s stellar parameters and Ramírez et al. (2007)’s NLTE correction scheme, whereas Ramírez et al. (2013) reports [O/H]=0.120.04 for HD3651A. Thus, while slightly higher, our [O/H] values still consistent with Ramírez et al.’s within errors; if we use [O/H] based on Ramírez et al.’s stellar parameters, and recalculate [C/H] with the same parameters, the resulting HD3651 C/O=0.66, in good agreement with our final value above (0.620.11).

Facilities: Magellan:Clay

## References

### Footnotes

1. slugcomment: Submitted to the Astrophysical Journal
2. affiliationtext: Correspondence to be directed to mrline@ucsc.edu
3. http://ssc.spitzer.caltech.edu/warmmission
/propkit/pet/magtojy
4. We note that uniform-in-log priors may cause issues for a larger parameter set. For larger parameter sets a Dirichlet prior should be used
5. van Leeuwen et al. 2007
6. Ramierez et al. (2013)
7. derived from our stellar abundance analysis described in the Appendix
8. the chemically derived bulk quantities are the “by hand” thermochemical model fits to the retrieved molecular abundances
9. Assumes an atmospheric metal depletion due to loss of O and silicates of 20% based on the chemical models
10. derived from our stellar abundance analysis described in the Appendix
11. assuming an O depletion from silicates of 28% for HD3651B and 26% for Gl570D based on the chemical model results
12. IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.
13. See http://kurucz.harvard.edu/grids.html.
14. Abundance derived through equivalent width analysis.
15. Abundance derived through synthesis analysis.
16.  LTE abundance.
17. , where is the standard deviation of the derived abundances and is the number of lines used to derive the abundance.
18. The C/O ratio – the ratio of the number of carbon atoms to oxygen atoms – is calculated in stellar abundance analysis as C/O= N/N=10 where log(N)=log(N/N)+12.

### References

1. Allard, F., Hauschildt, P. H., Baraffe, I., & Chabrier, G. 1996, ApJ, 465, L123
2. Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357
3. Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872
4. Benneke, B., & Seager, S. 2012, ApJ, 753, 100
5. Burgasser A. J., 2007, ApJ, 658, 617
6. Burgasser A. J., Burrows A., Kirkpatrick J. D., 2006b, ApJ, 639, 1095
7. Burgasser A. J., Kirkpatrick J. D., Cutri R. M., McCallon H., Kopan G., Gizis J. E., Liebert J., Reid I. N., Brown M. E., Monet D. G., Dahn C. C., Beichman C. A., Skrutskie M. F., 2000, ApJ, 531, L57
8. Burgasser, A. J., Geballe, T. R., Leggett, S. K., Kirkpatrick, J. D., & Golimowski, D. A. 2006a, ApJ, 637, 1067
9. Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856
10. Burrows, A., & Sharp, C. M. 1999, ApJ, 512, 843
11. Burrows, A., Hubbard, W. B., Lunine, J. I., & Liebert, J. 2001, Reviews of Modern Physics, 73, 719
12. Cohen, M., Wheaton, W. A., & Megeath, S. T. 2003, AJ, 126, 1090
13. Canty, J. I., Lucas, P. W., Yurchenko, S. N., et al. 2015, MNRAS, 450, 454
14. Czekala, I., Andrews, S. M., Mandel, K. S., Hogg, D. W., & Green, G. M. 2014, arXiv:1412.5177
15. Cushing, M. C., Marley, M. S., Saumon, D., et al. 2008, ApJ, 678, 1372
16. Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al. 2011, ApJ, 743, 50
17. Delgado Mena, E., Israelian, G., González Hernández, J. I., et al. 2010, ApJ, 725, 2349
18. Duquennoy A., Mayor M., 1988, A&A, 200, 135
19. Eilers, P. & Marx, B. 1996, Stat. Sci., 88
20. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
21. Forveille T., Beuzit J.-L., Delfosse X., Segransan D., Beck F., Mayor M., Perrier C., Tokovinin A., Udry S., 1999, A&A, 351, 619
22. Fortney, J. J. 2012, ApJ, 747, L27
23. Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013, ApJ, 775, 80
24. Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25
25. Geballe T. R., Saumon D., Leggett S. K., Knapp G. R., Marley M. S., Lodders K., 2001, ApJ, 556, 373
26. Ghezzi L., Cunha K., Smith V. V., de Araújo F. X., Schuler S. C., de la Reza R., 2010, ApJ, 720, 1290
27. Gliese W., 1969, Veroeffentlichungen des Astronomischen Rechen-Instituts Heidelberg, 22, 1
28. Gordon, S., & McBride, B. J.  1994, NASA Reference Publ. 1311
29. Hogg, D. W., Bovy, J., & Lang, D. 2010, arXiv:1008.4686
30. Irwin, P. G. J., Teanby, N. A., de Kok, R., et al. 2008, J. Quant. Spec. Radiat. Transf., 109, 1136
31. Kass, R. & Raftery, A.,  1995, JASA, 90, 430
32. Kirkpatrick, J. D. 2005, ARA&A, 43, 195
33. Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, ApJ, 793, L27
34. Lee, J.-M., Fletcher, L. N., & Irwin, P. G. J. 2012, MNRAS, 420, 170
35. Jullion, A., & Lambert, P.  2007, CSDA, 51
36. Lang, S., & Brezger, A.  2004, J. Comp. & Graph. Stats., 13.1
37. Leggett S. K., Marley M. S., Freedman R., Saumon D., Liu M. C., Geballe T. R., Golimowski D. A., Stephens D. C., 2007, ApJ, 667, 537
38. Line, M. R., Liang, M. C., & Yung, Y. L. 2010, ApJ, 717, 496
39. Line, M. R., Vasisht, G., Chen, P., Angerhausen, D., & Yung, Y. L. 2011a, ApJ, 738, 32
40. Line, M. R., Zhang, X., Vasisht, G., et al. 2012, ApJ, 749, 93
41. Line, M. R., Wolf, A. S., Zhang, X., et al. 2013, ApJ, 775, 137
42. Line, M. R., Knutson, H., Wolf, A. S., & Yung, Y. L. 2014a, ApJ, 783, 70
43. Line, M. R., Fortney, J. J., Marley, M. S., & Sorahana, S. 2014b, ApJ, 793, 33
44. Liu, M. C., Leggett, S. K., & Chiu, K. 2007, ApJ, 660, 1507
45. Luhman K. L., Patten B. M., Marengo M., Schuster M. T., Hora J. L., Ellis R. G., Stauffer J. R., Sonnett S. M., Winston E., Gutermuth R. A., Megeath S. T., Backman D. E., Henry T. J., Werner M. W., Fazio G. G., 2007, ApJ, 654, 570
46. Madhusudhan, N., & Seager, S. 2009, ApJ, 707, 24
47. Mariotti J.-M., Perrier C., Duquennoy A., Duhoux P., 1990, A&A, 230, 77
48. Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Science, 272, 1919
49. Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, ApJ, 756, 172
50. Moses, J. I., Visscher, C., Fortney, J. J., et al. 2011, ApJ, 737, 15
51. Mugrauer M., Seifahrt A., Neuhäuser R., Mazeh T., 2006, MNRAS, 373, L31
52. Nissen, P. E. 2013, A&A, 552, A73
53. Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16
54. Pinfield D. J., Jones H. R. A., Lucas P. W., Kendall T. R., Folkes S. L., Day-Jones A. C., Chappelle R. J., Steele I. A., 2006, MNRAS, 368, 1281
55. Rahman, M.M.  2005, Bayesian Curve Fitting with Roughness Penalty Prior Distributions, Diss. UBC
56. Rice, E. L., Barman, T., Mclean, I. S., Prato, L., & Kirkpatrick, J. D. 2010, ApJS, 186, 63
57. Rodgers, C. D. 2000, Inverse Methods for Atmospheric Sounding - Theory and Practice. Series: Series on Atmospheric Oceanic and Planetary Physics, ISBN: ¡ISBN¿9789812813718¡/ISBN¿. World Scientific Publishing Co. Pte. Ltd., Edited by Clive D. Rodgers, vol. 2, 2
58. Santos N. C., Israelian G., Mayor M., 2004, A&A, 415, 1153
59. Saumon, D., Marley, M. S., Cushing, M. C., et al. 2006, ApJ, 647, 552
60. Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327
61. Saumon D., Marley M. S., Abel M., Frommhold L., Freedman R. S., 2012, ApJ, 750, 74
62. Stephens, D. C., & Leggett, S. K. 2004, PASP, 116, 9
63. Stephens, D. C., Leggett, S. K., Cushing, M. C., et al. 2009, ApJ, 702, 154
64. Swain, M. R., Line, M. R., & Deroo, P. 2014, ApJ, 784, 133
65. Teske, J. K., Cunha, K., Schuler, S. C., Griffith, C. A., & Smith, V. V. 2013, ApJ, 778, 132
66. Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, ApJ, 574, 740
67. Trotta, R. 2008, Contemporary Physics, 49, 71
68. Tsuji, T., Ohnaka, K., Aoki, W., & Nakajima, T. 1996, A&A, 308, L29
69. Valenti J. A., Fischer D. A., 2005, ApJS, 159, 141
70. Yurchenko, S. N., Tennyson, J., Bailey, J., Hollis, M. D. J., & Tinetti, G. 2014, Proceedings of the National Academy of Science, 111, 9379
71. Allende Prieto, C., Lambert, D. L., & Asplund, M. 2001, ApJ, 556, L63
72. Allende Prieto, C., Barklem, P. S., Lambert, D. L., & Cunha, K. 2004, A&A, 420, 183
73. Asplund, M., Grevesse, N., Sauval, A. J., Allende Prieto, C., & Kiselman, D. 2004, A&A, 417, 751
74. Asplund, M. 2005, ARA&A, 43, 481
75. Asplund, M., Grevesse, N., Sauval, A. J., Allende Prieto, C., & Blomme, R. 2005, A&A, 431, 693
76. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
77. Bensby, T., Feltzing, S., & Lundström, I. 2004, A&A, 415, 155
78. Bernstein, R., Shectman, S. A., Gunnels, S. M., Mochnacki, S., & Athey, A. E. 2003, Proc. SPIE, 4841, 1694
79. Caffau, E., Ludwig, H.-G., Steffen, M., et al. 2008, A&A, 488, 1031
80. Delgado Mena, E., Israelian, G., González Hernández, J. I., Bond, J. C., Santos, N. C., Udry, S., Mayor, M., 2010, ApJ, 725, 2349
81. Ghezzi, L., Cunha, K., Smith, V. V., et al. 2010, ApJ, 720, 1290
82. Fabbian, D., Asplund, M., Barklem, P. S., Carlsson, M., & Kiselman, D. 2009, A&A, 500, 1221
83. Fabbian, D., Asplund, M., Carlsson, M., & Kiselman, D. 2006, A&A, 458, 899
84. Fischer, D. A., Butler, R. P., Marcy, G. W., Vogt, S. S., & Henry, G. W. 2003, ApJ, 590, 1081
85. Feltzing, S., & Gustafsson, B. 1998, A&AS, 129, 237
86. Fortney, J. J. 2012, ApJ, 747, L27
87. Gratton, R. G., Carretta, E., Eriksson, K., & Gustafsson, B. 1999, A&A, 350, 955
88. Gustafsson, B., Karlsson, T., Olsson, E., Edvardsson, B., & Ryde, N. 1999, A&A, 342, 426 Gustafsson 1999
89. Hinkel, N. R., Timmes, F. X., Young, P. A., Pagano, M. D., & Turnbull, M. C. 2014, AJ, 148, 54
90. Kelson, D. D. 2003, PASP, 115, 688
91. King, J. R., & Boesgaard, A. M. 1995, AJ, 109, 383
92. King, J. R., & Schuler, S. C. 2005, PASP, 117, 911
93. Kiselman, D. 1993, A&A, 275, 269
94. Kiselman, D. 2001, NewAR, 45, 559
95. Kupka, F., Piskunov, N., Ryabchikova, T. A., Stempels, H. C., & Weiss, W. W. 1999, A&AS, 138, 119
96. Lambert, D. L. 1978, MNRAS, 182, 249
97. Lambert, D. L., & Swings, J. P. 1967, Sol. Phys., 2, 34
98. Lee, Y. S., Beers, T. C., Allende Prieto, C., et al. 2011, AJ, 141, 90
99. Nissen, P. E., & Edvardsson, B. 1992, A&A, 261, 255
100. Nissen, P. E. 2013, A&A, 552, A73
101. Nissen, P. E., Chen, Y. Q., Carigi, L., Schuster, W. J., & Zhao, G. 2014, A&A, 568, A25
102. Petigura, E. A., & Marcy, G. W. 2011, ApJ, 735, 41
103. Ramírez, I., Allende Prieto, C., & Lambert, D. L. 2007, A&A, 465, 271
104. Ramírez, I., Allende Prieto, C., & Lambert, D. L. 2013, ApJ, 764, 78
105. Schuler, S. C., King, J. R., Hobbs, L. M., & Pinsonneault, M. H. 2004, ApJ, 602, L117
106. Schuler, S. C., Hatzes, A. P., King, J. R., Kürster, M., & The, L.-S. 2006a, AJ, 131, 1057
107. Schuler, S. C., King, J. R., Terndrup, D. M., et al. 2006b, ApJ, 636, 432
108. Storey, P. J., & Zeippen, C. J. 2000, MNRAS, 312, 813
109. Sneden, C. 1973, ApJ, 184, 839
110. Takeda, Y. 2003, A&A, 402, 343
111. Takeda, Y., & Honda, S. 2005, PASJ, 57, 65
112. Teske, J. K., Cunha, K., Schuler, S. C., Griffith, C. A., & Smith, V. V. 2013a, ApJ, 778, 132
113. Teske, J. K., Schuler, S. C., Cunha, K., Smith, V. V., & Griffith, C. A. 2013b, ApJ, 768, L12
114. Teske, J. K., Cunha, K., Smith, V. V., Schuler, S. C., & Griffith, C. A. 2014, ApJ, 788, 39
115. Thorén, P., & Feltzing, S. 2000, A&A, 363, 692
116. Tsantaki, M., Sousa, S. G., Adibekyan, V. Z., et al. 2013, A&A, 555, AA150
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters