© ESO, 2016

1. Introduction

The discovery and characterisation of extrasolar planets is progressing at a staggering pace, fueled by new instrumentation and data analysis methods. Main sequence low-mass stars, the M dwarfs, represent an important target sample because they are the most numerous stars in the Galaxy and host a large number of small planets (Bonfils et al. 2013; Dressing & Charbonneau 2013). Giant planets around M dwarfs, however, are found to have a low occurrence compared to their counterparts around Sun-like stars (e.g. Endl et al. 2006; Cumming et al. 2008; Bonfils et al. 2013), which is an expected outcome of the core accretion scenario for planet formation (Laughlin et al. 2004).

One example for a low-mass star harbouring giant planets is the M0 dwarf GJ 676A, which is part of a wide (~800 AU) binary system of M dwarfs located at a distance of ~17 pc from the Sun. Using radial-velocity monitoring, Forveille et al. (2011) discovered a giant planet (planet “b”) around GJ 676A with minimum mass Msini = 4.9 M_J and an orbital period of 2.9 yr. Forveille et al. (2011) also discovered an additional radial-velocity drift that could not be explained by the wide companion GJ 676B, but required the presence of a second outer companion to GJ 676A. Then Anglada-Escudé & Tuomi (2012, hereafter AT12) reported the detection of two additional super-Earth planets in short-period orbits (planets “d” and “e”) and confirmed the presence of the outer companion, probably a second gas giant (planet “c”).

GJ 676A thus represents a rare case of a planetary system with inner super-Earths and outer gas-giant planets around an M dwarf, a configuration reminiscent of our solar system. It appears that such systems are typically difficult to form (e.g. Raymond et al. 2008), yet other more compact examples have been found, e.g. around the M dwarf GJ  876 (Rivera et al. 2010) and the two Kepler transiting systems KIC 11442793 (Cabrera et al. 2014) and KOI 435 (Ofir & Dreizler 2013).

Here, we present new astrometric, radial velocity, and high-contrast imaging observations that allow us to better characterise the planetary system around GJ 676A.

2. Observations and data reduction

Upon the discovery of planet GJ 676A b, we initiated an astrometric programme to measure the star’s orbital reflex motion caused by the planet, which yields an accurate planet mass measurement by determining the orbital inclination. The minimum semi-major axis of GJ 676A’s barycentric orbit caused by planet “b” is ~0.7 milliarcsecond (mas), thus detectable with high-precision ground-based astrometry that reaches a per-epoch precision of approximately 0.1 mas (Lazorenko et al. 2011; Sahlmann et al. 2014).

2.1. FORS2/VLT astrometry

We obtained optical images of GJ 676A with the FORS2 camera (Appenzeller et al. 1998) installed at the Very Large Telescope (VLT) of the European Southern Observatory (ESO) between April 2010 and August 2012. The instrument setup and observation strategy, e.g. obtaining several dithered frames per epoch, the target position on CCD chip1, and constraints on airmass and atmospheric conditions, is very similar to that used in our exoplanet search survey (Sahlmann et al. 2014) and we reduced the data with the methods developed for that purpose (Lazorenko et al. 2009, 2014). A particularity of this programme is the choice of the OIII-6000 interference filter instead of the I-Bessel filter, which was necessary to decrease the contrast between the comparably bright GJ 676A and field stars.

Table 1 summarises the data we obtained and lists the epoch number, the mean date of the epoch exposures, the average airmass, and the average FWHM measured for star images. There are eight epochs spanning 865 d, and every epoch consists of 49 to 77 usable individual exposures (N_f) taken over Δt = 0.8 h on average, resulting in a total of 517 exposures. Figure 1 shows an example image. We used 272 reference stars located within a radius of 2.́1 of GJ 676A (I ≃ 8.6 mag, V ≃ 9.6 mag) to measure the motion of the star relative to the background field. The majority of the stars are faint and span a magnitude range of I ~ 16–19. Only four stars (including GJ 676B) are relatively bright and 2–5 mag fainter than GJ 676A.

The model for the astrometric reduction takes into account effects of various origin: instrumental (optical distortion, relative motion of the CCD chips), atmospheric (random image motion, differential chromatic refraction), and astrophysical, e.g. the displacement of reference stars due to proper motion and parallax. Optical distortion, in particular, is modelled by fitting the deformation of the reference star positions between frames with basic functions that are polynomials in x, y of powers 0,1,...,k/ 2 − 1. The even integer k is called the mode of the astrometric reduction and varies in the range of 6...12 which corresponds to polynomials of the power 2...5. The model is adjusted simultaneously to all available measurements and has thousands of free parameters, e.g. indexes of differential chromatic refraction and parallaxes of every star. The target object is not used for these reduction steps. Using the derived model, we then corrected the measured photocentre positions of the target object GJ 676A and extracted its astrometric parameters as described in Sect. 3.1. More details on the reduction principles can be found in Lazorenko (2006) and Lazorenko et al. (2009, 2011, 2014). Eventually, we obtained relative astrometric measurements of GJ 676A with an average per-epoch precision of 0.43 mas. The epoch astrometry is given in Table A.1.

GJ 676A is a relatively bright target and the exposure times have to be short to avoid saturation, reducing the number of measurable reference stars. Consequently, the nominal astrometric errors compare unfavourably with the precisions obtained for the ultracool dwarfs that we survey for orbiting planets (Sahlmann et al. 2014), mostly due to a lack of bright reference stars which resulted in higher reference-frame noise. One bright reference star is GJ 676B. To slightly improve precision, we set the parallax of GJ 676B to be equal to that for GJ 676A, noting that the difference of parallaxes in this binary system is small. Because the precision of the parallax measurement corresponds to ≈0.1 pc in distance, the uncertainty of the distance measurement between these two stars is ~20 000 AU. This greatly exceeds the expected value for the relative binary separation, which can be estimated from the sky-projected separation between the two stars. For the measured angular separation of ~47′′ at the distance of ~16.7 pc, the expected value is 800 AU, thus unresolved with our astrometry.

2.2. HARPS radial velocities

Forveille et al. (2011) used 69 high-precision radial-velocity measurements of GJ 676A obtained with the HARPS instrument (Mayor et al. 2003) between 2006 and 2010 to discover planet “b”. We have since continued to observe GJ 676A regularly with HARPS and collected 60 additional measurements between 2010 and 2016, bringing the total number of radial-velocity datapoints to 129. Some of these new observations were taken as part of the volume limited survey (Lo Curto et al. 2010) and of the follow-up on long-period planets (Moutou et al. 2015).

To extract radial velocities (RV), we constructed a high signal-to-noise spectrum by combining all GJ 676A spectra and computed the RV at each epoch with a chi-square minimisation relative to that master spectrum as described by Astudillo-Defru et al. (2015). Only for the two most recent measurements, we used the RV values given by the ESO standard instrument pipeline, which uses the cross-correlation function. The derived velocities used in this study are listed in Table A.2, available at the CDS. This is because those data were obtained after the HARPS upgrade, in which new octagonal fibres were installed (Lo Curto et al. 2015), which changed the line-spread function significantly and many more observations will be needed to generate a second master spectrum for observations taken after the HARPS upgrade. We thus treated the last two measurements as if they were taken with a different instrument and allowed for an offset in the model.

On 2016-03-05 (BJD 57 452.837858) we also obtained one HARPS observation of GJ 676B and measured its RV to − 39.3960 ± 0.0034 km s^-1, which was obtained with the standard pipeline and may include a small (≲1.5 m s^-1) zero-point offset from the measurements for GJ 676A.

3. Analysis of astrometric data

GJ 676A (HIP 85647) was observed 79 times over 1110 days by Hipparcos (ESA 1997) with a median astrometric uncertainty of σ_Λ = 6.1 mas (van Leeuwen 2007), hence covering the orbit of GJ 676A b 1.05 times. Using the methods of combining the radial-velocity orbital parameters with the Hipparcos astrometry described in Sahlmann et al. (2011b,a), we found no orbital signature in the Hipparcos astrometry data: both the permutation test and the F-test yield astrometric orbit significances below 1σ. However, we can use the Hipparcos observations to set an upper limit to the companion mass by determining the minimum detectable astrometric signal a_min of the individual target. When the data cover at least one complete orbit, Sahlmann et al. (2011b,a) showed that an astrometric signal-to-noise of S/N ≳ 6–7 is required to obtain a detection at the 3σ level, where and a is the semi-major axis of the detected barycentric stellar orbit. Using a conservative S/N-limit of 8, we derive the upper companion mass limit M_2,max as the companion mass which introduces the astrometric signal , where the factor 1 − e² accounts for the most unfavourable case of i = 90° and ω = 90° in which the astrometric signal is given by the semi-minor axis of the orbit. Using this criterion and a primary mass of M₁ = 0.71 M_⊙ (Forveille et al. 2011), we set an upper limit of 44 M_J to the mass of GJ 676A b, i.e. this companion must be a substellar object.

3.1. Analysis of FORS2 astrometry

We first analysed the FORS2 astrometry using the standard seven-parameter model without orbital motion as described in, e.g. Sect. 4.1 of Sahlmann et al. (2014). We obtained preliminary astrometric parameters of GJ 676A by obtaining the least-squares solution of Eq. (1) for the photocentre positions determined from the FORS2 images. When neglecting orbital motion, this is a linear model and the solution was obtained using matrix-inversion, taking into account the measurement uncertainties and covariances. (1)The parameters are given in Table 2, and Fig. 2 illustrates the results. In Table 2 we also compare these results to the final adopted solution (see next section) and to the Hipparcos catalogue values, which shows that FORS2 astrometry yields proper motions and a parallax that are compatible with Hipparcos.

As can be seen in the bottom panel of Fig. 2 the standard model does not fit the data well. The residual rms dispersion is 0.54 mas, which is larger than the average epoch precision of 0.43 mas and corresponds to a reduced . More importantly, the curved shape of the residuals hints towards a systematic effect rather than random noise.

To investigate whether the excess correlated signal is associated to the stellar reflex motion caused by planet “b”, we performed a combined, yet sequential, analysis of radial-velocity and astrometric data, which we adapted from Sahlmann et al. (2011b): using the spectroscopic orbital parameters of Forveille et al. (2011), we fitted the FORS2 epoch astrometry that was corrected for DCR with a seven-parameter model, where the free parameters are the inclination i, the longitude of the ascending node Ω, the parallax ϖ, and offsets to the coordinates (Δα^⋆, Δδ) and proper motions (Δμ_α^⋆, Δμ_δ). A two-dimensional grid in i and Ω was searched for its global χ²-minimum with a standard non-linear minimisation procedure. The statistical significance of the derived astrometric orbit was determined with a permutation test employing 1000 pseudo orbits. Uncertainties in the solution parameters were derived by Monte Carlo simulations that include propagation of RV parameter uncertainties. This method has proven to be reliable in detecting orbital signatures in the Hipparcos Intermediate Astrometric Data (e.g. Díaz et al. 2012; Sahlmann & Fekel 2013; Wilson et al. 2016).

This analysis yielded an orbit significance of 99.8% corresponding to better than 3σ on the basis of the permutation test. Therefore the astrometric orbit is clearly detected with the FORS2 astrometry. The preliminary parameters derived with this method are an orbital inclination of 42 ± 9°, corresponding to a planet mass of 7.3 ± 1.3 M_J, and an ascending node of Ω_seq = 200 ± 11°. The residual rms of the best solution is 0.28 mas (corresponding to a reduced χ² of 0.61) and significantly smaller than when employing the model without orbital motion. We repeated the same analysis with the updated orbital parameters of planet “b” obtained with new HARPS RV data (see Sect. 4), which agrees with the Forveille et al. (2011) solution within the uncertainties, and obtained essentially the same results. To derive more accurate model parameters and uncertainties, we performed a joint analysis of radial velocity and astrometry data, which is presented in Sect. 5.

3.2. Parallax correction

Because the astrometric reference stars are not located at infinity, a correction has usually to be applied to the relative parallax to convert it to absolute parallax that allows us to determine the distance to the system. As in Sahlmann et al. (2014), we used the Galaxy model of Robin et al. (2003) to obtain a large sample of pseudo-stars in the region around GJ 676A. The comparison between the model parallaxes and the measured relative parallaxes of stars covering the same magnitude range yields an average offset, which is the parallax correction Δϖ_galax. The absolute parallax ϖ_abs = ϖ − Δϖ_galax is larger than the relative parallax because the reference stars absorb a small portion of the parallactic motion, i.e. the parallax correction has to be negative.

Using N_s = 142 reference stars, we obtained a parallax correction of Δϖ_galax = + 0.12 ± 0.24 mas for GJ 676A. Because the reference stars are much fainter than GJ 676A, their parallaxes have large uncertainties, which translated into a large uncertainty of the parallax correction. The correction is smaller than its uncertainty, i.e. it is compatible with zero, and it has a positive value which is not allowed by definition. Therefore we did not apply the correction to the parallax of GJ 676A and we set ϖ_abs = ϖ.

In principle, a similar procedure should be applied to correct from relative to absolute proper motion. We refrain from doing so because proper motions are not critical parameters in the following analyses and their corrections will be small. In the future, the results of ESA’s Gaia mission will make it possible to determine model-independent parallax and proper motion corrections because Gaia will obtain accurate astrometry for many of the reference stars and for GJ 676A and GJ 676B themselves.

4. Analysis of radial velocities

We analysed the radial velocities by using the generalised Lomb-Scargle periodogram (Zechmeister & Kürster 2009) iteratively. We started by modelling the raw time series with one Keplerian (for planet “b”) plus a cubic drift (for planet “c”) and we computed the periodogram of the residuals of the best fit (top panel of Fig. 3). A powerful peak is seen around a period of 3.6 days. To give a false-alarm probability (FAP), we assumed that the residuals are caused by random noise and we generated virtual datasets by swapping the radial-velocity values randomly while retaining their dates. For every virtual set we computed a new periodogram and measured the power of the highest peak. In this way we obtained the statistical distribution of power maxima that is expected from a timeseries that contains solely noise. From that distribution, the power values corresponding to a FAP of 1, 10, and 50% are the power values found to be greater than 99, 90, and 50% of the distribution, respectively. In Fig. 3 those power levels are drawn with gray, dark-gray, and black lines, respectively, and the peak corresponding to planet “d” is recover with a FAP ≪1%. We next included an additional planet and applied a model composed of 2 Keplerians plus a cubic drift. The most prominent peak now has a period of ~36 days and a FAP marginally below 1%. After iterating once more with yet an additional planet, no significant-power periodicity is seen in the final residuals. The most powerful peak in the last residual periodogram is located at ~1600 days and has a FAP of ~10%. If we were to interpret this as the signature of a yet undiscovered planet “f” and modelled it accordingly, the corresponding mass is ~35 M_⊕, i.e. about two Neptune masses.

To derive the model parameters corresponding to the planetary signals, we performed a Markov chain Monte Carlo (MCMC) analysis of the HARPS radial velocities. The model incorporated three Keplerians (planets “b”+“d” and signal “e”), one circular Keplerian for planet “c” (we found that both a circular Keplerian and a cubic drift models reproduce the signal equally well), the systemic velocity, and an offset to account for the HARPS upgrade (see Sect. 5), for a total of 20 free parameters. The parameter values are reported in Table 3 and compared to the solution presented by AT12. Figure 4 shows the four Keplerian curves and the data. Due to the larger number of measurements, our parameters for planets “b” and “d” are generally more precise but in good agreement with the AT12 parameters. For planet “c” we could derive the first, yet preliminary, good constraints on period and minimum mass, which are further discussed in Sect. 5. In comparison to AT12 for the signal attributed to planet “e”, we find a lower eccentricity and a smaller signal amplitude, hence a smaller minimum planet mass of 8.1 ± 0.7 M_⊕.

We thus have recovered all signals previously reported by Forveille et al. (2011) and AT12. Before accepting the inner low-mass planets as genuine, however, we need to evaluate if the corresponding signals can alternatively be attributed to stellar activity. In particular, the period of planet “e” (~36 d) is close to the stellar rotation period of 41.2 ± 3.8 d measured by Suárez Mascareño et al. (2015). We also measured a chromospheric activity indicator that is very close to the one of GJ 205 (), which has a known rotation period of 33.6 d (Kiraga & Stepien 2007). The rotation period of GJ 676A is thus close to the orbital period of planet “e” and may induce a signal that is confused with the one of a planet. Furthermore, looking at the radial-velocity residuals when adjusting for all planets except for planet “e”, we noticed a strong correlation between the radial velocity and the activity index for epochs MJD 54 660–54 690, see Fig. 5. This range includes only 10 radial-velocity points, but when removing just those 10 points from the original time series, the periodogram power of signal “e” decreases by as much as 20% as shown in Fig. 6. We inspected several other activity indicators, e.g. the width of the cross-correlation function and its contrast, the S index, Hα and Sodium indices, and were not able to identify the stellar rotation period. In summary, we remain cautious with the interpretation of signal “e” as being caused by a planet because its period is close to the star’s rotation period. Additional data and analyses are required to undoubtably establish the planetary nature of signal “e”. On the contrary, the period of planet “d” is sufficiently short compared to the stellar rotation period to be accepted as a planet.

5. Joint analysis of radial velocities and astrometry

We applied an MCMC analysis to the individual radial-velocity data in Table A.2 and the astrometric measurements from our FORS2 observations. We used a global model with 21 free parameters that has five components:

• Radial velocity orbit of planet “b”: there are 5 orbital parameters:the period P, eccentricity e, argument of periastron ω, time of periastron passage T_P, companion mass M_b, and one offset γ₀ corresponding to the systemic velocity.

• Astrometric orbit of planet “b”: the model comprises seven free parameters, of which five are shared with the radial velocity model (P,e,ω,T_P,M_b) and two are uniquely constrained by astrometry: the inclination i and the ascending node Ω. We also included two nuisance parameters s_α and s_δ for the astrometry in RA and Dec, respectively, to account for the off-diagonal terms in the covariance matrix of the FORS2 astrometry in individual frames (Sahlmann et al. 2013). Because of the long period of planet “c”, its potential non-linear astrometric signature is much smaller than our uncertainties and the linear part will be absorbed by a bias in the proper motion values (see Sect. 5.4).

• Radial velocity signature of planet “c”: the orbital period of planet “c” is longer than the observation timespan (see Sect. 4). We modelled its signature with a circular Keplerian model (e = 0, ω = 0) that has three parameters: orbital period P_c, minimum mass M_csini, and time of ascending node T_c,P.

• Parallax and proper motion: the standard astrometric model has five free parameters (position offsets ,Δδ₀, parallax ϖ, proper motions μ_α^⋆,μ_δ) plus two parameters modelling differential chromatic refraction ρ and d.

• Radial velocity offset after HARPS upgrade: two RV measurement were taken with HARPS after the fibre upgrade (Lo Curto et al. 2015). To account for a potential instrumental offset, we included the free parameter Δγ₀ which is added to the RV data taken in the new HARPS configuration. With the help of 5 other stars with spectral types of M0–M4 observed with HARPS both before and after the upgrade, we found that the RV offset is compatible with zero with an uncertainty of 1.5 m s^-1.

We used the emcee package (Foreman-Mackey et al. 2013) to implement the MCMC and expressed the global model with the parameter vector θ composed of P, e, ω, T_P, γ₀, Δγ₀, log P_c, M_csini, T_P,c, , Δδ₀, ϖ, μ_α^⋆, μ_δ, ρ, d, s_α, s_δ, M_bsini, M_bcosi, and Ω, where we chose the pair M_bsini–M_bcosi instead of M_b–i to mitigate the effect of the strong correlation that naturally exists between those parameters. The host star mass M₁ was kept constant. The logarithm of the likelihood as a function of the parameter vector θ was expressed as (2)where the subscripts RV and AX denote the radial-velocity and astrometric part of the model, respectively. The radial velocity model ℳ_RV(θ) implements the Keplerian equations for planets “b” and “c” and the corresponding χ² is computed in the standard way (3)where denotes measured quantities and is the variance associated with the ith measurement. For astrometry, the model ℳ_AX(θ) implements Eq. (1) and we account for the nuisance parameters s_α and s_δ when computing the total variance of an individual measurement and in the lnℒ₀ term For some parameters we applied uniform priors (see Table 4) and imposed the following range limits: , M_bsini> 0, ϖ> 0, log 5000 < log P_c< log 9000, s_α> 0, and s_δ> 0. For the RV offset Δγ₀, we applied a Gaussian prior centred on zero with a width of 1.5 m s^-1. For a more general discussion on combined modelling of radial-velocity and astrometric data see Wright & Howard (2009) and Anglada-Escudé et al. (2012).

We did not include the astrometric orbit terms of planets “d” and “e” in the model because their signatures are approximately 0.05 and 0.5 micro-arcsecond, respectively, (estimated for edge-on orbits with the parameters of AT12), thus are negligible. Likewise, we did not include their radial-velocity terms (discussed in Sect. 4) because they have sufficiently small amplitudes and short periods that their omission does not affect the parameters of the large-amplitude and long-period signals of planets “b” and “c”.

Each of 160 walkers was initialised with a set of parameter values that was determined from the radial-velocity orbit parameters, from the standard astrometric fit in Sect. 3.1, and from the preliminary values for orbit inclination and planet mass derived from the sequential analysis in Sect. 3.1. Each walker was allowed to take 30 000 steps, of which we discarded the first 25%. The solution was therefore derived from distributions with 3.6 × 10⁶ samples. Table 5 lists the adopted solution parameters determined as the median of the posterior distributions with 1σ uncertainties.

5.1. Parallax and relative proper motions

Our parallax determination agrees with the Hipparcos parallax within the uncertainties, yet it is slightly smaller leading to a larger distance of 16.86 ± 0.07 pc to GJ 676A. This in turn would lead to a slightly higher mass estimate for the star, however, we assume a 10% uncertainty for the primary mass determination M₁ = 0.71 M_⊙ (Forveille et al. 2011, AT12), which renders this adjustment insignificant.

Our relative proper motions are discrepant from the absolute measurements derived by Hipparcos (van Leeuwen 2007) at the 5σ–8σ level (cf. Table 2). This can be explained by the intrinsically relative measurements accessible by FORS2, but also by the fact that Hipparcos and FORS2 proper motions are biased differently by the orbital motion of planet “c”. The high precision of our relative proper motion measurement allowed us to perform a detailed study of the relative orbit of the wide binary, which is presented in Sect. 6.

5.2. Radial velocity orbits

The radial velocity data and the best-fit model are shown in Fig. 7. The orbital parameters of planet “b” in Table 5 are compatible with the estimates of Forveille et al. (2011) and AT12, yet they are more precise because we have more data available. In particular, the minimum mass of planet “b” is M_bsini = 4.733 ± 0.011 M_J. The residual rms of our model that does not include any of the two inner planets “d” and “e” is 3.16 m s^-1.

Although our RV data do not cover one full revolution of planet “c”, its orbital period and RV amplitude is relatively well constrained. Assuming a circular orbit, we find a period of 7337 ± 95 days (~20 yr) with a semi-amplitude of ~89 m s^-1, which corresponds to a minimum planet mass of 6.8 ± 0.1 M_J, where the uncertainty does not include the 10% uncertainty on the mass of the star. The corresponding relative semi-major axis is 6.6 ± 0.1 AU, placing planet “c” between the orbital distances of Jupiter and Saturn in the solar system. Until the eccentricity of planet c’s orbit can be determined with additional measurements, these values should be considered preliminary.

5.3. Astrometric orbit and the mass of GJ 676A b

The MCMC chains converged towards stable solutions that produced quasi-Gaussian posterior distributions for most of the 21 free parameters, which indicated that the model is well-constrained and the astrometric orbit was detected with our FORS2 measurements. Figure A.1 displays the joint marginal distributions of all MCMC parameters and shows that correlations are generally weak, with the exception of the expected inter-dependencies between the periastron parameters T_P and ω and the chromatic refraction parameters ρ and d that are anti-correlated by design.

The joint marginal distributions of orbit parameters that are constrained by astrometry are shown in Fig. 8. All three (Ω, M_bsini, M_bcosi) are well constrained and are weakly correlated. After conversion to the M_b and i parameters, which are also displayed in Fig. 8, the lower mass limit imposed by the radial velocities becomes apparent as a pile-up at ~4.7 M_J, but a peak of the mass distribution at ~6.7 M_J can be clearly identified.

We conclude that the astrometric orbit of planet “b” was detected with our ground-based astrometry. The adopted parameter values are °, °, and the semi-major axis of GJ 676A’s reflex motion is mas. This corresponds to a mass of planet GJ 676A b of . These parameters (and the parallax and proper motions) are in good agreement with the preliminary values determined from the independent sequential analysis in Sect. 3.1. The planet mass uncertainty does not account for the uncertainty in the host star mass.

Figure 9 shows the astrometric orbital motion as a function of time. The epoch residuals of the best-fit orbit model have an rms dispersion of 0.28 mas (reduced ), which is almost twice smaller that the residuals of 0.54 mas for the standard model without orbit presented in Sect. 3.1. This is further strong evidence that the astrometric orbit was detected. The residual dispersion is smaller than the average uncertainty of 0.43 mas, which may indicate that the latter are slightly overestimated. Finally, we present the astrometric orbit of GJ 676A caused by planet “b” in the plane of the sky in Fig. 10.

5.4. FORS2 proper motion bias caused by planet “c”

Because our FORS2 measurements cover only ~12% of the orbital period of planet “c”, the corresponding reflex motion of GJ 676A results in a small bias of the measured proper motion reported in Tables 5 and 2. To quantify the bias we have to take into account the inclination i_c and orientation Ω_c of planet c’s orbit. The distribution of mutual inclinations for long-period giant planets in multi-planetary systems is essentially unknown because inclinations were constrained in only a few systems (McArthur et al. 2010; Fabrycky & Murray-Clay 2010; Correia et al. 2010; Pueyo et al. 2015). We conservatively assumed a circular orbit for planet “c” with an inclination that matches the one of planet “b”, which results in a planet mass of M_c = 9.5 M_J. Furthermore, we assumed that the ascending node Ω_c of planet “c” is randomly oriented, which translates into a random distribution of mutual inclinations between “b” and “c”.

We simulated astrometric data at the FORS2 observation epochs in a Monte Carlo fashion that take into account the orbit of planet “c” and fitted those data with the standard parallax + proper motion model. The difference between the proper motion determined from simulated data with and without planet “c” yields an estimate of the proper motion bias, which amounts to 0.0 ± 1.1 mas yr^-1 in both μ_α^⋆ and μ_δ. If instead we assume that the planets have both the same inclination and the same ascending node, the corresponding bias is − 0.5 ± 0.3 mas yr^-1 and − 1.4 ± 0.4 mas yr^-1 in μ_α^⋆ and μ_δ, respectively. These biases scale with the mass of planet “c” and may slightly increase if its orbit is eccentric.

6. Orbital motion of the GJ 676A – GJ 676B binary

We measured the proper motion of GJ 676B relative to GJ 676A, which is useful to characterise the dynamics of this wide binary system. The proper motion differences Δμ_α^⋆ = μ_α^⋆(GJ 676B) − μ_α^⋆(GJ 676A) and Δμ_δ = μ_δ(GJ 676B) − μ_δ(GJ 676A) are robust observables because they are independent of the zero-points of proper motion and of common systematic error components. We performed an astrometric reduction for GJ 676B, using the same procedures applied to GJ 676A (Sect. 2.1). Table 6 lists the values of and Δμ_δ obtained from our data and from three literature catalogues (PPMXL, NOMAD, and UCAC4) that present proper motion values for both stars, which means that those proper motions were determined in the same system, thus mitigating systematic error in the differences and Δμ_δ.

The literature relative motions are marginally significant and variable in their signs which leaves the direction of relative motion uncertain. In contrast, our astrometry clearly detects the relative motion of GJ 676B and GJ 676A. When investigating the relative orbital motion in this stellar system, we have to account for the influence of the giant planets around GJ 676A. The proper motion bias due to planet “b” is negligible because our FORS2 data cover almost a full orbit. However, the bias due to planet “c” that we estimated in Sect. 4 has the effect of increasing the formal uncertainties of GJ 676A’s proper motion by ± 1.1 mas yr^-1 both in RA and Dec.

To study the orbital configuration of the binary we used the constraints set by the projected distance between the two stars Δα^∗, Δδ and by the relative proper motion Δμ_α^⋆, Δμ_δ, which in the plane of the orbit is 8.87 ± 1.33 mas yr^-1 and which translates into a tangential velocity of V_t = 0.704 ± 0.106 km s^-1. We also know GJ 676A’s systemic velocity γ₀ = − 39.038 km s^-1  (Table 5) and we measured the radial velocity of GJ 676B as − 39.396 ± 0.003 km s^-1. The radial velocity of GJ 676B relative to GJ 676A is then ΔRV = − 0.358 ± 0.010 km s^-1, where we accounted for an additional RV uncertainty due to planet c’s incomplete orbit and to potential zero-point offsets. Thus, we have constrained all components of the GJ 676B’s relative velocity vector and of its relative spatial position, except for the distance z to GJ 676A along the line of sight.

The binary’s Keplerian motion is described by six parameters: the semi-major axis a of the relative orbit, eccentricity e, time of periastron passage, inclination i, argument of periastron ω, and ascending node Ω. We assumed a system mass of 1.02 M_⊙, with component masses of 0.71 M_⊙ and 0.29 M_⊙ for the stars (Forveille et al. 2011) and ~0.02 M_⊙ for the giant planets.

We ran Monte Carlo simulations to explore the allowed parameter values under the constraints discussed above. The result in the a − e plane is shown in Fig. 11. The allowed parameter space (at 3σ level) is split in two separated domains: for short line-of-sight distances z<z₀ between GJ 676A and GJ 676B, the allowed parameter space is the area inbetween the solid lines. For large distances z>z₀, the allowed area is delimited by dashed lines. The value z₀ ≈ + 360 AU is the particular point where the velocity vector points towards GJ 676A. The allowed range of z values is ± 7500 AU, which is unresolved with the FORS2 astrometry. However, Gaia may be sensitive to the cases of extreme z-values. For V = 7–12 stars, the parallax precision of Gaia¹ is ~0.01 mas (e.g. de Bruijne 2012) or ~800 AU if expressed in terms of the distance z. After solving for the astrometric orbit of GJ 676A due to planet “b”, Gaia will distinguish at 3σ confidence between solutions with large negative z< − 2400 AU (hatched area between solid lines in Fig. 11), large positive z> 2400 AU (filled area between dashed lines in Fig. 11), and intermediate z values.

The space of allowed a − e values is sensitive to errors in proper motion, especially in RA, whereas the uncertainties in RV are negligible. With the Gaia results, we expect that the allowed a − e values will be significantly better constrained, although they will neither resolve the correlation between a and e nor remove the ambiguity of solutions with z < z₀ and z > z₀.

Estimates of the angular orbital parameters are strongly affected by the uncertainty in the proper motion bias due to planet “c” and by the ambiguity in the z value. However, for moderately eccentric orbits with e < 0.9 we find that the inclination is constrained between 50 < i < 120°. The allowed ranges for the ascending node are 60 < Ω < 130° if z > z₀ and 260 < Ω < 300° if z < z₀.

The periastron distance a (1 − e) of the binary is of particular interest because of its influence on the dynamics of the planetary system around GJ 676A. We found that its upper limit is 1000–7000 AU almost independently of eccentricity, whereas its minimum value depends strongly on e. For nearly circular orbits, it is about 1000 AU and decreases to 100–200 AU for e ≈ 0.8 and z <z₀. For even more eccentric orbits, the periastron distance can be as small as 10 AU, and the limiting value is 5 AU (see Fig. 12), which also corresponds to the shortest orbital period of ~10 000 yr (Fig. 11). We conclude that the relative distance between GJ 676B and GJ 676A throughout their orbit most likely remains much larger than the extend of the known planetary system around GJ 676A, thus the wide binary is unlikely to have affected the formation and evolution of these planets. Additional RV measurements for GJ 676B and the measurements of the Gaia mission will make a better characterisation of the binary orbit possible.

6.1. Trigonometric parallax of GJ 676B

In Sect. 2.1 we explained that the trigonometric parallax of GJ 676B was derived under the assumption that it is indistinguishable from that of GJ 676A, which was necessary to mitigate the effects of the small number of reference stars. Therefore, the parallax 59.3 ± 0.3 mas in Table 5 applies equally to GJ 676B, which to our knowledge is the first parallax measurement for this star. Besides, the orbit modelling in this section yields an upper limit of ±7500 AU for the difference in distances z to these stars. Hence, their parallax difference cannot exceed ± 0.128 mas.

7. Prospects for directly imaging the gas giants around GJ 676A

Obtaining images of extrasolar planets makes it possible to directly measure their luminosities and spectra (e.g. Lagrange et al. 2010; Chilcote et al. 2015). With present-day instrumentation, this technique is mostly limited to young self-luminous giant planets in very long-period orbits (≳10 AU) around nearby young stars (e.g. Macintosh et al. 2015). The next generation of exoplanet imaging efforts on the ground and in space will target long-period, mature Jovian planets similar to GJ 676A b (Traub et al. 2016; Kasper 2015) and will pave the path towards observing the spectrum of a potential Earth-twin, which will require a large space observatory. Imaging the giant planets around GJ 676A is a stepping stone on this path.

The extreme technical requirements of high-contrast instrumentation, combined with the scientific motivation to prioritise atmospheric characterisation over blind imaging searches, creates a strong incentive to leverage radial-velocity and astrometry programmes to cull a target sample with well-constrained ephemeris and masses. Therefore, the task of assessing the observability of indirectly detected exoplanets is already underway, well in advance of the commissioning of such facilities (e.g. Howard & Fulton 2014; Crossfield 2013).

7.1. VLT/NaCo observations

GJ 676A was observed on 2010-10-15 with VLT/NaCo (Rousset et al. 2003; Lenzen et al. 2003) as part of a project aimed at imaging potential sub-stellar companions of red dwarfs that exhibit a radial-velocity drift (ESO Programme 086.C-0515(A), PI Montagnier). The 20-min K_S-band observing sequence was taken in saturated non-coronagraphic field-tracking mode, bracketed by two short non-saturated sequences for photometric calibration. The data were analysed using GRAPHIC (Hagelberg et al. 2016) and the reduced image is shown in Fig. 13. No additional point-source was detected in this snapshot observation. The detection limits shown in Fig. 14 were derived from the reduced image following a procedure similar to the one described in Chauvin et al. (2015), and the mass estimates are based on the BT-SETTL CIFIST2011 models (Baraffe et al. 2015) using a conservative age estimate of 5 Gyr for the system.

According to the best orbital solution determined in Sect. 5, the relative separation of planet “b” at the NaCo observation epoch was 70 mas. This is much smaller than the inner working angle (IWA) of ~200 mas achieved with this dataset, thus NaCo does not yield any constraint on planet “b”.

For planet “c” in the most conservative case where we assume that its orbit is seen edge-on and M_c = M_csini = 6.8 ± 0.1 M_J, the relative separation at the NaCo epoch was 170 mas. Figure 14 shows that at the smallest separation probed with NaCo we can exclude an object with mass ≳75 M_J with high confidence. On the basis of the NaCo data alone, we can thus exclude that the RV signature corresponding to planet “c” is caused by a stellar companion in an almost face-on orbit, which otherwise would be detectable in the images. Furthermore, the relative separation at the NaCo epoch increases with companion mass. If the orbital inclinations of planets “b” and “c” are equal, in which case M_c = 9.5 M_J, the relative separation at the NaCo epoch increases to 300 mas and the corresponding upper 5σ mass-limit derived from the images is ~52 M_J. Following this argument, we can set a lower limit of 10° to the inclination of planet “c” because this configuration with M_c = 39 M_J and separation 390 mas is ruled out by the NaCo images. The mass of planet “c” has therefore to be in the range of 6.8–39 M_J.

7.2. Thermal imaging of planets “b” and “c”

Even at its minimum mass of 4.9 M_J, Quanz et al. (2015) noted GJ 676A b as a potential target for thermal infrared imaging with the future METIS instrument on the European Extremely Large Telescope (E-ELT). Assuming an age of 5 Gyr, the planet is above the planned sensitivity limits of METIS in the thermal infrared L and M bandpasses. Since the true mass of GJ 676A b established here is yet higher, it is confirmed as an attractive target for ground-based infrared imaging efforts with the E-ELT or other extremely large telescopes with mid-infrared high-contrast instruments.

Moreover, we put forward planet “c” as a new promising target. Its minimum mass is M_csini = 6.8 ± 0.1 M_J, thus approximately the same as the mass of planet “b”, but GJ 676A c is on a wider orbit with larger relative separation from its host star. Using the estimates for the 6.8 M_J planet HD 111232 b given in Table 2 of Quanz et al. (2015), we extrapolated apparent L-band magnitudes of ~20.5 and ≲20.5 and maximum separations of 110 mas and 390 mas for GJ 676A b and GJ 676A c, respectively. Thus both planets can be observed with METIS given the expected limiting magnitude of L < 22.4 and IWA of 38 mas (Quanz et al. 2015).

We do not expect JWST to image the gas giants around GJ 676A, due to their small angular separations relative to the inner working angle limits of the various coronagraph modes. In principle, the 0.36′′ inner working angle of the four quadrant phase mask (FQPM) coronagraph at 11.4 microns enables searching for planets at the same angular separation as planet “c”, with apoastron angular separation 0.4′′ (Boccaletti et al. 2015). However, the predicted PSF subtraction residuals within 1′′ of the star will remain too high to reach the required contrast (~10^-5) to detect such a cool planet in the mid-infrared. The expected 3–5 μm planet-to-star contrast ratios of the planets (~2 × 10^-6) are also too extreme for the NIRCam coronagraph detection limits inside 0.5′′ (Beichman et al. 2010) and for NIRISS aperture masking interferometry (Artigau et al. 2014).

7.3. Reflected light observations with WFIRST

We investigated suitability of GJ 676A b and GJ 676A c for observation in reflected starlight with the visible wavelength, space-based coronagraph planned for NASA’s Wide-Field Infrared Survey Telescope (WFIRST, Noecker et al. 2016).

Due to their super-Jovian masses, we expect the planets’ internal energy to dominate the stellar irradiation in determining its effective temperature. In the giant-planet evolution models of Burrows et al. (2004), the effective temperatures (T_eff) of 5 Gyr-old gas giants of masses 6 M_J and 8 M_J are 216 K and 251 K, respectively. By comparison, even at planet b’s periastron (~1.24 AU), for an assumed GJ 676A bolometric luminosity of L_⋆ = 0.33 L_⊙, and a Bond albedo as low as 0.1, the equilibrium temperature would peak at only 185 K.

In the T_eff range of 200–300 K that evolution models predict for planets “b” and “c”, water condenses in the troposphere, resulting in a higher albedo than Jovian-type (cooler and NH₃ cloud-dominated) gas giants (Sudarsky et al. 2000). At a wavelength of 480 nm, the respective geometric albedos of Jupiter and Saturn are 0.46 and 0.39, respectively (Karkoschka 1998). In the atmosphere models computed by Sudarsky et al. (2005), a giant planet with T_eff similar to GJ 676A b would have a geometric albedo of 0.5 at 500 nm, rising towards 400 nm.

In Fig. 15 we show the estimated contrast of reflected light from planet GJ 676A b in the WFIRST blue channel as a function of time during a 6-yr window that may correspond to the WFIRST mission. We used the Keplerian orbital elements and a Lambert sphere scattering model with a classical phase function to compute the time-evolution of projected angular separation and the planet-to-star contrast ratio, and we assumed that the planets have radii equivalent to Jupiter. The angular separation plot captures the uncertainty in the ephemeris with a 1σ contour of an ensemble of draws from the posterior distributions of the orbital solution. The contrast curve is repeated for three geometric albedos, 0.4, 0.5, and 0.6.

In Fig. 16 we show contrast and separation in the sky plane. This representation allows us to define the orientations that corresponds to regions of optimal separation and contrast. In this way, both the epoch and the expected location of planet “b” are defined, which is crucial to optimise the efficiency of these resource-intensive observations. For instance, the instrument can be oriented in a way that the planet falls onto a preferred zone in the image plane.

At a wavelength of ~500 nm and maximum projected separation of 110 mas, the planet-to-star contrast is ~2.6 × 10^-9, above the nominal post-processed detection limits of the WFIRST Coronagraph Instrument (Krist et al. 2015). However, this angular separation is just inside the baseline 3.0 λ/D inner working angle of the instrument that corresponds to 120 mas at 465 nm, the central wavelength of the bluest imaging bandpass. Therefore, the feasibility of imaging GJ 676A b with WFIRST remains uncertain. Ultimately, it will depend on the final coronagraph mask specifications, which are subject to ongoing design studies trading between throughput, contrast, inner working angle, and bandwidth (Trauger et al. 2016; Zimmerman et al. 2016).

We performed a similar calculation for GJ 676A c using its preliminary ephemeris under the assumption of an edge-on circular orbit. We found that the planet-to-star contrast in reflected light reaches the range of 1 × 10^-9–2 × 10^-9 for separations between 200 mas and 390 mas. Because of its larger orbital separation planet “c” receives less insolation and is therefore as challenging as planet “b” at maximum separation in terms of contrast. At the same time its larger projected separation makes it a promising target for WFIRST coronagraphy.

8. Discussion

The astrometric measurement of the reflex stellar motion caused by orbiting planets is difficult with current instruments because the amplitude of the signal is approximately 1 mas or smaller. The detection of such small signals is significantly eased by targeting stars with giant planets previously characterised with radial velocities. So far, it has been predominantly achieved from space with the Hubble Space Telescope fine guidance sensor (Benedict et al. 2006, 2010; McArthur et al. 2014). In addition, upper limits on the masses of known exoplanets were set with the Hipparcos space mission (e.g. Perryman et al. 1996; Reffert & Quirrenbach 2011; Sahlmann et al. 2011a) and in rare cases from the ground (Anglada-Escudé et al. 2012).

Here, we demonstrate the detection of the astrometric orbit of an exoplanet host star with a ground-based instrument. Thanks to the outstanding astrometric precision achieved with FORS2/VLT and the detailed knowledge of the orbit’s spectroscopic orbital parameters from HARPS RV monitoring, the detection is made with high confidence and we can constrain the inclination and ascending node of GJ 676A’s orbit caused by planet “b”. As a direct consequence of measuring the orbital inclination, we determine the mass of planet “b” to .

We measured the minimum mass of the outer giant planet “c” M_csini = 6.8 ± 0.1 M_J under the preliminary assumption of a circular orbit and constrained its inclination to be larger then 10°. We do not set observational constraints on the orbital inclinations of the inner planets, but if the system is aligned, i.e. all planets share similar inclinations like the planets in the solar system, the masses of planets “c”, “d”, and “e” are ~40% higher than their minimum masses derived from RV. In the future, the determination of the mutual inclinations between the planets using astrometry or other techniques may hint on the dynamical history to the system.

9. Conclusions

We pursued the detailed characterisation of the planetary system around the M0 dwarf GJ 676A using a wide range of observational techniques (astrometry, radial velocity, high-contrast imaging) and instruments (FORS2, HARPS, NaCo, Hipparcos), which lead to the determination of the mass of planet “b” and the preliminary measurement of the minimum mass of planet “c”. We confirmed the presence of the inner planet “d” and of the periodic signal associated with planet “e”, whose period is close to the star’s rotation period. We also find tentative evidence for an additional periodic RV signal at ~1600 days.

We demonstrated how astrometry can leverage radial-velocity planet searches to identify the best targets for future high-contrast direct-imaging observations. The determination of the astrometric orbit yields the comprehensive ephemeris of the planet-star system, which is crucial for efficient planning of the observation timing and the setup/orientation of the instrument. Gaia’s astrometry will make it possible to extend these efforts to many other stars with already known or newly discovered giant planets. We showed that the outer giant planets around GJ 676A are promising targets for direct imaging of their thermal radiation and/or reflected light with future facilities, in particular extremely large ground-based telescopes like the E-ELT and space missions like WFIRST.

At a distance of 16.7 pc, GJ 676A – together with its rich planetary system with small inner planets, giant outer planets, and a wide binary companion – represents a fascinating outcome of the star and planet formation process.

---

Acknowledgments

J.S. is supported by an ESA Research Fellowship in Space Science. J.H. is supported by the Swiss National Science Foundation (SNSF). X.B., X.D. and T.F. acknowledge the support of the French Agence Nationale de la Recherche (ANR), under the programme ANR-12-BS05-0012 Exo-atmos, as well as funding from the European Research Council under the ERC Grant Agreement No. 337591-ExTrA. This research made use of the databases at the Centre de Données astronomiques de Strasbourg (http://cds.u-strasbg.fr), NASA’s Astrophysics Data System Service (http://adsabs.harvard.edu/abstract_service.html), the paper repositories at arXiv, of APLPY (Robitaille & Bressert 2012), and of ASTROPY, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013). Some figures were produced with corner.py (http://dx.doi.org/10.5281/zenodo.45906). The authors also made use of SCIPY (Jones et al. 2001), NUMPY (Oliphant 2007), IPYTHON (Pérez & Granger 2007), and MATPLOTLIB (Hunter 2007).

References

Appendix A: Figures and tables

All Tables

All Figures