Helios-r2: A New Bayesian, Open-source Retrieval Model for Brown Dwarfs and Exoplanet Atmospheres

For a given source function S_ν, the radiative transfer equation in a plane-parallel, semi-infinite atmosphere has a simple solution (e.g., Mihalas 1978), given by

$$\begin{eqnarray}&&{I}_{\nu }^{+}({\tau }_{\nu },\mu )={\int }_{{\tau }_{\nu }}^{\infty }{S}_{\nu }(t){e}^{-(t-{\tau }_{\nu })/\mu }{dt}/\mu ,\end{eqnarray} \tag{ 1 }$$

where ${I}_{\nu }^{+}$ is the outgoing intensity, i.e., for 0 < μ ≤ 1. In the following, we neglect scattering. Thus, the source function is simply given by

$$\begin{eqnarray}&&{S}_{\nu }({\tau }_{\nu })={B}_{\nu }(T({\tau }_{\nu })),\end{eqnarray} \tag{ 2 }$$

where ${B}_{\nu }(T({\tau }_{\nu }))$ is the Planck function with a temperature T at a given optical depth ${\tau }_{\nu }$. Equation (1) can be formally integrated with respect to the angular variable μ to yield the angular moments of the radiation field, such as the mean intensity or the flux. The resulting equations are known as the Schwarzschild–Milne equations and for the outgoing flux ${F}_{\nu }^{+}$ given by

$$\begin{eqnarray}&&{F}_{\nu }^{+}({\tau }_{\nu })=2\pi {\int }_{{\tau }_{\nu }}^{\infty }{S}_{\nu }({t}_{\nu }){E}_{2}({t}_{\nu }-{\tau }_{\nu }){{dt}}_{\nu },\end{eqnarray} \tag{ 3 }$$

where E₂ is the second exponential integral.

While it is possible to directly integrate Equation (3) to obtain the outgoing flux at the upper atmosphere (${\tau }_{\nu }=0$), this approach might lead to numerical and computational difficulties. Evaluating the exponential integrals can be costly in terms of computational time and becomes unstable at high optical depths. Additionally, numerically integrating the equation by using the trapezoidal rule can lead to rather large numerical errors that accumulate along the path of integration unless a high vertical resolution is used.

To circumvent these problems, we instead employ the so-called short characteristics method here, first described by Olson & Kunasz (1987). Essentially, this method solves the characteristic of the radiative transfer equation on a layer-by-layer basis. Additionally, to stabilize the integration, a weighting function of the form e^−τ/μ is introduced. For a single layer, Equation (1) can be written as

$$\begin{eqnarray}&&{I}_{\nu }^{+}({\tau }_{\nu ,i},\mu )={I}_{\nu }^{+}({\tau }_{\nu ,i-1},\mu ){e}^{-{\tau }_{\nu ,i}}+{\rm{\Delta }}{I}_{\nu ,i}^{+}({S}_{\nu },\mu ),\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu ,i}^{+}({S}_{\nu },\mu )={\alpha }_{i}{S}_{i+1}+{\beta }_{i}{S}_{i}+{\gamma }_{i}{S}_{i-1}.\end{eqnarray} \tag{ 5 }$$

Here, τ_ν,i refers to the optical depth in the ith layer, while α, β, and γ are coupling coefficients that connect the adjacent layers. Without scattering, the coefficient α is zero for outgoing rays. By assuming that the source function varies linearly within the layer, the coefficients are given by (Olson & Kunasz 1987)

$$\begin{eqnarray}&&{\beta }_{i}=1+\displaystyle \frac{{e}^{-{\rm{\Delta }}}-1}{{\rm{\Delta }}}\quad {\gamma }_{i}={e}^{-{\rm{\Delta }}}-\displaystyle \frac{{e}^{-{\rm{\Delta }}}-1}{{\rm{\Delta }}},\end{eqnarray} \tag{ 6 }$$

with

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{{\tau }_{i}-{\tau }_{i-1}}{\mu }.\end{eqnarray} \tag{ 7 }$$

It is also possible to use higher-order interpolants of the source function. For the parabolic case, the coefficients can be found in Olson & Kunasz (1987). While offering higher accuracy per layer, these coefficients are also computationally more expensive.

We solve Equation (4) for a set of angles μ, distributed according to a Gauß quadrature scheme (Gauß 1814) and then numerically integrate the intensities to obtain the outgoing flux

$$\begin{eqnarray}&&{F}_{\nu }^{+}=2\pi {\int }_{0}^{1}{I}_{\nu }(\mu )\mu \,d\mu .\end{eqnarray} \tag{ 8 }$$

For this study, we use two angles in the upward direction. This is equivalent to a four-stream discrete ordinate radiative transfer, which is usually sufficient for a problem without scattering and an atmosphere where the scale height is much smaller than the planet's radius. Compared to the direct solution of Equation (3), this method only requires the evaluation of two exponentials per layer.

The CPU part of Helios-r2 is also equipped with the radiative transfer library CDISORT (Stamnes et al. 1988; Hamre et al. 2013). This radiative transfer model uses a complex, multi-stream approach to solve the radiative transfer equation and provides the exact solution to the problem if the number of streams is large enough. While this library is usually too slow to run within a retrieval, we use it to verify the accuracy of our implementation of the short characteristic method.

For the spectral resolution, we usually use a constant step size of 1 cm⁻¹ in wavenumber space. This is equivalent to the one used by Line et al. (2015) who studied similar objects.

To relate the outgoing fluxes ${F}_{\nu }^{+}$ computed via Equation (8) to the actual ones measured by the observer (F_ν), the radius and distance of the brown dwarf need to be taken into account. We therefore scale the top-of-the atmosphere fluxes ${F}_{\nu }^{+}$ by the usual radius–distance relation

$$\begin{eqnarray}&&{F}_{\nu }={F}_{\nu }^{+}f{\left(\displaystyle \frac{{R}_{* }}{d}\right)}^{2},\end{eqnarray} \tag{ 9 }$$

where d is the distance to the brown dwarf, R_* is its radius, and f is a calibration factor.

For the retrievals performed in Section 5, we choose a radius of ${R}_{* }=1{R}_{{\rm{J}}}$ and use distances from the Gaia measurements. The calibration factor f is treated as a free parameter that describes the uncertainties in the flux calibration of the measured spectra, but it also includes the deviations of the actual brown dwarf radius from our assumed value of 1 R_J. Lastly, f also partially captures inadequacies of the forward model to describe the atmosphere of a brown dwarf in all its details, including the effect of a reduced emitting surface due to a potentially heterogeneous atmosphere.

Assuming that f only includes deviations with respect to the assumed a priori radius, it can be transformed into a derived radius (in the same units as R_*) via

$$\begin{eqnarray}&&R=\sqrt{f}.\end{eqnarray} \tag{ 10 }$$

It should, however, be noted that in full generality, Equation (10) does not provide a good radius estimate for the brown dwarf because f usually also includes other sources of uncertainties, as described above.

This work is focused on the wavelength range of the SpeX instrument, which extends from about 0.9 μm to roughly 2.4 μm. Within this range, we account for the major absorbers that are expected to be present (Line et al. 2015): CO₂, CO, CH₄, H₂O, H₂S, NH₃, H₂, He, as well as the alkali metals Na and K.

Calculations of the molecular line absorption cross-sections are done with our opacity calculator Helios-k (Grimm & Heng 2015). We use the ExoMol line lists where available (Barber et al. 2006; Yurchenko et al. 2011; Yurchenko & Tennyson 2014; Azzam et al. 2016), and the ones provided by HITEMP (Rothman et al. 2010), otherwise. The line wings are modeled by Voigt profiles (Voigt 1912) that describe the effects of thermal and pressure broadening. Additional details on the calculations can be found in S. L. Grimm et al. (2020, in preparation).

Collision-induced absorption of H₂–H₂ (Abel et al. 2011) and H₂–He (Abel et al. 2012) pairs are taken into account by the corresponding data provided within the HITRAN database (Karman et al. 2019).

2.3.1. Alkali Line Absorption Cross-sections

Within our forward model, we employ two different approximations for the description of the atmosphere's chemical composition. We either perform a free retrieval of the molecules' mixing ratios or use an equilibrium chemistry code to self-consistently calculate the molecular abundances.

For the equilibrium chemistry, we employ the FastChem model (Stock et al. 2018). More specifically, we here use version 2.0 of the model that features several enhancements over the previous version. Compared to FastChem 1.0, the new version does not require a pressure iteration and is valid for arbitrary elemental compositions (J. W. Stock et al. 2020, in preparation).

We use our standard set of about 550 chemical species that is included in the released version of FastChem. Due to the low temperatures expected for the brown dwarfs we aim to investigate, ions are, however, removed from the chemical network of FastChem. Their small abundances at low temperatures would increase the computational time of the chemistry by a factor of roughly two or more but, on the other hand, do not significantly change the number densities of the more important molecules.

A full chemistry calculation for a given temperature and pressure usually takes of the order of a few milliseconds down to one millisecond or less (in case of higher temperatures) on a single CPU. It is, thus, possible to run the chemistry model directly within the forward model.

The current version of FastChem does, however, not treat condensation. In equilibrium condensation chemistry models, the alkali metals are expected to condense into Na₂S and KCl for cool T-dwarf atmospheres (Marley et al. 2013). We simulate this effect by removing K and Na from the gas phase in the upper atmosphere once the temperature–pressure profile drops below 800 K.

One of the most important quantities that is required for the forward model is the temperature profile. In the past, several different approaches to this problem have been used. This includes using a profile described by a nine-parameter fitting function (Madhusudhan & Seager 2009), using an approximate solution of the radiative transfer equation under the condition of radiative equilibrium in a gray atmosphere (e.g., Lavie et al. 2017), or using a free layer-by-layer temperature retrieval (Irwin et al. 2008; Line et al. 2015). In this study, we will explore two different scenarios, one using the gray atmosphere approximation and a modified, free-temperature retrieval based on finite elements.

The temperature profile based on the solution of the radiative transfer equation in a gray, semi-infinite atmosphere under the condition of radiative equilibrium—usually referred to as Milne's problem—is given by (Mihalas 1978)

$$\begin{eqnarray}&&{T}^{4}({\tau }_{R})=\displaystyle \frac{3}{4}{\bar{T}}_{\mathrm{eff}}^{4}\left[{\tau }_{R}+{q}_{{\rm{H}}}({\tau }_{R})\right],\end{eqnarray} \tag{ 11 }$$

where τ_R is the Rosseland optical depth, ${\bar{T}}_{\mathrm{eff}}$ an effective temperature, and q_H(τ_R) the so-called Hopf function. If one assumes the Eddington approximation, i.e., the second (K) and the zeroth moment (J) of the specific intensity are related via $3K=J$, then q_H = 2/3. However, while the Eddington approximation is a good description for the deep interior of the atmosphere, its validity breaks down in the upper parts of the atmosphere, where the Eddington factor approaches unity. We, therefore, use the exact solution of the function q_H(τ_R) provided in Mihalas (1978), p. 72, to describe the temperature profile.

We note that ${\bar{T}}_{\mathrm{eff}}$ should not be interpreted as the actual effective temperature of the brown dwarf's atmosphere that we aim to retrieve. It is defined as the effective temperature of a gray atmosphere, where the spectral distribution of the radiation is a pure blackbody. The atmosphere of a brown dwarf, however, is far from being gray, and its spectrum does not resemble a blackbody curve at all because it is dominated by deep molecular absorption bands. Accordingly, we determine the actual effective temperature by a post-processing procedure (see Section 2.6 for details).

2.5.1. Finite Element Approach

For this approach, we partition the atmosphere into a number K_e of non-overlapping elements, distributed in log-pressure space (see Figure 1). These elements are not aligned to and in fact are fully independent of the previously mentioned discretization of the atmosphere into layers/levels.

Zoom In Zoom Out Reset image size

The discrete approximation ${T}_{h}^{k}(\mathrm{log}p)$ of the real (usually unknown) temperature ${T}^{k}(\mathrm{log}p)$ within each element is expressed such that

$$\begin{eqnarray}&&{T}_{h}^{k}(\mathrm{log}p):= \sum _{i=1}^{{N}_{p}}{T}_{h}^{k}(\mathrm{log}{p}_{i}^{k})\,{{\ell }}_{i}^{k}(\mathrm{log}p)\end{eqnarray} \tag{ 12 }$$

is a polynomial of order q on each element. Here, the functions ℓ_i^k are the Lagrange polynomials (Waring 1779; Lagrange 1795) through the grid points $\mathrm{log}{p}_{i}^{k}$, i.e.,

$$\begin{eqnarray}&&{{\ell }}_{i}^{k}(\mathrm{log}p):= \prod _{\displaystyle \genfrac{}{}{0em}{}{j=1}{j\ne i}}^{{N}_{p}}\displaystyle \frac{(x-{x}_{j})}{({x}_{i}-{x}_{j})}.\end{eqnarray} \tag{ 13 }$$

In a finite element approach, these would be the so-called trial functions. For a given order q, the number of local grid points N_p is given by

$$\begin{eqnarray}&&{N}_{p}=q+1.\end{eqnarray} \tag{ 14 }$$

Figure 1 shows an example of two second-order elements. In that case, every element has three nodes $\mathrm{log}{p}_{i}^{k}$, associated with three degrees of freedom (dof), which are the temperatures ${T}_{h}^{k}(\mathrm{log}{p}_{i}^{k})={T}_{i}^{k}$. We use the continuous formulation of the approximate solution in the following, i.e., the temperatures are continuous across element boundaries. This enforces a continuous temperature profile, which is expected in a brown dwarf atmosphere, and also reduces the overall number of dof.

The temperature at the element interfaces can, in principle, also be made disjoint, allowing the solution to have discontinuous jumps from one element to the next. In terms of finite element methods, this would lead to the discontinuous Galerkin method (e.g., Kitzmann et al. 2016). This allows for more flexibility in the temperature profile—especially in cases where the temperature profile is under-resolved—but might also lead to unphysical discontinuities in the retrieved temperature profile. In this work, we therefore focus on the continuous version.

In the following, we assume a polynomial order larger than zero ($q\gt 0$). The case of q = 0 refers to an atmosphere that is piecewise isothermal, which clearly does not satisfy the expected, continuous temperature profile.

Given the representation of the temperature in each element by Equation (12), the global, piecewise continuous solution can then be written as the direct sum of all K_e solutions

$$\begin{eqnarray}&&T(\mathrm{log}p)\backsimeq {T}_{h}(\mathrm{log}p)=\underset{k=1}{\overset{{K}_{{\rm{e}}}}{\displaystyle \bigoplus }}{T}_{h}^{k}(\mathrm{log}p).\end{eqnarray} \tag{ 15 }$$

This representation allows the evaluation of the temperature at any desired pressure p within the atmosphere.

The total number of free parameters N_T for a temperature profile with K_e elements of order q is given by

$$\begin{eqnarray}&&{N}_{{\rm{T}}}={K}_{{\rm{e}}}{N}_{p}-\left({K}_{{\rm{e}}}-1\right)={K}_{{\rm{e}}}\,q+1.\end{eqnarray} \tag{ 16 }$$

Instead of performing the calculations in Equations (12) and (13) in the log p space, all evaluations are done on a reference element, stretching from 0 to 1 and then transformed back into the pressure space. On the reference element, the points N_p are distributed according to a Gauß–Lobatto quadrature scheme (Lobatto 1852).

It should be noted that the approximation of T(log p) by Equation (15) does not assume any specific mathematical form of the global temperature profile. Any sufficiently continuous and smooth function can be approximated by a piecewise polynomial description. Very complex temperature profiles (e.g., temperature inversions, strong gradients) may require either a high-order polynomial or a larger number of elements. While Helios-r2 is designed to use any polynomial degree larger than unity, we normally restrict the retrieval to second-order elements. High-order polynomials are prone to suffer from so-called Runge's phenomenon (Runge 1901), i.e., they are susceptible to unphysical, local oscillations.

It would be natural to use the nodal values of the temperatures ${T}_{h}^{k}(\mathrm{log}{p}_{i}^{k})$ as free parameters for the retrieval. In many cases, though, this would give the temperature profile too much freedom, sometimes yielding unphysical temperature inversions in the lower atmosphere. Since one would not expect such inversions to occur based on the theory from brown dwarf atmospheres, we use a slightly different approach for the retrieval here.

The only, actually retrieved temperature is the one at the bottom of the modeled atmosphere, represented by ${T}_{1}^{1}$ in Figure 1. For all other temperatures, we retrieve a parameter b, such that, e.g., ${T}_{2}^{1}={b}_{2}{T}_{1}^{1}$. For first-order elements, b can be interpreted as the slope between two adjacent temperature nodes. By choosing a value for b of smaller or equal to unity, the temperature profile will be strictly monotonic. We found this approach to be much more stable than retrieving all temperatures individually. The same formulation can also be adapted to atmospheres with temperature inversions, by allowing b to exceed a value of unity.

An example for a Helios-r2 retrieval of an exoplanet atmosphere with a temperature inversion can be found in a separate publication (Bourrier et al. 2019). Due to the higher complexity of the temperature profile, four second-order elements are used in Bourrier et al. (2019), while the b values for the description of the profile are allowed to exceed unity.

Since neither of the two approaches we implemented in Helios-r to describe the temperature profile directly yield the effective temperature T_eff of the brown dwarf, we estimate this parameter in a post-processing procedure. Following Line et al. (2015), we calculate a high-resolution spectrum for every posterior parameter combination in the wavelength range from 1 to 20 μm. These spectra are then integrated over these wavelengths for the total outgoing flux, which is then converted into an effective temperature by using the Stefan-Boltzmann law (Stefan 1879; Boltzmann 1884).

The new version of Helios-r can also take an instrument profile function into account when calculating theoretical spectra. This function describes the spread of the flux at a certain wavelength across several pixels on the detector due to, for example, the finite slit width of the spectrograph. We here assume this profile to have a Gaussian (normal) distribution. The calculated high-resolution spectrum from the forward model is convolved with the instrument profile p to simulate the flux at a given wavenumber ν measured by each pixel on the detector

$$\begin{eqnarray}&&{F}_{\nu ,{\rm{d}}}={\int }_{0}^{\infty }{F}_{x}\,p(x-\nu ,\sigma ){dx},\end{eqnarray} \tag{ 17 }$$

where σ is the standard deviation of the profile that is a characteristic of the employed instrument. In practice, we do not perform the convolution over the entire wavenumber range, as this would be computationally extremely costly. Instead, the integration is stopped at a distance of $| x-\nu | =5\,\sigma$. Flux values outside of this range have a negligible impact on ${F}_{\nu ,{\rm{d}}}$ and, thus, can be safely neglected.

To make the model computationally as fast as possible, the numerically heavy part of the calculations is done on a GPU by using NVIDIA's CUDA language. This includes, for example, the interpolation of opacities, calculation of the optical depths, the solution of the radiative transfer equation, the convolution with an instrument profile, the binning of the high-resolution spectrum to the observational data, and the evaluation of the likelihood function.

The computationally most expensive part is usually the interpolation of opacities to the corresponding atmospheric temperatures and pressures. This operation, however, can be parallelized quite straightforwardly and runs much faster on a GPU than on a CPU, owing to the several thousands of cores available for simultaneous calculations on a graphics card.

The only major part of the model still running on a CPU is FastChem, as it cannot be efficiently parallelized for a GPU. Instead, several instances of FastChem are running in parallel on the CPU using OpenMP, depending on the number of available CPU cores.

In total, a typical evaluation of the forward model with 70 layers, 7100 wavelength points (corresponding to wavenumber steps of 1 cm⁻¹), and the equilibrium chemistry on a GeForce 2080 Ti requires about 20 ms of calculation time. Increasing the resolution to 0.3 cm⁻¹ (about 25,000 wavelength points) results in a calculation time of roughly 26.5 ms per model.

In principle, Helios-r2 can also be purely run on a CPU. The calculation times, however, can then be a factor of more than 100 times higher than those obtained employing a GPU.

Atmospheric retrieval tries to connect observational data ${\boldsymbol{D}}$ of an object with a probability distribution of a parameter set ${\boldsymbol{\Theta }}$, that are usually connected to physically relevant properties of the observed atmosphere. The data vector ${\boldsymbol{D}}$ is usually composed of a set of data points D_j taken at, e.g., different wavelengths. In addition to ${\boldsymbol{D}}$, the observational data is also characterized at each spectral point by a corresponding error σ_j.

Let ${{ \mathcal M }}_{i}$ be a model with a parameter vector ${\boldsymbol{\Theta }}\,=\left\{{{\rm{\Theta }}}_{1},{{\rm{\Theta }}}_{2},...,{{\rm{\Theta }}}_{N}\right\}$, containing N parameters Θ_n. In terms of atmospheric retrieval, ${{ \mathcal M }}_{i}$ is usually a simplified atmosphere model (forward model) that calculates a theoretical spectrum of the observed object based on a set of input parameters, such as molecular abundances, surface gravity, or the temperature profile. The joint prior probability distribution for ${\boldsymbol{\Theta }}$, $p\left({\boldsymbol{\Theta }}| {M}_{i}\right)$, describes the a priori knowledge or constraints we impose on the initial distribution of the parameters Θ_n.

The posterior probability distribution of ${\boldsymbol{\Theta }}$ for a specific model ${{ \mathcal M }}_{i}$ applied to the data ${\boldsymbol{D}}$ can then be expressed following Bayes' theorem (Bayes & Price 1763)

$$\begin{eqnarray}&&p\left({\boldsymbol{\Theta }}| {\boldsymbol{D}},{{ \mathcal M }}_{i}\right)=\displaystyle \frac{p\left({\boldsymbol{\Theta }}| {M}_{i}\right){ \mathcal L }\left({\boldsymbol{\Theta }}| {\boldsymbol{D}},{{ \mathcal M }}_{i}\right)}{{ \mathcal Z }\left({\boldsymbol{D}}| {{ \mathcal M }}_{i}\right)},\end{eqnarray} \tag{ 18 }$$

where ${ \mathcal L }\left({\boldsymbol{\Theta }}| {\boldsymbol{D}},{{ \mathcal M }}_{i}\right)$ is the likelihood and ${ \mathcal Z }\left({\boldsymbol{D}}| {{ \mathcal M }}_{i}\right)$ the so-called Bayesian evidence.

The likelihood function ${ \mathcal L }\left({\boldsymbol{\Theta }}| {\boldsymbol{D}},{{ \mathcal M }}_{i}\right)$ describes the probability of the model ${{ \mathcal M }}_{i}$ to match the data ${\boldsymbol{D}}$, given a set of parameters ${\boldsymbol{\Theta }}$. We here use the same likelihood function previously employed by, e.g., Benneke & Seager (2013), Line et al. (2015), or Lavie et al. (2017). Assuming that the observational points j each possess individual Gaussian errors, the (logarithm) of ${ \mathcal L }$ can be expressed by

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=\sum _{j=1}^{J}\left(-\displaystyle \frac{{\left[{D}_{j}-{D}_{j,{\rm{m}}}\left({{ \mathcal M }}_{i},{\boldsymbol{\Theta }}\right)\right]}^{2}}{2{\sigma }_{j}^{2}}-\displaystyle \frac{1}{2}\mathrm{ln}(2\pi {\sigma }_{j}^{2})\right),\end{eqnarray} \tag{ 19 }$$

with the theoretical observation ${D}_{j,{\rm{m}}}\left({{ \mathcal M }}_{i},{\boldsymbol{\Theta }}\right)$, calculated by the forward model M_i using the parameters ${\boldsymbol{\Theta }}$. The Bayesian evidence is formally given by the integral

$$\begin{eqnarray}&&{ \mathcal Z }\left({\boldsymbol{D}}| {{ \mathcal M }}_{i}\right)=\int p\left({\boldsymbol{\Theta }}| {M}_{i}\right){ \mathcal L }\left({\boldsymbol{\Theta }}| {\boldsymbol{D}},{{ \mathcal M }}_{i}\right)d{\boldsymbol{\Theta }}.\end{eqnarray} \tag{ 20 }$$

Evidently, since this is a multidimensional integral over the entire parameter space, the actual, direct evaluation of this integral is quite challenging when the number of parameters is large.

The Bayesian evidence can also be used to perform model comparisons. For two different models M_i and M_j applied to the same data ${\boldsymbol{D}}$, one can compute the so-called Bayes factor

$$\begin{eqnarray}&&{B}_{{ij}}=\displaystyle \frac{{ \mathcal Z }\left({\boldsymbol{D}}| {{ \mathcal M }}_{i}\right)}{{ \mathcal Z }\left({\boldsymbol{D}}| {{ \mathcal M }}_{j}\right)}.\end{eqnarray} \tag{ 21 }$$

This factor quantifies the strength of evidence in favor of model M_i over M_j to represent the measured data. On the Jeffreys scale (Kass & Raftery 1995), a value of 1, 3.2, and 10 correspond to no, substantial, and strong evidence in favor of model M_i over M_j, respectively. A decisive evidence is categorized by B_ij > 100.

Nested sampling is essentially a method that provides the possibility of evaluating the Bayesian evidence ${ \mathcal Z }\left({\boldsymbol{D}}| {{ \mathcal M }}_{i}\right)$ and the posterior distributions $p\left({\boldsymbol{\Theta }}| {\boldsymbol{D}},{{ \mathcal M }}_{i}\right)$. In this approach, the multidimensional integral is reduced to a one-dimensional one over the so-called prior mass. A full description of the mathematical procedure of the nested sampling method can be found in Skilling (2004).

As a specific implementation of the nested sampling method, we use the MultiNest code (Feroz & Hobson 2008; Feroz et al. 2009). We use the FORTRAN version of the library and couple it directly to Helios-r2, written in C++/CUDA.

The MultiNest code starts by drawing N_l samples from the parameter space. The points N_l are referred to as live points. As mentioned in Section 3.1, the parameter values are each subject to their own, individual prior distribution. Using the Helios-r2 forward model, a theoretical spectrum is generated for each of parameter sets of the live points. Finally, the likelihood is then evaluated via Equation (19).

At each iteration step, MultiNest replaces the live point with the smallest likelihood value with a new set of values from the parameter space. This new point is chosen such that the likelihood computed for this point is higher than for the one just discarded. By repeating this process, the nested sampling will localize the regions of highest likelihood in the parameter space. To efficiently sample the parameter space, MultiNest employs a simultaneous ellipsoidal nested sampling method developed by Feroz et al. (2009). We refer to Feroz et al. (2009) for a detailed description on this method.

After ${ \mathcal Z }$ is converged, the posterior distributions of the parameters ${\boldsymbol{\Theta }}$ are constructed by using all active live points, as well as those who have been previously removed during the iterative procedure.

The number of live points must be high enough to allow for a good coverage of the parameter space. Depending on the number of free parameters and types of priors, several hundred to thousands of live points are usually required. Ideally, convergence tests on the required number of live points should be performed to ensure that the retrieval has converged properly. For the simple test cases in Section 4, we use between 2000 and 4000 live points, while the retrievals of the actual brown dwarf data is done with 10,000 points.

For the retrieval of actual brown dwarf data, we use a slightly different version of the log-likelihood. Following Line et al. (2015), we include an inflation of the observational errors the calculation of $\mathrm{ln}{ \mathcal L }$. Effectively, the usual squared error ${\sigma }_{i}^{2}$ is replaced by a more general data error s_i², given by

$$\begin{eqnarray}&&{s}_{j}^{2}={\sigma }_{j}^{2}+{10}^{\epsilon }.\end{eqnarray} \tag{ 22 }$$

The last factor on the right-hand side is used to slightly inflate the original observational error ${\sigma }_{i}$. With these s_j², the log-likelihood function from Equation (19) is then given by

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}\sum _{j=1}^{J}\displaystyle \frac{{\left[{D}_{j}-{D}_{j,{\rm{m}}}\left({{ \mathcal M }}_{i},{\boldsymbol{\Theta }}\right)\right]}^{2}}{{s}_{j}^{2}}-\displaystyle \frac{1}{2}\mathrm{ln}(2\pi {s}_{j}^{2}).\end{eqnarray} \tag{ 23 }$$

The exponent in the error inflation term is added as an additional retrieval parameter. For the corresponding prior of , we employ the same assumptions as Line et al. (2015) and use a uniform prior with

$$\begin{eqnarray}&&0.01\cdot \min {\sigma }_{i}^{2}\leqslant {10}^{\epsilon }\leqslant 100\cdot \max {\sigma }_{i}^{2}.\end{eqnarray} \tag{ 24 }$$

This error inflation accounts for the fact that the simplified model physics are usually not able to describe all of the details in a measured brown dwarf spectrum. Without error inflation, the nested sampling would concentrate only on the data points with the smallest errors, neglecting other important wavelength regions in the process. By inflating the error bars to a certain degree, we, thus, give the retrieval model more freedom to fit the spectrum, usually resulting in retrieved parameters that are more comparable with those expected from the theory of brown dwarf physics.

For the first test, we use the Helios-r2 forward model to produce a high-resolution spectrum. For the temperature profile, we use the model output from the Sonora grid for the aforementioned parameters. The high-resolution spectrum is then binned to about 150 bins from 1 to 2.4 μm.

It should be noted that due to the differences in the two atmospheric models (e.g., chemistry or opacities), the resulting effective temperature differs from the Sonora model atmosphere. The effective temperature derived from integrating the high-resolution spectrum of the Helios-r2 forward model is 689 K. This value is slightly lower than the corresponding value of the original Sonora model (700 K).

We simulate point-wise, uncorrelated noise by shifting each point by a flux value randomly drawn from a normal distribution with a standard deviation equal to 0.2 times the median of all fluxes. The error of each bin flux is estimated by using the median fractional error of an actual SpeX observations of GJ 570 D (see Section 5).

4.1.1. Fixed Temperature Profiles

In the first test, we fix the temperature profile to the one from the Sonora output and set the f parameter to its predetermined value of 1. Thus, we are only retrieving the elemental abundances, the C/O ratio, the surface gravity and, via the aforementioned post-processing procedure, the effective temperature T_eff. This first, trivial test is purely testing the ability of the forward model to recover a mock spectrum generated by the same radiative transfer model using the same chemistry and opacities. The corresponding posteriors for these parameters are shown in Figure 2.

Zoom In Zoom Out Reset image size

Unsurprisingly, this first, simple test results in an almost perfect match of the estimated parameters to the actual values used to produce the input spectrum. All parameters are tightly constrained with quite small confidence intervals. As expected, we also find the well-known correlation between the elemental abundances and the surface gravity because of their direct influence on the atmospheric scale height.

In the second test, we now add the f scaling as a new free parameter and repeat the retrieval. The temperature profile is still kept fixed to the one used to produce the simulated observation. The resulting posteriors are shown in Figure 3.

Zoom In Zoom Out Reset image size

Again, all parameters are well constrained and compare extremely well to their actual values. An interesting feature of this retrieval that can noticed outright, is the strong degeneracy between the effective temperature and the surface gravity. As we explain later, this will have a large impact on the actual retrieval of brown dwarf atmospheres.

4.1.2. Free-temperature Retrieval

In a final test of retrieving the output of the forward model, we now also include the temperature profile in the retrieval. We study two different scenarios: a temperature profile following Milne's solution (free parameters: Rosseland opacity κ_R and temperature ${\bar{T}}_{\mathrm{eff}}$), as well as a free-temperature retrieval with three second-order elements, comprised of, in total, seven free parameters (T₁ and six coefficients b_i). Details on these temperature profiles can be found in Section 2.6. The posteriors of the two retrievals are shown in Figure 4.

Zoom In Zoom Out Reset image size

The left panel of Figure 4 implies that Milne's solution is only a simple approximation to the actual temperature profile. In the deep interior, the retrieved temperatures well match the actual ones, shown in blue. This is to be expected because the approximations made for Milne's solution are valid at high optical depths (see Mihalas 1978 for details). On the other hand, these approximations become less valid in the upper atmosphere. As a result, Milne's solution starts to deviate from the actual profile for pressures less than 1 bar. It is, thus, not surprising that other retrieved parameters deviate from their actual values, most notably the f factor that is predicted with larger value (1.10) than its actual value of 1. Apparently, the retrieval model uses this parameter to mitigate the shortcomings in the temperature profile. This emphasizes the fact that f should not be seen as a purely radius-related parameter but that it also includes deviations due to assumptions made for the forward model physics.

As already mentioned in Section 2.6, the free parameter ${\bar{T}}_{\mathrm{eff}}$ of the Milne profile should not be interpreted as the actual effective temperature. As shown in Figure 4, the value of ${\bar{T}}_{\mathrm{eff}}$ is almost 200 K smaller than the effective temperature derived by integrating the high-resolution posterior spectra.

The free-temperature retrieval, on the other hand, provides an excellent representation of the temperature profile. The retrieved profile matches the original input profile in the lower atmosphere. In the upper atmosphere, the confidence intervals for the temperature profile become larger because the spectrum is almost insensitive to this part of the atmosphere at the spectral resolutions and wavelength range we are using here. The original temperature profile is, however, included in the 1σ envelope of the retrieved one. It should be noted that in comparison to Line et al. (2015), our approach needs no additional smoothing and requires less free parameters.

In addition to the posterior distributions of the retrieved parameters, Figure 5 depicts the median and 1σ confidence levels of spectra that have been calculated for all points within the two posterior sets. As the figure suggests, there is almost no visible difference between the two distributions. Both median spectra fit the simulated data almost perfectly. Furthermore, the posteriors of the two retrievals are so tightly constrained that the 1σ confidence ranges of the median spectra are basically invisible in Figure 5. Thus, just based on the ability of the retrieved posterior values to fit the simulated spectrum, one would not be able to exclude the Milne approximation as a valid solution of the problem. The Bayesian evidence for the free-temperature retrieval is $\mathrm{ln}{ \mathcal Z }=14918.29$ while the one for Milne's solution is 14879.69, respectively. With a Bayes factor of ln B = 38.6, the free-temperature retrieval is decisively favored over the Milne one, even though it requires more free parameters to describe the temperature profile. This emphasizes the decisive role played by the non-gray opacities in controlling the temperature profile.

Zoom In Zoom Out Reset image size

4.1.3. Temperature Profile Retrieval

Finally, we explore the impact of the number of elements and polynomial orders on the free-temperature retrieval. As mentioned in Section 2.5.1, we normally restrict the polynomial order to a maximum of two. Figure 6 shows the results for first- and second-order elements.

Zoom In Zoom Out Reset image size

The results clearly suggest that three first-order elements are not enough to fully describe the temperature profile. It does not provide enough degrees of freedom to cover the detailed behavior of it. On the other hand, four or six elements seem to be already sufficient. However, since they are piecewise linear functions, they still show a small level of roughness at the element boundaries.

Second-order elements are, by construction, smoother than first-order ones. Thus, they can follow the original temperature profile much more closely. This is especially noticeable when comparing the first-order, six-element case with the second-order, three-element one. Both have the same number of degrees of freedom, but the second-order elements provide a much smoother fit. Based on the results of Section 2.5.1, we use either three second-order elements or six first-order ones in the actual retrievals of this study.

In this section, we test our retrieval model on model output from the Sonora grid of brown dwarf atmospheres (M. S. Marley et al. 2020, in preparation), which uses different implementations of opacities, chemistry, and radiative transfer. As previously mentioned, we choose a model with an effective temperature of 700 K, a surface gravity of 4.75, as well as solar elemental abundances.

The Sonora element abundances, however, are always given with respect to the bulk composition, i.e., they include the elements in the gas phase, as well as those present in condensates. The retrieval model, on the other hand, is only sensitive to the abundances in the gas phase. Thus, the retrieved metallicities and C/O ratios can differ from the original Sonora input values. In order to provide a more consistent comparison with the rainout chemistry of Sonora, we remove several heavier elements from the FastChem equilibrium chemistry, such as Fe, Mg, or Si.

We also again assume a stellar radius of 1 Jupiter radius, a distance of 10 pc, and an f factor of 1. The simulated observation is created the same way as for the Helios-r2 test case. The comparison is performed for two different cases: with and without the calibration factor f. The temperature profile is freely retrieved in both cases.

The posteriors for the retrieval without the f parameter, shown in the upper panel of Figure 7, are well constrained, with only small standard deviations. The retrieved values for log g and the effective temperature are a bit less accurate than in the previous Helios-r2 test retrieval, owing to differences in the two atmospheric models. This is most likely caused by different opacity line lists, differences in the chemical networks, or radiative transfer methods. Nonetheless, the retrieved values are quite close to the ones from the Sonora grid.

Zoom In Zoom Out Reset image size

The metallicity derived by Helios-r2 is slightly sub-solar, while the C/O shows an enrichment in carbon compared to solar element abundances used by Sonora. This increased C/O ratio is expected because the latter model considers the removal of chemical species via condensation. Since this includes also oxygen-bearing condensates, oxygen will have a smaller than solar elemental abundance in the gas phase, which results in a super-solar C/O ratio.

In a second test, we now add the f factor to the retrieval (Figure 7, lower panel). The posteriors imply that the f parameter has a very strong impact on the other retrieval parameters. Instead of the expected value of 1, we obtain the much higher value of 1.39. This clearly also influences the other posterior distributions compared to the previous case without f. The effective temperature decreased slightly, while the surface gravity and the metallicity are affected more strongly. The retrieval seems to be partially misled by the f parameter and uses it to mitigate differences in the two atmospheric models. As already mentioned in Section 4.1.2, f should not be viewed as a parameter that is only related to the brown dwarf radius, as it also includes contributions due to, e.g., choices of molecular line lists or model physics assumptions.

The temperature profile is retrieved quite accurately in both cases. It follows the Sonora one in the lower atmosphere but is slightly shifted to lower temperatures in the case that includes the calibration factor as a free parameter. As expected, the profile has a wider confidence interval in the upper atmosphere where the spectrum becomes insensitive to the temperature. The Sonora profile, however, is still included within this interval in both cases.

The Bayes factor of these two retrievals is 1.37. Evaluating this factor using the Jeffreys scale (see Section 3.1) implies that there is very weak to no evidence of favoring either of the two different models.

The median spectra of both retrieval tests are shown in Figure 8. The figure clearly implies that both retrievals worked—in the sense that they provide a solution to fit the simulated observational data. Like in the previous case of the Helios-r2 retrieval test, both cases—despite their different retrieved parameters—are almost indistinguishable. Larger differences between the two cases can be identified in the region near the 1 μm peak. Here, the retrieval including the f factor seems to fit the Sonora model spectrum better. Thus, the f parameters try to provide a more accurate fit in this region.

Zoom In Zoom Out Reset image size

After comparing the various opacity sources with those in the Sonora model (D. Saumon 2020, private communication), we found the cause of these differences in the resonance ling wings of the alkali metals K and Na. The Sonora model currently uses older versions of the profiles from N. Allard that were limited to a perturber density of 10¹⁹ cm⁻³. The alkali opacity of Helios-r2, on the other hand, is based on newer calculations (Allard et al. 2016, 2019) that are now valid for H₂ densities up to 10²¹ cm⁻³. Furthermore, the new versions of the alkali line wings fall off much faster at large distances from the line center than the older ones used in the Sonora model. Thus, even the impact of Na can still be seen at 1.1 μm in the Sonora model, whereas the new sodium far-wing line profiles used in Helios-r2 have already decayed to undetectable values.

GJ 570 D—Our spectrum of GJ 570 D has been taken by the SpeX instrument, spanning the wavelength range from 0.8 to about 2.5 μm (Burgasser et al. 2006a). The spectral resolution varies between about 85 and 300 throughout the spectrum. The spectrograph's slit width for this observation is 0''5. With an image scale of 015 per pixel, the spectral flux at a given wavelength in the spectrum is, thus, sampled onto 3.3 pixels on the CCD. This oversampling is accounted for in our retrieval by using an appropriate instrument profile (see Section 2.7 for details). In order to obtain the wavelength-dependent standard deviation for the instrument profile, we estimate the local width of a pixel in wavelengths by using the values of neighboring wavelength points in the spectrum. The result multiplied by 3.3 is then approximately the FWHM of the Gaussian instrument profile.

As mentioned by Line et al. (2015), the oversampling also results in the neighboring pixel being not statistically independent, because the flux information is partly duplicated. Therefore, we follow the approach of Line et al. (2015) and use only every third pixel in our retrieval.

We use 2MASS photometric data to flux-calibrate the spectrum of GJ 570 D. The spectrum is scaled by a multiplicative factor that is separately computed for the J (15.32 ± 0.05 mag), H (15.27 ± 0.09 mag), and K_S (15.24 ± 0.16 mag) bandpasses (see Cushing et al. 2005, for details). Uncertainties in the scale factor take into account spectral measurement errors and photometric uncertainties. We adopt the weighted mean of these three values for our final flux calibration scale factor. The calibrated spectrum of GJ 570 D is shown in Figure 12.

Epsilon Indi Ba & Bb—The two brown dwarfs in the Indi system have been previously observed by King et al. (2010) using the FORS2 instrument at the Very Large Telescope (VLT) in the optical wavelength range and the ISAAC spectrograph (VLT) in the near-infrared and infrared. Details on the calibration of the spectrum can be found in King et al. (2010). In the following, we focus on the near-infrared ISAAC measurement that covers about the same wavelength range as the SpeX measurement of GJ 570 D. With about 20,000 wavelength points, the spectral resolution of the ISAAC measurement is much higher than the one provided by the SpeX prism.

Using the full resolution of 20,000 wavelength points, however, is problematic in the framework of the nested sampling Bayesian retrieval used in this study. The definition of the likelihood function in Equation (19) is essentially the χ² distance between the measured flux values and those predicted from the forward model. At high dimensions, this distance is known to lose its mathematical meaning (Beyer et al. 1999), which is commonly referred to as the curse of dimensionality. It should be noted that due to its lower spectral resolution, spectra obtained with the JWST will not have this issue.

To avoid this issue, we integrate the 20,000 wavelength bins to a lower resolution. In total, we use about 400 bins, which is equivalent to the full resolution of the SpeX instrument. Due to the binning to a resolution that is much lower than the original one, no instrument profile has to be accounted for, and this effect is, therefore, neglected for Ind Ba and Bb. The spectra of Indi Ba & Bb are shown in Figure 16.

For each of the three brown dwarfs, we split the retrieval calculations into two different categories, a first one assuming equilibrium chemistry and a second using a free chemistry approach. A summary of all retrieval parameters is given in Table 1. For the equilibrium chemistry model, we use the overall metallicity M/H (assuming solar element abundance ratios) and the C/O ratio as free parameters. In the free chemistry approach, we retrieve for the mixing ratios x_i of H₂O, NH₃, CH₄, CO, CO₂, H₂S, and K. We note that in contrast to Line et al. (2015), we do not use the mixing ratio of sodium as a free parameter. The Na abundance is obtained from the potassium mixing ratio by using their solar elemental abundance ratio. As mentioned in Section 4.2, by using the new far-wing resonance line-wing profiles, our calculated spectra are insensitive to sodium beyond a wavelength of 1 μm.

Table 1.  Summary of Retrieval Parameters and Prior Distributions Used for the Free and Equilibrium Chemistry Models

Parameter	Prior
	Type	Values
log g	uniform	3.5–6.0
d	Gaussian	measureda
f	uniform	0.1–5.0
T1	uniform	1000–3000
bi	uniform	0.3–0.95
	uniform	min(σj)–max(σj)
Equilibrium chemistry
M/H	uniform	0.1–5.0
C/O	uniform	0.1–4.0
Free chemistry
xi	log-uniform	10−12–10−1

Note.

^aWe use distances inferred from Gaia parallax measurements (Gaia Collaboration 2018). For the GJ 570 system, the measured distance is 5.8819 ± 0.0029 pc, for Ind 3.6389 ± 0.0033 pc.

Download table as:  ASCIITypeset image

The species' mixing ratios as a function of pressure are assumed to be isoprofiles, i.e., the retrieved abundances should be interpreted as the mean abundances within the visible atmosphere. The abundances of H₂ and He are obtained by assuming that they make up the rest of the atmosphere, using the solar ratio of their elemental abundances.

The metallicity and C/O ratio are obtained in a post-process procedure. The C/O ratio is computed by dividing the sum of the carbon-bearing species by the sum of oxygen-bearing species, each weighted by the number of carbon or oxygen atoms present in the molecule. The metallicity [M/H] is approximated by summing up the constant mixing ratios for each species weighted by the number of metal atoms and divided by the abundance of hydrogen. The result is then compared to the sum of solar metals relative to hydrogen.

In addition, we introduce the distance d as a retrieval parameter. We use the measured distances and the corresponding errors with a Gaussian prior (see Table 1). This procedure will propagate the error in the measured distances through all other retrieval parameters.

The temperature profile is assumed to be characterized by either six first-order or three second-order elements. In both cases, the profile is described by seven free parameters. Thus, for a free chemistry retrieval, we have in total 18 free parameters, while the equilibrium chemistry approach requires 13.

A summary of the retrieval results for all three brown dwarfs is given in Table 2. Additionally, the table also lists the corresponding parameters obtained by other studies for comparison. The posterior distributions are shown and discussed within the next subsections.

Table 2.  Summary of Retrieved Parameters for the Three Brown Dwarfs and Comparison with Previously Reported Values

	Parameter	This Work		Previous Work
		Equilibrium Chemistry	Free Chemistry
GJ 570 D	Teff (K)	${730}_{-17}^{+18}$	${703}_{-20}^{+17}$	${714}_{-23}^{+20}$ (Line et al. 2015)
				759 ± 63 (Filippazzo et al. 2015)
				780–820 (Burgasser et al. 2006a)
				800–820 (Saumon et al. 2006)
				900 (Testi 2009)
				948 ± 53 (Del Burgo et al. 2009)
	log g	${4.61}_{-0.08}^{+0.08}$	${5.01}_{-0.19}^{+0.13}$	${4.5}_{-0.5}^{+0.5}$ (Del Burgo et al. 2009)
				${4.76}_{-0.28}^{+0.27}$ (Line et al. 2015)
				4.90 ± 0.5 (Filippazzo et al. 2015)
				5.0 (Testi 2009)
				5.09–5.23 (Saumon et al. 2006)
				5.1 (Burgasser et al. 2006a)
	R (RJ)	${1.00}_{-0.09}^{+0.10}$	${1.13}_{-0.06}^{+0.05}$	${1.14}_{-0.09}^{+0.10}$ (Line et al. 2015)
	C/O	${0.83}_{-0.08}^{+0.09}$	${1.11}_{-0.09}^{+0.09}$	0.95–1.25a, 0.70b (Line et al. 2015)
	[M/H]	$-{0.15}_{-0.04}^{+0.05}$	$-{0.13}_{-0.08}^{+0.06}$	−0.29–0.04a, −0.15b (Line et al. 2015)
	$\mathrm{ln}{ \mathcal Z }$	4775.73	4775.06
${\boldsymbol{\epsilon }}$ Indi Ba	Teff (K)	${1339}_{-19}^{+19}$	${1420}_{-16}^{+16}$	1250–1300 (Kasper et al. 2009)
				1300–1400 (King et al. 2010)
	log g	${5.49}_{-0.10}^{+0.06}$	${5.62}_{-0.07}^{+0.07}$	5.2–5.3 (Kasper et al. 2009)
				5.27 ± 0.09 (Dieterich et al. 2018)c
				5.50 (King et al. 2010)
	R (RJ)	${0.73}_{-0.02}^{+0.02}$	${0.55}_{-0.01}^{+0.01}$
	C/O	${0.44}_{-0.03}^{+0.04}$	${0.95}_{-0.03}^{+0.02}$
	[M/H]	$-{0.70}_{-0.07}^{+0.06}$	${0.89}_{-0.23}^{+0.17}$	−0.2 (King et al. 2010)
	M (${M}_{{\rm{J}}})$	${67}_{-12}^{+8.4}$	${50}_{-6.8}^{+7.8}$	$75\pm 0.82$ (Dieterich et al. 2018)
	$\mathrm{ln}{ \mathcal Z }$	4319.21	4352.01
${\boldsymbol{\epsilon }}$ Indi Bb	Teff (K)	${768}_{-25}^{+26}$	${992}_{-21}^{+22}$	875–925 (Kasper et al. 2009)
				880–940 (King et al. 2010)
	log g	${5.11}_{-0.05}^{+0.05}$	${4.85}_{-0.19}^{+0.17}$	4.9–5.1 (Kasper et al. 2009)
				5.24 ± 0.09 (Dieterich et al. 2018)c
				5.25 (King et al. 2010)
	R (RJ)	${0.73}_{-0.02}^{+0.02}$	${0.71}_{-0.03}^{+0.04}$
	C/O	${1.21}_{-0.08}^{+0.09}$	${0.84}_{-0.07}^{+0.07}$
	[M/H]	$-{0.30}_{-0.06}^{+0.06}$	$-{0.34}_{-0.11}^{+0.12}$	−0.2 (King et al. 2010)
	M (${M}_{{\rm{J}}})$	${77}_{-4.2}^{+2.2}$	${15}_{-4.6}^{+6.0}$	70.1 ± 0.68 (Dieterich et al. 2018)
	$\mathrm{ln}{ \mathcal Z }$	4366.03	4417.71

Notes.

^aRetrieved parameter. ^bDerived from post-process chemistry model. ^cDerived parameter, based on the measured dynamical mass and assuming $R=1\pm 0.1{R}_{{\rm{J}}}$.

Download table as:  ASCIITypeset image

The results for Helios-r2 with equilibrium chemistry are shown in Figure 9, while the ones with the free chemistry retrieval are depicted in Figures 10 and 11. The corresponding spectra for both retrievals are presented in Figure 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Overall, most parameters seem to be quite well constrained in both cases. For the equilibrium chemistry forward model, we obtain a sub-solar metallicity of $-{0.15}_{-0.04}^{+0.05}$ and a super-solar C/O ratio of 0.83. The retrieved metallicity and C/O compares well with that of the host star reported by Line et al. (2015): −0.22–0.12 and 0.65–0.97, respectively. It should, however, be noted that our reported metallicity and C/O ratio corresponds to the pure gas phase and neglect the losses of elements due to condensation. The atmosphere's intrinsic, bulk element abundances, thus, might differ from our reported values.

Our retrieved value for the surface gravity of ${4.61}_{-0.08}^{+0.08}$ also matches the one reported by Line et al. (2015) within their confidence intervals. When converting the retrieved f parameter into a stellar radius via Equation (10), we obtain a value of about 1 Jupiter radius. Combined with our log g value, we estimate a substellar mass of about ${17}_{-3.0}^{+3.8}$ M_J, which is smaller than the one reported by Line et al. (2015) (${31}_{-16}^{+24}$ M_J) but still contained within their confidence interval.

The free chemistry retrieval, on the other hand, yields slightly different results. With ${5.01}_{-0.19}^{+0.13}$, the retrieved log g is higher than for the equilibrium chemistry case, while the derived C/O ratio is about 1.11 ± 0.09, indicating an atmosphere that is enriched in carbon or, via condensation of oxygen-bearing species, depleted in O. The directly retrieved C/O ratio found by Line et al. (2015) is 0.95–1.25 and, thus, well within our 1σ confidence interval.

As suggested by the posterior distributions in Figure 10, we can constrain the abundances of water, methane, ammonia, and potassium. Upper bounds of roughly 10⁻⁵ are obtained for CO₂ and H₂S, while an upper limit of 10⁻² is found for CO. When comparing the retrieved molecules' abundances with those from the free chemistry retrieval of Line et al. (2015), our model yields very similar median values.

Even though the value of the surface gravity is now higher than the one of the Line et al. (2015) study, it is still consistent with those from other publications. As Table 2 indicates, the spread in reported surface gravities for GJ 570 D extends from 4.5 to 5.23. Additionally, the f factor, and thus the inferred radius, is larger than for the equilibrium chemistry case, which also results in a derived mass of about 53 Jupiter masses—more than a factor of three higher than in the previous case.

Both retrievals also have similar posterior distributions for the effective temperature, both of which result in T_eff values of slightly larger than 700 K. This, again, is consistent with previous estimates by Line et al. (2015) and Filippazzo et al. (2015), even though higher temperatures have also been obtained by, e.g., Del Burgo et al. (2009) and Testi (2009; see Table 2).

The temperature profiles obtained in both cases are quite similar, follow the form expected from the theory of brown dwarf atmospheres, and compare well to the one retrieved by Line et al. (2015). The smallest confidence intervals are found around pressures of 1 bar where most of the spectrum originates from. Larger intervals are obtained at lower and higher pressures. As mentioned in Section 2.5.1, our profiles are continuous and smooth by construction without any additional smoothing parameters required by other representations of the retrieved temperature profile.

Even though the results of the two different retrievals differ in terms of, e.g., surface gravity or metallicity, the spectra generated from the posterior distributions are surprisingly similar. As Figure 12 suggests, the median spectra only differ in details. For a larger part of the wavelength range, they are almost indistinguishable. Overall, the fit of the theoretical spectra to the actual observed one of GJ 570 D is quite good. The largest differences between the two are found at the 1.05 μm peak, where both of our retrievals have smaller values than the measured brown dwarf spectrum. This effect is also noticeable in the corresponding spectra shown in Line et al. (2015). It is possible that the description of the potassium far line wings from Allard et al. (2016), which have a large impact in this region, is still not satisfactorily representing the actual line shapes encountered in the atmospheres of brown dwarfs.

The resulting Bayesian evidence $\mathrm{ln}{{ \mathcal Z }}_{\mathrm{ec}}$ of the equilibrium chemistry retrieval is 4775.73, while the free chemistry forward model yields a value of $\mathrm{ln}{{ \mathcal Z }}_{\mathrm{fc}}=4775.06$. The corresponding Bayes factor $B={Z}_{\mathrm{ec}}/{Z}_{\mathrm{fc}}$ is 1.95. On the Jeffreys scale (Kass & Raftery 1995), this indicates that there is no evidence to favor either one of the two different chemistry approaches. Both are equally likely to explain the data.

In the following, we present our retrieval analysis of the two brown dwarfs in the Indi system, using the observational data from King et al. (2010). In contrast to our retrieval of GJ 570 D from the previous subsection, we here have to impose an upper limit on the derived substellar masses of the brown dwarfs to obtain more realistic values for the retrieval parameters. Based on estimates of the hydrogen-burning main-sequence edge (Burrows et al. 2001), we employ an upper mass limit of 80 M_J. The same approach has also been used by, for example, Line et al. (2015).

5.4.1. Equilibrium Chemistry Retrieval

The resulting posterior distributions for the equilibrium chemistry forward model are presented in Figure 13, while those for free chemistry retrieval are shown in Figures 14 and 15. The spectra and comparison to observations are given in Figure 16. A summary of the results and a comparison to the corresponding values from other studies is again presented in Table 2. We note that the potential impact of clouds is not accounted for in the retrieval calculations presented in this section. A corresponding retrieval for Eps Ind Ba with a gray cloud layer is shown in Appendix B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the case of equilibrium chemistry, we obtain quite different results for the overall metallicity of the two brown dwarfs. For Indi Ba, we retrieve a value of −0.70 ± 0.07, while Indi Bb yields −0.30 ± 0.06, respectively. With a C/O ratio of 0.44 ± 0.04 for Indi Ba and 0.79 ± 0.09 for Indi Bb, both are predicted to be enriched in oxygen compared to the solar value. However, just like the metallicity, the C/O ratios also differ by almost a factor of two.

Another striking difference is the retrieved calibration parameter f. For the T1 dwarf, we obtain a value that is much smaller than unity, while for Indi Bb, the retrieved parameter is about 1.47 ± 0.21. Since both brown dwarfs are part of the same system and, therefore, have the same distance to the observer, the distinctively different f factor cannot originate from an erroneous distance estimate. The very different predicted values for f also result in inferred radii that differ by quite a wide margin. The radius of Indi Ba is estimated to be 0.73 Jupiter radii, which is smaller than one would expect from a brown dwarf, while our results for Indi Bb infer a radius of 1.21 R_J. As mentioned earlier, the very low, retrieved radius for the early T dwarf could also be the result of a heterogeneous atmosphere that has a smaller effective emitting area. Alternatively, f might also again compensate for missing or oversimplified model physics.

Our retrieved surface gravities for both brown dwarfs are also quite high. In particular, the value of ${5.49}_{-0.10}^{+0.06}$ for Indi Ba is higher than those reported by most previous studies (see Table 2). For Indi Bb, we obtain a value of 5.11 ± 0.06, which, on the other hand, agrees very well with previous estimates. It should, however, be noted that these values are influenced by the restriction of our retrievals to a total derived mass of 80 Jupiter masses, which is thought to be the upper mass limit for brown dwarfs before the star becomes heavy enough to ignite the hydrogen burning in its core. As can be clearly noted in the correlation plots of Figure 13, the posteriors for f and log g are cut at higher values. A fully free retrieval would probably have resulted in even higher values of the surface gravity. Dynamical masses for the two companions have been reported by Dieterich et al. (2018). With 75 M_J ( Indi Ba) and 70.1 M_J ( Indi Bb), respectively, these masses also quite high. At least for Indi Ba, the measured mass is contained within the retrieved 1σ confidence interval.

The same also applies to the inferred masses (see bottom panel of Figure 13). The posterior for Indi Bb is clearly cut at the upper mass limit, while the one for Indi Ba is also skewed toward higher values of M.

Our derived equilibrium temperature of about 1339 K in the case of the T1.5 dwarf Indi Ba falls within the predicted range of 1250–1400 K published in earlier studies. For Indi Bb, a T6 brown dwarf, we obtain a value of 768 K that is cooler than the lower bound (875 K) estimated by Kasper et al. (2009).

One striking difference between the two brown dwarfs is the retrieved temperature profiles. The T1 dwarf Indi Ba shows a peculiar, shallow lapse rate in the lower atmosphere, which is absent in both later T dwarfs, GJ 570 D and Indi Bb, of our study. Such shallow profiles are normally not expected from the standard brown dwarf atmosphere models. They usually predict much steeper lapse rates in the lower atmosphere that are either given by the dry/moist adiabates or a temperature profile in radiative equilibrium (Marley & Robinson 2015). This behavior might be caused by the lack of an isothermal cloud layer in the lower atmosphere in the current retrieval model.

5.4.2. Free Chemistry Retrieval

5.4.3. Comparison of Spectra