The HADES RV Programme with HARPSN at TNG
VIII. GJ15A: a multiple wide planetary system sculpted by binary interaction^{★}
^{1}
INAF – Osservatorio Astrofisico di Torino,
Via Osservatorio 20,
10025
Pino Torinese, Italy
email: m.pinamonti.astro@gmail.com
^{2}
Dipartimento di Fisica e Astronomia, Università di Padova,
Via Marzolo 8,
35131
Padova, Italy
^{3}
INAF – Osservatorio Astronomico di Padova,
Vicolo dell’Osservatorio 5,
35122
Padova, Italy
^{4}
INAF – Osservatorio Astronomico di Palermo,
Piazza del Parlamento 1,
90134
Palermo, Italy
^{5}
INAF – Osservatorio Astrofisico di Catania,
Via S. Sofia 78,
95123
Catania, Italy
^{6}
INAF – Osservatorio Astronomico di Trieste,
Via G. B. Tiepolo 11,
34143
Trieste, Italy
^{7}
Fundación Galileo Galilei – INAF,
Ramble José Ana Fernandez Pérez 7,
38712
Breña Baja,
TF, Spain
^{8}
Dipartimento di Fisica e Chimica, Università di Palermo,
Piazza del Parlamento 1,
90134
Palermo, Italy
^{9}
INAF – Osservatorio Astronomico di Brera,
Via E. Bianchi 46,
23807
Merate, Italy
^{10}
Instituto de Astrofísica de Canarias (IAC),
38205
La Laguna,
Tenerife, Spain
^{11}
Departamento de Astrofísica, Universidad de La Laguna,
38206
La Laguna,
Tenerife, Spain
^{12}
Institut de Ciències de l’Espai (ICE, CSIC), Campus UAB,
C/ de Can Magrans s/n,
08193
Cerdanyola del Vallès, Spain
^{13}
Observatoire Astronomique de l’Université de Genève,
1290
Versoix, Switzerland
^{14}
Institut d’Estudis Espacials de Catalunya (IEEC),
C/ Gran Capità 24,
08034
Barcelona, Spain
Received:
22
December
2017
Accepted:
10
April
2018
We present 20 yr of radial velocity (RV) measurements of the M1 dwarf Gl15A, combining five years of intensive RV monitoring with the HARPSN spectrograph with 15 yr of archival HIRES/Keck RV data. We have carried out an MCMCbased analysis of the RV time series, inclusive of Gaussian Process (GP) approach to the description of stellar activity induced RV variations. Our analysis confirms the Keplerian nature and refines the orbital solution for the 11.44day period super Earth, Gl15A b, reducing its amplitude to 1.68_{−0.18}^{+0.17} m s^{−1} (M sin i = 3.03_{−0.44}^{+0.46} M_{⊕}), and successfully models a longterm trend in the combined RV dataset in terms of a Keplerian orbit with a period around 7600 days and an amplitude of 2.5_{−1.0}^{+1.3} m s^{−1}, corresponding to a superNeptune mass (M sin i = 36_{−18}^{+25} M_{⊕}) planetary companion. We also discuss the present orbital configuration of Gl15A planetary system in terms of the possible outcomes of Lidov–Kozai interactions with the wideseparation companion Gl15B in a suite of detailed numerical simulations. In order to improve the results of the dynamical analysis, we have derived a new orbital solution for the binary system, combining our RV measurements with astrometric data from the WDS catalogue. The eccentric Lidov–Kozai analysis shows the strong influence of Gl15B on the Gl15A planetary system, which can produce orbits compatible with the observed configuration for initial inclinations of the planetary system between 75° and 90°, and can also enhance the eccentricity of the outer planet well above the observed value, even resulting in orbital instability, for inclinations around 0° and 15°−30°. The Gl15A system is the multiplanet system closest to Earth, at 3.56 pc, and hosts the longest period RV subJovian mass planet discovered so far. Its orbital architecture constitutes a very important laboratory for the investigation of formation and orbital evolution scenarios for planetary systems in binary stellar systems.
Key words: techniques: radial velocities / stars: individual: GJ15A / binaries: visual / instrumentation: spectrographs / planets and satellites: detection / planets and satellites: dynamical evolution and stability
Based on observations made with the Italian Telescopio Nazionale Galileo (TNG), operated on the island of La Palma by the INAF – Fundación Galileo Galilei at the Roche de Los Muchachos Observatory of the Instituto de Astrofísica de Canarias (IAC); photometric observations made with the robotic telescope APT2 (within the EXORAP programme) located at Serra La Nave on Mt. Etna.
© ESO 2018
1 Introduction
Extrasolar planetary systems have always shown a huge variety of orbital architectures, ever since the very first exoplanet orbiting a main sequence star was discovered (Mayor & Queloz 1995). In the following two decades the number of known exoplanets has increased to more than 3500^{1}, mainly discovered from transit (e.g. Kepler) and radialvelocity (RV; e.g., HARPS, HARPSN) observations. The detected systems encompass a wide range of masses, radii, and orbital separations, including small mass rocky Earthtwins (e.g. Pepe et al. 2011).
Lowmass M dwarfs are the most common main sequence stars, comprising ~ 70% of the stars in the solar neighbourhood (Henry et al. 2006). Furthermore, M dwarfs have become the most promising ground for the hunt for lowmass, rocky planets (e.g. Dressing & Charbonneau 2013; Sozzetti et al. 2013; AstudilloDefru et al. 2017b), due to their more advantageous mass and radius ratios compared to solartype stars. There is a solid evidence, arising both from HARPS and Kepler observations, that super Earths and Neptunes are commonly found in multiple systems (e.g. Udry et al. 2017; Rowe et al. 2014, and references therein).
The nearby M1 dwarf Gl15A was studied by Howard et al. (2014), who find a shortperiod superEarth orbiting the star: they measure a period of 11.44 d and an amplitude of 2.94 m s^{−1}. They also study the activity signals of the host star, identifying the rotational period of the star to be 44 d, both from the S_{HK} index analysis and also from the precise photometric lightcurve, collected with the automatic photometric telescope (APT) at Fairborn Observatory (Eaton et al. 2003). Gl15A has also a known M3.5 binary companion, Gl15B, identified astrometrically by Lippincott (1972) from a small fragment of its orbit, with a measured orbital separation of 146 AU, and an orbital period of 2600 yr.
In this framework we present the clear detection of a longperiod eccentric superNeptune planet around Gl15A. We have used five years of highprecision Doppler monitoring with the High Accuracy Radial velocity Planet Searcher for the Northern emisphere (HARPSN) highresolution (resolving power R ~ 115 000) optical echelle spectrograph (Cosentino et al. 2012) at the Telescopio Nazionale Galileo (TNG), combined with 15 yr of archival RVs data from the Lick–Carnegie Exoplanet Survey (LCES) HIRES/Keck Precision Radial Velocity survey (Butler et al. 2017). The HARPSN data were collected as part of the Harpsn red Dwarf Exoplanet Survey (HADES) programme, a collaboration between the Italian Global Architecture of Planetary Systems (GAPS; Covino et al. 2013; Desidera et al. 2013; Poretti et al. 2016) Consortium^{2}, the Institut de Ciències de l’Espai de Catalunya (ICE) and the Instituto de Astrofísica de Canarias (IAC). We also confirm the presence and update the amplitude of the Gl15A b RV signal, significantly reducing its minimum mass, while confirming the orbital period measured by Howard et al. (2014). Our findings are also discussed in the light of the recent nondetection of Gl15A b in the CARMENES visual RV time series by Trifonov et al. (2018).
With tens of exoplanets orbiting binary stars observed to date (e.g. Eggenberger 2010), including a confirmed wide binary system with planetary companions orbiting both components (Desidera et al. 2014; Damasso et al. 2015), it is debated how the presence of stellar companions has influenced the distribution of planetary orbital parameters. For this reason, we derived new orbital parameters for the stellar companion Gl15B, using HARPSN RV measurements along with astrometric measurements from the WDS (Washington Double Star; Mason et al. 2001) catalogue, and performed several numerical simulations to test the dynamical influence of Gl15B on the planetary system.
We describe the Doppler and photometric measurements collected for the analysis in Sect. 2, and then in Sect. 3 we describe the host star and binary companion updated properties. The complete analysis of the RVs and activity indices of the system is presented in Sect. 4. We then analyse the binary orbit and the perturbations produced on the two planets orbits in Sect. 5. We summarize and discuss our findings in Sect. 6.
2 Observations and HIRES catalogue data
As part of the HADES RV programme, Gl15A has been observed from BJD = 2456166.7 (27th August 2012) to BJD = 2457772.4 (18th January 2017). The total number of data points acquired was 115 over a time span of 1605 days. The HARPSN spectra were obtained using an exposure time of 15 minutes, and achieving an average signaltonoise ratio (S/N) of 150 at 5500 Å. Of the 115 epochs, 49 were obtained within the GAPS time and 67 within the Spanish time. Since the simultaneous ThAr calibration could contaminate the CaII H&K lines, which are crucial in the analysis of stellar activity of M dwarfs (Forveille et al. 2009; Lovis et al. 2011a), the observations were gathered without it. However, to correct for instrumental drift during the night, we used other GAPS targets spectra, gathered by the Italian team during the same nights as the Gl15A observations using the simultaneous ThAr calibration.
The data reduction and RV extraction were performed using the TERRA pipeline (TemplateEnhanced Radial velocity Reanalysis Application; AngladaEscudé & Butler 2012), which is considered to be more accurate when applied to Mdwarfs, with respectto the HARPSN Data Reduction Software (DRS; Lovis & Pepe 2007). For a more thorough discussion of the DRS and TERRA performances on the HADES targets, see Perger et al. (2017). The rms of the TERRA RVs is 2.69 m s^{−1}, while the mean internal error is 0.62 m s^{−1}. The TERRA pipeline also corrected the RV data for the perspective acceleration of Gl15A, dv_{r} ∕t = 0.69 m s^{−1} yr^{−1}.
We also used the HIRESKeck binned data for Gl15A in our analysis, downloaded from the LCES HIRES/Keck Precision Radial Velocity Exoplanet Survey. The HIRES time series spans 6541 days, from BJD = 2450461.8 (13th January 1997) to BJD = 2457002.7 (11th December 2014). These data were newly reduced with respect of the original RVs used by Howard et al. (2014), and also new data have been observed with the HIRES spectrograph since the paper came out. For a full description of the catalogue and data reduction, see Butler et al. (2017). We discarded the last data point of the HIRES time series (BJD = 2457002.7), since it is almost ~10 m s^{−1} off with respect of the rest of the data. In our analysis we thus used 169 HIRES RVs, showing a variation of 3.26 m s^{−1} and an average internal error of 0.84 m s^{−1}. The RV dataset from the Keck archive was already corrected for the perspective acceleration of the host star.
We combined the HARPSN and HIRES RV datasets, obtaining a 290 points time series, spanning 7310 days. The complete RV time series has rms = 3.08 m s^{−1}, mean error = 0.71 m s^{−1}. Figure 1 shows the combined RV time series, with the respective mean RV subtracted from each dataset to visually compensate for the offset expected between the measurements of the two instruments.
As a part of the analysis of the binary orbit described in Sect. 5, we obtained 5 RV data points of the companion star Gl15B with HARPSN, collected from BJD = 2457753.9 (31 December 2016) to BJD = 2457771.9 (18 January 2017). The rms of the time series is 1.26 m s^{−1}, while the mean internal error is 2.42 m s^{−1}.
As for all targets in the HADES survey, Gl15A was also monitored photometrically by the EXORAP (EXOplanetary systems Robotic APT2 Photometry) programme, to estimate the stellar rotation from the periodic modulation in the differential light curve. The observations were performed by INAFCatania Astrophysical Observatory at Serra la Nave on Mt. Etna (+14.973°E, +37.692°N, 1725 m a.s.l.) with the APT2 telescope, which is an 80 cm f/8 Ritchey–Chretien robotic telescope. The detector is a 2 k × 2 k e2v CCD 230–42, operated with standard JohnsonCousins UBV RI filters. The IDL routine aper.pro was used to implement the aperture photometry. The data were collected from the 13th of August 2013 to the 15th of June 2017, for a total 242 and 233 data points, for B and V photometry, respectively.
Fig. 1 Combined HARPSN and HIRES RV time series. As described in Sect. 4.3, we considered the HIRES data as splitinto two separate datasets for the purpose of the RV analysis. 

Open with DEXTER 
3 Stellar properties of Gl15A and Gl15B
The star Gl15A is a high proper motion nearby (π = 280.74 ± 0.31 mas) earlyM dwarf, of type M1. We used the stellar parameters published by Maldonado et al. (2017), which were calculated applying the empirical relations by Maldonado et al. (2015) on the same HARPSN spectra from which we derived the RV time series. This technique calculates stellar temperatures from ratios of pseudoequivalent widths of spectral features, and calibrate the metallicity on combinations and ratios of different features. Although such techniques are mainly used for solartype stars, Maldonado et al. (2015) proved them to be equally effective for lowmass stars.
Howard et al. (2014) derived a rotational period of 43.82 ± 0.56 d, both from their KeckHIRES measurements of the S_{HK} index and their APT photometric observations at Fairborn Observatory. Recently, Suárez Mascareño et al. (2018) analysed the potential signatures of magnetic activity in the CaII H&K and Hα activity indicators of the HADES Mdwarfs sample. For Gl15A they computed a mean level of chromospheric emission , and a rotation period P_{rot} = 45.0 ± 4.4 d, fully consistent with the value found by Howard et al. (2014). The definition of can be extended for application on Mdwarfs spectra (Suárez Mascareño et al. 2018, and references therein).
The stellar parameters for Gl15A used in this work are summarized in the left column of Table 1. We can see that most are fully consistent with the values used by Howard et al. (2014).
The stellar companion of Gl15A, Gl15B, is a type M3.5 dwarf whose orbit was measured astrometrically by Lippincott (1972). For the purpose of our orbital analysis in Sect. 5, we took five HARPSN spectra of Gl15B during the last observing season at TNG (see next section) and we applied on them the Maldonado et al. (2015) techniques to calculate updated stellar properties. The derived values are listed in the right column of Table 1.
Stellar parameters for the two stars Gl15A and Gl15B.
Photometric analysis
In order to identify the potential modulation in the stellar photometry due to the presence of active regions, we used the GLS periodogram (Zechmeister & Kürster 2009) to analyse the B and V time series collected within the EXORAP survey, composed of 242 and 233 points, respectively, taken over five seasons from 13 August 2013 to 15 June 2017. No evident periodicity is found in either time series, in contrast with the findings of Howard et al. (2014) who found a clear signal at 43.82 days (see their Fig. 4), identified as the rotation period of the star, with a corresponding signal seen in the CaII activity index. The discrepancy between the two analysis may be due to the different photometricprecision of the observations: the amplitude of the signal found by Howard et al. (2014) was only of 6 mmag, which is below the sensitivity of the EXORAP survey for targets in the magnitude range of Gl15A.
4 Spectroscopic data analysis
In Fig. 1 an evident longperiod signal can be seen. Howard et al. (2014) identified in their HIRES data a hint of a long period decreasing trend, which they included in their model as a constant negative acceleration of − 0.26 ± 0.09 m s^{−1} yr^{−1}, and concluded it to be consistent with the orbit of the Gl15AB system as calculated by Lippincott (1972). But considering also the new HARPSN data, it becomes clear how the longterm RV variation cannot be modelled as a linear trend.
Even if it could no longer be treated as a constant acceleration, the longperiod signal could still be due to the binary reflex motion of Gl15A, so we investigate this possibility. To model the possible RV signal due to the stellar companion, we follow the procedure used by Kipping et al. (2011) to study the presence of a longperiod companion in the HATP31 system: they modelled the longperiod signal as a quadratic trend, and then derived a range of orbital parameters from the quadratic coefficients. The ratio between the semiamplitude, K_{B}, and period, P_{B}, of the companion signal is given by the secondorder term of the trend, , as (their Eq. (4)) (1)
and, since K_{B} depends on the orbital period and on the mass of the two stars, the orbital period can be derived as a function of the secondorder term and the masses .
From the fit of the complete RV time series, we obtain m s^{−1} day^{−2}. Using the stellar masses listed in Table 1, the resulting period is P_{B} ≃ 680 yr, which would correspond to a semimajor axis a_{B} ≃ 63 AU. This solution is clearly unphysical, since the presumed semimajor axis is less than half the observed orbital separation of the two objects.
Moreover, several highcontrast imaging surveys ruled out the presence of additional comoving stellar objects close to Gl15A: van Buren et al. (1998) excluded the possibility of objects with M_{⋆} > 0.084M_{⊙} at separations of 9−36 AU, while Tanner et al. (2010) ruled out objects up to a magnitude contrast of ΔK_{s} ≃ 6.95 mag within 1^{″} (3.57 AU) and ΔK_{s} ≃ 10.24 mag within 5^{″} (17.8 AU).
The longperiod signal observed in the RV time series is therefore very unlikely to be caused either by Gl15B or by additional stellar companions. Instead, a longperiod planetarymass companion orbiting around Gl15A could be a possible explanation of the observed signal. We now investigate this hypothesis, by analysing both the potential presence of a Keplerian signal in the RVs and the stellar activity signals in the chromospheric indicators, which could also cause longperiod variations due to the star magnetic cycle.
4.1 The MCMC model
Bearing in mind that a signal tightly linked to the stellar rotation period is clearly present in the RV data, as shown by Howard et al. (2014), we have selected the Gaussian process (GP) regression as a useful, and physically robust, tool both to investigate the presence of periodicities in the chromospheric activity indicators and to mitigate the stellar activity contribution to the RV variability. The GP regression is becoming a commonly used method to suppress the stellar activity correlated “noise” in RV time series (e.g. Dumusque et al. 2017, and references therein), especially when adopting a quasiperiodic covariance function. This function is described by four parameters, called hyperparameters, and it can model some of the physical phenomena underlying the stellar noise through a simple but efficient analytical representation. The quasiperiodic kernel is described by the covariance matrix: (2)
where t and t′ indicate twodifferent epochs.
This kernel is composed of a periodic and an exponential decay term, so that it can model a recurrent signal linked to stellar rotation, also taking into account the finitelifetime of the active regions. Such approach is therefore particularly suitable when modelling signals of shortterm timescales, as those modulated by the stellar rotation period.
About the covariance function hyperparameters, h is the amplitude of the correlations; θ represents the rotation period of the star; w is the length scale of the periodic component, linked to the size evolution of the active regions; and λ is the correlation decay timescale, that can be related to the active regions lifetime. In Eq. (2), σ_{instr, data}(t) is the data internal error at time t for each instrument; σ_{instr, jit} is the additional uncorrelated “jitter” term, one for each instrument, that we added in quadrature to the internal errors in the analysis of the RV datasets to take into account additional instrumental effects and noise sources neither included in σ_{instr, data}(t) nor modelled by the quasiperiodic kernel; is the Kronecker delta function.
In the GP framework, the loglikelihood function to be maximized by the Markov chain Monte Carlo (MCMC) procedure is (3)
where n is the number of the data points, K is the covariance matrix built from the covariance function in Eq. (2), and is the data vector.
We use the publicly available emcee algorithm (ForemanMackey et al. 2013) to perform the MCMC analysis, and the publicly available GEORGE Python library to perform the GP fitting within the MCMC framework (Ambikasaran et al. 2015). We used 150 random walkers to sample the parameter space. The posterior distributions have been derived after applying a burnin as explained in Eastman et al. (2013, and references therein). To evaluate the convergence of the different MCMC analyses, we calculated the integrated correlation time for each of the parameters, and ran the code a number of steps around 100–1000 times the autocorrelation times of all the parameters, depending on the complexity of each analysis. This ensured that the emcee walkers thoroughly sample the parameter space (ForemanMackey et al. 2013).
4.2 Analysis of the activity indexes
We first investigated both the HIRES and HARPSN CaII H&K and Hα index time series, in order to test the potential stellar origin of the long period variation seen in the combined RV time series shown in Fig. 1. No longperiod trend is found in either the HIRES or HARPSN datasets. Suárez Mascareño et al. (2018) found a magneticcycle type periodicity at P_{cycle} = 2.8 ± 0.5 yr in the HARPSN CaII H&K and Hα indexes, even if we find no similar signal in the respective HIRES time series. Nevertheless, the period of this cycle is far from the time span of the longperiod signal we investigate, so it could not be the origin of it. Another clue of the stellar origin of thelongperiod modulation of the RVs could be the correlation between the activity indexes and RV time series, which we computed via the Spearman’s rank correlation coefficients. No significant correlation was identified ( for all the indexes). These facts suggest that the longperiod signal is not due to the star magnetic cycle or chromospheric activity, and reinforces the hypothesis of it be due to a wideorbit planetary companion.
To further test the effect of the activity of Gl15A, we then investigated the stellar activity behaviour over the four seasons covered by HARPSN observations by analysing chromospheric activity indicators based on the CaII H&K, and Hα spectroscopic lines. They were extracted using the definition of the line cores and of the reference intervals given in Gomes da Silva et al. (2011). We analysed measurements obtained only from HARPSN spectra because they represent an homogeneous dataset. By spanning more than 1600 days, the HARPSN data alone can provide robust insights into the long and shortterm stellar activity variability. The average S/N was 18 for the CaII H&K index and 213 for the Hα.
We performed an analysis of the activity indicators based on a GP regression, as detailed in the previous section. By adopting the same covariance function (Eq. (2)), we used to model the stellar contribution present in the RV variations in the following section, our primary goal is to investigate some properties of the stellar activity and use them to constrain some parameter priors in the analysis of the radial velocities. This represents a reasonable expectation, because neither the activity indicators nor the RVs have been prewhitened, and the variability patterns in the former could be present with similar properties in the latter, as was noticed by Affer et al. (2016) during the analysis of another HADES target. In this sense, results from the analysis of the activity indicators can be used to train the GP regression of the RVs, by keeping unchanged the way the stellar activity contribution is modelled. Here we have used Eq. (2) to describe the variability correlated with the stellar rotation period P_{rot}, by adopting a uniform prior for the corresponding hyperparameter θ which is constrained between 40 and 60 days (the list of priors is shown in Table 2).
The MCMC analysis was stopped after running for around 1000 times the highest autocorrelation time. The posterior distributions of the model parameters are shown in Fig. 2, and the bestfit estimates are listed in Table 2. We note that the estimates of the stellar rotation periods are well constrained but slightly different for the two activity indicators. The value found from the CaII H&K index is very similar to that found by Howard et al. (2014) for the CaII H&K Sindex derived from HIRES spectra. Thehyperparameter θ found in theanalyses of the two time series is in very good agreement with the rotation period found by Suárez Mascareño et al. (2018) with a GLS analysis of the activity indexes time series of Gl15A.
Priors and bestfit results for the Gaussian process regression analysis of the chromospheric activity indicators extracted from the HARPSN spectra of Gl15A.
Fig. 2 Posterior distributions of the fitted (hyper)parameters of the GP quasiperiodic model applied to the time series of CaII H&K (upper panel) and Hα (lower panel) activity indexes. The vertical lines denote the median (solid) and the 16th–84th percentiles (dashed). 

Open with DEXTER 
4.3 Analysis of the combined HIRES and HARPSN RV time series
First, we studied the HARPSN RV time series for confirmation of the presence of the Gl15A b signal found by Howard et al. (2014). We performed a GLS analysis to identify any periodicity in our data in P ∈[2, 100] d. In the periodogram, shown in Fig. 3, we can see that the P = 11.44 d period of planet Gl15A b is clearly recovered. As a confirmation of the coherence of the signal throughout the entire HARPSN time series, Fig. 4a shows the evolution of the period, amplitude and phase recovered by the GLS periodogram as a function of the number of observations: we can clearly see how from ~ 70 forward the period and phase of the signal remain stable, with small oscillations of the order of the final error on the parameters, even if the remainder of the observations covers almost one and a half years. We also studied the periodogram of the CaII H&K and Hα HARPSN time series, and found no significant signal at periods close to the inner planet period P_{b}, thus confirming its planetary nature as announced by Howard et al. (2014).
We note that, as can be seen in Fig. 3, the final amplitude recovered by the GLS periodogram on the HARPSN dataset, K_{b} = 1.82 ± 0.31 m s^{−1}, is smaller than the one published by Howard et al. (2014), K_{b,How} = 2.94 ± 0.28 m s^{−1}. The decreasing behaviour of the amplitude recovered by GLS with increasing number of observation is not unexpected, as the sampling of the signal strongly influence the periodogram structure and fit, and we tested that a similar behaviour can be observed also in simulated datasets with the same epochs as our HARPSN time series but containing only the planetary signal and white noise, as shown in Fig. 4b. The lower amplitude value than that found by Howard et al. (2014) can be explained similarly, since the HARPSN time series is composed of five years of continuous highcadence observations of the target, while the dataset used by Howard et al. (2014) consisted in mostly sparse measurements with an intensive highcadence monitoring only in the last seasons, which could affect the signal recovery. Moreover, it is worth noticing that the HIRES data from the now public LCES HIRES/Keck archive have been reduced with a different and more effective technique (Butler et al. 2017) than the data used by Howard et al. (2014). Performing a GLS analysis of the two time series we obtain the same peak periodicity but an amplitude K_{b} ~ 23% smaller in the new archive data. Since no other shortperiod signal emerges from the GLS analysis of the RVs, we focussed our attention on the study of the longperiod signal, which we assume to be due to a planetarymass wideorbit companion.
To estimate the orbital and physical parameters of the known planet Gl15A b and the candidate companion, hereafter Gl15A c, we have performed an MCMC analysis of the RVs. Following the prescription of Howard et al. (2014), we modelled the RVs dividing the HIRES dataset in a preupgrade and a postupgrade sublist, due to the HIRES CCD upgrade occurred in August 2004. Each subsample was then treated as an independent dataset, with its own zeropoint offset and additional jitter term. In this case, the in Eq. (3) represents the RV residuals, obtained by subtracting the Keplerian signal(s) from the original RV dataset.
The general form for the models that we tested in this work is given by the following equation: (4)
where n_{planet} = 1, 2; ν is a function of time t, time of the inferior conjuntion T_{0,j}, orbital period P_{j}, eccentricity e and argument of periastron ω_{j}; γ_{instr} is the RV offset, one for each instrument; GP is the stellar noise modelled with the Gaussian Process. The term is the residual acceleration of the system, which can include also the RV acceleration due to the orbital motion within the Gl15AB binary system, which is difficult to predict due to the large uncertainties in the orbital parameters (as we will discuss in Sect. 5.1). Instead of fitting separately e_{j} and ω_{j}, we use the auxiliary parameters and to reduce the covariance between e_{j} and ω_{j}. Howard et al. (2014) find the eccentricity of planet Gl15A b to be consistent with zero, and concluded a circular orbit to be the best fit for the data. To confirm this conclusion or to evaluate the real value of the eccentricity e_{b}, we assumed in our analysis a Keplerian orbit for both planets. The eventual eccentricity orbit for Gl15A b would also be of interest as a potential consequence of the dynamical interaction of planet b with Gl15B, as we will discuss in Sect. 5.
With the exception of very few parameters, we assumed uniform priors in our analysis. Our choice for the range of λ is justified by the results obtained from the GP analysis of the CaII and Hα spectroscopic activity indexes (see Table 2), and from a preliminary, quick MCMC test on the data, which showed that the chains were converging towards values not far from the expected stellar rotation period. For the semiamplitudes K of the Doppler signals, we used the modified invariant scale prior: (5)
following the prescription in Gregory (2010). For the orbital period of Gl15A b there was a strong constraint from the results of Howard et al. (2014), confirmed by our GLS analysis previously discussed. Nevertheless, since the same dataset from Howard et al. (2014) is used in our analysis, it would be improper to use this results as an informative prior. We instead adopted an uninformative prior over a small period interval centred on the value found by Howard et al. (2014), e.g. P_{b} ∈ [10.9, 12.0]. For the orbital period of the candidate planet Gl15A c, we used a logarithmic prior to assume all the orders of magnitude equally likely.
The MCMC analysis was stopped after running for around 300 times the highest autocorrelation time. Table 3summarizesthe bestfit values for each of the 19 free parameters of our model, together with the details of the prior distributions used to draw the samples. The posterior distributions of the model parameters are shown in Fig. 5.
About the bestfit values for the free parameters, we first note that the semiamplitude of the Doppler signal induced by Gl15A b is even lower than the estimate obtained from the GLS periodogram, differing from the estimate of Howard et al. (2014) (2.94 ± 0.28 m s^{−1}) for more than 3.5σ, resulting in a lower minimum mass. In addition to the previously discussed effect, this is a consequence of our global model, where planetary signals are fitted jointly with a term describing the stellar correlated noise, and also explains the lower values of σ_{HIRES, jit} required by our fit.
Our dataset does not cover a complete orbit of the outer planet Gl15A c, then we cannot reliably constrain the eccentricity, which appears to be significant within less than 1.5σ (e < 0.40 at a 68% level of confidence). Our estimate for the minimum mass of Gl15A c places this companion in the group of the socalled superNeptunes (as planets in the [20, 80]M_{⊕} mass range are generally referred to).
The eccentricity of Gl15A b also resulted consistent with zero within 1.5σ, with a much lower value (e < 0.13 at a 68% level of confidence). Our analysis thus agrees with the results by Howard et al. (2014), who concluded that there was no evidence for an eccentric orbit to be preferable to a circular one.
Figure 6 shows the RV curves folded at the bestfit orbital periods for the known planet Gl15A b and the detected outer companion Gl15A c. In Fig. 7 we show the RV residuals, after the two Keplerians have been removed, with superposed the bestfit stellar correlated, quasiperiodic signal.
As for the GP hyperparameters, the stellar rotation period θ assumes a value which is in agreement with that expected, and the higher uncertainties than those for the activity indicators (Table 2) should be due to the longer time span covered by the RVs, which include HIRES data. The longer timespan should also explain why the activeregion evolutionary timescale λ sets on a value less than those found for the activity indicators. Without any ancillary data available as photometry, covering the same time span of the RVs, we cannot conclude if the value for λ is actually physically robust, but a value not far from the rotation period seems not unreliable, as already observed in studies of different systems (e.g. Affer et al. 2016). Also, the h hyperparameter, corresponding to the mean amplitude of the stellar signal, is fully compatible with the expected value derived by Suárez Mascareño et al. (2018) from the mean activity level of the star, , that is K_{exp} = 1.9 ± 0.4 m s^{−1}.
Fig. 3 GLS periodogram of the complete HARPSN RV dataset. 

Open with DEXTER 
Fig. 4 Panel a: evolution of the GLS orbital parameters as a function of the number of HARPSN RV measurements; panel b: evolution of the GLS orbital parameters as a function of the number of simulated data. The error bars on the final points are shown as references. 

Open with DEXTER 
Planetary parameter bestfit values obtained through a joint modelling of Keplerian signals and correlated stellar noise, using Gaussian process regression.
Fig. 5 Posterior distributions of the fitted (hyper)parameters of the twoplanet model, where the stellar correlated noise has been modelled with a GP regression using a quasiperiodic kernel. On the yaxis is shown the logarithm of the product between the likelihood and the prior. The vertical lines denote the median (solid) and the 16th–84th percentiles (dashed). 

Open with DEXTER 
4.4 Comparison with CARMENES results
In a very recent work, Trifonov et al. (2018) publish an RV analysis at optical wavelengths of the Gl15A system as part of the CARMENES survey. The RVs derived from the visible arm of CARMENES show no evidence of the presence of Gl15A b, thus the authors conclude it to be an artefact of the stellar noise. The authors also identify a longterm trend in the combined Keck+CARMENES dataset, which they propose to be due to a distant lowmass companion. As discussed in Sect. 4.3, our HARPSN data alone clearly confirm the presence of the 11.44 d period due to Gl15A b, albeit with a reduced amplitude and mass, and the analysis of the combined HARPSN+Keck datasets makes a decisive case for the existence of the longperiod lowmass giant Gl15A c based on a fiveyear time span of HARPSN observations, partly overlapping with the Keck timeseries. The orbital elements for Gl15A c derived in this work have larger formal uncertainties than those reported in Trifonov et al. (2018) due to the likely much longer orbital period inferred for the planet, but the inferred companion mass is fully in agreement with the preliminary CARMENES estimate.
In Fig. 8a we can see the last season of HARPSN observations, that overlaps with most of the CARMENES data, compared toour twoKeplerian plus correlated stellar noise model. This season is clearly dominated by activityrelated variations in the optical HARPSN spectra, while the internal errors for HARPSN are typically three times smaller than those of CARMENES (0.6 m s^{−1} vs. 1.8 m s^{−1}). We nonetheless tested the effect of combining the HARPSN and CARMENES datasets: the GLS periodogram, shown in Fig. 8b, clearly peaks on the 11.44 d period of Gl15A b, and remarkably resembles Fig. 3, recovering the same amplitude as that recovered on the HARPSN dataset alone.
We also tested our twoplanet + stellar noise model described in Sect. 4.3 on the complete dataset obtained combining the HIRES, HARPSN, and CARMENES time series (running the MCMC code for around 100 times the maximum autocorrelation time), and we obtained values of the system parameters entirely in line with those presented in Table 3. In particular, the recovered Doppler semiamplitude and minimum mass for the inner planet are m s^{−1} and , respectively. These values are well within 1σ of the results from the analysis of only the HIRES and HARPSN datasets: no significant drop in the amplitude of Gl15A b signal is observed. We also note that the residual jitter found in the analysis is small ( m s^{−1}) signifying that again the GP is able to model the stellar noise for this third dataset as well. The GP parameters are almost unchanged with respect to the values presented in Table 3 ( m s^{−1}, d, and d), meaning that the activity signal description derived from the first two datasets is robust also in modelling the CARMENES data. Given the intrinsically higher quality of the HARPSN data that drive the GP regression modelling and in which the coherence of the signal from the inner planet Gl15A b is clearly present, we decide to stick to the results in Table 3 for the purpose of the dynamical interaction analysis described in Sect. 5.2.
Trifonov et al. (2018) also stated, as an argument in favour of the nonKeplerian interpretation of Gl15A b, that the 11.44 d signal disappeared when analysing the last two years of HIRES observations, subsequent to the time series used by Howard et al. (2014). Studying the same time span of data, we observed that, when ignoring the outlier described in Sect. 2, the 11.44 days signal is still clearly visible in the GLS periodogram.
Fig. 6 Phasefolded RV curves for panel a Gl15A b and panel b Gl15A c. Each curve shows the residuals after the subtraction of the other planet and the stellar correlated signal. The red curve represents the bestfit Keplerian orbit, while the red dots and error bars represent the binned averages and standard deviations of the RVs. 

Open with DEXTER 
Fig. 7 Upper panel: best fit stellar quasiperiodic signal (blue line) compared to the RV residuals. Lower panels: magnification of the highcadence HIRES/KECK observations (left panel) and of the last HARPSN observing season (right panel). 

Open with DEXTER 
Fig. 8 Panel a: comparison between the last season HARPSN data (black triangles) and the overlapping CARMENES data (grey dots). The red solid line represents our twoplanet plus stellar noise model, with the pink shaded area representing the 1σ uncertainties on the GP hyperparameters; panel b: GLS periodogram of the combined HARPSN and CARMENES RV datasets. 

Open with DEXTER 
5 Binary orbital interaction
The orbital eccentricity of the outer planet, Gl15A c, is quite uncertain: as shown in Fig. 5, the posterior distribution has no clear peak, only constraining the eccentricity towards small values, with a 68% probability to be <0.40 and 95% probability to be <0.72. This would normally point towards the adoption of a circular orbit as best fit solution for the system, lacking a significant evidence of eccentricity, but this would be to reckon without the stellar companion. The presence of Gl15B cannot be ignored, especially when studying a wide orbit such as that of Gl15A c.
We thus investigated how the dynamical interaction with the companion star could influence the Gl15A system. One of the main mechanisms to excite exoplanets eccentricities is the Lidov–Kozai effect, in which the presence of an external perturber causes oscillations of the eccentricity, e, and the inclination, i, with the same period but opposite phase. It was originally studied by Lidov (1962) and Kozai (1962) to compute the orbits of high inclination small solar system bodies, like asteroids and artificial satellites, and it is strongly dependent on the eccentricities and mutual inclination of the involved objects. Thus, in order to better estimate the strength of Lidov–Kozai oscillations in the Gl15A planetary system, we need to understand as precisely as possible the orbit of the stellar companion Gl15B.
5.1 Orbital modelling from astrometry and RV data
The first obstacle in the dynamical analysis of the system was the poorly constrained orbital parameters of the companion. Lippincott (1972) estimated a period of 2600 yr from ~ 100 yr of astrometric measurements, which coversless than 4% of the orbit.
Dealing with a similar case of an eccentric planet hosted by a wide binary system, Hauser & Marcy (1999) developed a technique for constraining longperiod binary orbital parameters, combining astrometric and radial velocity measurements. The method is based on the fact that, being Newtonianmechanics deterministic, knowing the instantaneous full position and velocity vector, [x, y, z, V_{x}, V_{y}, V_{z}], of one mass with respect to the other, you can compute the exact orbit of the system. Therefore, even by observing a small fragment of the orbit, we should be able to gather all the orbital parameters of the stellar companion. Of course astrometry, which is restrained in the plane of the sky, cannot provide the complete 3D information needed for this analysis, so additional data from radial velocity, to compute the third component of the velocity vector, and parallax distance, to convert astrometric positions in Cartesian coordinates, is necessary.
Following the procedure of Hauser & Marcy (1999), we downloaded 122 astrometric observations of the Gl15 system from the WDS catalogue, spanning from 1860 to 2015. The variations of the position angle θ of Gl15B relative to Gl15A and the angular separation ρ are shown in Fig. 9. To derive θ and ρ at a specific time, along with their derivatives dθ∕dt and dρ∕dt, which we need to calculate the velocity component in the plane of the sky, we fitted the data with a secondorder polynomial:
From these we can easily derive dθ_{fit}∕dt and dρ_{fit}∕dt as
Adopting the parallax value from Table 1, π_{P} = 280.74 ± 0.31 mas, we can derive the Cartesian position and velocity using their Eqs. (1)–(4):
Thus, we gained 4 of the desired physical components, [x, y, V_{x}, V_{y}], as a function of time t, through θ_{fit} and ρ_{fit} (Eqs. (6) and (7)).
To obtain the third component of the velocity vector, V_{z}, we used the combined Doppler information of the two stars: from the same HARPSN spectra we used to obtain the RV time series for the planet detection, collected as illustrated in Sect. 2, we extracted the absolute radial velocities with the DRS pipeline. We used the DRS pipeline instead of TERRA since the latter only produce relative RV, which cannot be used in the comparison of two different objects. For Gl15A we take all the 115 HARPSN epochs, and subtract from the absolute RV the planetary and stellar signals, as derived in Sect. 4.3. For Gl15B we used the five spectra taken on January 2017, as described in Sect. 2. The two datasets are illustrated in Fig. 10. To derive the binary orbit, we need to know all the position and velocity components at the same instant. Therefore, we fitted the two RV time series with firstorder polynomials. From the difference of the two linear fits, we obtained the relative lineofsight velocity V_{z} of Gl15B with respect to Gl15A. We were then able to then select an epoch, and compute the values of [x, y, V_{x}, V_{y}, V_{z}] for that time. The selected epoch to derive the binary orbit is BJD = 2457754.5, that is 1st January 2017, which is well inbetween all the datasets, and close to the Gl15B RV time series, which is the shortest and most uncertain.
The last missing piece of the puzzle is, of course, the lineofsight separation z, which cannot be measured, but can be constrained imposing the condition of a bound orbit for the binary system. This is expected due to the similar spectral type and proximity of the two stars. The condition of bound orbit translate into a condition on the total energy of the system: (14)
where μ = G(M_{A} + M_{B}), and, of course, .
From this we get a range of acceptable z values − 400 ≲ z ≲ 400 AU. From every value of z is possible to compute all the orbital parameters [P, e, i, ω, Ω, T_{P}] for the binary system (see Hauser & Marcy 1999, Eqs. (7)–(15)). The results are shown in Fig. 11. As we can see, there is a wide variety of possible orbits, with completely different eccentricities and orientations, even with the bound orbit constraint, and this is still insufficient for any meaningful dynamical analysis.
The procedure by Hauser & Marcy (1999) provides the orbital solution for every single value of z, but does not in any way distinguish between the more probable configurations. But those orbital configurations are not all equally likely: from theory and observations of binary systems we know the expected distributions for different orbital parameters. From this information we can extract some priors to help us identify the most probable orbital configuration for the system, that is the best fit value of the lineofsight separation z.
To do this we performed a Monte Carlo simulation in which the priors on the orbital parameters are injected via rejection sampling. Not all the a priori distributions of the orbital parameters are known, so we restricted the prior selection to the parameters that have a strong impact on the outcome and for which a good information is available. Due to the central role played by the eccentricity of the perturber in the Lidov–Kozai perturbation, we applied a prior on e, to have the best possible constraint on it. We selected the eccentricity distribution from Tokovinin & Kiyaeva (2016), who studieda sample of 477 wide binaries within 67 pc from the Sun with median projected separation of ~ 120 AU, very close to that of the Gl15AB system. The aforementioned prior is (15)
The second choice is a prior on the orbital period, in order to penalize long period poorly bound orbits. The chosen prior is the one suggested by Duchêne & Kraus (2013) for lowmass binary stars, that is a lognormal distribution, with ā ≈ 5.2 AU and σ_{logP} ≈ 1.3.
The results of the Monte Carlo simulation are illustrated in Fig. 12. As we can see, there are three main peaks in the distribution. The third one, on the far right, represents orbits that are almost unbound, so it can be ignored, both because the two stars are expected to be in a binary system, and because even if bound, such wide orbits would have no influence whatsoever on the dynamics of the planetary system we intendedto study. The latter can also be said of the peak on the left, which corresponds to a period of yr, and a semimajor axis of AU, again too distant from the planetary system to have a significant influence.
We thus used the central peak of Fig. 12 to derive the orbital parameters best solutions and error bars. To do this we fitted a truncated normal distribution in the range z ∈ [−200, 200] AU, and use the mean value μ_{z} to calculate the corresponding best solution orbital parameters as described before. The upper and lower errors on the orbital parameters were calculated by taking μ_{z} + σ_{z} and μ_{z} − σ_{z} and deriving the corresponding values of [P, e, i, ω, Ω, T_{P}]. The orbital parameters solutions and errors are listed in Table 4.
Assuming the best fit orbital parameters listed in Table 4, we were able to calculate the RV signal caused by Gl15B on Gl15A, and compare it with the value of the residual acceleration found by our MCMC analysis. Gl15B RV signal is approximately linear, as is to be expected due to the small fraction of the orbit covered by the time series, and correspond to an acceleration m s^{−1} d^{−1}, which is fully compatible with the result from our MCMC listed in Table 3.
Fig. 9 Position angle (panel a) and angular separation (panel b) of Gl15B with respect to Gl15A. Position angle is measured from north towards east. The red lines in both panels represent the second order polynomial fits. 

Open with DEXTER 
Fig. 10 Radial velocities of (panel a) Gl15A (after subtracting the planetary and stellar signals) and (panel b) Gl15B. The solid lines show the respective first order polynomial fits. 

Open with DEXTER 
Fig. 11 Orbital parameters of Gl15B as a function of the lineofsight separation z with Gl15A. 

Open with DEXTER 
Fig. 12 Distribution of z resulting from the Monte Carlo simulation with the e and P priors injected. 

Open with DEXTER 
Best orbital parameters for the Gl15 binary system from the MC simulation with priors on period and eccentricity.
5.2 Lidov–Kozai Interaction modelling
The results of the previous section show that the most likely orbit of the stellar companion Gl15B has a high eccentricity, e = 0.53. This is a further clue that strong orbital perturbation could affect the planetary system.
The interaction mechanism studied is the Eccentric Lidov–Kozai effect (commonly reffered to with the literaturecoined acronym EKL). This mechanism applies to hierarchical triplebody systems, and consists in eccentricity and inclination oscillations on timescales much larger than the orbital period of the influenced body. We can safely ignore the mutual interaction between the two planets of the Gl15A system, except in the event of close encounters, and thus treat their interaction with the binary separately, as threebody systems.
The EKL mechanism is very sensitive to the mutual orientations of the orbits of the perturber, (the companion star), and of the influenced body (the planet). We have derived i, ω, Ω of Gl15B, but we do not knoweither i or Ω of the two planets Gl15A b and Gl15A c. Thus, some assumptions are to be made about their orbital orientation. To compare the results of the dynamical interaction model with the posterior distribution found by the MCMC analysis in Sect. 4.3, we calculated the fraction of time spent by Gl15A c below e = 0.40 and e = 0.72 (f(e < 0.40) and f(e < 0.72)), which can be considered a proxy of the probabilities (68% and 95%, respectively) to observe the system with eccentricities under those thresholds (e.g. Anderson & Lai 2017). The same can be done to study the interaction of Gl15B with Gl15A b, in which case the 68% and 95% probability levels correspond to e = 0.13 and e = 0.25.
We performed some preliminary tests to verify the various mechanisms which could be involved in the dynamical evolution of the system: we proved that the inner planet, Gl15A b, is too distant from Gl15B for any significant interaction. This is in good agreement with the extremely low level of eccentricity found in our MCMC analysis. We thus focussed our efforts on studying the orbit of the newly discovered Gl15A c; we also checked for the influence of dissipative tides, which could lessen the orbital eccentricity, and find that they act on much longer timescales than the EKL mechanism, so we neglected them in our analysis.
We denote the orbital parameters of Gl15B with the subscript_{B} and the ones of Gl15A c with_{c}. All the EKL integrations were performed with a timescale of ~10 Myr. Since we cannot rule out planet–planet scattering events to have occurred in the early phases of the system’s evolution, we consider different values of the initial eccentricity e_{c,0} = 0.0, 0.1, 0.2, 0.3, 0.4. The other unknown is the longitude of the ascending node of the outer planet Ω_{c}. However, the longitude of the ascending node influences the EKL interaction mainly by changing the relative inclination the two orbits, θ, which derives from the parameters of the two orbits as
and thus can be ignored so long as θ is conveniently sampled, which can be obtained by changing i_{c}. For our analysis we fixed Ω_{c} at 0° and varied the planet orbital inclination in the interval i_{c} ∈ [−5°, 90°].
The results are shown in Fig. 13. In some cases the EKL oscillations are so extreme that the numerical integration has to be stopped due to the planet passing too close to the host star; the corresping systems are clearly unstable, and so we considered f(e < 0.40) = f(e < 0.72) = 0 to represent the incompatibility with the observed case. We also considered as unstable the cases in which the outer planet’s orbit becomes too close to the inner planet’s orbit, possibly producing planet–planet scattering events. In other terms, when (17)
where e_{c,max} is the maximum eccentricity reached during the EKL oscillations. In these cases we also considered f(e < 0.40) = f(e < 0.72) = 0.
As we can see in both the panels of Fig. 13, there are regions in the i_{c} space that are unstable regardless of the initial eccentricity of the planet e_{c,0}: between 15° and 30° and around 0°. Only the solid black line, corresponding to e_{c,0} = 0, behaves differently, showing the system to be stable between 15° and 30° and unstable at larger inclinations. We can also see that the Lidov–Kozai interaction is weak for i_{c} ~ 90° and strengthen as i_{c} decreases. The top panel of Fig. 13 shows that for i_{c} ∈ [75°, 90°] the resulting eccentricity ranges are compatible with the observed values. The constraints for the e < 0.72 threshold are somewhat looser, but pointing in the same direction.
As previously mentioned, we do not know the initial eccentricity of the system, so a more robust way to consider the dynamical evolution is to consider the average of the results of the single integrations. This can be seen in Fig. 14, which confirms the trends discussed above, with initial inclinations around 0° and between 15° and 30° leading to an unstable system, and inclinations higher than 75° producing a good agreement with the observed orbital configuration.
Fig. 13 Fraction of time spent below the e = 0.40 (panel a) and e = 0.72 (panel b) thresholds. The different linestyles and symbols correspond to different values of e_{c,0}: solid and plus signs – e_{c,0} = 0; dotted and asterisks – e_{c,0} = 0.1; dashed and diamonds – e_{c,0} = 0.2; dash dot andtriangles – e_{c,0} = 0.3; dash dot dot and squares e_{c,0} = 0.4. The horizontal dashed grey lines indicate the corresponding probability from the MCMC – panel a: 68%; panel b: 95%). 

Open with DEXTER 
Fig. 14 Fraction of time spent below the e = 0.40 (solid grey) and e = 0.72 (dotted black) thresholds, averaged on the initial eccentricity e_{c,0}. The errorbars show the standard deviation. 

Open with DEXTER 
6 Discussion and conclusions
We present in this paper the fifth planet detected by the HADES programme conducted with HARPSN at TNG. This longperiod planet was found orbiting the planethost M1 star Gl15A, from the analysis of high precision, highresolution RV measurements collected as part of the survey in conjunction with archive RV data from the HIRES/Keck spectrograph.
The different trends observed in the two datasets suggest the presence of a longperiod companion, which is confirmed by the homogeneous Bayesian analysis of the combined RV time series. The known inner planet Gl15A b is also recovered. The minimum masses derived from our analysis are and for the inner and outer planets, respectively. The mass we find for Gl15A b is much smaller than the value found by Howard et al. (2014), which was almost double due to the higher signals amplitude. The smaller value that we find can be easily explained by the additional information brought by the highprecision HARPSN RVs, along with the new calibration of the archival HIRES data published by Butler et al. (2017). The combined dataset is almost twice as large that the one analysed by Howard et al. (2014), stretched on a significantly longer time span, with better sampling and precision. This, together with the simultaneous modelling of the stellar activity signal, can explain the much smaller uncertainty on the minimum mass of Gl15A b. It also highlights the importance of taking into account the chromospheric stellar activity for the correct identification of planetary signals. It is worth noticing that instead the orbital period P_{b} is almost unchanged from the previous estimate. Our fit also places an even lower upper limit to the eccentricity of Gl15A b than that found by Howard et al. (2014;e < 0.13 at a 68% level of confidence), thus confirming their conclusion that the planet’s motion is best described by a circular orbit.
With its period of ≃21 yr, Gl15A c is the longestperiod subJovian planet detected up to date with the RV method^{3}, the second being HD 10180 h (Lovis et al. 2011b) with a period of approximately six years, and a minimum mass of 65.74M_{⊕} (Kane & Gelino 2014). With the confirmed presence of two widely spaced planetary mass companions, Gl15A is now the multiplanet system closest to our Sun, at a distance of only 3.5620 ± 0.0039 pc, slightly smaller than the distance to YZ Cet (AstudilloDefru et al. 2017a), 3.69 ± 0.11 pc.
We also compared the results of our analysis of the HIRES and HARPSN RV data with the recent results by Trifonov et al. (2018) based on CARMENES highcadence monitoring of the target. Trifonov et al. (2018) found no evidence for the presence of planet b, while we showed in Sect. 4.4 that the signal can be clearly detected when combining the HIRES, HARPSN and CARMENES data, without any loss in significance. The reason why CARMENES does not see the signal of Gl15A b is not fully clear, but the higher quality, and much longer time span, of the HARPSN data, combined with our modelling of the stellar activity quasiperiodic signal could be a possible explanations for the nondetection based on the CARMENES data alone. Another factor contributing to this nondetection could be the sampling of the CARMENES RV data, which appears not to be homogeneously distributed over the 11.44 d period (see Fig. 10 from Trifonov et al. 2018).
The CARMENES visual arm contains a spectral region that extends all the way to 0.95μm; in other words, a significantly redder spectral range than that covered by HIRES and HARPSN. As the amplitude of activity induced RV variations is known to be chromatic (Reiners et al. 2010), the nondetection of the 11.44 d signal in the CARMENES time series could be an indication of a wavelength dependentamplitude of the signal, which would clearly indicate its stellar origin. However, the CARMENES visual arm spectral range still significantly overlaps with that of HIRES and HARPSN and thus it must be affected by stellar activity in a similar way. Based on the higher RV precision of HARPSN, allowing a detailed modelling of quasiperiodic stellar signal (as shown in Fig. 7), and on the coherence of the period and phase of the signal over the 20 yr time span covered by the combined HIRES and HARPSN time series, the Keplerian origin of the signal still seems the most straightforward explanation for the observed data. It would, however, be interesting to carry out a systematic study on both CARMENES and HARPSN time series of this target, adopting the same strategy outlined by Feng et al. (2017) for Tau Ceti, that is to study separately the RVs derived from different regions of the spectra, for a sistematic investigation of potential differences in the amplitude of the 11.44 d signal, as a function of the wavelength, but this lies well beyond the scope of this paper.
Dwarf stars are known to turn up much more frequently in multiple systems than they do in isolation, with a binary fraction as high as ≃57% for Sunlike stars (Duquennoy & Mayor 1991) and somewhat lower for M dwarfs (Bergfors et al. 2012). Many young binaries possess either circumstellar or circumbinary disks (e.g. Monin et al. 2007), and the existence of stable planetary orbits in binary systems was postulated well in advance of the first exoplanets discoveries (Dvorak 1982).
Early studies proposed different massperiod relations for planets around binaries and single stars (Zucker & Mazeh 2002) but in the following years the evidence of such diversity decreased (e.g. Desidera & Barbieri 2007; Eggenberger 2010), until most recently Ngo et al. (2017) claimed there not to be any difference of planetary properties between the two kind of systems. On the other hand, recent works such as Moutou et al. (2017) find statistical evidence for a much higher binary fraction in extrasolar systems hosting eccentric exoplanets than in the ones hosting only circular planets: this points towards the confirmation of the role of stellar multiplicity in orbital excitation of planetary systems, as predicted by theoretical studies which suggested a strong orbital influence of stellar companions on planetary systems, via mechanisms such as the eccentric Lidov–Kozai (EKL) oscillations.
Our numerical analysis of the EKL effect proved the strong influence of the Gl15B on the planetary system. We show that for a narrow range of initial inclination, 75°−90°, the outer planet maintains a low eccentricity orbit, regardless of the initial status in which the system was due to possible planet–planet scattering events. We also pointed out the presence of a forbidden ranges of inclination, 15°− 30° and ~ 0°, in which the Lidov–Kozai interaction become so strong that no stable orbit can be achieved, regardless of the initial eccentricity of Gl15A c.
The orbital parameters of Gl15A c have still large uncertainties due to the observation time span shorter than the orbital periods, and the semiamplitude K_{c} is significant only at a 3σ level, although the strong combined observational evidence from RV and imaging leaves no doubt as to the presence of a longperiod planetarymass companion. Additional RV observations in the years to come will, however, be very helpful in better constrainingthe orbit, and thus the mass of Gl15A c.
Our knowledge of this system will be also greatly improved by the results published in future Gaia data releases. For a circular orbit and assuming the minimum mass value for Gl15Ac, the expected astrometric signature on the primary is 570 μas. Gaia astrometry will only cover ~ 20%−25% of the full orbit. However, based on the Torres (1999) formalism curvature effects in the stellar motion should typically amount to 20− 30μas yr^{−2}, thus they should be easily revealed by Gaia, that for such a bright star as Gl15A will be able to deliver endofmission proper motion accuracies ≲10μas yr^{−1} (e.g. Gaia Collaboration 2016b). Even with the orbital solutions now available, our analysis shows how interesting dynamical studies can be performed on the system, which, due to the presence of the eccentric binary companion, is an excellent playground to test orbital interaction mechanisms and their influence on the evolution of planetary systems.
Acknowledgements
GAPS acknowledges support from INAF through the “Progetti Premiali” funding scheme of the Italian Ministry of Education, University, and Research. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/20072013) under Grant Agreement No. 313014 (ETAEARTH). The HARPSN Project is a collaboration between the Astronomical Observatory of the Geneva University (lead), the CfA in Cambridge, the Universities of St. Andrews and Edinburgh, the Queen’s University of Belfast, and the TNGINAF Observatory. M.P. thanks E.T. Russo for the helpful discussions in the final stages of the analysis. J.I.G.H., R.R.L., A.S.M., and B.T.P. acknowledge financial support from the Spanish Ministry project MINECO AYA201456359P. J.I.G.H. also acknowledges financial support from the Spanish MINECO under the 2013 Ramón Cajal programme MINECO RYC201314875. A.S.M also acknowledges financial support from the Swiss National Science Foundation (SNSF). We acknowledge the computing centres of INAF – Osservatorio Astronomico di Trieste/Osservatorio Astrofisico di Catania, under the coordination of the CHIPP project, for the availability of computing resources and support. This research has made use of the Washington Double Star Catalog maintained at the U.S. Naval Observatory. We gratefully acknowledge an anonymous referee for her/his insightful comments that materially improved an earlier version of this manuscript.
Appendix A Observation log for Gl15A
In this section we report the observational data collected with the HARPSN spectrograph as part the HADES project and used in the present study. We list the observation dates (barycentric Julian date or BJD), the radial velocities (RVs), and the Hα and S_{HK} indices, calculated by the TERRA pipeline. The RV values are corrected for perspective acceleration. The RV errors reported are the formal ones, not including the jitter term, while the Hα and S_{HK} errors are due to photon noise through error propagation.
Data of the 115 observed HARPSN spectra of Gl15A.
References
 Affer, L., Micela, G., Damasso, M., et al. 2016, A&A, 593, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ambikasaran, S., ForemanMackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 2 [Google Scholar]
 Anderson, K. R., & Lai, D. 2017, MNRAS, 472, 3692 [NASA ADS] [CrossRef] [Google Scholar]
 AngladaEscudé, G., & Butler, R. P. 2012, ApJS, 200, 15 [NASA ADS] [CrossRef] [Google Scholar]
 AstudilloDefru, N., Díaz, R. F., Bonfils, X., et al. 2017a, A&A, 605, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AstudilloDefru, N., Forveille, T., Bonfils, X., et al. 2017b, A&A, 602, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bergfors, C., Brandner, W., Hippler, S., et al. 2012, in From Interacting Binaries to Exoplanets: Essential Modeling Tools, eds. M. T. Richards & I. Hubeny, IAU Symp., 282, 460 [NASA ADS] [Google Scholar]
 Butler, R. P., Vogt, S. S., Laughlin, G., et al. 2017, AJ, 153, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Cosentino, R., Lovis, C., Pepe, F., et al. 2012, in Groundbased and Airborne Instrumentation for Astronomy IV, Proc. SPIE, Vol. 8446, 84461V [CrossRef] [Google Scholar]
 Covino, E., Esposito, M., Barbieri, M., et al. 2013, A&A, 554, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
 Damasso, M., Biazzo, K., Bonomo, A. S., et al. 2015, A&A, 575, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desidera, S., & Barbieri, M. 2007, A&A, 462, 345 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desidera, S., Sozzetti, A., Bonomo, A. S., et al. 2013, A&A, 554, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Desidera, S., Bonomo, A. S., Claudi, R. U., et al. 2014, A&A, 567, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dressing, C. D., & Charbonneau, D. 2013, ApJ, 767, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Borsa, F., Damasso, M., et al. 2017, A&A, 598, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 [NASA ADS] [Google Scholar]
 Dvorak, R. 1982, Oesterreichische Akademie Wissenschaften Mathematisch naturwissenschaftliche Klasse Sitzungsberichte Abteilung, 191, 423 [NASA ADS] [Google Scholar]
 Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Eaton, J. A., Henry, G. W., & Fekel, F. C. 2003, Astrophys. Space Sci. Lib., 288, 189 [NASA ADS] [Google Scholar]
 Eggenberger, A. 2010, EAS Pub. Ser., 42, 19 [CrossRef] [Google Scholar]
 Feng, F., Tuomi, M., Jones, H. R. A., et al. 2017, AJ, 154, 135 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Forveille, T., Bonfils, X., Delfosse, X., et al. 2009, A&A, 493, 645 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2016a, A&A, 595, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016b, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gomes da Silva, J., Santos, N. C., Bonfils, X., et al. 2011, A&A, 534, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gregory, P. 2010, Bayesian Logical Data Analysis for the Physical Sciences (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Hauser, H. M., & Marcy, G. W. 1999, PASP, 111, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Henry, T. J., Jao, W.C., Subasavage, J. P., et al. 2006, AJ, 132, 2360 [NASA ADS] [CrossRef] [Google Scholar]
 Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27 [Google Scholar]
 Howard, A. W., Marcy, G. W., Fischer, D. A., et al. 2014, ApJ, 794, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Kane, S. R., & Gelino, D. M. 2014, ApJ, 792, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M., Hartman, J., Bakos, G. Á., et al. 2011, AJ, 142, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1962, AJ, 67, 591 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Leggett, S. K. 1992, ApJS, 82, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Lidov, M. L. 1962, Planet. Space Sci., 9, 719 [NASA ADS] [CrossRef] [Google Scholar]
 Lippincott, S. L. 1972, AJ, 77, 165 [NASA ADS] [CrossRef] [Google Scholar]
 Lovis, C., & Pepe, F. 2007, A&A, 468, 1115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lovis, C., Dumusque, X., Santos, N. C., et al. 2011a, ArXiv eprints [1107.5325] [Google Scholar]
 Lovis, C., Ségransan, D., Mayor, M., et al. 2011b, A&A, 528, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maldonado, J., Affer, L., Micela, G., et al. 2015, A&A, 577, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maldonado, J., Scandariato, G., Stelzer, B., et al. 2017, A&A, 598, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mason, B. D., Wycoff, G. L., Hartkopf, W. I., Douglass, G. G., & Worley, C. E. 2001, AJ, 122, 3466 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Monin, J.L., Clarke, C. J., Prato, L., & McCabe, C. 2007, Protostars and Planets V, 395 [Google Scholar]
 Moutou, C., Vigan, A., Mesa, D., et al. 2017, A&A, 602, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ngo, H., Knutson, H. A., Bryan, M. L., et al. 2017, AJ, 153, 242 [NASA ADS] [CrossRef] [Google Scholar]
 Pepe, F., Lovis, C., Ségransan, D., et al. 2011, A&A, 534, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perger, M., GarcíaPiquer, A., Ribas, I., et al. 2017, A&A, 598, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poretti, E., Boccato, C., Claudi, R., et al. 2016, Mem. Soc. Astron. It., 87, 141 [NASA ADS] [Google Scholar]
 Reiners, A., Bean, J. L., Huber, K. F., et al. 2010, ApJ, 710, 432 [NASA ADS] [CrossRef] [Google Scholar]
 Rowe, J. F., Bryson, S. T., Marcy, G. W., et al. 2014, ApJ, 784, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Sozzetti, A., Bernagozzi, A., Bertolini, E., et al. 2013, Eur. Phys. J. Web Conf., 47, 03006 [CrossRef] [Google Scholar]
 Suárez Mascareño, A., Rebolo, R., González Hernández, J. I., et al. 2018, A&A, 612, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tanner, A. M., Gelino, C. R., & Law, N. M. 2010, PASP, 122, 1195 [NASA ADS] [CrossRef] [Google Scholar]
 Tokovinin, A., & Kiyaeva, O. 2016, MNRAS, 456, 2070 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, G. 1999, PASP, 111, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Trifonov, T., Kürster, M., Zechmeister, M., et al. 2018, A&A, 609, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Udry, S., Dumusque, X., Lovis, C., et al. 2017, A&A, submitted [arXiv:1705.05153] [Google Scholar]
 van Buren, D., Brundage, M., Ressler, M., & Terebey, S. 1998, AJ, 116, 1992 [NASA ADS] [CrossRef] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zucker, S., & Mazeh, T. 2002, ApJ, 568, L113 [NASA ADS] [CrossRef] [Google Scholar]
https://exoplanetarchive.ipac.caltech.edu/  18/09/2017
https://exoplanetarchive.ipac.caltech.edu/ – 28/09/2017
All Tables
Priors and bestfit results for the Gaussian process regression analysis of the chromospheric activity indicators extracted from the HARPSN spectra of Gl15A.
Planetary parameter bestfit values obtained through a joint modelling of Keplerian signals and correlated stellar noise, using Gaussian process regression.
Best orbital parameters for the Gl15 binary system from the MC simulation with priors on period and eccentricity.
All Figures
Fig. 1 Combined HARPSN and HIRES RV time series. As described in Sect. 4.3, we considered the HIRES data as splitinto two separate datasets for the purpose of the RV analysis. 

Open with DEXTER  
In the text 
Fig. 2 Posterior distributions of the fitted (hyper)parameters of the GP quasiperiodic model applied to the time series of CaII H&K (upper panel) and Hα (lower panel) activity indexes. The vertical lines denote the median (solid) and the 16th–84th percentiles (dashed). 

Open with DEXTER  
In the text 
Fig. 3 GLS periodogram of the complete HARPSN RV dataset. 

Open with DEXTER  
In the text 
Fig. 4 Panel a: evolution of the GLS orbital parameters as a function of the number of HARPSN RV measurements; panel b: evolution of the GLS orbital parameters as a function of the number of simulated data. The error bars on the final points are shown as references. 

Open with DEXTER  
In the text 
Fig. 5 Posterior distributions of the fitted (hyper)parameters of the twoplanet model, where the stellar correlated noise has been modelled with a GP regression using a quasiperiodic kernel. On the yaxis is shown the logarithm of the product between the likelihood and the prior. The vertical lines denote the median (solid) and the 16th–84th percentiles (dashed). 

Open with DEXTER  
In the text 
Fig. 6 Phasefolded RV curves for panel a Gl15A b and panel b Gl15A c. Each curve shows the residuals after the subtraction of the other planet and the stellar correlated signal. The red curve represents the bestfit Keplerian orbit, while the red dots and error bars represent the binned averages and standard deviations of the RVs. 

Open with DEXTER  
In the text 
Fig. 7 Upper panel: best fit stellar quasiperiodic signal (blue line) compared to the RV residuals. Lower panels: magnification of the highcadence HIRES/KECK observations (left panel) and of the last HARPSN observing season (right panel). 

Open with DEXTER  
In the text 
Fig. 8 Panel a: comparison between the last season HARPSN data (black triangles) and the overlapping CARMENES data (grey dots). The red solid line represents our twoplanet plus stellar noise model, with the pink shaded area representing the 1σ uncertainties on the GP hyperparameters; panel b: GLS periodogram of the combined HARPSN and CARMENES RV datasets. 

Open with DEXTER  
In the text 
Fig. 9 Position angle (panel a) and angular separation (panel b) of Gl15B with respect to Gl15A. Position angle is measured from north towards east. The red lines in both panels represent the second order polynomial fits. 

Open with DEXTER  
In the text 
Fig. 10 Radial velocities of (panel a) Gl15A (after subtracting the planetary and stellar signals) and (panel b) Gl15B. The solid lines show the respective first order polynomial fits. 

Open with DEXTER  
In the text 
Fig. 11 Orbital parameters of Gl15B as a function of the lineofsight separation z with Gl15A. 

Open with DEXTER  
In the text 
Fig. 12 Distribution of z resulting from the Monte Carlo simulation with the e and P priors injected. 

Open with DEXTER  
In the text 
Fig. 13 Fraction of time spent below the e = 0.40 (panel a) and e = 0.72 (panel b) thresholds. The different linestyles and symbols correspond to different values of e_{c,0}: solid and plus signs – e_{c,0} = 0; dotted and asterisks – e_{c,0} = 0.1; dashed and diamonds – e_{c,0} = 0.2; dash dot andtriangles – e_{c,0} = 0.3; dash dot dot and squares e_{c,0} = 0.4. The horizontal dashed grey lines indicate the corresponding probability from the MCMC – panel a: 68%; panel b: 95%). 

Open with DEXTER  
In the text 
Fig. 14 Fraction of time spent below the e = 0.40 (solid grey) and e = 0.72 (dotted black) thresholds, averaged on the initial eccentricity e_{c,0}. The errorbars show the standard deviation. 

Open with DEXTER  
In the text 