UNIFORM ATMOSPHERIC RETRIEVAL ANALYSIS OF ULTRACOOL DWARFS. I. CHARACTERIZING BENCHMARKS, Gl 570D AND HD 3651B

The retrieval model is the Bayesian engine that optimizes the forward model to fit the data. We use the Markov chain Monte Carlo (MCMC) 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. (2014a), 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's theorem,

$$\begin{eqnarray}&&p({\boldsymbol{x}}| {\boldsymbol{y}})=\displaystyle \frac{{\mathcal{L}}({\boldsymbol{y}}| {\boldsymbol{x}})p({\boldsymbol{x}})}{E},\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{x}}$ is the parameter vector described in Section 2.3, ${\boldsymbol{y}}$ is the data vector (in our case the spectrum), $p({\boldsymbol{x}}| {\boldsymbol{y}})$ is the posterior probability distribution, ${\mathcal{L}}({\boldsymbol{y}}| {\boldsymbol{x}})$ is the likelihood distribution that penalizes poor fits to the data, p(${\boldsymbol{x}}$) is the prior that represents any external constraints, and E 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:

$$\begin{eqnarray}&&\mathrm{ln}{\mathcal{L}}({\boldsymbol{y}}| {\boldsymbol{x}})=-\displaystyle \frac{1}{2}\displaystyle \underset{i=1}{\overset{n}{\displaystyle \sum }}\displaystyle \frac{{({y}_{i}-{F}_{i}({\boldsymbol{x}}))}^{2}}{{s}_{i}^{2}}-\displaystyle \frac{1}{2}\mathrm{ln}(2\pi {s}_{i}^{2}).\end{eqnarray} \tag{ 2 }$$

Here the index i denotes the ith data point, in our case some property at a single wavelength bin, y is the measured flux, F(${\boldsymbol{x}}$) is the modeled flux that comes out of the forward model (Section 2.3), and s is the data error given by

$$\begin{eqnarray}&&{s}_{i}^{2}={\sigma }_{i}^{2}+{10}^{b},\end{eqnarray} \tag{ 3 }$$

where σ is the measured error for the ith data point and b is a free parameter. Differing from Line et al. (2014a), we modify the standard error on the data point by the factor 10^b to account for underestimated uncertainties and/or unknown missing forward model physics (Tremaine et al. 2002; Hogg et al. 2010; Foreman-Mackey et al. 2013), 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 10^b 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, which is normally excluded from standard fitting routines owing to the unchanging data errors. Because the data errors include the free parameter b, 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, p(${\boldsymbol{x}}$), can be broken up into several pieces as

$$\begin{eqnarray}&&p({\boldsymbol{x}}| {\boldsymbol{y}})=p({\boldsymbol{T}})p({\boldsymbol{x}}\prime )p(\gamma ),\end{eqnarray} \tag{ 4 }$$

where p(${\boldsymbol{x}}\prime$) is the prior on the log of the gas mixing ratios, the instrumental parameters, gravity, and the radius-to-distance scaling, while p(${\boldsymbol{T}}$) and p(γ) 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.

Table 2.  Summary of Priors for Each of the Parameters

Parameter	Prior
$\mathrm{log}{f}_{i}$	Uniform-in-loga with $\mathrm{log}{f}_{i}$ $\geqslant -12$, ${{\rm{\Sigma }}}_{i}{f}_{i}\leqslant 1$
(R/D)${}^{2}$, $\mathrm{log}g$	Uniform, constrained by $1\leqslant {{gR}}^{2}/G\leqslant 80\;{M}_{{\rm{J}}}$
Ti	See Equation (5)
Δ λ	Uniform (−10–10 nm)
b	Uniform, $0.01\times \mathrm{min}({\sigma }_{i}^{2})\leqslant {10}^{b}\leqslant 100\times \mathrm{max}({\sigma }_{i}^{2})$
γ	Inverse Gamma ($\tilde{{\rm{\Gamma }}}(\gamma ;\alpha ,\beta )$), $\alpha =1,\beta =5\times {10}^{-5}$

Note.

^aWe note that uniform-in-log priors may cause issues for a larger parameter set. For larger parameter sets a Dirichlet prior should be used.

Download table as:  ASCIITypeset image

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 gravity, we can then constrain the mass. We know for brown dwarfs that the mass cannot exceed ∼75–80 M_J (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 ${M}_{{\rm{J}}}$. This constraint prevents the retrieved radii and gravities from entering an unphysical region of parameter space.

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 are 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; Benneke & Seager 2012; Line et al. 2012, 2013). 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 parameterization 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 oftentimes 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 and 100 independent TP points. This was the approach used in Line et al. (2014a). Because many of the atmospheric levels are degenerate, the retrieved profiles often result in unphysical oscillations, or ringing (Rodgers 2000). The standard remedy to this problem is to implement some a priori covariance matrix, or a smoothing kernel (or Tikhonov regularization), given some set smoothing length scale and a priori width (Irwin et al. 2008). While this reduces wild oscillations, this smoothing length scale and width must be chosen a priori 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 nonparametric regression literature (Lang & Brezger 2004; Rahman 2005; Jullion & Lambert 2007). The goal is to allow flexibility to fit for each of the independent TP points while preserving smoothness, without having to set a priori 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 nonparametric 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

$$\begin{eqnarray}&&\mathrm{ln}p({\boldsymbol{T}})=-\displaystyle \frac{1}{2\gamma }\underset{i=1}{\overset{N}{\displaystyle \sum }}{({T}_{i+1}-2{T}_{i}+{T}_{i-1})}^{2}-\displaystyle \frac{1}{2}\mathrm{ln}(2\pi \gamma ).\end{eqnarray} \tag{ 5 }$$

Inside the sum is the discrete second derivative of the temperature profile at each level, i, 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 warrant 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 do not justify rough TP profiles, γ will be small, resulting in a larger penalty to rough profiles. As mentioned in Section 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 as long as there are enough evenly spaced knots to sufficiently resolve potential structure (Eilers & Marx 1996). We have tested higher numbers of knots (30 versus 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 eight 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 ∼5%–10% (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.

The observations of the K4 dwarf Gl 570A were conducted on 2014 July 13 (UT) with the Magellan Inamori Kyocera Echelle (MIKE) spectrograph (Bernstein et al. 2003) on the 6.5 m Landon Clay (Magellan II) Telescope at Las Campanas Observatory. Three frames of 100 s each were taken of the target with the 05 × 5'' slit and 1 × 1 binning. On 2014 July 12 (UT), three frames of 500 s each with the 035 × 5'' slit and 1 × 1 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 that 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 signal-to-noise ratio (S/N) in the Gl 570A spectrum was ∼160 and ∼200 at 6300 Å. Continuum normalization, order stitching, and Doppler-shift correction were performed with standard packages in IRAF.⁹

Though the stellar parameters of Gl 570A have been measured previously (e.g., Feltzing & Gustafsson 1998; Thorén & Feltzing 2000; Ghezzi et al. 2010; Lee et al. 2011), for consistency we rederived the T_eff, log g, microturbulence (ξ), and [Fe/H] values from the MIKE spectra, based on the methods from our previous work (e.g., Teske et al. 2014). Briefly, Gl 570A's stellar parameters were derived from EW measurements of Fe i and Fe ii. We used the iron line list of Tsantaki et al. (2013), optimized for cool stars (T_eff < 5000 K) by matching spectroscopic and infrared flux method temperatures. We forced zero correlation between [Fe i/H] and lower excitation potential (χ) to set T_eff, zero correlation between [Fe i/H] and reduced EW [log(EW/λ)] to set ξ, and zero difference (within two decimal places) between [Fe i/H] and [Fe ii/H] to set log g. 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 grids.¹⁰ EWs 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 N(Fe) values for the Sun were determined with our Vesta spectrum and a solar Kurucz model with T_eff = 5777 K, log g = 4.44 dex, [Fe/H] = 0.00 dex, and ξ = 1.38 km s⁻¹. In Gl 570A, 89 Fe i and 10 Fe ii lines were measured; the line properties and EWs are provided in Table 6.

Table 6.  Lines Measured, Equivalent Widths, and Abundances

Ion	λ	χ	log g f	EW${}_{\odot }$	log ${N}_{\odot }$	EW${}_{\mathrm{GL}570{\rm{A}}}$	log ${N}_{\mathrm{GL}570{\rm{A}}}$	EW${}_{\mathrm{HD}3651{\rm{A}}}$	log ${N}_{\mathrm{HD}3651{\rm{A}}}$
	(Å)	(eV)		(mÅ)		(mÅ)		(mÅ)
C i	5052.17	7.685	−1.24	33.3	8.38	...	...	25.0	8.75
C i	5380.34	7.695	−1.57	20.9	8.44	...	...	15.1	8.81
[C i]	8727.13	1.26	−8.165	5.3	8.43	6.8	8.568	...	...
[O i]	6300.30	0.00	−9.717	5.4 (5.0)	8.67a (8.62)a	7.9	8.59a	6.6	8.82a
					8.68b (8.61)b	⋯	8.54b	⋯	8.82b
O i	7771.94	9.15	0.369	71.5 (69.7)	8.86c (8.83)c	12.3	8.85c	38.4	9.14c
	7774.17	9.15	0.220	61.5 (63.2)	8.86c (8.88)c	11.2	8.93c	38.3	9.28c
	7775.39	9.15	0.001	49.0 (44.8)	8.86c (8.78)c	7.30	8.87c	25.1	9.15c

Notes. The number abundances (log N) for HD 3651A listed in this table are calculated as an example with Allende Prieto et al. (2004) stellar parameters. Any value in parentheses refers to a HIRES solar measurement; all solar-normalized HD 3651A abundances in the text are relative to HIRES solar measurements.

^aAbundance derived through equivalent width analysis. ^bAbundance derived through synthesis analysis. ^c LTE abundance.

Download table as:  ASCIITypeset image

Our final parameters and errors for Gl 570A 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. 2013a, 2013b, 2014, 2015)—the change in ${T}_{\mathrm{eff}}$ (ξ) required to cause a correlation coefficient r 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 g was calculated differently, through an iterative process described in detail in Bubar & 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 T_eff, log g, and ξ, as well as the uncertainty in the mean (σ_μ)¹¹ of each abundance. The uncertainty due to the stellar parameters is measured from the sensitivity of the abundance to each parameter for changes of ±150 K in T_eff, ±0.25 dex in log g, and ±0.30 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.

Table 7.  Gl 570A Stellar Parameters

Parameter	This Work	Feltzing	Thorén	Ghezzi et al. (10)	Lee et al. (11)
		& Gustafsson (98)	& Feltzing (00)
T${}_{\mathrm{eff}}$ (K)	4686 ± 47	4585	4585	4799 ± 72	4615
log g (cgs)	4.37 ± 0.27	4.70	4.58	4.60 ± 0.16	4.36
ξ(km s−1)	1.03 ± 0.16	1.0	1.0	0.77 ± 0.08	...
[Fe/H] (dex)	−0.05 ± 0.17	0.04	0.04	0.03 ± 0.03	−0.05

Download table as:  ASCIITypeset image

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

The bright (V = 5.88) K0 dwarf HD 3651A hosts a 0.2 M_J 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/N ∼ 200 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 (Howard et al. 2014), taken in the same configuration as the HD 3651A spectrum.

Table 8.  HD 3651A Previously Measured Stellar Parameters

Parameter	Allende Prieto et al. (04)	Delgado Mena et al. (10)	Petigura & Marcy (11)	Ramírez et al. (13)
T${}_{\mathrm{eff}}$ (K)	5117 ± 94	5173	5221	5303 ± 63
log g (cgs)	4.58 ± 0.04	4.37	4.45	4.56 ± 0.03
ξ(km s−1)	0.93	0.74	...	0.64 ± 0.12
[Fe/H] (dex)	0.13	0.12	0.16	0.18 ± 0.07

Download table as:  ASCIITypeset image

The determination of [C/H] in Gl 570A required a different technique than our previous work, owing to the low temperature of the star. The carbon abundances for Gl 570A 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]_avg = 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_eff and lower log g (e.g., Asplund et al. 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 decreasing them, 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 Gl 570A.

Thus, the [C/H] for Gl 570A 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 & Swings 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 and 8730 Å, gathered from the Vienna Atomic Line Database (VALD; Kupka et al. 1999), we derive A(C)_{Gl 570A} = 8.57. This also agrees with our derivation using the IRAF-measured EW (5.30 mÅ) and the "abfind" driver in MOOG with the Gl 570A model derived from the same MIKE spectra. Via the same procedure, we measure A(C)_solar = 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 (T_eff, log g, ξ) and σ_μ as described above for the case of iron.

Zoom In Zoom Out Reset image size

The [C i] line used to measure the carbon abundance of Gl 570A was not available in our HD 3651A 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 HD 3651A (∼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 HD 3651A. We use our EW 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 we 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 Gl 570A does not apply.

As discussed in Teske et al. (2013, 2014) and by many others in the past (e.g., Nissen & Edvardsson 1992; King & Boesgaard 1995; Asplund et al. 2004; Schuler et al. 2006a, 2006b; 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 an 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 the ⁶⁰Ni i +⁵⁸Ni i feature at 6300.335 Å with log g f(⁶⁰Ni) = −2.695 and log g f(⁵⁸Ni) = −2.275 as derived by Bensby et al. (2004). This resulted in [O/H]_6300,blends = −0.08 for Gl 570A and 0.19–0.23 for HD 3651A, 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 nearby unblended lines, to represent the instrument point-spread function, 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]_6300,synth = −0.13 for Gl 570A and 0.17–0.27 for HD 3651A, depending on the set of stellar parameters. For our final [O/H]₆₃₀₀ we took the mean of these measurements, −0.11 ± 0.19 dex for Gl 570A and 0.18–0.25 dex for HD 3651A.

The O i triplet suffers NLTE effects due to the dilution of each line's source function with respect to the Planck function (e.g., Kiselman 1993, 2001; Gratton et al. 1999), 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 Gl 570A and HD 3651A. We measured the EWs of the three O i triplet lines and derived an A(O) for each line in Gl 570A, HD 3651A, and their respective solar standards (see Table 6); the average [O/H]_LTE of the line-by-line differences with the Sun is 0.05 ± 13 for Gl 570A (line 3 of Table 9) and 0.17–0.38 for HD 3651 (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 Gl 570A, HD 3651, 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 are in the range 0.13–0.21, for Gl 570A in the range 0.04–0.07, and for HD 3651 in the range 0.52–0.85. The resulting [O/H]_NLTE, averaged over all three lines and all three sources of NLTE corrections, is 0.14 dex for Gl 570A (line 4 of Table 9) and 0.25–0.46 dex for HD 3651, depending on the stellar parameters. In the case of Gl 570A, based on previous NLTE abundance uncertainty calculations, we assume the same uncertainty as the LTE abundance; for HD 3651, again the σ (standard deviations) across the three lines and three sources are listed in Table 10.

Our [O/H]_avg differs from the [O/H] values for Gl 570A 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 Gl 570A), would lower the C/O ratio to 0.54, which matches the solar value (C/O_⊙ = 0.55 ± 0.10; Asplund et al. 2009; Caffau et al. 2008). 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's (2011) derived parameters for Gl 570A (T_eff of 4744 K, log g of 4.76 dex, [Fe/H] of 0.10 dex, 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/O = 10${}^{8.43+.37}/{10}^{8.69+(0.12+0.20)/2}=0.89$. Performing the same exercise with Feltzing & Gustafsson's (1998) derived stellar parameters (T_eff of 4585 K, log g of 4.70 dex, [Fe/H] of 0.04 dex, ξ of 1.00 km s⁻¹) results in [Ni/H] = 0.14 and C/O = ${10}^{8.43+.33}/{10}^{8.69+(0.40+0.14)/2}=0.63$, where in both formulae the average [O/H] = ([O/H]${}_{\mathrm{triplet},\mathrm{LTE}}+$[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.

The [C/H] values we derive for HD 3651A 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 EW 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.'s (2004) [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 owing 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; Storey & Zeippen 2000 for [O i]) differ from those used by Delgado Mena et al. (Allende Prieto et al. 2001, for Ni i; Lambert 1978 for [O i]), but this is not enough to account for the 0.12 dex difference. Without published EW measurements or synthesis fits, it is unclear why our [O/H]₆₃₀₀ for HD 3651A 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 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 HD 3651A (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 HD 3651A from Ramírez et al. (who do not list their measured EWs) results in [O/H]${}_{\mathrm{triplet},\mathrm{LTE}}=0.11$, σ = 0.05, which is slightly higher than, but overlaps within errors of, Ramírez et al.'s [O/H]${}_{\mathrm{triplet},\mathrm{LTE}}$ = 0.05 ± 0.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]${}_{\mathrm{triplet},\mathrm{NLTE}}=0.19$, σ = 0.04, using our triplet EWs, Ramírez et al.'s stellar parameters, and Ramírez et al.'s (2007) NLTE correction scheme, whereas Ramírez et al. (2013) report [O/H]${}_{\mathrm{triplet},\mathrm{LTE}}$ = 0.12 ± 0.04 for HD 3651A. Thus, while slightly higher, our [O/H] values are still consistent with Ramírez et al.'s within errors; if we use [O/H]${}_{\mathrm{triplet},\mathrm{LTE}}$ based on Ramírez et al.'s stellar parameters and recalculate [C/H] with the same parameters, the resulting HD 3651 C/O = 0.66, in good agreement with our final value above (0.62 ± 0.11).