An Atmospheric Retrieval of the Brown Dwarf Gliese 229B

The methane-bearing object that confirmed the existence of brown dwarfs in 1995 has been the subject of intense study and debate since its discovery. Separated from its primary at a distance of 778 ± 01, Gl 229B was initially found as a proper motion companion to an M1V star, Gl 229A (Oppenheimer et al. 1995; Nakajima et al. 1995). Initial observations of Gl 229B were done by Oppenheimer et al. (1995) who obtained a low-resolution near-infrared spectrum (1.0–2.5 μm) on the Hale 200 inch telescope, later followed up by higher resolution spectroscopic observations that would cover 0.8–5.0 μm (Geballe et al. 1996; Schultz et al. 1998; Saumon et al. 2000). These spectroscopic studies confirmed the presence of H₂O and CH₄ as well as CO in excess of the abundance predicted by chemical equilibrium. Broadband photometric observations obtained by Matthews et al. (1996), Golimowski et al. (1998), Leggett et al. (1999), Golimowski et al. (2004), combined with a Hipparcos parallax for Gl 229A (Perryman et al. 1997), since updated in Gaia Collaboration et al. (2021), resulted in a robust determination of its bolometric luminosity.

Several works have attempted to reproduce Gl 229B's spectra using thermo-chemical equilibrium models and subsequently derive its fundamental parameters to varied results. One initial spectroscopic study used PHOENIX grid models (Allard et al. 1996) to place upper limits on the effective temperature at 1000 K in an attempt to constrain the mass. However, uncertainties in the age of this system cause a significant challenge in determining fundamental parameters, such as mass, across several forward modeling attempts. As spectroscopic and photometric observations improved, Saumon et al. (2000) employed evolutionary models of Burrows et al. (1997) but were still not able to constrain log(g) better than 5.0 ± 0.5, although they predicted an atmosphere depleted in heavy metals, particularly oxygen. More recent work by Nakajima et al. (2015) attempted to use measured abundances and an age estimate of Gl 229A to determine an approximately solar metallicity for Gl 229B, assuming coevality of the pair. Filippazzo et al. (2015) derives fundamental parameters based on the evolutionary models of Baraffe et al. (2003) and Saumon & Marley (2008) to create a distance-calibrated SED and directly integrated bolometric luminosity (L_Bol). A summary of fundamental parameters for Gl 229B can be found in Table 6.

While uncertainty in age estimates for the system posed a challenge in spectroscopic modeling carried out shortly after its discovery, a dynamical mass measurement for Gl 229B highlights this tension even further. Brandt et al. (2021) used astrometry from Gaia Collaboration et al. (2021) to constrain a dynamical mass at 71.4 ± 0.6 M_Jup, placing Gl 229B at the edge of the stellar mass boundary. All evolutionary model predictions prior to this work estimate an upper limit on age to be <5 Gyr based on its bolometric luminosity and derived effective temperature, with one work even suggesting Gl 229B is as young as 16 Myr (Leggett et al. 2002a). A predicted mass as high as that reported in Brandt et al. (2020, 2021) acutely challenges this age estimate, since an object this massive would take much longer to cool to its reported luminosity and approximate temperature of 900 K. In fact, Gl 229B is on a growing list of T dwarfs whose high masses conflict with what one might expect from evolutionary models (for example, Indi AB; Dieterich et al. 2018). It is unclear at this point whether this is due to unresolved issues in the models or observational biases. However, in this work, we attempt to derive best-fit fundamental parameters for Gl 229B as well as investigate the plausibility of an anomalously high mass.

In 1995, Nakajima et al. (1995) observed Gl 229A, a known M1V dwarf (Kirkpatrick et al. 1991; Cushing et al. 2006), with the Adaptive Optics Coronograph at the Palomar 60 inch telescope in search of a proper motion companion. Gl 229A was initially chosen as a target of interest in a search for brown dwarf companions to stars within 15 pc of the Sun. Its known low space motion (Leggett 1992) as well as a precisely measured distance in Hipparcos (Perryman et al. 1997) made it a promising candidate. Since the discovery of Gl 229B, Gl 229A has also been a target of several observational studies, in an attempt to further our understanding of the physics of low mass stars and the fundamentals of this system in particular. However, there are still seemingly discordant conclusions on its age and chemical composition which have added to the mystery of this system as a whole.

Based on its kinematics and coronal activity, it was originally believed to be younger than the Sun, but at least as old as the Hyades (∼0.8 Gyr), about 0.5–5 Gyr (Nakajima et al. 1995). However, M dwarf ages are notoriously difficult to determine due to the long activity lifetimes of fully convective stars (West et al. 2008). Gl 229A is a prime example of this dilemma, as later studies are in complete disagreement over its age. Leggett et al. (2002a) used evolutionary models from Allard et al. (2001) to argue that its kinematics and flare-star designation support a much younger age of 16–45 Myr. While, more recently, Brandt et al. (2020) used chromospheric and coronal activity, along with gyrochronology relations from Angus et al. (2019), to estimate an age of 2.6 ± 0.5 Gyr. This is further complicated as Brandt et al. (2020) reports a system age of 7–10 Gyr based on the dynamical mass measurement of Gl 229B. They also calculate an age based on Gl 229A's membership in the kinematic thin disk, but note that the resulting young age is strongly disfavored by the low levels of chromospheric and coronal activity.

In better agreement among different models are metallicity measurements for Gl 229A. Using medium and high-resolution spectra, Nakajima et al. (2015), Gaidos & Mann (2014), Neves et al. (2014), Mould (1978) and Schiavon et al. (1997) find Gl 229A to be consistent with solar values. There is one exception among the literature which found a best fitting model with [M/H] = −0.5 (Leggett et al. 2002a). However, spectroscopic measurements do not seem to support such a low metallicity. Reported metallicity measurements for Gl 229A are listed in Table 7. Additionally, Nakajima et al. (2015) reported a C/O ratio for Gl 229A based upon inferred bulk carbon and oxygen abundances published in Tsuji & Nakajima (2014) and Tsuji et al. (2015) and found a slightly supersolar value of 0.68 ± 0.12 (as compared to a solar value of 0.55 from Asplund et al. 2009), making it the only determined C/O ratio for Gl 229A.

For the purposes of this retrieval work, we began using near-infrared data obtained using CGS4 at the United Kingdom 3.8 m Infrared Telescope (UKIRT) (Geballe et al. 1996) whose spectra was updated with improved photometry by Leggett et al. (1999). Table 1 lists the photometry used in this work.

Table 1. Properties of Gliese 229B

Parameter	Value	Reference
Spectral Type	T7p	2
Astrometry
R.A.	06h10m3461	1
decl.	− 21h51m5266	1
π (mas)	173.57 ± 0.017	1
μα (mas yr−1)	−135.692 ± 0.011	1
μδ (mas yr−1)	−719.178 ± 0.017	1
Photometry
YMKO (mag)	15.17 ± 0.1	3
JMKO (mag)	14.01 ± 0.05	4
HMKO (mag)	14.36 ± 0.05	4
KMKO (mag)	14.36 ± 0.05	4
$L{{\prime} }_{\mathrm{MKO}}$ (mag)	12.24 ± 0.05	5
$M{{\prime} }_{\mathrm{MKO}}$ (mag)	11.74 ± 0.11	5
SED-Derived Parameters
Radius (RJup )	0.94 ± 0.15	6
Mass (MJup )	41 ± 24	6
log(g)	4.96 ± 0.46	6
Teff (K)	927 ± 79	6
log(LBol /LSun)	−5.21 ± 0.05	6

References. (1) Gaia Collaboration et al. (2021), (2) Burgasser et al. (2006), (3) Hewett et al. (2006), (4) Leggett et al. (2002b), (5) Golimowski et al. (2004), (6) Filippazzo et al. (2015).

Download table as:  ASCIITypeset image

To further constrain the retrieved chemical abundances and other fundamental parameters, L band spectral data from Oppenheimer et al. (1998) and M band spectral data from Noll et al. (1997) were included, making a combined spectral wavelength coverage 1.0–5.0 μm. Table 2 lists the spectra used in this work.

The forward model consists of the radiative transfer solver, thermal profile, and opacity and scattering properties as a function of wavelength. The forward model solves for emergent flux from radiative transfer using the two stream technique of Toon et al. (1989). This includes scattering, first introduced by McKay et al. (1989) and later used by Marley et al. (1996); Saumon & Marley (2008) and Morley et al. (2012). We use a 64 pressure layer atmosphere (65 levels) with geometric mean pressures in range −4 < log P < 2.3 in bars, spaced at 0.1 dex intervals.

4.1.1. Thermal Profile

4.1.2. Gas Opacities

4.1.3. Gas Abundances

4.1.4. Cloud Modeling

The cloud parameterizations we utilize in this work closely follow the methodology described in Burningham et al. (2017) and Gonzales et al. (2020, 2021). As in Burningham et al. (2017, 2021), we define two categories of clouds, "slab" and "deck." Both slab and deck clouds have an opacity distributed among layers in pressure space and an optical depth determined by cloud designation as either gray or non-gray. Total optical depth (τ_cloud) for a gray cloud model is calculated at 1 μm. For a non-gray cloud application, we use a power-law distribution to describe the optical depth, τ = τ₀ λ^α , where τ₀ is the optical depth at 1 μm. For the case of a non-gray cloud, we designate another model parameter, the power (α) in the optical depth. We discuss a cloud's optical depth in terms of extinction and assume an absorbing cloud by setting the single scattering albedo to zero, as done in Gonzales et al. (2020, 2021).

Beyond the gray and non-gray cloud parameterizations, we follow the work of Burningham et al. (2021) by testing different condensate species under the assumption of Mie scattering. In particular, we investigate the impact of zinc sulfide (ZnS) and potassium chloride (KCl) condensates using refractive indices from Wakeford & Sing (2015) and pre-tabulated Mie coefficients as a function of particle radius and wavelength. Wavelength-dependent optical depths and phase angles in each pressure layer are calculated by integrating the cross-sections and Mie efficiencies over the particle size distribution in that layer for a given cloud species. Particle number density in a layer is calibrated to the optical depth at 1 μm as determined by the cloud type (slab or deck). We assume either a Hansen (Hansen 1971) or lognormal distribution of particle sizes. The Hansen distribution for particle number n with radius r is given as:

$$\begin{eqnarray}&&n(r)\propto {r}^{\tfrac{1-3b}{b}}{e}^{-\tfrac{r}{{ab}}},\end{eqnarray} \tag{ 1 }$$

where a and b are the effective radius and spread of the distribution, respectively. These parameters are defined as:

$$\begin{eqnarray}&&a=\frac{{\int }_{0}^{\infty }r\pi {r}^{2}n(r){dr}}{{\int }_{0}^{\infty }\pi {r}^{2}n(r){dr}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&b=\displaystyle \frac{{\int }_{0}^{\infty }{(r-a)}^{2}\pi {r}^{2}n(r){dr}}{{a}^{2}{\int }_{0}^{\infty }\pi {r}^{2}n(r){dr}}.\end{eqnarray} \tag{ 3 }$$

A deck cloud is parameterized by the cloud top pressure, P_top, the decay height, ${\rm{\Delta }}\mathrm{log}P$ and the cloud particle single scattering albedo, which we set to zero. The cloud top pressure is defined as the point in pressure space at which the optical depth passes unity, or τ = 1 (looking down). The decay scale pressure describes how the optical depth changes with changing pressure from the cloud deck and is defined as $d\tau /{dP}\propto \exp ((P-{P}_{\mathrm{deck}})/{\rm{\Phi }}$), where P_deck is the height at which the cloud is optically thick and Φ = P_top(10${}^{{\rm{\Delta }}\mathrm{log}P}$-1)/10${}^{{\rm{\Delta }}\mathrm{log}P}$ is the decay scale of the cloud in bars. The deck cloud becomes optically thick at P_top, so for P > P_top, optical depth increases following the decay function until it reaches Δτ_layer = 100. The deck cloud becomes opaque with increasing pressure relatively quickly as a result of the decay function, so we obtain little atmospheric information from deep below the cloud top. To account for this, the pressure-temperature (P-T) profile below the cloud deck is an extension of the gradient and spread at the cloud top pressure.

The distinguishing marker for a slab cloud is the addition of another parameter to determine the total optical depth at 1 μm (τ_cloud), since we can "see" the bottom of this type of cloud. The optical depth is distributed throughout the cloud as d τ/dP ∝ P (looking down), where its maximum value is reached at the bottom of the slab (highest pressure). In principle, the slab can have any optical depth, but we restrict the prior to 0.0 ≤ τ_cloud ≤ 100.0. Instead of considering the decay scale pressure, as we did for a deck cloud, we consider the cloud thickness in dlogP and parameterize for cloud top pressure, P_top.

The retrieval process depends on the chosen elements of the parameter set. Changing the elements that are passed to the forward model have effects on the resultant spectrum. Optimizing the forward model's fit to the data by varying the parameter set, or state-vector, takes place within a Bayesian framework. A detailed explanation of this framework can be found in Burningham et al. (2017). To summarize, Brewster applies Bayes' theorem to calculate the "posterior probability," p( x ∣ y ), the probability of a set of parameters' ( x ) truth value given some data ( y ), in the following way:

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

where ${\mathscr{L}}({\boldsymbol{x}}| {\boldsymbol{y}})$ is the likelihood that quantifies how well the data match the model, p( x ) is the prior probability on the parameter set and p( y ) is the probability of the data marginalized over all parameter values, also known as the Bayesian evidence. As detailed in Burningham et al. (2017), the original version of Brewster uses EMCEE (Foreman-Mackey et al. 2013) to sample posterior probabilities and inflate errors using a tolerance parameter to allow for unaccounted sources of uncertainty. The EMCEE chain requires tens of thousands of iterations over hundreds of parallel walkers in order to converge and can often require several days of computing to complete a single round of parameter estimation, depending on wavelength range and model complexity. Additionally, the use of EMCEE as a sampler can result in models with degenerate parameter solutions or local, rather than global, maximum likelihoods.

In this new instance of Brewster, the posterior probability space is explored using the PyMultiNest sampler (Buchner 2014; Feroz et al. 2011) which utilizes nested sampling to discover the set of parameters with the maximum likelihood given the data. PyMultiNest is a Bayesian inference tool that explores the parameter space to maximize the likelihood of the forward model fit to the data. The samples in an n-dimensional hypercube, or state-vector, are translated into parameter values via the "prior-map." The prior-map function is how prior probabilities are set for each parameter and are then transformed into appropriate parameter values to be used in the forward model. This algorithm is equipped to handle a parameter space that may contain multiple posterior modes and/or degeneracies in moderately high dimensions.

Unlike the procedure in Burningham et al. (2017), we retrieve radius and mass directly which are then used to numerically determine a value for gravity. The radius is determined from the scaling factor required to match the absolute flux from the forward model to the data and the measured parallax. The radius is restricted to be within the range 0.5–2.0 R_Jup and the mass between 1 and 80 M_Jup as suggested by Saumon & Marley (2008), COND (Baraffe et al. 2003) and DUSTY (Chabrier et al. 2000; Baraffe et al. 2002) substellar evolutionary models. In order to investigate a dynamical mass measurement of 71.4 ± 0.6 M_Jup (Brandt et al. 2021), we restrict the probability space for mass to be 70–72 M_Jup in just one of our model investigations. Table 3 lists the priors used in our modeling, which are based on those defined in Burningham et al. (2017).

Table 3. Priors for Gl 229B Retrieval Models

Parameter	Prior
gas volume mixing ratio	uniform, log fgas ≥ −14, Σgas fgas ≤ 1
thermal profile (Tbottom, Ttop, Tmiddle, Tq1, Tq3)	uniform, 0 K < T < 4000 K
radius	uniform, 0.5RJup ≤ R ≤ 2.0RJup
mass a	uniform, 1MJup ≤ M ≤ 80MJup
cloud top b	uniform, −4 ≤ logPCT ≤ +2.3
cloud decay scale c	uniform, 0 < logΔPdecay < 7
cloud thickness d	uniform, logPCT ≤ log(PCT+ΔP) ≤ 2.3
cloud total optical depth (extinction)	uniform, 0 ≤ τcloud ≤ 100
single scattering albedo	constant, ω0 = 0
wavelength shift	uniform, −0.01 < Δλ < 0.01 μm
tolerance factor	uniform, log(0.01 x $\min ({\sigma }_{i}^{2})$) ≤ b ≤ log(100 x $\max ({\sigma }_{i}^{2}$))

Notes.

^a This mass prior range was constrained to 70–72 M_Jup for a single, cloudless model. See Section 5.2. ^b For a deck cloud, this is the pressure where τ_cloud = 1, for a slab cloud this is the top of the slab. ^c Decay height for cloud deck above the τ_cloud = 1.0 level. ^d Thickness and τ_cloud retrieved only for slab cloud.

Download table as:  ASCIITypeset image

We began this modeling by investigating the impact of running a retrieval model on near-infrared data only (1.0–2.5 μm) as opposed to the entire combined infrared spectrum. We are unable to constrain CO when using only NIR data, which is known to be in excess of thermo-chemical equilibrium abundance in the photosphere (Oppenheimer et al. 1998). This result is consistent with the work of Line et al. (2015, 2017) which could not constrain CO, CO₂ or H₂S with NIR data, alone. The first observable CO feature in a T dwarf spectrum is expected to be the first vibration-rotation band (1 – 0) at 4.7 μm which was observationally confirmed in Gl 229B by Noll et al. (1997) and then again by Oppenheimer et al. (1998). Therefore, with the addition of extended infrared data (2.98–5.0 μm), we are able to constrain CO abundance across models.

For all models we use the distance-calibrated (to 10 pc) SED beginning at 1.0 μm and out to 5.0 μm. This spectrum calibration differs from the method used in Burningham et al. (2017) in which they calibrated spectra to the 2MASS J-band photometry and used the object's true distance in their initialization.

We retrieve the following gases known to sculpt T dwarf spectra: H₂O, CO, CH₄, NH₃, Na, and K. We tie Na and K together as a single element in the state vector, assuming a Solar ratio taken from Asplund et al. (2009) (Line et al. 2015; Burningham et al. 2017; Gonzales et al. 2020, 2021). We are consistent with Gonzales et al. (2020) in excluding CO₂ and H₂S from our gas list, as we are still unable to constrain these even with the inclusion of longer wavelength data. We test multiple cloud parameterizations, beginning with the most simple, cloudless model and building to a 5 parameter Mie scattering cloud slab model.

Across models, we focused on making changes in our approach to cloud parameterization while holding fixed the gases included in each model, as well as gas abundance method and alkali opacity. In order to compare these different models, we use a calculation of the Bayesian evidence, specifically the logEvidence or logEv, where the highest logEv is preferred. We use the following selection criterion from Kass & Raftery (1995) to distinguish between two models, with evidence against the lower logEv as:

• 1.  
  0 < ΔlogEv < 0.5: no preference worth mentioning
• 2.  
  0.5 < ΔlogEv < 1: positive
• 3.  
  1 < ΔlogEv < 2: strong
• 4.  
  ΔlogEv > 2: very strong

We began by building from the least complex model (cloud-free) to the most complex slab cloud model. We compare Burrows and Allard alkali opacities for the best fit model as there is conflicting evidence over which is preferred for T dwarf retrievals (Line et al. 2017; Gonzales et al. 2020).

5.1.1. Retrieved Gas Abundances and Fundamental Parameters

Figure 2 shows the posterior probability distributions for retrieved gas abundances, mass and radius as well as log(g), T_eff, L_Bol, C/O ratio, [M/H], [C/H] and [O/H] from our best fit model, which are calculated based on retrieved quantities. We list the values from Figure 2 in Table 5 for ease of reading.

Zoom In Zoom Out Reset image size

Table 5. Cloudless, Allard Alkali Model Parameters

Parameter	Value
Model	Best Fit a	Mass Constrained b
Retrieved
H2O	−3.53 ± 0.04	−3.49 ± 0.03
CO	$-{4.59}_{-0.24}^{+0.21}$	$-{4.56}_{-0.25}^{+0.22}$
CH4	−3.38 ± 0.05	−3.32 ± 0.03
NH3	−4.58 ± 0.05	−4.51 ± 0.03
Na+K	−5.69 ± 0.03	−5.67 ± 0.03
Radius (RJup)	1.10 ± 0.04	1.12 ± 0.04
Mass (MJup)	${50}_{-9}^{+12}$	71 ± 1
Derived
log g (dex)	${5.01}_{-0.09}^{+0.10}$	5.15 ± 0.03
Teff (K)	${834}_{-15}^{+17}$	${829}_{-17}^{+19}$
log(LBol /LSun)	$-{5.25}_{-0.01}^{+0.02}$	$-{5.25}_{-0.01}^{+0.02}$
C/O	${1.38}_{-0.07}^{+0.08}$	1.45 ± 0.07
C/O c	${1.06}_{-0.05}^{+0.06}$	${1.11}_{-0.05}^{+0.06}$
[M/H]	−0.24 ± 0.05	−0.19 ± 0.03
[C/H]	−0.01 ± 0.05	0.05 ± 0.03
[O/H]	$-{0.41}_{-0.04}^{+0.05}$	−0.37 ± 0.04

Notes.

^a Best fit model without any prior constraints placed on retrieved parameters. ^b Best fit model with a prior constraint placed on the mass. ^c Oxygen-corrected C/O ratio (see Section 8).

Download table as:  ASCIITypeset image

The derived T_eff and log(g) are calculated from the retrieved radius and mass along with the parallax measurement. The scale factor (R²/D²) is calculated from the retrieved radius and parallax. T_eff is then determined using this derived scale factor and by integrating the flux in the resultant forward model spectra between 0.6 and 20 μm. We find that our retrieval-based T_eff is cooler than the semi-empirical value found from the SED, although the large uncertainty on the SED generated value allows for 1.2σ agreement. However, this ∼100 K temperature difference is due to our retrieved radius being ∼0.16 R_Jup larger than its SED value. Our retrieved radius and mass both agree within 1σ to the values determined from evolutionary models when generating the SED. Our derived value for log(g) is also within 1σ agreement with its SED value. Due to the uncertainty on mass being relatively large for both our retrieval model and the SED generated value, as well as inconsistency with the reported dynamical mass, we discuss this fundamental parameter in greater detail in Section 6.1.

The C/O ratio is calculated under the assumption that all of the oxygen exists in H₂O and CO and all of the carbon exists in CO and CH₄. The following equations were used to derive a value for metallicity, [M/H]:

$$\begin{eqnarray}&&{{f}}_{{{\rm{H}}}_{2}}=0.84\times (1-{{f}}_{\mathrm{gas}})\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{N}_{{\rm{H}}}=2{{f}}_{{{\rm{H}}}_{2}}{N}_{\mathrm{tot}}\end{eqnarray} \tag{ 5b }$$

$$\begin{eqnarray}&&{N}_{\mathrm{element}}=\sum _{\mathrm{molecules}}{n}_{\mathrm{atoms}}{{f}}_{\mathrm{molecule}}{N}_{\mathrm{tot}}\end{eqnarray} \tag{ 5c }$$

$$\begin{eqnarray}&&{N}_{M}=\sum _{\mathrm{element}}\displaystyle \frac{{N}_{\mathrm{element}}}{{N}_{{\rm{H}}}}\end{eqnarray} \tag{ 5d }$$

$$\begin{eqnarray}&&[{\rm{M}}/{\rm{H}}]=\mathrm{log}\left(\displaystyle \frac{{N}_{M}}{{N}_{\mathrm{sun}}}\right)\end{eqnarray} \tag{ 5e }$$

where ${f}_{{{\rm{H}}}_{2}}$ is the fraction of H₂, f_gas is the total gas fraction, N_H is the number of neutral hydrogen atoms, N_element is the number of atoms for each element, n_atoms is the number of atoms for a given element contained in a single molecule and N_tot is the total number of gas molecules. In our calculation of metallicity, N_Sun is determined using the same formula as N_M , using the sum of the solar abundances compared to H. Both [C/H] and [O/H] were calculated using the same procedure for total metallicity using just single element abundances.

As we do not have an SED-based comparison for C/O ratio or metallicity, we discuss these parameters in greater detail in Section 8.

5.1.2. Temperature–Pressure Profile, Contribution Function and Abundances

Figure 3(a) shows the retrieved temperature versus pressure (T-P) profile for the best fit cloudless model for Gl 229B. Overplotted in this figure are the Sonora grid models (Marley et al. 2021) using Solar-metallicity. The slight thermal inversion apparent in the very top of the atmosphere (∼0.01–0.001 bar) in our retrieved profile is a computational consequence of our initialized five point profile and not thought to actually exist in that way. As a result, we have much larger 1σ and 3σ confidence ranges in that pressure space. Our retrieved profile is in best agreement with the Sonora Solar-metallicity log(g) = 5.0, 850 K model throughout the photosphere (∼0.5 – 20 bar). While the 1σ and 3σ bounds are relatively small in this pressure range, the log(g) = 5.0, 850 K model fits the profile almost entirely within these bounds. This model agrees to within 3σ in the upper atmosphere, as well. Comparing our retrieved profile to the Solar-metallicity log(g) = 5.0750 K model, our profile is ∼50 K warmer throughout most of the photosphere, although we do see 1σ agreement at the top of the atmosphere (≤0.3 bar) as well at as the bottom of the photosphere. For the log(g) = 5.0, 950 K model, it is clear that our profile is at least 100 K cooler throughout the photosphere and deeper into the atmosphere.

Zoom In Zoom Out Reset image size

Figure 3(b) shows the contribution function indicating the location of an optical depth of τ = 1 for the gas opacity. The contribution function is calculated in each pressure layer in the following way:

$$\begin{eqnarray}&&C(\lambda ,P)=\displaystyle \frac{B(\lambda ,T(P)){\int }_{{P}_{1}}^{{P}_{2}}d\tau }{{e}^{{\int }_{{P}_{0}}^{{P}_{2}}d\tau }},\end{eqnarray} \tag{ 6 }$$

where B(λ, T(P)) is the Planck function, P₀ is the pressure at the top of the atmosphere, P₂ is the pressure at the bottom of the layer and P₁ is the pressure at the top of the layer. The majority of the flux in the near-infrared spans a relatively large portion of the atmosphere from ∼0.5 to 80 bar compared to the flux contribution in the early mid infrared that originates higher in the photosphere around ∼0.2–15 bar.

Figure 3(c) shows retrieved gas abundances as compared to thermo-chemical equilibrium grid model abundances. These chemical grids were calculated using the NASA Gibbs minimization CEA code (McBride & Gordon 1994) based on prior thermo-chemical models (Fegley & Lodders 1994, 1996; Lodders 1999, 2002, 2004, 2010; Lodders & Fegley 2002, 2006; Visscher et al. 2006, 2010; Visscher 2012; Moses et al. 2012, 2013). These grids are then used to determine thermo-chemical equilibrium abundances of various atmospheric species for pressure in the range 1 microbar – 300 bar and temperatures in the range 300–4000 K. In this case, we compare against chemical grids of Solar metallicity and supersolar C/O ratio of 1.0 relative to the Asplund et al. (2009) Solar value of 0.55.

Since our model calculates uniform-with-altitude abundances, as opposed to the more physically plausible varying-with-altitude mixing ratios (see Section 4), our retrieved values for CO and Na+K show the median retrieved abundance despite their equilibrium abundances trailing off toward lower values at the top of the photosphere. This retrieved abundance can be considered an average of abundances probed at different pressure levels in the photosphere. Therefore, we expect to find a median retrieved abundance close to the expected value in the middle of the photosphere. More generally, we would at least expect the retrieved abundance to be within the range of top-of-atmosphere and bottom-of-atmosphere abundances set by the thermo-chemical equilibrium grids. Here, we can see photospheric agreement around 10 bar for CO and 3 bar for Na+K. A similar result is shown for NH₃ where the median retrieved abundance is in range of its equilibrium abundance in the photosphere.

There is a distinction for our retrieved H₂O and CH₄ which show a slightly depleted abundance as compared to their equilibrium predictions. H₂O, CH₄ and CO are plotted together in Figure 3(c) as their chemical mixing influences the photospheric abundances we observe. We know from Oppenheimer et al. (1998) and Noll et al. (1997) that there are unexpectedly high amounts of CO in the atmosphere of Gl 229B based on spectral features, particularly at 4.7 μm, which suggest disequilibrium chemistry. A slightly depleted abundance of H₂O and CH₄ could be a result of strong vertical mixing that also causes CO abundance in excess of thermo-chemical equilibrium predictions in the photosphere. This is due to mixing timescales being shorter than the CO → CH₄ chemical timescale (Fegley & Lodders 1996; Marley & Robinson 2015), where the net reaction for this conversion in Gl 229B is

$$\begin{eqnarray}&&\mathrm{CO}+3{{\rm{H}}}_{2}\to {\mathrm{CH}}_{4}+{{\rm{H}}}_{2}{\rm{O}}\end{eqnarray} \tag{ 7 }$$

(Visscher & Moses 2011). We discuss further causes for depleted oxygen abundance in Section 8. In general, since we expect chemical disequilibrium in the atmosphere of Gl 229B and we are freely retrieving gas abundances (as opposed to retrieving thermo-chemical equilibrium abundances), it is not surprising that our abundances for H₂O and CH₄ differ slightly from these chemical grids.

5.1.3. Retrieved Spectrum versus Forward Model Fit

Figures 4(a) and (b) show our retrieved spectrum as compared to the observed data. The retrieved spectrum fits the near-infrared portion of the observed spectrum well, even though it struggles to reach the top of the flux peaks in the J and H spectral bands. The retrieved spectrum is able to fit the methane features at 1.63, 1.67 and 1.71 μm particularly well. It also fits the narrow absorption features due to water on the shortward side of the H-band peak. There is also a good fit to the top of the flux peak in the K-band. The choice of alkalies had a significant impact on how our model fit the K-I doublet in the J spectral band. In particular, we find the use of Allard alkalies fits this K-I doublet much better than a cloudless model with Burrows alkalies. This difference in the retrieved spectrum J-band between the Allard and Burrows cloud models is likely a result of how pressure broadening for these lines is calculated for objects in this temperature regime. The retrieved spectrum does a good job of fitting the fundamental methane absorption feature of the L spectral band, particularly in the region from 3.6 to 4.1 μm. The retrieved spectrum also fits the general shape of the M spectral band data and attempts to fit the CO feature at 4.67 μm, providing a generally good fit within uncertainty bounds. This could be due to difficulty fitting disequilibrium abundances of CO, as CO is also the most poorly constrained of all the gases in our model. It is important to note that the narrow feature around ∼4.8 μm is due to an incorrect removal of a telluric line (Noll et al. 1997) and therefore we do not expect our model to fit this.

Zoom In Zoom Out Reset image size

Figures 4(c) and (d) show our retrieved spectrum and observed data in comparison to Sonora model synthetic spectra that bracket our retrieved T_eff . We find that the observed spectrum is best fit by the 850 K Solar metallicity model, overall, but note that it does not reach the flux peak in the L-band. We also point out that none of the Sonora models seem to adequately fit the observed M-band data. While model spectra can provide good spectral fits in certain bands, no single model can simultaneously fit J, H, K, L and M spectral bands as well as our best fit retrieval model spectrum.

In this section, we present results from constraining the posterior probability range for the mass from 10–80 M_Jup to 70–72 M_Jup on our best fit model. The motivation for placing a prior on this model is the dynamical mass reported in Brandt et al. (2021) that predicts a companion mass for Gl 229A of 71.4 ± 0.6 M_Jup . Placing a prior on the model inherently increases model confidence, so we do not use our Bayesian evidence parameter to compare with our best fit model (i.e., retrieval model without prior knowledge placed on retrieved parameters). However, we intend to use this model constraint to compare retrieved gas abundances and fundamental parameters to our best fit model, in order to understand whether the photospheric chemistry can be explained by a 70 M_Jup object as well as a 50 M_Jup object in Section 6.1.

5.2.1. Retrieved Gas Abundances and Fundamental Parameters

Figure 5 shows the posterior probability distributions for retrieved gas abundances and radius in addition to the derived log(g), T_eff, L_Bol, C/O ratio, [M/H], [C/H] and [O/H] for the mass constrained model. We list the values from Figure 5 in Table 5 for ease of reading.

Zoom In Zoom Out Reset image size

The procedure to calculate T_eff, log(g), C/O ratio and metallicity are described in Section 5.1. Again, we have 1σ agreement with the retrieved radius and derived log(g) to the values determined from the SED. The mass constrained model found a T_eff slightly cooler than the best fit model, however, we still have agreement within 1.2σ to the SED value.

As compared to parameters of the best fit model, the T_eff, radius, and log(L_Bol/L_Sun) of the constrained model are all within 1σ agreement. We have 1.5σ agreement with log(g) values which is expected considering the higher mass of the constrained model. Additionally, we have retrieved abundances of H₂O, CO, Na and K in 1σ agreement and CH₄ and NH₃ within 1.5σ agreement between models. For an object of this temperature, we expect to see increased abundance in both CH₄ and NH₃ as a result of an increased log(g) as compared to CO, for example, which is thought to be less sensitive to gravity (Zahnle & Marley 2014). Overall, we find the chemistry and temperature to be consistent between our best fit, cloudless model and a cloudless model with a constraint at the dynamical mass value.

5.2.2. Temperature–Pressure Profile

Figure 6(a) shows the retrieved T-P profile for the mass constrained cloudless model as compared to our best fit cloudless profile as well as Sonora grid models. As our derived T_eff is only slightly cooler (∼5 K) for this mass constrained model than our best fit model, we find the Sonora Solar-metallicity log(g) = 5.0, 850 K model fits well throughout the top of the photosphere (1–5 bar). At pressures >10 bar and above the photosphere we have better agreement with the 750 K model. The 1σ and 3σ bounds, which were previously noted to be small for the best fit model, are shown to almost disappear below 1 bar. This is due to increased model confidence as a result of the prior constraint placed on mass. However, we still see good agreement to the log g = 5.0, 850 K model throughout the top of the photosphere until about P = 5 bar where the Sonora model diverges from our profile to slightly warmer temperatures. Above the photosphere, we have agreement within 3σ to both 750 K and 850 K models. In contrast, the mass constrained model profile is approximately 100 K cooler throughout the photosphere than the Solar-metallicity log g = 5.0, 950 K model.

Zoom In Zoom Out Reset image size

As shown in Figure 6(a), the best fit profile is nearly identical to the mass constrained profile, only diverging to slightly warmer temperatures (∼10 K) toward the bottom of the photosphere (around 5 bar).

5.2.3. Retrieved Spectrum versus Best Fit Model

Of all fundamental parameters to consider, constraining the mass of Gl 229B has been a top priority as an updated dynamical mass of 71.4 ± 0.6 M_Jup reported in Brandt et al. (2021) is at odds with several evolutionary model predictions that placed Gl 229B in the 30–55 M_Jup range (Nakajima et al. 1995, 2015; Allard et al. 1996; Marley et al. 1996). While it is difficult to determine a precise age for Gl 229A, gyrochronology, thin disk kinematics and coronal and chromospheric activity all disfavor an old age, placing the Gl 229 system at an intermediate age of 2–6 Gyr (Brandt et al. 2020). Evolutionary models (ex: Burrows et al. 1997; Allard et al. 2001; Saumon & Marley 2008; Phillips et al. 2020) predict that a T dwarf with a mass of 71.4 ± 0.6 M_Jup would have to be 4–7 Gyr older than the predicted age of Gl 229A to have cooled to an observable log(L_Bol/L_Sun) ≈ −5.2. As a result, best fit evolutionary models place Gl 229B at or below 55 M_Jup in order to be consistent with this age-mass–luminosity relation. While our best fit model independently deduced a mass of ${50}_{-9}^{+12}$ M_Jup, we also present results using our best fit model (cloudless, Allard alkalies) with a constraint on the posterior probability space for mass from 1–80 to 70–72 M_Jup.

From our results in Section 5.2 with the constrained mass prior, we see an increase in all molecular abundances as compared to our best fit model. However, the abundances of the mass constrained model are all within or close to the 1σ bounds of the best fit model. We find no difference in any retrieved abundance worth mentioning. There is an increase in our calculated log(g) from ${5.01}_{-0.09}^{+0.10}$ to 5.15 ± 0.03 which is expected, due to the increased mass. However, we do see a relative increase in retrieved radius from 1.10 ± 0.04 to 1.12 ± 0.04 R_Jup which is unexpected since brown dwarfs are predicted to contract with age. However, these retrieved radii agree within 1σ so while this result is unexpected it is not unreasonable. Retrieved and calculated values for both models are listed in Table 5.

While our model places no constraint on the age of Gl 229B, we present results showing the plausibility that the spectrum and atmospheric chemistry can be fit consistently as well by a 70 M_Jup model as it can by a 50 M_Jup model.

In Table 6 we list our SED derived and retrieved fundamental parameters as compared to values from the literature for Gl 229B. We find that our semi-empirical SED derived L_Bol, radius, mass and log(g) agree within 1σ to our best fit model retrieved parameters, while T_eff agrees within 1.2σ.

Table 6. Comparison of Fundamental Parameters for Gliese 229B

Source	Radius (RJup)	Mass (MJup)	log g	Teff (K)	log(LBol/LSun)	C/O	[M/H]	Age (Gyr)
This Paper a	0.94 ± 0.15	41 ± 24	4.96 ± 0.46	927 ± 79	−5.21 ± 0.05	⋯	0.0	0.5–10
This Paper b	1.10 ± 0.04	${50}_{-10}^{+12}$	${5.01}_{-0.09}^{+0.10}$	${834}_{-15}^{+17}$	$-{5.25}_{-0.01}^{+0.02}$	${1.38}_{-0.07}^{+0.08}$	−0.24 ± 0.05	⋯
This Paper c	1.12 ± 0.04	71 ± 1	5.15 ± 0.03	${829}_{-17}^{+19}$	$-{5.25}_{-0.01}^{+0.02}$	1.45 ± 0.07	−0.19 ± 0.03	⋯
Naka95	⋯	20–50	⋯	<1200	−5.398	⋯	⋯	0.5–5
Naka15	⋯	30–38	4.87 ± 0.12	825 ± 25	⋯	⋯	0.13 ± 0.07	1–2.5
Legg02	⋯	>7	3.5 ± 0.5	1000 ± 100	−5.21 ± 0.02	⋯	−0.5	0.016-0.045
Legg99	⋯	25–35	⋯	900	−5.18 ± 0.04	⋯	⋯	0.5–1
Alla96	⋯	40–58	5.3 ± 0.2	1000	−5.21 ± 0.1	⋯	⋯	0.5–5
Matt96	⋯	⋯	⋯	913	−5.194	⋯	⋯	⋯
Saum00	⋯	15–73	5.0 ± 0.5	950 ± 80	−5.21 ± 0.04	⋯	−0.3 ± 0.2	>0.2
Bran20	⋯	70 ± 5 d	5.433 ± 0.033	1025 ± 15	−5.208 ± 0.007 e	⋯	⋯	7–10
Howe22	1.11 ± 0.03	41.6 ± 3.3	${4.93}_{-0.03}^{+0.02}$	${869}_{-7}^{+5}$	⋯	1.13 ± 0.03	-0.07 ± 0.03	⋯

Notes.

^a SED. ^b Best Fit Model. ^c Mass Constrained Best Fit Model. ^d This dynamical mass has since been updated in Brandt et al. (2021) with improved astrometry from Gaia EDR3 (Gaia Collaboration et al. 2021) to be 71.4 ± 0.6. ^e This value is taken from Filippazzo et al. (2015) which we have updated in this work.

References. Naka95: Nakajima et al. (1995), Naka15: Nakajima et al. (2015), Legg02: Leggett et al. (2002a), Legg99: Leggett et al. (1999), Alla96: Allard et al. (1996), Matt96: Matthews et al. (1996), Saum00: Saumon et al. (2000), Bran20: Brandt et al. (2020), Howe22: Howe et al. (2022).

Download table as:  ASCIITypeset image

Comparing our retrieval-derived log(g) to values published in the literature, we find agreement within 1σ between our best fit model and all literature predictions except Leggett et al. (2002a) and Brandt et al. (2020), however, we note that the published log(g) from Leggett et al. (2002a) is in 2σ disagreement with all other published values due to the young age estimate of the system. Our mass constrained model log(g) is in agreement within 1σ with Allard et al. (1996) and Saumon et al. (2000) and within 2.5σ to Nakajima et al. (2015). Our constrained model is not in agreement with the reported log(g) from Brandt et al. (2020), however, the reported uncertainties for both values are notably small.

The T_eff we derive for both the best fit and mass constrained model is cooler than all literature values except Nakajima et al. (2015), which agrees within 1σ. We retrieved the same L_Bol for both constrained and best fit models and find its value is slightly smaller than those reported in the literature but still see 2σ agreement to nearly all models, despite small uncertainties.

In Table 7, we list reported metallicities for Gl 229A, which, similar to reported values for Gl 229B, range from subsolar to supersolar. For both best fit and mass constrained retrieval models, we find a subsolar metallicity which is within 1σ and 2σ of the value reported for Gl 229A in Schiavon et al. (1997) and Neves et al. (2014), respectively. However, our derived metallicity skews subsolar due to retrieved oxygen abundances that are lower than expected for both models. Subsequently, we use carbon as a tracer for overall object metallicity (J. Gaarn et al. 2022, in preparation) and calculate [C/H] = −0.01 ± 0.05 (best fit model) and [C/H] = 0.05 ± 0.03 (mass constrained model). A roughly Solar metallicity is in agreement with all published values for Gl 229A, except Leggett et al. (2002a) which reports a distinctly subsolar metallicity.

We list in Table 6 the retrieval-based C/O ratio for Gl 229B and in Table 7 the C/O ratio reported for Gl 229A from Nakajima et al. (2015). Nakajima et al. (2015) determines a C/O ratio using carbon and oxygen abundances from Tsuji & Nakajima (2014) and Tsuji et al. (2015), a notoriously difficult task for M dwarfs due to molecular absorption features throughout their spectra. Tsuji & Nakajima (2014) and Tsuji et al. (2015) use CO as an indicator of bulk carbon abundance and H₂O as an indicator of bulk oxygen abundance, finding logA_C = −3.27 ± 0.07 and logA_O = −3.10 ± 0.02, respectively. The C/O ratio we find in this work is calculated based on the abundances of all retrieved carbon-bearing molecules (CH₄, CO) as compared to all retrieved oxygen-bearing molecules (H₂O, CO). We apply a correction to the oxygen abundance (30% increase) to account for oxygen sequestered in silicate grains deeper in the atmosphere (Line et al. 2015, 2017; Zalesky et al. 2019, 2022). For this correction, we assume the primary oxygen sink is enstatite (Mg₂Si₂O₆) and therefore account for the removal of three oxygen atoms for every magnesium silicate. We do this in an attempt to probe at a bulk, as opposed to atmospheric, C/O ratio for Gl 229B, but note that this correction is still an approximation of the chemistry happening in the interior. The carbon and oxygen abundances of Gl 229A, B are listed in Table 8 for ease of comparison.

Table 8. Carbon and Oxygen Abundances

	log AC	log AO	C/O Ratio
Gl 229A	−3.27 ± 0.07 a	−3.10 ± 0.02 b	0.68 ± 0.12 c
Gl 229B (BF)	−3.35 ± 0.05	$-{3.49}_{-0.04}^{+0.05}$	${1.38}_{-0.07}^{+0.08}$
Gl 229B (BF) d	−3.35 ± 0.05	$-{3.38}_{-0.04}^{+0.05}$	${1.06}_{-0.05}^{+0.06}$
Gl 229B (MC)	−3.29 ± 0.03	−3.45 ± 0.04	1.45 ± 0.07
Gl 229B (MC) d	−3.29 ± 0.03	−3.34 ± 0.04	${1.11}_{-0.05}^{+0.06}$

Notes. (BF) indicates results from the best fit retrieval model while (MC) indicates results from the best fit model with a mass constraint.

^a Tsuji & Nakajima (2014). ^b Tsuji et al. (2015). ^c Nakajima et al. (2015). ^d Corrected oxygen abundance and subsequent C/O ratio (see Section 8).

Download table as:  ASCIITypeset image

It is evident for both best fit and mass constrained models that our C/O ratios are not in agreement with that of the primary. Our best fit retrieved C/O is in agreement within 3σ and our mass constrained, 4σ. Despite the median reported value for Gl 229A being supersolar, our C/O ratios are both >1. C/O ratio is theorized as a tracer of formation mechanism such that stellar companions with elevated C/O ratios as compared to their primary likely formed due to core accretion in the disk (e.g., Lodders 2004; Madhusudhan et al. 2011; Öberg et al. 2011; Madhusudhan 2012; Konopacky et al. 2013). While we do find higher C/O ratios in our models for Gl 229B, it is extremely unlikely that an object ≥50 M_Jup formed by core accretion in a disk around an M1V star (Bowler et al. 2015; Schlaufman 2018; Mercer & Stamatellos 2020). Furthermore, we compare the abundances of carbon and oxygen and find that the carbon values for both the best fit and mass constrained model agree with that of the primary within 1.2σ and 1σ, respectively. We find this as evidence of coevality and consider the disagreement in oxygen abundance as a result of unexplained chemistry in the atmosphere of Gl 229B. As in Line et al. (2017), we consider that a better understanding of the thermo-chemical mechanisms producing oxygen-bearing condensates is required to estimate a true oxygen abundance in cooler atmospheres. While magnesium silicates are the dominating oxygen sink in brown dwarf atmospheres, we hypothesize additional oxygen sinks, such as iron silicates, could contribute to an oxygen-depleted atmosphere. However, an origin for such influential abundances of alternative metals is unknown.

To place our work in context with another retrieval analysis on Gl 229B, we turn again to the results from Howe et al. (2022). As a result of the increased abundance of oxygen-bearing absorbers as compared to our best fit model (see Section 7), it is not surprising that they find a relatively smaller C/O ratio of 1.13 ± 0.03. It is also important to note that this work does not make an oxygen correction as we do, which would likely decrease this ratio even further. However, our work does align with Howe et al. (2022) in finding an overall supersolar C/O ratio, which, in their work, was speculated to be a result of unresolved disequilibrium chemistry in the atmosphere of Gl 229B.

There has been a noticeable trend developing toward supersolar C/O ratios in retrieval models of late T dwarfs, to which our work on Gl 229B now contributes (Figure 7). It is unclear whether this is a nod toward formation pathways, unresolved atmospheric chemistry in these types of objects or a bias in retrieval codes. This is certainly an active area of research that highlights the importance of brown dwarf companions, where studies can be anchored by the chemistry of the primary star (e.g., J. Gaarn et al. 2022, in preparation). In this section, we return to the retrievals of Zalesky et al. (2019, 2022) and Line et al. (2017) to add our work to the growing chemical trends previously reported in these works.

Zoom In Zoom Out Reset image size

Of the 11 T dwarfs retrieved in Line et al. (2017), seven are isolated brown dwarfs with no potential for comparative chemistry to a Solar-type star. However, whether companion or not, a similar trend emerges for the entire sample of T dwarfs toward supersolar C/O ratio, even after the implemented oxygen correction factor. Of these seven isolated objects, six were reported to have supersolar C/O ratios in the range 0.63– 1.14. It is important to make a distinction here that, as a result of unconstrained CO and CO₂, Line et al. (2017) calculated both C/O ratio and metallicity using CH₄ and H₂O abundances only. We note this because non-equilibrium chemistry, particularly with CO, could be a factor in the oxygen depletion seen in these retrievals. As a result, unconstrained abundances for CO and CO₂ could possibly be driving supersolar C/O ratios rather than unresolved chemistry or code biases.

Next, we consider our work against Zalesky et al. (2019, 2022) which initially presented retrieval models of a sample of 14 isolated late-type T dwarfs and Y dwarfs (T8-Y1), later expanded to a sample of 50 late-type T dwarfs (T7-T9). Of the six T dwarfs in the Zalesky et al. (2019) sample (T8-T9.5), they again find a trend toward supersolar C/O ratio in the late T regime with derived values in the range 0.7– 1.47(with exception for one T9 object with a retrieved C/O of 0.57 ± 0.07). In this work, the C/O ratio is calculated using H₂O, CH₄, CO and CO₂, despite only having upper limits on CO and CO₂. Consistent with Line et al. (2017), Zalesky et al. (2019) also includes a C/O correction factor of 30%. Unlike our work, Zalesky et al. (2019) reports slightly enhanced metallicities for this subset of T dwarfs. However, they report no apparent trend between high C/O and metallicity for late-T objects. Building on this initial study, Zalesky et al. (2022) also reports supersolar C/O ratios in range 0.65–1.28 for approximately 75% of their sample. For the comparative objects we discussed in 7, WISEJ024124.73 − 365328.0 has a retrieved C/O ratio of 0.82 ± 0.07, WISEJ112438.12 − 042149.7 a value of ${0.83}_{-0.08}^{+0.07}$, and WISEPCJ221354.69 + 091139.4, ${0.63}_{-0.08}^{+0.12}$. As our chemical abundances for Gl 229B were most similar to WISEJ024124.73 − 365328.0 and WISEJ112438.12 − 042149.7, it is not surprising that we find our C/O ratio closer in value to these objects than WISEPCJ221354.69 + 091139.4. While our model for Gl 229B still reports a higher C/O ratio than all three of these temperature comparisons, we can still place our model alongside this apparent trend of supersolar C/O ratios in T dwarf retrievals. We visualize this trend in 7.

The appeal of studying companion objects lies largely in the draw of understanding formation pathways. For brown dwarfs specifically, we hope that by comparing the atmospheric chemistry to that of its primary star (or brown dwarf companion) we can probe at the system age and make a determination as to whether the two objects formed together and if they formed in the same way (i.e., via gravitational fragmentation). As brown dwarf formation mechanisms are still unclear, studying co-moving systems provides one possible path to make advances in this topic.

In this section, we compare our findings on the Gl 229 system to the comparative chemistry of co-moving systems Gl 570 and HD 3651 (Line et al. 2015, 2017), Indi (Kitzmann et al. 2020), HR 8799 (Konopacky et al. 2013; Ruffio et al. 2021) and HR 7672 (Wang et al. 2022). We also plot the C/O ratios of these brown dwarf companion systems in Figure 7. These systems are of particular interest in relation to this work since they are all subjects of retrieval studies. Since Gl 570D, HD 3651B and Indi Bb are all late-T type objects, we would expect similar trends to emerge chemically (i.e., a high C/O ratio due to sequestered oxygen). However, we use the HR 8799 and HR 7672 systems as potential contrast, as we would not necessarily expect the same thermodynamic behavior in objects significantly hotter or cooler than Gl 229B. By placing the Gl 229 system in context with other retrieved co-moving systems, we can investigate trends in chemical abundance across companion objects with the ultimate goal of understanding underlying brown dwarf chemistry and formation.

Gl 570 (K4V+T7.5):Gl 570D is an ideal companion candidate, as it is widely separated from a well-studied main sequence star, Gl 570A. Line et al. (2015) reports metallicity and C/O ratio for the primary as 0.05 ± 0.17 and 0.81 ± 0.16, respectively. As previously discussed, Line et al. (2015, 2017) applies a correction to the retrieved atmospheric oxygen abundance to probe at bulk abundance for Gl 570D, giving metallicity and C/O ratio values of $-{0.15}_{-0.09}^{+0.07}$ and ${0.79}_{-0.23}^{+0.28}$, respectively. While the spread for C/O spans the range of solar to supersolar and metallicity spans subsolar to solar, there is agreement to 1σ for these chemical parameters between Gl 570D and its primary, which is in contrast to our model.

HD 3651 (K0V+T7.5):HD 3651B is another widely separated companion to a well-studied main sequence star, HD 3651A. Line et al. (2015) reports metallicity and C/O ratio for the primary as 0.18 ± 0.07 and 0.62 ± 0.11, respectively. The oxygen-corrected metallicity and C/O ratio of HD 3651B are ${0.08}_{-0.06}^{+0.05}$ and ${0.89}_{-0.17}^{+0.21}$, respectively. Clearly, this system shares a similarity to our work in that the primary star is of solar C/O and its brown dwarf companion is distinctly supersolar.

Indi (K4.5V + T1.5 + T6.5): The Indi system is another useful example where the main sequence primary can help put constraints on fundamental parameters. King et al. (2010) reports a metallicity measurement of −0.2 for the primary which is in agreement only with Indi Bb, the T6.5 dwarf, at a retrieved value of $-{0.34}_{-0.11}^{+0.12}$. The subsolar metallicity is predicted by the chemistry of the primary star. However, this work does not correct for condensation of oxygen-bearing condensates, so we would expect better agreement after this assumption. A supersolar C/O ratio reported for both Indi Bb at 0.84 ± 0.07 and Indi Ba at 0.95 ± 0.03 is also predicted to be closer to solar after a correction. Similar to Line et al. (2015, 2017) this work was not able to place constraints on CO, CO₂, or H₂S, but does, in fact, calculate a C/O with all oxygen- and carbon-bearing molecules. Another interesting feature of this work, similar to ours, is that dynamical measurements predicted both brown dwarf companions to be ∼50 - 70 M_Jup (Dieterich et al. 2018; Chen et al. 2022), while retrieval models place Ba and Bb at ${50}_{-6.6}^{+7.8}$ and ${15}_{-4.6}^{+6}$ M_Jup, respectively.

HD 7672 (G0+L4): HD 7672 is another system in which the solar-type primary is well studied, with reported abundances [Fe/H] = −0.04 ± 0.07, [C/H] = −0.08 ± 0.05, [O/H] = 0.15 ± 0.06 and C/O = 0.56 ± 0.11 (Wang et al. 2022). Retrieval-derived abundance for the brown dwarf companion are reported as [C/H] = −0.24 ± 0.05, [O/H] = −0.19 ± 0.06 and C/O = 0.52 ± 0.02. This is an interesting divergence from our work where the C/O ratio is within 1σ agreement to its primary, but [C/H], which we use as a proxy for overall object metallicity, is subsolar. Since this brown dwarf companion is an L-type, we should not need to correct for oxygen-bearing condensates, although Wang et al. (2022) does state that retrieved carbon and oxygen abundances are atmospheric (as opposed to bulk) and suggests that inefficient mixing could be a cause of these metallicity differences.

HR 8799 (A5/F0+bcde): This system provides an interesting contrast to ours, as this young star has four substellar companions with masses in the range 5–13 M_Jup orbiting within 15–70 AU, which is beyond the H₂O ice line in the disk (Konopacky et al. 2013). Initial studies from Konopacky et al. (2013), found enhanced C/O ratio for the b, c companions as compared to the primary, which is approximately solar at ${0.54}_{-0.09}^{+0.12}$, suggesting a core accretion scenario. However, Ruffio et al. (2021) concluded that all four planets have stellar C/O ratios, which could imply either core accretion or gravitational instability as a formation mechanism. This result is in agreement with retrieval results on HR 8799c,e (Mollière et al. 2020; Wang et al. 2020) but in disagreement with the retrieval study of Lavie et al. (2017), which found all four companions to be oxygen-enriched relative to their star. This opens up a particularly interesting discussion on the usefulness of the C/O ratio as a formation diagnostic for objects near or above the deuterium burning limit. We include HR 8799c in Figure 7 as a contrast to higher mass, brown dwarf companions that show distinct supersolar C/O trends, as well as chemical non-uniformity to their primary star.

While the C/O ratio has the potential to be a powerful formation diagnostic, the complications that arise in our attempt to determine bulk C/O ratios in substellar objects currently prevent us from using it as conclusive evidence. While we expect oxygen depletion in T-type objects, we still see unexpectedly incongruous C/O ratios for the Gl 229, HD 3651 and Indi systems. However, retrieval studies on Gl 570 and HD 7672 were able to match C/O ratio across companions within 1σ. Additionally, in the HR 8799 system where we initially predicted superstellar C/O ratios, retrieval results found substellar or stellar ratios. From these varied retrieval results so far, we do not see a unique chemical trend emerging (i.e., anomalously high C/O ratios for all T-type objects, indicating unresolved chemistry in that temperature regime or anomalously high C/O ratios for all brown dwarfs, suggesting a retrieval code bias). However, we also consider the possibility that both unresolved chemistry and code bias are contributing to the retrieval results we present.