Precise Dynamical Masses of ε Indi Ba and Bb: Evidence of Slowed Cooling at the L/T Transition

We measure the relative positions of ε Indi Ba and Bb using nine years of monitoring by the Nasmyth Adaptive Optics System (NAOS) + Near-Infrared Imager and Spectrograph (CONICA), NACO for short (Lenzen et al. 2003; Rousset et al. 2003). We use images taken by the S13 Camera on NACO in the J, H, and K_s passbands. Our images come from Program IDs 072.C-0689(F), 073.C-0582(A), 074.C-0088(A), 075.C-0376(A), 076.C-0472(A), 077.C-0798(A), 078.C-0308(A), 079.C-0461(A), 380.C-0449(A), 381.C-0417(A), 382.C-0483(A), 383.C-0895(A), 384.C-0657(A), 385.C-0994(A), 386.C-0376(A), 087.C-0532(A), 088.C-0525(A), 089.C-0807(A), and 091.C-0899(A), all PI McCaughrean, and 381.C-0860(A), PI Kasper.

The S13 camera on NACO has a field of view (FOV) of 14'' × 14'' and a plate scale of ≈13.2 mas pix⁻¹. Most observing sequences consisted of about five dithered images in each filter. The binary system HD 208371/2 was usually observed on the same nights and in the same mode to serve as an astrometric calibrator. We use a total of 939 images of ε Indi Ba and Bb, taken over 56 nights of observations from 2004 to 2013 for which we have contemporaneous imaging of HD 208371/2.

We perform basic calibrations on all of these images. For each night, we use contemporaneous dark images to identify bad pixels and to remove static backgrounds. We construct and use a single, master flat field for all images. We mask pixels for which the flat-field correction deviates by more than 20% from its median or for which the standard deviation of the dark frames is more than five times its median standard deviation. We then subtract the median dark image and divide by the flat-field image.

The data quality varies depending on the observing conditions and the performance of the adaptive optics (AO) system. Therefore, we apply a selection criterion to exclude poor quality data. We first extract the sources in the images using the DAOPHOT program as implemented in the photutils python package (Stetson 1987; Bradley et al. 2021). We obtain estimates of the following parameters for ε Indi Ba and Bb: centroid, sharpness (a DAOPHOT parameter that characterizes the width of the source), roundness (a DAOPHOT parameter that characterizes the symmetry of the source), and signal to noise ratio (SNR). We discard images where one or both of the two targets fall outside the field of view and, for the remaining images, we apply the following cutoffs in DAOPHOT detection parameters to exclude highly extended, highly elongated and low signal to noise images: sharpness ≥ 0.3, −0.5 ≤ roundness ≤ 0.5 and SNR ≥ 25. We then visually inspect the remaining images to remove ones with bad pixels (cosmic rays or optical defects) landing on or near the target objects, and ones with AO correction artifacts that survived our DAOPHOT cut. Table 1 summarizes the final data set selected for relative astrometry measurements.

The long-term absolute position of ε Indi B was monitored with the FOcal Reducer and low dispersion Spectrograph (FORS, Appenzeller et al. 1998) installed on ESO's UT1 telescope at the Very Large Telescope (VLT). The FORS system consists of twin imagers and spectrographs FORS1 and FORS2, collectively covering the visual and near UV wavelength. The absolute astrometry monitoring was done with the FORS2 imager coupled with two mosaic MIT CCDs; the camera has a pixel scale of 0126/pixel in its unbinned mode and a field of view (FOV) of ≈86 × 86.

The FORS2 monitoring of ε Indi B covers a long temporal baseline beginning in 2005 and ending in 2016. Our images come from Program IDs 072.C-0689(D), 075.C-0376(B), 076.C-0472(B), 077.C-0798(B), 078.C-0308(B), 079.C-0461(B), 380.C-0449(B), 381.C-0417(B), 382.C-0483(B), 383.C-0895(B), 384.C-0657(B), 385.C-0994(B), 386.C-0376(B), 087.C-0532(B), 088.C-0525(B), 089.C-0807(B), and 091.C-0899(B), all PI McCaughrean. The FORS2 focal plane consists of two CCDs, chip1 and chip2. We only consider the data taken with the chip1 CCD. Over the 12 yr of absolute position monitoring, 940 images were taken with chip1 over 88 epochs. For the majority of the epochs, 10 dithered images in the I_BESSEL filter were obtained, with a 20 s exposure time for each image. We exclude 36 blank image frames over four epochs between 2009 August 21 and November 3, resulting in a final total of 904 image frames for our analysis. A summary of the FORS2 data is given in Table 2. These 904 science frames are bias-corrected and flat-fielded using the normalized master values generated from median combination of the flat and bias frames obtained in the same set of observing programs.

To measure the relative separations of the two brown dwarfs in the NACO data, we need to fit their PSFs. Cardoso (2012) has demonstrated that Moffat is the best analytical profile for the NACO data compared to Lorentzian and Gaussian. During the epochs when the projected separations of the two brown dwarfs are small, the two PSFs are only separated by one or two full widths at half maximum (FWHM). As a result, the flux near the center of one source has non-negligible contributions from the wings of the other source. This could introduce significant biases in the measured positions if fitting a PSF profile to each source separately. Therefore, we implement a joint fit of the two PSFs using a sum of two elliptical Moffat profiles:

$$\begin{eqnarray}&&\mathrm{Counts}(x,y)={f}_{1}{\psi }_{1}(x,y)+{f}_{2}{\psi }_{2}(x,y)\end{eqnarray} \tag{ 1 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{\psi }_{i}(x,y) & = & {f}_{i}\left(1+{{\boldsymbol{c}}}_{{\bf{1}}}{\left(x-{x}_{i}\right)}^{2}+2{{\boldsymbol{c}}}_{{\bf{2}}}(x-{x}_{i})(y-{y}_{i})\right.\\ & & +\,{\left.{{\boldsymbol{c}}}_{{\bf{3}}}{\left(y-{y}_{i}\right)}^{2}\right)}^{-\beta }\end{array}\end{eqnarray} \tag{ 2 }$$

where ψ_i is a general elliptical 2D Moffat profile centered at {x_i , y_i } with peak intensity f_i . Our model is the sum of two such profiles with different fluxes at different locations, sharing the same morphology, i.e., the same { c ₁ , c ₂ , c ₃ }. Instead of fitting for { c ₁ , c ₂ , c ₃ } directly, we fit for three equivalent parameters: {fwhm_x , fwhm_y , ϕ }, which are the FWHMs of the elliptical Moffat profile along the x and y axes, and the counter-clockwise rotation angle of the PSF, respectively. These physical parameters are related to { c ₁ , c ₂ , c ₃ , β} through the following equations:

$$\begin{eqnarray}&&{{\boldsymbol{c}}}_{{\bf{1}}}=\displaystyle \frac{{\cos }^{2}{\boldsymbol{\phi }}}{{\sigma }_{x}^{2}}+\displaystyle \frac{{\sin }^{2}{\boldsymbol{\phi }}}{{\sigma }_{y}^{2}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\boldsymbol{c}}}_{{\bf{2}}}=\displaystyle \frac{\sin 2{\boldsymbol{\phi }}}{2{\sigma }_{x}^{2}}-\displaystyle \frac{\sin 2{\boldsymbol{\phi }}}{2{\sigma }_{y}^{2}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\boldsymbol{c}}}_{{\bf{3}}}=\displaystyle \frac{{\sin }^{2}{\boldsymbol{\phi }}}{{\sigma }_{x}^{2}}+\displaystyle \frac{{\cos }^{2}{\boldsymbol{\phi }}}{{\sigma }_{y}^{2}}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\mathrm{fwhm}}_{x,y}=2{\sigma }_{x,y}\sqrt{({2}^{1/\beta }-1)}.\end{eqnarray} \tag{ 6 }$$

For each background subtracted image, we fit for the sum of two PSFs by minimizing χ² over 10 parameters: {x₁, y₁, x₂, y₂, f₁, f₂, fwhm_x , fwhm_y , ϕ , β}. In this case, χ² is defined by:

$$\begin{eqnarray}&&{\chi }^{2}=\sum _{i}^{{n}_{{pix}}}\displaystyle \frac{{\left({D}_{i}-{f}_{1}\,{\psi }_{1,i}-{f}_{2}\,{\psi }_{2,i}\right)}^{2}}{{\sigma }_{i}^{2}}.\end{eqnarray} \tag{ 7 }$$

We use scipy's nonlinear optimization routines (Virtanen et al. 2020) to minimize χ² over the eight nonlinear parameters {x₁, y₁, x₂, y₂, fwhm_x , fwhm_y , ϕ , β} and, for each trial set of nonlinear parameters, we solve for the best-fit linear parameters {f₁, f₂} analytically and marginalize over them.

In order to measure precise relative astrometry, we must measure and correct various instrumental properties and atmospheric effects that can alter the apparent separation and position angle (PA) of ε Indi Ba and Bb. In this section, we describe our calibrations for the instrument plate scale and orientation, distortion correction, and differential atmospheric refraction.

4.2.1. Plate Scale, Orientation, and Distortion Correction

We calibrate the plate scale and the north pointing of the NACO S13 camera using NACO's observations of a nearby wide separation binary, HD 208371/2, observed concurrently with the science data over the ∼10 yr relative orbit monitoring period. We calibrate the separation and PA of the binary in NACO data against the high precision measurements from Gaia EDR3 for HD 208371/2:

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{sep}}{\mathrm{arcsec}}=8.90612+0.00011\left(\mathrm{Jyear}-2016.0\right)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{PA}}{\mathrm{degree}}=348.10345-0.00040\left(\mathrm{Jyear}-2016.0\right).\end{eqnarray} \tag{ 9 }$$

The uncertainties on these do depend on the epoch, but with proper-motion uncertainties ≲40 μas yr⁻¹, positional uncertainties are only ≈0.5 mas even extrapolated 10 years before Gaia. This represents a fractional uncertainty in separation below 10⁻⁴ and contributes negligibly to our error budget.

To measure the separation and PA of the calibration binary, we use the Moffat PSF fitting algorithm described in Section 4.1. Since the binary is widely separated, a joint PSF fit in this case is effectively equivalent to fitting a single 2D Moffat profile for each star separately (albeit with the same structure parameters for each star's Moffat function). The calibration results are shown in Figure 2. We measure an overall average plate scale of 13.260 ± 0.001, but we also note that the plate scale seems to increase slightly with time from 2004 to 2010. Both the plate scale and the increasing trend agree with other measurements in the literature (Chauvin et al. 2010; Cardoso 2012). In Cardoso (2012), the same calibration binary was used to derive the plate scales but a different reference measurement for the binary was used. Adjusting their results to the more precise Gaia measurement of the binary brings their plate scale into agreement with ours. The PA zero-point of the instrument varies from observation to observation, and has a long-term trend as well. This is in agreement with the analysis in Plewa (2018).

Zoom In Zoom Out Reset image size

The distortion correction was shown to be of little significance for the NACO S13 camera (Trippe et al. 2008) due to the small field of view. For completeness, we still apply the distortion correction derived by Plewa et al. (2018).

4.2.2. Differential Atmospheric Refraction and Annual Aberration

The dominant atmospheric effect that needs to be corrected for is differential atmospheric refraction (Gubler & Tytler 1998). When a light ray travels from vacuum into Earth's atmosphere, it is refracted along the zenith direction and changes the observed zenith angle of the source, making the apparent zenith angle, z, deviate from the true zenith angle in the absence of an atmosphere, z₀:

$$\begin{eqnarray}&&z={z}_{0}+R\end{eqnarray} \tag{ 10 }$$

where R is the total refraction angle experienced by the light ray. The amount of this refraction depends on atmospheric conditions, the wavelength of the incoming light, and the zenith angle of the object. Therefore, for two objects at different positions in the sky and with different spectral types, the total refraction angles are different and can alter the apparent separation and PA of the objects. We can write this differential refraction, ΔR, in terms of two components, one due to their difference in color, and one due to the difference of their true zenith angles (Gubler & Tytler 1998):

$$\begin{eqnarray}&&{\rm{\Delta }}R={\rm{\Delta }}{R}_{\mathrm{color}}+{\rm{\Delta }}{R}_{{\rm{\Delta }}{z}_{0}}.\end{eqnarray} \tag{ 11 }$$

For ε Indi Ba and Bb, the second term is much smaller as they are separated by only <1'', and produced negligible effects on the final results compared to the first term. We included both effects for completeness. The total differential refraction can be calculated with:

$$\begin{eqnarray}&&{\rm{\Delta }}R={R}_{2}({n}_{2},{z}_{2})-{R}_{1}({n}_{1},{z}_{1})\end{eqnarray} \tag{ 12 }$$

where the n_i 's are the effective refractive indices of the Earth's atmosphere for the target sources. The parameter n_i depends on the effective central wavelength (λ_i ) of the target in the observed passband, and on observing conditions, most commonly pressure (P), temperature (T), humidity (H), and altitude (z). Cardoso (2012) calculated the effective central wavelengths for ε Indi Ba and Bb in the J, H, and K_s bands by integrating high-resolution spectra of the two brown dwarfs. To calculate the refractive index, n_i (λ_i , P, T, H, z), we use the models in Mathar (2007) covering a wavelength range of 1.3 μm to 24 μm. Then, the total refraction can be approximately expressed as (Smart 1977):

$$\begin{eqnarray}&&R(n,z)\approx (n-1)\tan (z).\end{eqnarray} \tag{ 13 }$$

A comparison of the separations along the zenith direction of ε Indi Ba and Bb are shown in Figure 3. We can clearly see the systematic differences between the J, H, and K_s bands due to differential refraction before the correction. After applying the correction, the three bands are brought to much better agreement as well as having a smaller total scatter around the mean.

Zoom In Zoom Out Reset image size

Annual aberration is a phenomenon that describes a change in the apparent position of a light source caused by the observer's changing reference frame due to the orbital motion of the Earth (Bradley 1727; Phipps 1989). We correct for the differential annual aberration, the difference in aberration between ε Indi Ba and Bb, in relative astrometry by transforming the measured positions of ε Indi Ba and Bb to a geocentric reference frame using astropy. The effect is generally a small fraction of the relative astrometry error bars and has negligible impact on the relative orbit fit. For absolute astrometry, the aberration is absorbed by the linear component of the distortion correction.

4.2.3. PSF Fitting Performance and Systematics

In order to understand how well our PSF fitting algorithm described in Section 4.1 performs, we investigate the systematic errors and potential biases of the algorithm in this section, and adjust the errors of our results accordingly.

To do this, we crop out boxes around the stars in the calibration binary, HD 208371/2, and use them as empirical PSFs. We build a collection of such PSF stamps from the images of the calibration binary based on AO quality and SNR. We use these PSFs stamps and empty background regions of the NACO data to generate mock data sets containing overlapping PSFs. For each such mock image, we randomly select one empirical PSF from the collection and place two copies of this PSF onto the background of an ε Indi B image. We scale the fluxes of the two PSF copies to be similar to those of ε Ba and Bb in a typical image. We then generate a large sample of these mock images at various separations and PAs. Since the calibration binary stars are widely separated, these empirical PSFs are effectively free of nearby star contamination. We then perform the PSF fitting described in Section 4.1 on the mock images and compare the measurements to the true, known separations and PAs.

The results for this test are shown in Figure 4. Each data point is the root mean square residual from fitting 400 mock images at the same separation but with various PAs. The errors we find from these mock data sets are slightly larger but within the same order of magnitude as the scatter in our ε Indi B measurements. We also find that the residuals of these mock data measurements do increase as the PSF overlap becomes significant, but they remain at the milliarcsecond level even at a separation equal to the closest separation in the ε Indi B data set. The performance for the K_s band is slightly worse due to the large flux ratio of the system in K_s . Overall, our joint PSF fitting algorithm has sub-milliarcsec errors across all three bands for widely separated sources, and has within a few milliarcseconds errors for overlapping sources. For our final relative astrometry results for ε Indi B, we add the systematic errors shown in Figure 4 to the measurement errors of the relative astrometry in quadrature.

Zoom In Zoom Out Reset image size

The final relative astrometry results are shown in Table 3. These are measured by applying the calibrations described in Section 4.2 and jointly fitting for the positions of ε Indi Ba and Bb for every selected image, using the PSF fitting method described in Section 4.1. We take the mean and the error on the mean for every epoch, and add the systematic errors quadratically to the measurement errors as described in Section 4.2.3. In Figure 5, we show a PSF fit for the simple case where the two PSFs are effectively isolated, as well as a PSF fit from the epoch with the closest projected separation and hence maximum PSF overlap. For each case, we take the fit with the median squared residual from that epoch to demonstrate the typical residual level.

Zoom In Zoom Out Reset image size

We now seek to measure the position of ε Indi Ba relative to a set of reference stars in the FORS2 images with known absolute astrometry. We approach the problem in stages. First, we fit for the pixel positions of all stars in the frame. We then use Gaia EDR3 astrometry of a subsample of these stars to construct a distortion map. Next, we use our fit to the NACO data (Section 6) to fix the relative positions of ε Indi Ba and Bb. Finally, we use the PSFs of nearby reference stars to model the combined PSF of ε Indi Ba and Bb and measure their position in a frame anchored by Gaia EDR3.

We begin by measuring stellar positions in pixel coordinates and using them to derive a conversion between pixel coordinates (x, y) and sky coordinates (α, δ), i.e., a distortion correction. We identify 46 Gaia sources in the field of view of the FORS2 images; these will serve as reference stars to calibrate and derive the distortion corrections.

We fit elliptical Moffat profiles to retrieve each individual reference star's pixel location (x, y) on the detector. These are Gaia sources with known α and δ measurements propagated backwards from Gaia EDR3's single star astrometry in epoch 2016.0. We adopt the same module used for relative astrometry described in Section 4.1 to fit for the reference stars' positions. For each star, we fit for three additional parameters: the FWHMs along two directions and a rotation angle in between.

We assume a polynomial distortion solution of order N for FORS2:

$$\begin{eqnarray}&&{\alpha }_{\mathrm{model}}^{* }\equiv \alpha \cos \delta =\sum _{i=0}^{N}\sum _{j=0}^{N-i}{a}_{{ij}}{x}^{i}{y}^{j}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\delta }_{\mathrm{model}}=\displaystyle \sum _{i=0}^{N}\displaystyle \sum _{j=0}^{N-i}{b}_{{ij}}{x}^{i}{y}^{j}\end{eqnarray} \tag{ 15 }$$

minimizing

$$\begin{eqnarray}&&{\chi }^{2}=\sum _{k=1}^{{n}_{\mathrm{ref}}}\left[{\left(\displaystyle \frac{{\alpha }_{k}^{* }-{\alpha }_{\mathrm{model},k}^{* }}{{\sigma }_{\alpha * ,k}}\right)}^{2}+{\left(\displaystyle \frac{{\delta }_{k}-{\delta }_{\mathrm{model},k}}{{\sigma }_{{\delta }_{k}}}\right)}^{2}\right].\end{eqnarray} \tag{ 16 }$$

This defines a linear least-squares problem because each of the a_ij appears linearly in the data model. To avoid numerical problems, we define x = y = 0 at the center of the image and subtract α_ref = 181327, δ_ref = −56789 from all Gaia coordinates. To determine the best model for distortion correction, we compare second-, third-, and fourth-order polynomial models. We derive distortion corrections excluding one Gaia reference star at a time. We then measure the excluded star's positions, and use the distortion correction built without using this star to derive its absolute astrometry. The consistency of the best-fit astrometric parameters with the Gaia measurements, and the scatter of the individual astrometric measurements about this best-fit sky path, both act as a cross-validation test of the distortion correction. For most stars, a second-order correction outperforms a third-order correction on both metrics. This also holds true dramatically for ε Indi B itself, with a second-order distortion correction providing substantially smaller scatter about the best-fit sky path.

Once we have a list of pixel coordinates (x, y) and sky coordinates (α_*, δ) for all of our reference stars, we derive second-order distortion corrections for each image. To avoid having poorly fit stars drive the results, we clip reference stars that are ≥10σ outliers. Figure 6 shows an example of an image frame indicating the displacement of the distortion-corrected centroids according to Gaia with respect to their original "uncorrected" centroid locations on the detector. The empty red circles are Gaia stars that were discarded as outliers.

Zoom In Zoom Out Reset image size

We now seek to measure the position of ε Indi Ba on the distortion-corrected frame defined by the astrometric reference stars. We cannot fit the brown dwarfs' positions in the same way as the reference stars: their light is blended in most images. Instead, we first fix their relative position using an orbital fit to the relative astrometry (Sections 4.2 and 6). We then model the two-dimensional image around ε Indi Ba and Bb as a linear combination of the interpolated PSFs of the five nearest field stars.

With the relative astrometry fixed, our fit to the image around ε Indi Ba and Bb has nine free parameters: five for the normalization of each reference PSF, one for the background intensity, one for the flux ratio between Ba and Bb, and two for the position of ε Indi Ba. The fit is linear in the first six of these parameters. We solve this linear system for each set of positions and flux ratios, and use nonlinear optimization to find their best-fit values in each image. We then fix the flux ratio to its median best-fit value of 0.195 and perform the fits again, optimizing the position of ε Indi Ba in each FORS2 image.

Figure 7 shows two examples of the residuals to this fit. The residual intensity exhibits little structure, whether the two components are strongly blended (bottom panel) or clearly resolved (top panel).

Zoom In Zoom Out Reset image size

Our fit produces pixel coordinates of ε Indi Ba in each frame. Our use of self-calibration ensures that these pixel coordinates are in the same reference system as the astrometric standard stars. We then apply the distortion correction derived from these reference stars to convert from pixel coordinates to absolute positions in R.A. and decl.. Another important calibration for absolute astrometry is the correction for atmospheric dispersion. However, our data were taken with an atmospheric dispersion corrector (ADC) in place, which has not been sufficiently well-characterized to model and remove residual dispersion (Avila et al. 1997). We therefore use only the azimuthal projection of the absolute astrometry in the orbital fit. The effects and implications of the ADC and residual atmospheric dispersion are discussed in Section 6.2.2.

We use the relative astrometry measurements summarized in Table 3 to fit for a relative orbit and obtain the orbital parameters. For this, we use an adaptation of the open-source orbital-fitting python package, orvara (Brandt et al. 2021b), and fit for seven orbital parameters: period, the angular extent of the semimajor axis ( a _ang, hereby referred to as the angular semimajor axis), eccentricity (e), argument of periastron (ω), time of periastron (T₀), longitude of ascending node (Ω), and inclination (i). The corner plot for the MCMC chain is shown in Figure 9. The best-fit relative orbit is shown in Figure 10, and the best-fit orbital parameters are summarized in Table 4. The reduced χ² is 0.77 which suggests that we may be slightly overestimating the errors, especially for the earlier epochs. This is possibly because the earlier epochs in general have higher quality data, while we used empirical PSFs of a wider range of qualities in order to generate a large enough sample for the error inflation estimate described in Section 4.2.3. Nevertheless, we are able to produce an excellent fit and obtain very tight constraints on the orbital parameters thanks to high quality direct imaging data and a long monitoring baseline that almost covers an entire period.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Orbital Fit of the ε Indi B System

Fitted Parameters	Posterior mean ± 1σ
Period (yr)	11.0197 ± 0.0076
Angular semimajor axis (mas)	661.58 ± 0.37
Eccentricity	0.54042 ± 0.00063
ω (deg)	328.27 ± 0.12
Ω (deg)	147.959 ± 0.023
Inclination (deg)	77.082 ± 0.032
μα* (mas yr−1)	3987.41 ± 0.12
μδ (mas yr−1)	−2505.35 ± 0.10
$\left(\tfrac{{M}_{\mathrm{Bb}}}{{M}_{\mathrm{Ba}}+{M}_{\mathrm{Bb}}}\right)$	0.4431 ± 0.0008
ϖ (mas)	274.99 ± 0.43
Reduced χ2	1.00
Derived Parameters	Posterior mean ± 1σ
a (AU)	2.4058 ± 0.0040
System mass (MJup)	120.17 ± 0.62
MassBa (MJup)	66.92 ± 0.36
MassBb (MJup)	53.25 ± 0.29

Download table as:  ASCIITypeset image

We have derived optical geometric distortion corrections for all the FORS2 images in Section 4.4. We describe here our approach to fit for astrometric models to the reference stars, field stars, and most importantly the ε Indi B system. We fit standard five-parameter astrometric models, with position, proper motion, and parallax, to the reference stars and field stars in the field of view of the FORS2 images. The results from the fits for reference stars match within 20% from the proper motions and parallaxes Gaia provided. For the binary system ε Indi B, we fit a six-parameter astrometric solution, adding an extra parameter, which is the ratio between the semimajor axes of the orbits of the two components. We also review, and ultimately project out, the effects of atmospheric dispersion. The wavelength-dependent index of refraction of air causes an apparent, air-mass-dependent displacement between the redder brown dwarfs ε Indi Ba and Bb and the bluer field stars along the zenith direction.

The results from absolute astrometry give proper motions, parallax, and a ratio between the semimajor axes, which can then be converted into a mass ratio and individual masses. In conjuction with our previous relative astrometry results, full Keplerian solutions can be derived that completely characterize the orbits of both ε Indi Ba and ε Indi Bb.

6.2.1. Astrometric Solution

The astrometric solution for a single and isolated reference star or background star is a five-parameter linear model in terms of reference pixel coordinates in R.A. and decl., proper motions in R.A. and decl., and the parallax. A star's instantaneous position (α*, δ) would be its position $({\alpha }_{\mathrm{ref}}^{* },{\delta }_{\mathrm{ref}})$ at a reference epoch, plus proper motion (μ_α*, μ_δ ) multiplied by the time since the reference epoch t_ref, and parallax ϖ times the so-called parallax factors Δπ_α* and Δπ_δ :

$$\begin{eqnarray}\left[\begin{array}{ccccc}1 & 0 & t-{t}_{\mathrm{ref}} & 0 & {\rm{\Delta }}{\pi }_{{\alpha }_{* }}\\ 0 & 1 & 0 & t-{t}_{\mathrm{ref}} & {\rm{\Delta }}{\pi }_{\delta }\end{array}\right]\left[\begin{array}{c}{\alpha }_{\mathrm{ref}}^{* }\\ {\delta }_{\mathrm{ref}}\\ {\mu }_{\alpha * }\\ {\mu }_{\delta }\\ \varpi \end{array}\right]=\left[\begin{array}{c}\alpha * \\ \delta \end{array}\right].\end{eqnarray} \tag{ 17 }$$

To test the robustness of the distortion corrections in R.A. and decl. that we have derived for each image, we "reverse engineer" by excluding a particular reference star from the fit and solve for the astrometric solution of that star based on the discussion above for comparison to the Gaia parameters. In particular, we focus on the reference stars close to ε Indi B.

For the binary system ε Indi B, the astrometric solution demands an additional parameter r_Ba: the ratio between the semimajor axis of ε Indi Ba about the barycenter to the total semimajor axis a. The parameter r_Ba is related to the binary mass ratio by

$$\begin{eqnarray}&&{r}_{\mathrm{Ba}}=\displaystyle \frac{{M}_{\mathrm{Bb}}}{{M}_{\mathrm{Ba}}+{M}_{\mathrm{Bb}}}.\end{eqnarray} \tag{ 18 }$$

The model becomes

$$\begin{eqnarray}\begin{array}{l}\left[\begin{array}{cccccc}1 & 0 & t-{t}_{\mathrm{ref}} & 0 & {\rm{\Delta }}{\pi }_{\alpha * } & {a}_{\alpha * }\\ 0 & 1 & 0 & t-{t}_{\mathrm{ref}} & {\rm{\Delta }}{\pi }_{\delta } & {a}_{\delta }\end{array}\right]\left[\begin{array}{c}{\alpha }_{\mathrm{ref}}^{* }\\ {\delta }_{\mathrm{ref}}\\ {\mu }_{\alpha * }\\ {\mu }_{\delta }\\ \varpi \\ {r}_{\mathrm{Ba}}\end{array}\right]\\ \quad =\,\left[\begin{array}{c}\alpha {}_{\mathrm{Ba}}^{* }-0.5{\dot{\mu }}_{\alpha * }{\left(t-{t}_{\mathrm{ref}}\right)}^{2}\\ {\delta }_{\mathrm{Ba}}-0.5{\dot{\mu }}_{\delta }{\left(t-{t}_{\mathrm{ref}}\right)}^{2}\end{array}\right].\end{array}\end{eqnarray} \tag{ 19 }$$

We also take into account the perspective acceleration that occurs when a star passes by the observer and its proper motion is exchanged for radial velocity. This effect is more significant for ε Indi than for remote stars. We employ constant perspective accelerations of ${\dot{\mu }}_{\alpha * }$ = 0.165 mas yr⁻² in R.A. and ${\dot{\mu }}_{\delta }$ = 0.078 mas yr⁻² in decl. for the ε Indi B system based on Gaia EDR3 measurements and assuming the radial velocity measured for ε Indi A. We adopt a reference epoch t_ref = 2010. With an astrometric baseline of ∼10 yr, this gives a displacement of $0.5\dot{\mu }{(t-{t}_{\mathrm{ref}})}^{2}\approx 2$ mas at the edges of the observing window, where $\dot{\mu }$ is the acceleration. The perspective acceleration, because it is known, is included in the right-hand side of Equation (19).

6.2.2. Residual Atmospheric Dispersion

The FORS2 imaging covers twelve years, with data taken over a wide range of air masses. This makes it essential to correct for atmospheric dispersion caused by the differential refraction of light of different colors as it passes through the atmosphere. The degree of dispersion is related to the wavelength of light, the filter used, and the air mass, but is always along the zenith direction. Many of the FORS2 images were taken at very high air mass. The typical air mass will vary over the course of the year, because the time of observation will vary depending on what part of night the target is up.

All of the FORS2 images were taken with an atmospheric dispersion corrector, or ADC. The residual dispersion depends on filter, air mass, and position on the FOV, but is typically tens of mas (Avila et al. 1997). This is smaller than the system's parallax and angular semimajor axis, but only by a factor of ≈10. Further, the ADC is only intended to provide a full correction to a zenith angle of 50° (Avila et al. 1997). At lower elevations, it is parked at its maximum extent; Cardoso (2012) applied an additional correction to these data. Because the residual dispersion is only in the zenith direction, we perform two fits to the absolute astrometry of ε Indi Ba. First, we use our measurements in R.A. and decl. directly. Second, we use the parallactic angle θ to take only the component of our measurement along the azimuth direction, which is immune to the effects of differential atmospheric refraction.

We project the data into the altitude-azimuth frame by left-multiplying both sides of Equation (19) by the rotation matrix

$$\begin{eqnarray}{\boldsymbol{R}}=\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right].\end{eqnarray} \tag{ 20 }$$

The top row of Equation (20) corresponds to the azimuth direction, while the bottom row corresponds to the altitude direction.

Fitting in both the altitude and azimuth directions produces a parallax of 263 mas, in agreement with Cardoso (2012) but much lower than both the Hipparcos and Gaia values for ε Indi A. We then perform a fit only in the azimuth direction: we multiply both sides of Equation (19) by the top row of Equation (20). We exclude eight 3σ outliers and assume a uniform per-epoch uncertainty of 8.01 mas to give a normalization factor that gives a reduced χ² value of 1.00. This procedure results in a parallax of 274.99 ± 0.43 mas, in good agreement with both the Hipparcos and Gaia measurements. We note that the 25 yr time baseline between Hipparcos and Gaia causes a small parallax difference. The star ε Indi A has a radial velocity of 40.5 km s⁻¹ (Gaia Collaboration et al. 2018), which translates to a fractional change of 3 × 10⁻⁴ or a decrease in parallax of about 0.08 mas over 25 yr. This difference is much smaller than the uncertainties of any of these parallax measurements.

Figure 11 shows the residual to the best-fit model using only azimuthal measurements: the top panel shows the residuals in altitude, while the bottom panel shows the residuals in azimuth. An upward trend and nonzero mean are seen in the altitude component of the parallax as a function of altitude, but no dependence on altitude was seen in the azimuth-based parallax. This confirms that the altitude component of the position measurements is corrupted by residual atmospheric dispersion of a magnitude consistent with expectations (Avila et al. 1997).

Zoom In Zoom Out Reset image size

The six-parameter azimuth-component-only astrometric solution gives a mass ratio of 0.4431 ± 0.0008 between the binary brown dwarf ε Indi Ba and Bb. This mass ratio is consistent with Cardoso (2012). The only differences in our approaches arise from our usage of Gaia EDR3 to anchor the distortion correction and our account of atmospheric dispersion by only taking the azimuthal projection of the motion of the system. Our parallax of 274.99 ± 0.43 mas agrees with the parallax value from Hipparcos and that of Dieterich et al. (2018).

The relative astrometry orbital fit provides a precise period and angular semimajor axis. With a parallax from absolute astrometry, we convert the angular semimajor axis to distance units: 2.4058 ± 0.0040 au. We then use Kepler's third law to calculate a total system mass of 120.17 ± 0.62 M_Jup. Finally, the mass ratio derived from absolute astrometry provides individual dynamical masses of 66.92 ± 0.36 and 53.25 ±0.29 M_Jup for ε Indi Ba and Bb, respectively. Table 4 shows the results of each component of the orbital fit, with the final individual mass measurements in the bottom panel. The uncertainty on these masses is dominated by uncertainty in the parallax: the mass ratio is constrained significantly better than the total mass. In the following section, we use both our individual mass constraints and our measurement of the mass ratio to test models of substellar evolution.

In general, tests of substellar luminosity as a function of time are either dominated by the uncertainty in the age or in the mass. In the case of the ε Indi B system, with highly precise masses having 0.5% errors, the uncertainty in the system age ($t={3.5}_{-1.0}^{+0.8}$ Gyr) is by far the dominant source of uncertainty.

Figure 12 shows the measured joint confidence intervals on luminosities and age of ε Indi Ba and Bb compared to evolutionary model predictions given their measured masses. The measured luminosity–age contours overlap all model predictions to within ≈1σ or less for both components. To quantitatively test models and observations, we compared our model-derived substellar cooling ages for ε Indi Ba and Bb to ε Indi A's age posterior, finding that they are all statistically consistent with the stellar age.

Zoom In Zoom Out Reset image size

Our results for ε Indi Ba and Bb are comparable to other relatively massive (50–75 M_Jup) brown dwarfs of intermediate age (1–5 Gyr) that also broadly agree with evolutionary model predictions of luminosity as a function of age (Brandt et al. 2021a). These include objects such as HR 7672 B (Brandt et al.2019), HD 4747 B (Crepp et al. 2018), HD 72946 B (Maire et al. 2020), and HD 33632 Ab (Currie et al. 2020).

However, despite agreeing with models in an absolute sense, it is evident in Figure 12 that the ATMO-2020 and Burrows et al. (1997) models prefer a younger age for ε Indi Ba than for ε Indi Bb. To examine the statistical significance of this difference in model-derived ages between the two components, we now consider only their measured masses and luminosities, excluding the rather uncertain stellar age.

Evolutionary models of brown dwarfs, from some of the earliest theoretical calculations up to modern work (e.g., Burrows & Liebert 1993; Phillips et al. 2020), typically predict a power-law relationship between mass and luminosity with a slope of ${\rm{\Delta }}\mathrm{log}L/{\rm{\Delta }}\mathrm{log}M=2.5$–3.0. This general agreement between models with very different assumptions, and that vary greatly in other predictions such as the mass of the hydrogen-fusion boundary, can be seen in the slopes of isochrones for 40–60 M_Jup brown dwarfs in Figure 13.

Zoom In Zoom Out Reset image size

One set of models, from Saumon & Marley (2008), that substantially alters the atmospheric boundary condition as objects cool from T_eff = 1400–1200 K predicts a much shallower slope from the M–L relation during that phase of evolution (Figure 13). These so-called hybrid models provide the best match to the M–L relation as measured in binaries composed of late-L dwarf primaries (T_eff ≈ 1400 K) and early-T dwarf secondaries (T_eff ≈ 1200 K), objects that straddle this evolutionary phase (Dupuy et al. 2015; Dupuy & Liu 2017). A fundamental prediction of these models is that during the L/T transition, objects of similar luminosity can have wider-ranging masses than in other models. The other chief prediction is that luminosity fades more slowly during the L/T transition, so that the brown dwarfs emerging from this phase are more luminous than in other models.

The star ε Indi B is the only example of a binary with precise individual masses where one component is an L/T transition brown dwarf and the other is a cooler T dwarf. This provides a unique test of the M–L relation, where the cooler brown dwarf is well past the L/T transition and the other is in the middle of it. According to hybrid models, the brown dwarf within the L/T transition will be experiencing slower cooling, so it would be more luminous than in other models. On the other hand, with the immediate removal of cloud opacity in hybrid models below 1200 K, a brown dwarf will cool even faster than predicted by other, non-hybrid models. These two effects predict that a system like ε Indi B will, in fact, have an especially steep M–L relation.

Our measured masses give a particularly steep slope for the M–L relation of ${\rm{\Delta }}\mathrm{log}L/{\rm{\Delta }}\mathrm{log}M=5.37\pm 0.08$ between the L/T-transition primary ε Indi Ba and the cooler secondary ε Indi Bb. The only evolutionary models that predict such a steep slope are the hybrid models of Saumon & Marley (2008).

To quantitatively test models, we compared the model-derived cooling ages of ε Indi Ba and Bb, given their measured masses and luminosities (Figure 14). Models like ATMO-2020 that assume a single, cloud-free atmospheric boundary condition are 3.9σ inconsistent with our measurements. At 6.9σ, models from Burrows et al. (1997) are even more inconsistent because the bunching up of isochrones around the end of the main sequence, which has a similar effect as the bunching up of isochrones due to slowed cooling in the hybrid models, occurs at higher masses than in ATMO-2020.

Zoom In Zoom Out Reset image size

The ε Indi B system therefore provides further validation of hybrid evolutionary models, where the atmosphere boundary condition is changed drastically over the narrow range of T_eff corresponding to late-L and early-T dwarfs. No longer just within the L/T transition, but affirming the consequences of slowed cooling during the L/T transition to cooler brown dwarfs (T_eff < 1000 K).

Brown dwarfs that have both directly measured masses and individually measured spectra have long been used in another type of benchmark test that tests for consistency between evolutionary models and the atmosphere models that they use as their surface boundary condition. Comparison of model atmospheres to observed spectra allows for determinations of T_eff, log(g), and metallicity. Evolutionary models predict brown dwarf radii as a function of mass, age, and metallicity. Combining these radii with empirically determined luminosities produce mostly independent estimates of ${T}_{\mathrm{eff}}\,={({L}_{\mathrm{bol}}/4\pi {R}^{2}{\sigma }_{\mathrm{SB}})}^{1/4}$, and with masses gives estimates of $\mathrm{log}({\text{}}g)=\mathrm{log}({GM}/{R}^{2})$. (Evolutionary model radii have a small dependence on the model atmospheres and thus estimates of T_eff and log(g) from their radii are not strictly, completely independent.) There are many examples of such benchmark tests, ranging from late-M dwarfs (e.g., Kenworthy et al. 2001; Zapatero Osorio et al. 2004; Dupuy et al. 2010), L dwarfs (e.g., Bouy et al. 2004; Dupuy et al. 2009; Konopacky et al. 2010), and T dwarfs (e.g., Liu et al. 2008; Dupuy & Liu 2017).

From King et al. (2010), ε Indi Ba and Bb have perhaps the most extensive and detailed spectroscopic observations (0.6–5.1 μm at up to R ∼ 5000) of any brown dwarfs with dynamical mass measurements. They found that BT-Settl atmosphere models (Allard et al. 2012) with parameters of T_eff = 1300–1340 K and log(g) = 5.50 dex best matched ε Indi Ba. For ε Indi Bb, they found T_eff = 880–940 K and $\mathrm{log}({\text{}}g)=5.25$ dex.

We computed evolutionary model-derived values for T_eff and log(g) to compare to the model atmosphere results of King et al. (2010). The most precise estimates result from using the mass and luminosity to derive a substellar cooling age and then interpolating T_eff and log(g) from the same evolutionary model grid using the measured mass and the cooling age. The SM08 hybrid models gave T_eff = 1312 ± 13 K and 972 ± 13 K for ε Indi Ba and Bb, respectively, and $\mathrm{log}({\text{}}g)=5.365\pm 0.006$ dex and 5.288 ± 0.003 dex. These evolutionary model-derived values agree remarkably well with the model atmosphere results, which were based on an atmosphere grid with discrete steps of 20 K in T_eff and 0.25 dex in log(g).

ATMO-2020 models are only strictly appropriate for ε Indi Bb, and they give T_eff = 992 ± 13 K and $\mathrm{log}({\text{}}g)\,=5.311\pm 0.003$ dex. This effective temperature is ≈4σ higher than the BT-Settl model atmosphere temperature. ATMO-2020 models are actually based on this family of model atmospheres (BT-Cond and BT-Settl should be effectively equivalent at this T_eff), so this suggests a genuine ≈50 K discrepancy between atmosphere model-derived T_eff (too low) and evolutionary model-derived T_eff (too high). If so, this could be due a to a combination of systematics in atmosphere models (e.g., non-equilibrium chemistry, inaccurate opacities) and/or ATMO-2020 evolutionary model radii (10%–20% too high).