Free Access
Issue
A&A
Volume 618, October 2018
Article Number A111
Number of page(s) 11
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201833626
Published online 18 October 2018

© ESO 2018

1. Introduction

After its discovery by Luhman (2013), the binary brown dwarf WISE J104915.57–531906.1 (Luhman 16; hereafter LUH 16) located at 2 pc from the Sun has been observed extensively. The astrometric follow-up by Boffin et al. (2014) used observations with the FORS2 camera on the Very Large Telescope (VLT) between April and June of 2013 to refine the parallax value to 2.020 ± 0.019 pc and to claim indications for the presence of a massive substellar companion around one of the binary components.

Including additional epochs obtained in 2014 with FORS2, (Sahlmann & Lazorenko 2015) derived the relative (500.23 mas) and absolute parallax ϖabs = 500.51 ± 0.11 mas of this system, obtained an upper mass limit for a potential third body of 2 M Jup, and determined a mass ratio q = 0.78 ± 0.10 for the LUH 16 binary. After publication of Gaia Data Release 1 (DR1 Gaia Collaboration 2016a, 2016b; Lindegren et al. 2016), we updated the FORS2 distortion correction, which led to the updated values of 501.139 mas and ϖabs = 501.419 ± 0.11 mas for the relative and absolute parallax, respectively (Lazorenko & Sahlmann 2017).

Using Hubble Space Telescope (HST) observations in 2014–2016, Bedin et al. (2017) derived a relative 501.118 mas and absolute parallax 501.398 ± 0.093 mas for LUH 16. Because of a longer observation time span, they were able to obtain estimates of the orbital elements of the system, in particular, 31.3 ± 7.9 years for the orbital period and 0.463 ± 0.064 for eccentricity. Bedin et al. (2017) found inconsistencies between their HST astrometry and the Sahlmann & Lazorenko (2015) astrometry at the level of 10–20 mas and therefore did not use FORS2 measurements for the orbit fitting.

Garcia et al. (2017) performed an independent analysis of the FORS2 observations obtained in 2013–2014 and added five more epochs in 2015, using Gaia DR1 for the transformation from CCD positions to the International celestial reference frame (ICRF). They added relative astrometry from the Gemini South Multiconjugate Adaptive Optics System (GeMS, Ammons et al. 2016), ESO photographic plates, CRIRES radial velocity, and archival observations from the Deep Near-Infrared Survey of the Southern Sky (DENIS). The HST astrometry was not published at that time and was not used.

Here we report new results from the combination of all available datasets, which increases the time span of high-precision astrometry from ∼2 years covered by the Bedin et al. (2017) study to 3.5 years. In comparison with our first study (Sahlmann & Lazorenko 2015), we present improved astrometric calibrations of the FORS2 data. We then update the orbital parameters and dynamical masses of the binary on the basis of a combination of FORS2, HST, GeMS, ESO archive, Gaia DR1, and Gaia DR2 (Lindegren et al. 2018; Gaia Collaboration 2018) astrometry, and the CRIRES relative radial velocity to constrain the inclination.

2. FORS2 astrometry

2.1. Observational data

We used imaging observations obtained with three large telescopes. This dataset contains 22 FORS2 frame series in 2013–2014 obtained by Boffin et al. (2014) with the high-resolution mode and a 0.126′′/px pixel scale. The dates of observations and the average full width at half-maximum (FWHM) are given in Table 1. The 2015 observations were obtained when the relative binary separation was ∼0.7′′, which is too small for very precise astrometry, and we therefore did not use these data.

We used the 36 HST measurements published in Bedin et al. (2017) that were obtained with the Wide Field Camera 3 in 2014–2016, which increases the total time span of high-precision observations to 3.5 years with approximately even sampling. This dataset is supplemented with six series of images obtained by Ammons et al. (2016) in 2014–2015 with Gemini/GeMS. Because GeMS uses adaptive optics, it produces well-separated and therefore well-measured images of the components. In addition, these infrared observations have low sensitivity to differential chromatic refraction (DCR) effects. Some GeMS images were taken close in time to FORS2 epochs, which is useful to independently control the quality of calibrations in the FORS2 astrometry. The ESO archive contains some more archived images, for example, photographic R-band images obtained with the Schmidt telescope in 1984. We used these data, reduced to ICRF by Garcia et al. (2017), to improve the quality of our reduction. In total, we used 64 frame series covering 3.5 years, and one more distant epoch of ESO-R in 1984.

Gaia DR2 includes two-parameter solutions, that is, positions only, for the two components of LUH 16. The identifiers are Gaia DR2 5353626573562355584 for LUH 16 A and Gaia DR2 5353626573555863424 for LUH 16 B. Because of the large astrometric uncertainties (2–47 mas) and the potential biases in the DR2 positions due to unaccounted orbital motion and blending, we did not use these data in our analysis.

Table 1.

Dates of FORS2 observations, average FWHM ε (in arcseconds), and DCR-related functions used in Eq. (1).

2.2. Layout of investigation

Several studies (Boffin et al. 2014; Sahlmann & Lazorenko 2015; Garcia et al. 2017) used the same FORS2 dataset (Table 1) but arrived at discrepant astrometric results. This illustrates the difficulties of ground-based sub-milliarcsecond astrometry, which are related to DCR, biases in measuring positions of close binary components, and other effects analyzed below.

Here we aim at mitigating these effects. We begin with a basic astrometric reduction applied to the FORS2 observations, which mitigates atmospheric image motion, geometric field distortion, DCR, and other effects to a level of about 0.1 mas. We started from the measured positions of stars x′, y′ obtained in our previous study (Sahlmann & Lazorenko 2015). The reduction was performed in a sequence of repeating steps, which is necessary to estimate the optimal radius of the reference field and remove systematic errors (Sahlmann et al. 2014; Lazorenko et al. 2011, 2014).

The following analysis is specific to LUH 16, and we used FORS2, HST, and GeMS observations to apply a model that includes new calibrations of the FORS2 measurements. Solutions were derived iteratively by splitting the full model into independent blocks, finding approximate solutions for model parameters within these blocks, and repeating computations until converging to the final solution. This expanded model includes

  • five astrometric parameters xc, yc, µx, ϖ for the barycenter motion of the binary;

  • two DCR parameters ρ, d, which affect the barycenter motion as seen through the atmosphere (FORS2 only).

  • seven orbital elements;

  • the mass ratio q;

  • 12 parameters , for each binary component and each annual period, to calibrate for the bias in positions caused by the seeing and flux variations (FORS2 only);

  • offset ε0 for the effective average seeing (FORS2 only);

  • two calibration parameters to remove the offsets between the equatorial positions of the binary barycenter of FORS2 and HST.

In Sect. 2.3 we describe the transformation of the measured positions x′, y′ to raw xraw, yraw positions of each LUH 16 component in the reference frame of FORS2, free of the optical distortions. The positional residuals Δ′ of the least-squares fit derived at this step provide information on the relative displacement of CCD chips (Lazorenko et al. 2014), the residual correlations in the CCD space (Sect. 2.4), and on variations with seeing (Sect. 2.5). Calibration of these effects produces “clean” astrometric positions , , that are still affected by DCR, however (Sect. 2.5). Some of the above calibrations require knowledge of the orbital motion of the LUH 16 components on the sky, however, which was derived by iteration starting from a rough initial approximation.

The FORS2 and HST/GeMS positions were obtained in different reference systems. Whereas the FORS2 reference frame is local and set by the field stars, the HST/GeMS positions are given in ICRF. Therefore, adopting ICRF as a common system for this investigation, we transformed , to positions xA, yA, xB, and yB for each A and B component in ICRF (Sect. 2.6). These positions, combined with HST and GeMS data, were used to derive the solution for the orbital and barycenter motion (Sect. 3).

2.3. Basic astrometric reductions

To reduce the FORS2 observations, we employed our usual method, which ensures an astrometric precision of about 0.1 mas for isolated sources (Sahlmann et al. 2014; Lazorenko et al. 2011, 2014). The measured photocenter positions x′ and y′ for every star in the field (reference or target) imaged at time t are represented by the model(1)

which describes the displacement of the stellar photocenter in the three-dimensional space formed by the coordinate axes x, y oriented along RA, Dec, and time t. Here x0, y0 are the approximate CCD reference star positions, which are fixed and define the directions for the coordinate axes x, y of the reference frame at the adopted reference time t0 = 56 500 MJD (=J2013.56742) close to the average epoch. Equation (1) is solved by the least-squares fit, using all available measurements of reference stars. This yields seven astrometric parameters for each star: the coordinate offsets 0 and 0, the proper motion μ̂x and μ̂y, the relative parallax ϖ̂, the DCR index ρ̂ describing atmospheric displacement of the stellar image depending on the star’s color, and the reverse displacement d, produced by the longitudinal atmospheric dispersion compensator (LADC; Avila et al. 1997). Equation (1) also contains the parallax factors px, py and the DCR related functions f1, x = tan z sin γ, f1, y = tan z cos γ, f2, x = tan zLsin γ, and f2, y = tan zLcos γ, where z is the zenith distance, zL is the LADC parameter, and γ is an angle between direction to zenith and north.

The key component of Eq. (1) is the function Φ(x′,y′), which is derived from reference stars only and models the sum of atmospheric image motion and geometric distortion variations for each individual frame. It is a full polynomial of order k/2 − 1, where the model parameter k is an even integer between 4 and 16. Here we set k = 10 as the value that yields the most stable results. Using solution for Φ(x′,y′), we determined raw positions(2)

of each A and B component in the reference frame, which are the measured positions free of geometric distortion variations. Then from Eq. (1) we derived all astrometric parameters of the target objects (coordinate offsets, proper motion, parallax, and DCR parameters).

The LUH 16 components perform a nonlinear motion in the sky that is due to the curvature of the orbit, and the positional residuals xraw, yraw display the measurement of the orbital segment ψ expressed as ψx = BX + GY, ψy = AX + FY with the Thiele-Innes parameters A, B, F, and G. The measured segment represents the residual of ψ affected by the least-squares fit (Eq. (1)). Because of the correlation between the parameters of the system Eq. (1) with ψ, their estimates are significantly biased. To emphasize this feature, we assigned a “hat” to the parameters of Eq. (1), in contrast to the “actual” corresponding parameters introduced later. The effect is illustrated by the difference of 3.53 mas in parallax values derived for components A and B with this preliminary approach, which is larger than the measurement uncertainty of 0.07 mas and is unphysical.

2.4. Correlations across CCD

The raw positions xraw, yraw were corrected for space-dependent systematic errors as described below. The residuals Δ of the least-squares fit of Eq. (1) for reference stars are nearly random values with a scatter that corresponds to a typical FORS2 astrometric precision of 0.1–0.2 mas. However, they quite often display a correlation pattern of different type across CCD space. These correlations are traced and used to calibrate the raw positions in Eq. (2) of LUH 16 for similar systematic errors. In this way, we detected and removed the relative displacement between CCD chips 1 and 2, which can affect the positions by 1–2 mas. Another type of systematic errors is seen as an oscillating autocorrelation function across CCD space. Correlations occur as a natural consequence of the polynomial fit of any measurements limited in space and therefore are not specific for FORS2. For stars without orbital motion, these errors were mitigated as described in Lazorenko et al. (2014), but for LUH 16, the effect is masked by the dominating signal ψ, which should first be subtracted from Δ′. As an initial approximation for ψ, we used polynomials of time that smoothen the individual residuals Δ′ within each annual period, allowing us to form the corrected residuals Δ′−ψ, which were then investigated and calibrated for errors correlated across the CCD. Later on, we used the actual signal ψ to refine the calibrations iteratively.

2.5. Correlation of FORS2 positions with seeing

The images of the components of LUH 16A and LUH 16B partly overlap because the separation between the sources (1.4′′ on 2013 and 1.0′′ and 2014) was comparable to the FWHM (Table 1), which was about 0.69′′ in RA and 0.76′′ in Dec on average for both periods. We therefore expect that the measurements of CCD positions can be biased depending on seeing. This effect should cause a correlation between the positional residuals and the seeing for individual images.

We found that this correlation is of low significance when determined for individual frame series, therefore we split the data into the years 2013 and 2014, during which the separation between the A and B components changed only little. To fit the measured positional residuals Δ′x − ψx and Δ′y − ψy separately for 2013 and 2014, we used the linear expressions(3)

where εx and εy stand for the FWHM along x and y axes, ε̄ is the average FWHM for the annual period, ε0 is an offset discussed below, and k0, k1, and k2 are the free model parameters. In addition to the linear dependence on ε, the model Eq. (3) includes a small but statistically significant component proportional to the change in flux I in the image relative to its annual average flux I0.

The coefficients k0, k1, and k2 determined by the least-squares fit of the residuals Δ′−ψ with the functions (3) were used to correct the measured values x′, y′, xraw, and yraw, thus obtaining the positions(4)

which correspond to the average observational conditions set by ε̄ and I0. In this way, we decreased the internal scatter of the positional residuals within each frame series and within each annual period, but did not change their average values in 2013 and 2014 significantly, thus the orbital parameters are hardly affected by this correction.

However, this does not yet solve the problem because the measurements that are not biased by seeing should correspond not to ε̄ but to ε = ε̄ + ε0 offset to a better seeing by some unknown value ε0. The key parameters derived at this phase of the reduction are the linear coefficients k1, x and k1, y (for components A and B each) used later to determine ε0 in Sect. 3 as a free parameter of the relative orbital motion.

The derived coefficients of Eq. (3) except for the offsets k0 are given in Table 2. These coefficients were initially obtained with the residuals Δ′−ψ using a polynomial approximation of ψ. Then we iteratively improved this estimate by applying ψ derived for the full solution of the orbital motion. We note that the coefficients k1 are larger for component B, as expected, because it is fainter, and they are approximately equal for both annual periods because the distance between the components did not change much. The values have opposite signs for components A and B. Considering the sign of k1, this means that for poor seeing, the measured positions of the components are offset in the direction away from the other component. The amplitude of the corrections in the positions for the maximum variations in seeing of about ±0.2′′ can reach ±1 mas, which significantly exceeds the measurement uncertainty.

Table 2.

Coefficients of Eq. (3) of the systematic bias in the photocenter position of components A and B.

2.6. Transformation into ICRF

It is convenient to perform the following analysis in the observational plane, where we have derived positions Eq. (4) related to the FORS2 reference frame. The metric of this frame, however, is deformed by the optical distortion and other effects that are significant even at the small separation between the LUH 16 components. Therefore we transformed these positions into an absolute system that is now well represented by the Gaia source catalogs of DR1 (Lindegren et al. 2016) or DR2 (Lindegren et al. 2018). As a suitable epoch for the transformation we adopted the epoch TDR1=J2015.0. The Gaia DR1 catalog was also used by Bedin et al. (2017) and Garcia et al. (2017) as reference to measure the positions of LUH 16. In this way, we obtained all positions in the ICRF at the epoch TDR1.

2.6.1. Transformation to the DR1 epoch

As explained in Lazorenko & Sahlmann (2017), the transformation was made only with stars imaged on the upper chip because the bottom CCD chip2 of FORS2 is rotated and shifted relative to the upper chip1, and the use of both chips leads to large residuals of about 50 mas. We tangent-plane projected the Gaia equatorial coordinates to the CCD plane at the point X0 = 1025.739, Y0 = 109.640 px near the position of LUH 16 in chip1 at the average time of observations, which corresponds to α0 = 162.31073, δ0 = −53.31815, and compared the projected Gaia DR1 and DR2 source positions with those in the FORS2 reference system. Positions of DR2 were computed at the epoch TDR1 using the DR2 proper motions. To bridge the epoch difference between the reference epoch t0 of FORS2 and the comparison epoch TDR1, we corrected the FORS2 positions with proper motions μ̂x, μ̂y of Eq. (1) known with about twice better internally consistent precision as compared to DR2 (but systematically offset from the absolute, see Sect. 2.6.2).

The following least-squares fit of positional residuals was made with functions for RA, and for Dec, which are a sum of two-dimensional polynomials of the maximum power n formed on coordinates X0 and Y0 of reference stars only. We tested all polynomial orders, starting from a linear model with n = 1. The optimal order n = 4 was found by using an F-test for the rejection of non-significant coefficients as described in Lazorenko & Sahlmann (2017). Thus, using 79 stars identified as common between FORS2 and DR2 on the upper chip, we obtained the average rms of the residuals of 1.10 mas, which is approximately consistent with an average uncertainty 0.67 mas for Gaia DR2 stars and 0.38 mas for FORS2 stars at this epoch. A similar transformation was also made with DR1. In this case, with 73 common stars, the average rms of the residuals DR1 – FORS2 increased to 2.03 mas, which is explained by the larger uncertainty of the DR1 positions. In the following, we proceed with the results based on Gaia DR2.

With above procedure, we derived coefficients of functions . Then, for each epoch of observations, we computed rectangular positions of LUH 16 A and B linked to the ICRF, which along RA are computed as(5)

The expressions for yA, yB along Dec are similar. The transformation uncertainty, which depends on the covariation function of Fn, is 0.21 mas for the stars in the center of the reference field, but it degrades to 0.43 mas for our targets located at the edge of this field (we recall that this reference field is not centered on the target because all stars in chip2 were rejected). When the transformation to the common system is made with DR1, these uncertainties are about twice larger. The computed positions Eq. (5) may differ, within the above uncertainty limits, by some offset Δ from the “real” positions in ICRF computed with error-free function Fn. Considering that the function Fn changes slowly across the CCD, we can assume that this offset is approximately equal for components A and B and for all their measurements at different epochs. The offset Δ, if exists, is canceled out in the distances between the components, but affects the barycenter position by Δ. Because its value can be non-negligible, we included Δ in the model Eq. (6) for the motion of LUH 16.

For similar reasons, the HST positions of Bedin et al. (2017) can be offset from the ICRF by about 0.3 mas, thus we expect to find systematic discrepancies between FORS2 and HST measurements by ∼0.5 mas.

In the left columns of Table 3 we give the positions Eq. (5), linked to the ICRF, that are affected by the DCR and reduced to the average seeing and flux with Eq. (3) and the coefficients from Table 2, assuming that ε0 = 0.

Table 3.

FORS2 Cartesian measured positions xA, yA of the A and xB, yB of the B component, corrected for variations of seeing and flux (assuming that ε0 = 0), but including the DCR displacements and the measurement uncertainties σA, σB.

2.6.2. Conversion of proper motions and parallaxes to absolute motions and parallaxes

By definition, the differential proper motions μ̂ and parallaxes ϖ̂ of FORS2 systematically differ from their absolute values μ(abs) and ϖ(abs) by some constant offsets Δμ and Δϖ, thus μ(abs) = μ̂ + Δμ and ϖ(abs) = ϖ̂ + Δϖ (Lazorenko et al. 2009, 2014). While this is not important for deriving the relative angular motion of the binary, these parameters of the barycenter derived in Sect. 3.2 need to be corrected to absolute values. In our previous studies of nearby sources (e.g., Sahlmann et al. 2014), we derived statistical corrections for the FORS2 relative parallaxes ϖ by comparing our data with a galaxy model (Robin et al. 2003). For LUH 16 we thus derived the correction term Δϖ=0.28 ± 0.01 mas (Sahlmann & Lazorenko 2015).

With the availability of Gaia DR2, it becomes possible to determine the parallax correction by directly comparing the FORS2 and DR2 parallaxes of field stars. We used only those stars that define the system of FORS2 astrometric parameters (Lazorenko et al. 2014, Sect. 3.3.1); they are typically brighter than G = 20. Because we used equal weights Pi in (Lazorenko et al. 2014, Eq. (2)) for these stars, the estimate of Δϖ is the simple arithmetic average offset. There were 65 stars in chip1 and 32 stars in chip2 that were also reduced to ICRF specifically for the purpose of deriving Δϖ and Δμ. This direct comparison (Fig. 1, upper panel) yields a parallax correction value of ϖDR2 − ϖFORS2 = Δϖ=0.342 ± 0.056 mas, which within error bars agrees with the value reported by Sahlmann & Lazorenko (2015), thereby indirectly validating the model of Robin et al. (2003) and our previous statistical method of deriving Δϖ at a level of ±0.1 mas in this field.

This assessment does not account for the DR2 parallax zero-point of −29 μas discussed in Lindegren et al. (2018), which was determined for quasars and may be variable in amplitude as a function of magnitude, position, and other parameters. Since Gaia parallaxes are therefore too small, the DR2 zero-point offset increases the value of the parallax correction we derive by the equivalent amount.

In a similar way, we compared the proper motions in RA and obtained a linear dependence μx(DR2)=Δμx + cxμx(FORS2) with an offset Δμx = −6.618 ± 0.111 mas yr−1 and a nearly unity coefficient cx = 0.953 ± 0.037 mas yr−1 (Fig. 1). The rms of the fit residuals of 1.1 mas yr−1 is small but exceeds the uncertainties of 0.18 and 0.47 mas yr−1 for proper motions of FORS2 and DR2, respectively. In Dec, using the similar relation μy(DR2)=Δμy + cyμx(FORS2), we derived Δμx = 2.103 ± 0.120 mas yr−1, cy = 0.957 ± 0.052 mas yr−1, and a residual scatter of 1.2 mas yr−1. The calibration terms Δμx, Δμy, and Δϖ determined here were used later in Sect. 3.2 to convert the relative barycenter proper motion and parallax into the absolute values, neglecting potential Gaia DR2 proper motion zero-point offsets.

thumbnail Fig. 1.

Comparison between relative FORS2 and absolute DR2 parallaxes (upper panel) and proper motions in RA (lower panel) for brighter (filled circles; they were used to determine the correction) and fainter stars (open circles), with linear fit functions (solid lines).

Open with DEXTER

3. Barycenter and orbital motion

The measured positions xA, yA xB, yB derived with Eq. (5) contain directly measured target positions, all calibrations for the effects described above, and although it is not explicitly evident, contain a full information on all measurements of the reference stars. The reference data were not discarded because they were integrated into the basic astrometric model and then were extracted again at any exposure moment and any CCD pixel in a refined form, for instance, filtered from the atmospheric image motion.

The measured positions from Eq. (5) were modeled by the sum of the barycenter motion of the system and its orbital motion in RA,(6)

and similar for the positions in Dec yA, yB. The model includes coordinate offsets xc, yc, proper motion μx(=μα cos δ), μy, relative parallax ϖ, Thiele-Innes parameters A, B, F, and G, and the mass-ratio q. The DCR parameters ρA, ρB, dA, dB, and the offset ε0 to the average FWHM from Eq. (3) are applicable to FORS2 only, which we flag by λ = 1 for FORS2, and λ = 0 in the other cases. The constants Δx, Δy are defined in the previous Sect. 2.6.1 and are the systematic offsets to the positions xA, yA, xB, yB, which depend on the instrument. These equations can therefore contain Δx(HST), Δx(FORS2), Δy(HST), or Δy(FORS2). Because the observational information for GeMS is insufficient, we assumed that Δ(GeMS)=0.

From Eq. (6) we can obtain an expression for the barycenter motion in the sky that is independent of the orbital parameters, and an expression for the orbital motion that is independent of the barycenter motion. For the orbital motion, we have(7)

and for the barycenter motion, we find(8)

where ρ = (qρB + ρA)/(1 + q) and d = (qdB + dA)/(1 + q) are the effective DCR parameters that apply to the barycenter, and are the effective values of and . Four individual parameter values ρA, ρB, dA, and dB are obtained by combining the effective parameters ρ, d and their differences ρB − ρA, dB − dA derived from Eq. (7).

3.1. Deriving orbital parameters

3.1.1. Deriving ε0

The orbital parameters were obtained by fitting FORS2, HST, GeMS, and ESO-R measurements with the model Eq. (7). Because these equations are nonlinear, we initially used approximate values of the nonlinear parameters e, T0, and P spread over a wide range with a small increment and computed the elliptical orbital coordinates X, Y of component B relative to A at the observation time. Thus, for each set of e, T0, and P, we linearized Eq. (7) in the parameters A, B, F, and G, and ρ, d, initially not including the term ε0. The least-squares solution then yielded the rms of the fit residuals, which indicates the likelihood of the tested group of e, T0, and P. As a final solution, we adopted the parameter set that yielded the best fit. The positional residuals obtained in this way (Fig. 2) show a clear systematic discrepancy between FORS2 and GeMS, and the rms of the residuals σ = 0.57 mas is large. Therefore we complemented the model Eq. (7) with the term ε0 and obtained ε0 = −0.192 ± 0.014′′. Thus we expect that at seeing ε̄ + ε0 = 0.50″, the measurements are not biased by the light from the nearby companion. The new solution yielded better σ = 0.357 mas, and the systematic offset between FORS2 and GeMS measurements obtained at comparable epochs was mitigated (Fig. 3).

3.1.2. Initial and final solutions

The barycenter and orbital motion were initially computed with weights according to the measurement uncertainties σcat associated with the data. However, the measured rms of the positional residuals exceeds the model expectation by a factor of K that is different for each data set and the barycenter or the orbital model. We therefore assumed that Kσcat better represents the real uncertainties, and accordingly updated the weights in Eqs. (7) and (8). Then all computations in Sects. 3.1 and 3.2 were repeated until convergence for K. Table 4 represents the summary for the measured rms and K values for each instrument and separately for the barycenter and orbital motion. This table shows that the residual rms is roughly equal for HST and FORS2, but the excess K is more significant in the barycenter motion for both space- and ground-based astrometry, which can be caused by some common problem, for example, related to the conversion from CCD space into ICRF. In the following, we refer to the solutions with K = 1 for all data sets and with K taken from Table 4 as the initial and final solution, respectively. In this way, the sum of squares of the normalized residuals χ2 decreased the high initial value of ∼305 to χ2 = 130.00, as expected for 2 × 65 epoch measurements.

Because the observed orbital segment is short, the range of solutions with nearly equal rms is relatively large. We find that practically the same fit quality can be obtained when the parameters e, P, and T0 change by about ±0.01, ±1 yr, and ±0.1 yr, respectively. We applied the approach of Lucy (2014), which is applicable to the fit with a model that includes the subsets of a standard nonlinear (e, P, T0) and linear (A, B, F, G) parameters1. We modified this method because for the particular case of Eq. (7), a linear group has to include ( ρB − ρA), (dB − dA), and ε0, and thus contains seven parameters. In contrast to the original algorithm of Lucy (2014), which takes into account the correlation between two pairs of parameters (B, G) and (A, F) only, we deal with the correlations between seven parameters of a linear group as defined by their covariance matrix D. This difference was removed by Jacobi transformation (Press et al. 1986) of the matrix D into its diagonalized form U = VTDV, where V is a matrix whose columns are seven eigenvectors of D, and the diagonal elements of U are the eigenvalues of D. This allowed us to closely follow the idea of the Lucy (2014) algorithm, using for transformation of the orthogonal variables z (Lucy 2014, Sect. A.4), to the linear subset of parameters (Lucy 2014, Eqs. (A.11) and (A.15)) an expression VU1/2z where U1/2 is the diagonal matrix whose diagonal elements are .

The resulting median values with 68.3% confidence intervals are given in Table 5. The differences between the initial (K = 1 for all data sets) and final (K given in Table 4) are smaller than the uncertainties of the latter.

The single-epoch ESO-R measurement from Garcia et al. (2017), which was obtained approximately one orbital period before the main body of high-precision but short-term observations, deserves particular attention. When including that data point, the observations cover the full orbit and the uncertainties in the orbital parameters and masses improve by approximately a factor of three, in spite of its large nominal σcat ≈ 150 mas (and final Kσcat ≈ 85 mas) uncertainty. This agrees with the finding of Garcia et al. (2017) that ESO-R data help to break the degeneracy between their model parameters. In our case, we have sufficient data to avoid the degeneracy even without the ESO-R data point, but its incorporation improves the precision of our solution.

In Table 5 we give the total system mass Mtot = (a/ϖabs)3/P2 and individual A and B masses computed with q = 0.8519 and absolute ϖabs derived in the next section. The uncertainty of these dynamical masses is about 1%, an approximately fourfold improvement compared to Garcia et al. (2017). Since the biases of the ESO-R measurement are difficult to quantify, we performed several tests to asses potential systematic errors in the dynamical masses: We varied the ESO-R measurement by Garcia et al. (2017) within the given uncertainty shifting it in both coordinates by uσcat, towards (u >  0) or away from (u <  0) the model predicted position and repeated our analysis, see Fig. 4. We find that the precision of the orbital and mass parameters are hardly affected by the exact value of the ESO-R measurement, and the median values can vary within σM. We therefore concluded that our formal parameter uncertainties are correct and three to four times smaller than reported by Garcia et al. (2017) if the ESO-R measurement is affected by random errors only. However, the accuracy of the individual masses depends on the systematic error/bias in the ESO-R measurement, which is unknown and may exceed the random errors (a case discussed in Sect. 2.5 in the context of FORS2). This problem will likely disappear when the next orbital inflection point will have been covered (in ∼2021), which will reveal the true bias in the ESO-R measurement.

If the ESO-R measurement is not used, we obtain slightly higher masses of Mtot = (62.6  ±  1.8) MJup, MA = (33.8  ±  0.9) MJup, and MB = (28.8  ±  0.8) MJup, which are still within the confidence intervals for the final solution in Table 5, but the uncertainties increased to 3% (Fig. 4).

A minor source of bias in the masses may stem from the seeing correction of the FORS2 measurements (Sect. 2.5), which partially relies on the external HST and GeMS data.

thumbnail Fig. 2.

Residuals in RA (lower panel) and Dec (upper panel) between the FORS2 (open circles), HST (triangles), and GeMS (filled circles) for the relative position of the secondary component of LUH 16 with ε0 = 0 in Eq. (7).

Open with DEXTER
thumbnail Fig. 3.

Same as Fig. 2, but with ε0 = −0.192′′ in Eq. (7).

Open with DEXTER
Table 4.

rms of positional residuals and the excess K in the uncertainties σcat for each data set.

Table 5.

Median orbital parameters of LUH 16AB, with 68.3% confidence intervals for initial and final solutions (with σcat and Kσcat uncertainties of the measurements, respectively), and the corresponding data derived by Garcia et al. (2017).

thumbnail Fig. 4.

Median system mass with 1σ uncertainty ranges derived when applying an artificial shift to the ESO-R measurement. The shifts were applied in both coordinates by uσcat toward (u >  0) or away from (u <  0) the model-predicted position, and the parameters KESO-R = 0.57 (open square boxes) or KESO-R = 1 (open circles) were used. The adopted value of Mtot given in Table 5 is shown as a filled square. For comparison, the 1σ uncertainty range (gray area) and the median value (solid horizontal line) computed without the ESO-R measurements are shown.

Open with DEXTER
Table 6.

Radial relative velocity VA − VB measured by CRIRES (Garcia et al. 2017) and the computed minimum/maximum values for our orbital solutions within the 68.3% confidence interval.

3.1.3. Comparison to Bedin et al. (2017) and Garcia et al. (2017)

Except for eccentricity, our parameter values agree with Bedin et al. (2017) within the error bars, but our estimates have smaller confidence intervals because we have access to a longer time span (∼3.5 instead of 2 years). Our estimates also agree with the median parameter values obtained by Garcia et al. (2017) with the Markov chain Monte Calro (MCMC) analysis. Because of a larger volume of observations and improved reduction of FORS2 data, our confidence intervals are 2–10 fold narrower for most parameters, except for the orbital period P, for which the interval 27.1–27.6 yr given by Garcia et al. (2017) is smaller than the interval 27.1–27.9 yr derived by us. In contrast to Garcia et al. (2017), we did not include radial velocity measurements in the fit, which may explain the slightly less well constrained period.

The inclination i >  90 that we obtained corresponds to the observed clockwise motion, whereas the value of i <  90 given by Bedin et al. (2017, Table 8), and Garcia et al. (2017) corresponds to a counterclockwise motion, assuming identical parameterization. In particular, we define ω as the argument of periastron for the barycentric orbit of the primary. The use of their orbital elements leads to equatorial positions with rearranged coordinate axes, that is, the computed RA is in fact Dec. From the definition of the Thiele-Innes parameters it follows that this inconsistency can be removed by changing i to 180 − i and Ω to 270 − Ω. For instance, the parameters i = 79.21 and i = 79.5, and Ω = 130.3 and Ω = 130.12 given by Bedin et al. (2017), Table 8) and Garcia et al. (2017), respectively, are converted into i = 100.79, i = 100.5, and Ω = 139.7, Ω = 139.88 which matches our estimates well. The relative motion in the binary system is shown in Fig. 5, which essentially agrees with Garcia et al. (2017) also in terms of the families of allowed orbital configurations.

We also computed the radial velocity of component LUH16 A relative to B for all possible orbital solutions within the 68.3% confidence interval and compared these estimates (Table 6) with the corresponding CRIRES measurements presented by Garcia et al. (2017). These velocity measurements were not included in the fit because they do not increase the time span of the data set and have a small weight compared to the orbital astrometry. However, we confirmed that our solution is compatible and that the inclination is 0 <  i <  180. We found good agreement between our model-predicted minimum/maximum velocity and the data (Table 6) within the uncertainty range of the measurements.

We searched for a potential seeing dependence of the residuals ΔRA in RA and ΔDec in Dec of the fit Eq. (7) that might be due to errors in the photocenter measurements caused by overlapping images. Errors of this kind are expected to occur along the line connecting the components of LUH 16. Therefore, projecting the residuals ΔRA and ΔDec onto this line, we formed the equivalent residuals Δ||, which we can correlate with ε. Because the binary system was oriented approximately at 45 relative to the coordinate axes, the projected residuals are found simply as , with the positive values corresponding to the displacement of LUH 16 B toward A. Figure 6 shows that there is no significant dependence of Δ|| on seeing.

For comparison, we present Δ|| computed based on HST, GeMS, and astrometric measurements by Garcia et al. (2017, Table 2). With these data, we computed the orbital motion with Eq. (6), neglecting the DCR terms because these measurements are free from DCR, and obtained the residuals ΔRA and ΔDec nearly equal to those shown by Garcia et al. (2017, Fig. 10). Then we formed the values Δ|| shown in Fig. 6 (open circles), which reveal a clear dependence on seeing. Their rms of individual values is 1.2 mas, which significantly exceeds the value of 0.35 mas for our reduction. This comparison demonstrates the efficiency of our mitigation of seeing-dependent errors in the epoch-average positions using statistics of residuals Δ′−ψ within individual frame series.

thumbnail Fig. 5.

Motion of component B relative to A for the direct data fit (solid line) for all solutions expected within the 68.3% confidence interval (indicated by the space between the dashed lines), and the measured positions with FORS2 (open circles), HST (triangles), GeMS (squares), and ESO-R obtained in 1984 (the cross with error bars).

Open with DEXTER
thumbnail Fig. 6.

Residuals along the projected distance between components B and A for FORS2 measurements for our reduction (filled circles, σ = 0.35 mas) and the reduction by Garcia et al. (2017; open circles, σ = 1.2 mas).

Open with DEXTER
Table 7.

Barycenter parameters for two versions of model Eq. (8), without (Δ = 0), and with separate offsets (Δ ≠ 0) in FORS2 and HST measurements, and the solution by Garcia et al. (2017).

3.2. Deriving the mass ratio and astrometric parameters

By fitting all HST, FORS2, and GeMS observations with the model for the barycenter motion of Eq. (8) in its simple form with Δ = 0 and ε0 derived in Sect. 3.1, we derived the first solution (left column of Table 7), yielding q = 0.865. This value of q gives the smallest rms of 0.359 mas (0.352 and 0.366 mas along RA and Dec) for the residuals of all observed and model positions shown in the lower panel of Fig. 7 in RA.

To set the uncertainty limits for q, we used the Monte Carlo method, creating random samples of FORS2 and HST measured positions, and processed them with different values of q. For each sample, we determined the best q value as that which corresponds to the minimum rms of the positional residuals. The scatter of these values around the average is characterized by the 1σ interval of 0.0024 and corresponds to the statistical uncertainty of the q value determination. With q fixed to 0.865, Eq. (8) is solved by the least-squares fit as a linear system of equations, and the uncertainty σfit of the parameters derived thus is given in Table 7.

Equation (8), however, is a nonlinear system if q is considered as a free model parameter. The direct solution of this expanded model is same as the solution of the linear model with q = 0.865, but the uncertainties of parameters are quite different. They were estimated by the above Monte Carlo simulation when, together with q, for each random sample we also derived all parameters of the model of Eq. (8). Thus, for each parameter we obtained a set of solutions that corresponds to the expected statistical variations of q. Thus computed 1σ deviations σstat for each parameter are given in Table 7. Unlike the linear model, some parameters of the expanded model (e.g., xc and μx) are highly correlated, which is indicated by the inequality σstat ≫ σfit. In parameter space, the values of these parameters are concentrated along a line with a small scatter corresponding to σfit.

The residuals for GeMS shown in Fig. 7 are qualitatively equal to those obtained by Garcia et al. (2017, Fig. 12). In comparison to FORS2 and HST residuals, they are more scattered, possibly because of insufficiently precise transformation into ICRF. Because of a large scatter σ, the barycenter parameters of LUH 16 and the rms of the residuals were computed with zero weights for the GeMS data. However, we expect that the zero-point errors do not affect the distances between components A and B, and the scatter in the distances measured with GeMS and shown in Fig. 3 are indeed much smaller than those in Fig. 7. Indirectly, these features support our assumption on the origin of Δ in Sect. 2.6, and that it is instrument dependent and should be included in the model of Eq. (6).

The large uncertainties of the chromatic parameters ρ and d in Table 7 are due to the large and expected correlation of about −0.9996, which means that the DCR displacement correction by the LADC is modeled by Eq. (1) with an uncertainty smaller than 0.1 mas.

The high precision in the xc, yc offset determination (±0.066 and ±0.097 mas), which is much better than the uncertainty of the transformation Eq. (5) into the ICRF of about 0.5 mas, motivated us to try out a model function with different offsets Δ(HST) for HST and Δ(FORS2) for FORS2. This produced the second solution with better rms of 0.340 mas achieved at q = 0.8519 (Table 7). In this model we cannot distinguish the offsets, however, and instead obtained xc + Δx(HST)=28.433 ± 0.162, xc + Δx(FORS2)=28.243 ± 0.063, yc + Δy(HST)=329.843 ± 0.159 and yc + Δy(FORS2)=330.449 ± 0.104 mas. Each of these sums corresponds to the zero-point of the barycenter in the FORS2 or HST system, but they are significantly, about 3–4 mas, different from those derived for the first solution. We verified that this is due to the difference in q but not because different offsets Δ were introduced.

The difference between the FORS2 and HST zero-points Δx(FORS2)−Δx(HST)= − 0.190 and Δy(FORS2)−Δy(HST)= + 0.606 mas is comparable to the uncertainty in the zero-points of the transformation into the ICRF, giving a potential explanation for the discrepancy. To confirm this, we derived the corresponding values for a transformation (Sect. 2.6) using Gaia DR1 and obtained significantly larger differences Δx(FORS2)−Δx(HST)= − 0.816 in RA and Δy(FORS2)−Δy(HST)= + 1.199 mas in Dec. This demonstrates the difference between Gaia DR1 and DR2 and supports our assumption on the origin of this discrepancy.

Figure 7 shows the difference between the two solutions: with use of common xc, yc (bottom panel) and offsets different for each data set (middle panel), which we illustrate here in RA. Visually, the gain in overall rms is small (σ decreased from 0.36 to 0.33 mas in RA), but for the subset of FORS2 residuals, the improvement is more significant (0.31 mas versus 0.24 mas). We applied the F-test of additional parameters to find that the improvement is statistically significant. With all FORS2 and HST measurements along RA and Dec, we derived F ≈ 10 with 2 and 114 degrees of freedom, which means that the simpler model is true with a probability < 10−4.

Finally, we converted the barycenter parameters obtained with separate offsets in Table 7 into the absolute system with the correction terms determined in Sect. 2.6.2. We obtained absolute proper motions and parallaxes of μx = −2767.502 ± 0.147 mas yr−1, μy = 356.856 ± 0.150 mas yr−1, and ϖabs = 501.557 ± 0.082 mas.

thumbnail Fig. 7.

Residuals between the FORS2 (open circles), HST (triangles), and GeMS (filled circles) measured barycenter position of LUH 16 and the model of Eq. (8). The middle and upper panels are obtained assuming that offsets Δ are different for each telescope, while the lower panel shows the residuals in RA for the standard model of Eq. (8) with no offset (Δ = 0). The σ values refer to the average rms of residuals in RA and Dec for both HST and FORS2.

Open with DEXTER

3.3. Corrected FORS2 positions

Now we can apply corrections to the FORS2 measurements to remove the difference in offsets Δ. Since we have no handle on the proportion with which the discrepancy occurs, the errors were assumed equal in magnitude but opposite in sign, that is, Δ(FORS2)= − Δ(HST)=Δ. The barycenter zero-point xc, yc and the offsets Δ to the system of equatorial positions at the location of LUH 16 on the CCD obtained in this assumption are given in Table 7. The corrections x(FORS2)  −  Δx, y(FORS2)  −  Δy, x(HST)+Δx, y(HST)+Δy reduce these two data sets of positions to a system anchored at the ICRF. In Table 3 we give the positions(9)

corrected also for DCR and other calibrations, where the expressions for LUH 16 B are similar.

4. Conclusion and discussion

We presented an improved reduction of the FORS2 astrometric measurements of the LUH 16 binary and showed that they are consistent with the HST and GeMS positions in the literature. In our previous reduction (Sahlmann & Lazorenko 2015), the positions of LUH 16 were linked to the ICRF using USNO-B, and as noted by Bedin et al. (2017), are inconsistent with the HST measurements at a level of 10–20 mas. This was due to a flawed transformation into the ICRF that used stars in both chips of FORS2. Before Gaia DR1, this effect could not be diagnosed and lead to a significant bias in RA, Dec of about ±50 mas (Lazorenko & Sahlmann 2017). It had only a negligible effect on the relative positions of the LUH 16 components via the pixel scale and rotation of the coordinate axes, however. The main reason for the inconsistency noted by Bedin et al. (2017) was the incorrect transformation into the ICRF with Eq. (5), where as the argument of functions Fn we used the positions of LUH 16 , extrapolated to the USNO-B epoch J2000.0 instead of using the positions at the actual epochs.

Gaia DR2 provides the reference frame for determining absolute proper motions and parallaxes. However, we argue that very large ground-based telescopes provide us with competitive differential astrometry for faint 16–21 mag stars. In the LUH 16 field analyzed in this paper, the relative proper motions and parallaxes were determined with FORS2 with an average precision of 0.18 mas yr−1 and 0.13 mas, respectively, for 16–21 mag stars in common with DR2. In spite of the short one-year duration of FORS2 observations, this precision is about twice better than the corresponding values 0.48 mas yr−1 and 0.31 mas in the DR2 catalog. The good astrometric performance of FORS2 is due to the higher signal-to-noise ratio in the position estimation.

We refined the individual dynamical masses of LUH 16 to 33.5 ± 0.3 MJup (component A) and 28.6 ± 0.3 MJup (component B), which corresponds to a relative precision of ∼1% and is three to four times more precise than the estimates of Garcia et al. (2017). We found that a minor bias with 1σ of the ESO-R resolved astrometry based on a 1984 photographic plate leads to a change in the dynamical masses that does not exceed the random error of our determination. The exact characterization of any residual bias in mass will be resolved in the future when the high-precision astrometry will cover a larger portion of the orbit.

Because of the importance of LUH 16 as an extremely well-studied system of nearby brown dwarfs, the refinement of the parallax, orbital, and dynamical mass parameters will continue, with the purpose of increasing the knowledge of the physical parameters of the system. We expect that the astrometric data set presented here will contribute significantly to this process.


1

Alternatively, the orbital parameters and their confidence intervals can be derived using a Markov chain Monte Carlo method (Gregory 2005) or permutation (Wilson et al. 2016; Cubillos et al. 2017).

Acknowledgments

This work has made use of data from the ESA space mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

All Tables

Table 1.

Dates of FORS2 observations, average FWHM ε (in arcseconds), and DCR-related functions used in Eq. (1).

Table 2.

Coefficients of Eq. (3) of the systematic bias in the photocenter position of components A and B.

Table 3.

FORS2 Cartesian measured positions xA, yA of the A and xB, yB of the B component, corrected for variations of seeing and flux (assuming that ε0 = 0), but including the DCR displacements and the measurement uncertainties σA, σB.

Table 4.

rms of positional residuals and the excess K in the uncertainties σcat for each data set.

Table 5.

Median orbital parameters of LUH 16AB, with 68.3% confidence intervals for initial and final solutions (with σcat and Kσcat uncertainties of the measurements, respectively), and the corresponding data derived by Garcia et al. (2017).

Table 6.

Radial relative velocity VA − VB measured by CRIRES (Garcia et al. 2017) and the computed minimum/maximum values for our orbital solutions within the 68.3% confidence interval.

Table 7.

Barycenter parameters for two versions of model Eq. (8), without (Δ = 0), and with separate offsets (Δ ≠ 0) in FORS2 and HST measurements, and the solution by Garcia et al. (2017).

All Figures

thumbnail Fig. 1.

Comparison between relative FORS2 and absolute DR2 parallaxes (upper panel) and proper motions in RA (lower panel) for brighter (filled circles; they were used to determine the correction) and fainter stars (open circles), with linear fit functions (solid lines).

Open with DEXTER
In the text
thumbnail Fig. 2.

Residuals in RA (lower panel) and Dec (upper panel) between the FORS2 (open circles), HST (triangles), and GeMS (filled circles) for the relative position of the secondary component of LUH 16 with ε0 = 0 in Eq. (7).

Open with DEXTER
In the text
thumbnail Fig. 3.

Same as Fig. 2, but with ε0 = −0.192′′ in Eq. (7).

Open with DEXTER
In the text
thumbnail Fig. 4.

Median system mass with 1σ uncertainty ranges derived when applying an artificial shift to the ESO-R measurement. The shifts were applied in both coordinates by uσcat toward (u >  0) or away from (u <  0) the model-predicted position, and the parameters KESO-R = 0.57 (open square boxes) or KESO-R = 1 (open circles) were used. The adopted value of Mtot given in Table 5 is shown as a filled square. For comparison, the 1σ uncertainty range (gray area) and the median value (solid horizontal line) computed without the ESO-R measurements are shown.

Open with DEXTER
In the text
thumbnail Fig. 5.

Motion of component B relative to A for the direct data fit (solid line) for all solutions expected within the 68.3% confidence interval (indicated by the space between the dashed lines), and the measured positions with FORS2 (open circles), HST (triangles), GeMS (squares), and ESO-R obtained in 1984 (the cross with error bars).

Open with DEXTER
In the text
thumbnail Fig. 6.

Residuals along the projected distance between components B and A for FORS2 measurements for our reduction (filled circles, σ = 0.35 mas) and the reduction by Garcia et al. (2017; open circles, σ = 1.2 mas).

Open with DEXTER
In the text
thumbnail Fig. 7.

Residuals between the FORS2 (open circles), HST (triangles), and GeMS (filled circles) measured barycenter position of LUH 16 and the model of Eq. (8). The middle and upper panels are obtained assuming that offsets Δ are different for each telescope, while the lower panel shows the residuals in RA for the standard model of Eq. (8) with no offset (Δ = 0). The σ values refer to the average rms of residuals in RA and Dec for both HST and FORS2.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.