The HADES RV Programme with HARPS-N@TNG††thanks: 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 program) located at Serra La Nave on Mt. Etna.
Key Words.:techniques: radial velocities - stars: individual: Gl15A - stars: binaries: visual - instrumentation: spectrographs - planets and satellites: detection - planets and satellites: dynamical evolution and stability
We present 20 years of radial velocity (RV) measurements of the M1 dwarf Gl15A, combining 5 years of intensive RV monitoring with the HARPS-N spectrograph with 15 years of archival HIRES/Keck RV data. We carry out an MCMC-based 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.44-day period super Earth, Gl15A b, reducing its amplitude to m s ( ), and successfully models a long-term trend in the combined RV dataset in terms of a Keplerian orbit with a period around 7600 days and an amplitude of m s, corresponding to a super-Neptune mass ( ) 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 wide-separation companion Gl15B in a suite of detailed numerical simulations. In order to improve the results of the dynamical analysis, we derive 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 and , and can also enhance the eccentricity of the outer planet well above the observed value, even resulting in orbital instability, for inclinations around and .
The Gl15A system is the multi-planet system closest to Earth, at pc, and hosts the longest period RV sub-jovian 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.
Extrasolar planetary systems always showed a huge variety of orbital architectures, ever since the very first exoplanet orbiting a main sequence star was discovered (mayorqueloz95). In the following two decades the number of known exoplanets grew up to more than 3500111https://exoplanetarchive.ipac.caltech.edu/ - 18/09/2017, mainly discovered from transit (e.g., Kepler) and radial-velocity (e.g, HARPS, HARPS-N) observations, over a wide spread range of masses, radii and orbital separations, aiming down to the identification of small mass rocky Earth-twins (e.g. pepeetal11).
Low mass M dwarfs are the most common main sequence stars, comprising of the stars in the Solar Neighbourhood (henryetal06). Furthermore, M dwarfs have become the most promising ground for the hunt for low-mass, rocky planets (e.g. dreschar2013; sozzettietal13; astudillodefruetal2017), due to their more advantageous mass and radius ratios compared to solar-type 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. udryetal2017; roweetal2014, and references therein).
The nearby M1 dwarf Gl15A was studied by howardetal2014, who found a short-period super-Earth orbiting the star: they measured a period of d and an amplitude of m s. They also studied the activity signals of the host star, identifying the rotational period of the star to be d, both from the index analysis and also from the precise photometric light-curve, collected with the automatic photometric telescope (APT) at Fairborn Observatory (eatonetal2003). Gl15A has also a known M3.5 binary companion, Gl15B, identified astrometrically by lippincott1972 from a small fragment of its orbit, with a measured orbital separation of AU, and an orbital period of yr.
In this framework we present the clear detection of a long-period eccentric super-Neptune planet around Gl15A, from 5 years of high-precision Doppler monitoring with the HARPS-N high-resolution (resolving power ) optical echelle spectrograph (cosentinoetal2012) at the Telescopio Nazionale Galileo (TNG), combined with 15 years of archival RV data from the LCES HIRES/Keck Precision Radial Velocity survey (butleretal2017). The HARPS-N data were collected as part of the HADES (Harps-n red Dwarf Exoplanet Survey) programme, a collaboration between the Italian GAPS (Global Architecture of Planetary Systems, covinoetal2013; desideraetal2013; porettietal2016) Consortium222http://www.oact.inaf.it/exoit/EXO-IT/Projects/Entries/2011/12/27_GAPS.html, 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 howardetal2014. Our findings are also discussed in the light of the recent non-detection of Gl15A b in the CARMENES visual RV time series by trifonovetal2018.
With tens of exoplanets orbiting binary stars observed to date (e.g. eggenberger2010), including a confirmed wide binary system with planetary companions orbiting both components (desideraetal2014; damassoetal2015), it is debated how the presence of stellar companions influenced the distribution of planetary orbital parameters. For this reason, we derived new orbital parameters for the stellar companion Gl15B, using HARPS-N RV measurements along with astrometric measurements from the WDS (Washington Double Star, masonetal2001) 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 catalog data
As part of the HADES RV programme, Gl15A has been observed from BJD (27th August 2012) to BJD (18th January 2017). The total number of data points acquired was 115, over a time span of days. The HARPS-N spectra were obtained using an exposure time of minutes, and achieving an average signal-to-noise ratio () of at Å. Of the 115 epochs, were obtained within the GAPS time and within the Spanish time. Since the simultaneous Th-Ar calibration could contaminate the Ca II H & K lines, which are crucial in the analysis of stellar activity of M dwarfs (forveilleetal2009; lovisetal2011), 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 Th-Ar calibration.
The data reduction and RV extraction were performed using the TERRA pipeline (Template-Enhanced Radial velocity Re-analysis Application, anglada-escudebutler2012), which is considered to be more accurate when applied to M-dwarfs, with respect to the HARPS-N Data Reduction Software (DRS, lovispepe2007). For a more thorough discussion of the DRS and TERRA performances on the HADES targets see pergeretal2017. The rms of the TERRA RVs is m s, while the mean internal error is m s . The TERRA pipeline also corrected the RV data for the perspective acceleration of Gl15A, m s yr.
We also use in our analysis the HIRES-Keck binned data for Gl15A, downloaded from the LCES HIRES/Keck Precision Radial Velocity Exoplanet Survey. The HIRES time series spans days, from BJD (13th January 1997) to BJD (11th December 2014). These data were newly reduced with respect of the original RVs used by howardetal2014, and also new data have been observed with the HIRES spectrograph since the paper came out. For a full description of the catalog and data reduction see butleretal2017. We discarded the last data point of the HIRES time series (BJD ), since it is almost m s off with respect of the rest of the data. We thus use in our analysis HIRES RVs, showing a variation of m s and an average internal error of m s. The RV dataset from the Keck archive was already corrected for the perspective acceleration of the host star.
We combine the HARPS-N and HIRES RV datasets, obtaining a points time series, spanning days. The complete RV time series has rms m s, mean error m s . 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 Section 5, we obtained with HARPS-N 5 RV data points of the companion star Gl15B, collected from BJD (31st December 2016) to BJD (18th January 2017). The rms of the time series is m s, while the mean internal error is m s .
As for every target in the HADES survey, Gl15A was also monitored photometrically by the EXORAP (EXOplanetary systems Robotic APT2 Photometry) program, to estimate the stellar rotation from the periodic modulation in the differential light curve.
The observations were performed by INAF-Catania Astrophysical Observatory at Serra la Nave on Mt. Etna (+14.973E, +37.692N, 1725m a.s.l.) with the APT2 telescope, which is an 80cm f/8 Ritchey-Chretien robotic telescope. The detector is a 2k2k e2v CCD 230-42, operated with standard Johnson-Cousins 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 the B and V photometry respectively.
3 Stellar Properties of Gl15A and Gl15B
Gl15A is a high proper motion nearby ( mas) early-M dwarf, of type M1. We used the stellar parameters published by maldonadoetal2017, which were calculated applying the empirical relations by maldonadoetal2015 on the same HARPS-N spectra from which we derived the RV time series. This technique calculates stellar temperatures from ratios of pseudo-equivalent widths of spectral features, and calibrate the metallicity on combinations and ratios of different features. Although such techniques are mainly used for solar-type stars, maldonadoetal2015 proved them to be just as effective on low-mass stars.
howardetal2014 derived a rotational period of d, both from their Keck-HIRES measurements of the index and their APT photometric observations at Fairborn Observatory. Recently suarezmascarenoetal2017b analysed the potential signatures of magnetic activity in the CaII HK and H activity indicators of the HADES M-dwarfs sample. For Gl15A they computed a mean level of chromospheric emission , and a rotation period d, fully consistent with the value found by howardetal2014. The definition of can be extended for application on M-dwarfs spectra (suarezmascarenoetal2017b, and references therein).
All the stellar parameters for Gl15A used in this work are summarized in the left column of Table 1. We can see that most of them are fully consistent with the values used by howardetal2014.
|Spectral Type||M1 maldonadoetal2017||M3.5 This work|
|FeH dex||maldonadoetal2017||This work|
|Mass M||maldonadoetal2017||This work|
|Radius R||maldonadoetal2017||This work|
The stellar companion of Gl15A, Gl15B, is a type M3.5 dwarf whose orbit was measured astrometrically by lippincott1972. For the purpose of our orbital analysis in Sect. 5, we took 5 HARPS-N spectra of Gl15B during the last observing season at TNG (see next Section) and we applied on them the maldonadoetal2015’s techniques to calculate updated stellar properties. The derived values are listed in the right column of Tab. 1.
3.1 Photometric analysis
In order to identify the potential modulation in the stellar photometry due to the presence of active regions, we analyse with the GLS periodogram (zechkur2009) the B and V time series collected within the EXORAP survey, composed of 242 and 233 points respectively, taken over five seasons from the 13th of August 2013 to the 15th of June 2017. No evident periodicity is found in either time series, in contrast with the findings of howardetal2014 who found a clear signal at 43.82 days (see their Figure 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 photometric precision of the observations: the amplitude of the signal found by howardetal2014 was only of mmag, which is below the sensitivity of the EXORAP survey for targets in the magnitude range of Gl15A.
4 Spectroscopic data analysis
In Figure 1 an evident long period signal can be seen. howardetal2014 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 m s yr, and concluded it to be consistent with the orbit of the Gl15AB system as calculated by lippincott1972. But considering also the new HARPS-N data, it becomes clear how the long term RV variation cannot be modeled as a linear trend.
Even if it could no longer be treated as a constant acceleration, the long period 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 kippingetal2011 to study the presence of a long-period companion in the HAT-P-31 system: they modeled the long-period signal as a quadratic trend, and then derived a range of orbital parameters from the quadratic coefficients. The ratio between the semi-amplitude, , and period, , of the companion signal is given by the second-order term of the trend, , as (their Equation (4)):
and, since depends on the orbital period and on the mass of the two stars, the orbital period can be derived as a function of the second order term and the masses .
From the fit of the complete RV time series we obtain: m s day. Using the stellar masses listed in Table 1, the resulting period is yr, which would correspond to a semi-major axis AU. This solution is clearly unphysical, since the presumed semi-major axis is less than half the observed orbital separation of the two objects.
Moreover, several high-contrast imaging surveys ruled out the presence of additional co-moving stellar objects close to Gl15A: vanburenetal1998 excluded the possibility of objects with M at separations of AU, while tanneretal2010 ruled out objects up to a magnitude contrast of mag within (3.57 AU) and mag within (17.8 AU).
The long-period signal observed in the RV time series is therefore very unlikely to be caused either by Gl15B or by additional stellar companions. Instead a long-period planetary-mass 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 long-period 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 howardetal2014, 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 timeseries (e.g. dumusque2017, and references therein), especially when adopting a quasi-periodic 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 quasi-periodic kernel is described by the covariance matrix
where and indicate two different epochs.
This kernel is composed by a periodic and an exponential decay term, so that it can model a recurrent signal linked to stellar rotation, also taking into account the finite-lifetime of the active regions. Such approach is therefore particularly suitable when modeling signals of short-term timescales, as those modulated by the stellar rotation period.
About the covariance function hyperparameters, is the amplitude of the correlations; represents the rotation period of the star; 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. 4.1, is the data internal error at time for each instrument; is the additional uncorrelated ’jitter’ term, one for each instrument, that we add 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 nor modeled by the quasi-periodic kernel; is the Kronecker delta function.
In the GP framework, the log-likelihood function to be maximized by the Markov chain Monte Carlo (MCMC) procedure is:
where is the number of the data points, K is the covariance matrix built from the covariance function in Equation (4.1), and is the data vector.
We use the publicly available emcee algorithm (foreman13) to perform the MCMC analysis, and the publicly available GEORGE Python library to perform the GP fitting within the MCMC framework (ambi14). We used 150 random walkers to sample the parameter space. The posterior distributions have been derived after applying a burn-in as explained in eastman13 (and references therein). To evaluate the convergence of the different MCMC analyses we calculate the integrated correlation time for each of the parameters, and run 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 ensures that the emcee walkers thoroughly sample the parameter space (foreman13).
4.2 Analysis of the activity indexes
We first investigate both the HIRES and HARPS-N CaII HK and H index time series, in order to test the potential stellar origin of the long period variation seen in the combined RV times series shown in Figure 1. No long period trend is found in either the HIRES or HARPS-N datasets. suarezmascarenoetal2017b found a magnetic-cycle type periodicity at yr in the HARPS-N CaII HK and H indexes, even if we don’t find any similar signal in the respective HIRES time series. Nevertheless the period of this cycle is far from the time span of the long-period signal we care to investigate, so it could not be the origin of it. Another clue of the stellar origin of the long-period 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 suggest that the long-period signal is not due to the star magnetic cycle or chromospheric activity, and reinforces the hypothesis of it be due to a wide-orbit planetary companion.
To further test the effect of the activity of Gl15A, we then investigate the stellar activity behaviour over the four seasons covered by HARPS-N observations by analysing chromospheric activity indicators based on the Ca II H and K, and H spectroscopic lines. They were extracted using the definition of the line cores and of the reference intervals given in gomesdasilva11. We analyse here only measurements obtained from HARPS-N spectra because they represent an homogeneous dataset. By spanning more than 1600 days, the HARPS-N data alone can provide robust insights into the long- and short-term stellar activity variability. The average S/N was 18 for the CaII H&K index and for the H.
We performed an analysis of the activity indicators based on a Gaussian process (GP) regression, as detailed in the previous Section. By adopting the same covariance function (Eq. 4.1) we use 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 pre-whitened, and the variability patterns in the former could be present with similar properties in the latter, as was noticed by afferetal2016 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 modeled. Here we use Eq. 4.1 to describe the variability correlated with the stellar rotation period , 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 best-fit 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 HK index is very similar to that found by howardetal2014 for the CaII HK S-index derived from HIRES spectra. The hyper-parameter found in the analyses of the two time series is in very good agreement with the rotation period found by suarezmascarenoetal2017b with a GLS analysis of the activity indexes time series of Gl15A.
|Jump parameter||prior||Best-fit value|
4.3 Analysis of the combined HIRES and HARPS-N RV time series
First we study the HARPS-N RV time series for confirmation of the presence of the Gl15A b signal found by howardetal2014. We perform a GLS analysis to identify any periodicity in our data in d. In the periodogram, shown in Fig. 3, we can see that the d period of planet Gl15A b is clearly recovered. As a confirmation of the coherence of the signal throughout the entire HARPS-N 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 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 HARPS-N time series, and found no significant signal at periods close to the inner planet period , thus confirming its planetary nature as announced by howardetal2014.
It is worth noticing that, as can be seen in Fig. 3, the final amplitude recovered by the GLS periodogram on the HARPS-N dataset, m s, is smaller than the one published by howardetal2014, m s. 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 HARPS-N time series but containing only the planetary signal and white noise, as shown in 4b. The lower amplitude value than that found by howardetal2014 can also be explained similarly, since the HARPS-N time series is composed of five years of continuous high-cadence observations of the target, while the dataset used by howardetal2014 consisted in mostly sparse measurements with an intensive high-cadence 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 (butleretal2017) than the data used by howardetal2014, and performing a GLS analysis of the two time series we obtain the same peak periodicity but an amplitude smaller in the new archive data.
Since no other short period signal emerges from the GLS analysis of the RVs, we focus our attention on the study of the long-period signal, which we assume to be due to a planetary-mass wide-orbit companion.
To estimate the orbital and physical parameters of the known planet Gl15A b and the candidate companion, hereafter Gl15A c, we have performed a Markov chain Monte Carlo analysis of the RVs. Following the prescription of howardetal2014, we model the RVs dividing the HIRES dataset in a pre-upgrade and a post-upgrade sublist, due to the HIRES CCD upgrade occurred on August 2004. Each subsample is then treated as an independent dataset, with its own zero-point 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 equation
where ; is a function of time t, time of the inferior conjuntion , orbital period , eccentricity e and argument of periastron ; is the RV offset, one for each instrument; is the stellar noise modeled 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 Sec. 5.1). Instead of fitting separately and , we use the auxiliary parameters and to reduce the covariance between and . howardetal2014 found 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 , we assume 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 with Gl15B, as we will discuss in Sec. 5.
Except for very few parameters, for our analysis we assumed uniform priors. Our choice for the range of is justified by the results obtained from the GP analysis of the Ca II 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 semi-amplitudes of the Doppler signals we used the modified invariant scale prior:
following the prescription in gregory10. For the orbital period of Gl15A b there was a strong constraint from the results of howardetal2014, confirmed by our GLS analysis previously discussed. Nevertheless, since the same dataset from howardetal2014 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 centered on the value found by howardetal2014, e.g. . 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 3 summarizes the best-fit 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 best-fit values for the free parameters, we first note that the semi-amplitude of the Doppler signal induced by Gl15A b is even lower than the estimate obtained from the GLS periodogram, differing from the estimate of howardetal2014 (2.940.28 m s) 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 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- ( at a 68 level of confidence). Our estimate for the minimum mass of Gl15A c places this companion in the group of the so-called super-Neptunes, as are generally referred planets in the M mass range.
The eccentricity of Gl15A b also resulted consistent with zero within 1.5-, with a much lower value ( at a 68 level of confidence). Our analysis thus agrees with the results by howardetal2014, who concluded that there was no evidence for an eccentric orbit to be preferable to a circular one.
Fig. 6 shows the RV curves folded at the best-fit 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 best-fit stellar correlated, quasi-periodic signal.
As for the GP hyper-parameters, the stellar rotation period assumes a value which is in agreement with that expected, and the higher uncertainties than those for the activity indicators (Tab. 2) should be due to the longer timespan covered by the RVs, which include HIRES data. The longer timespan should also explain why the active-region evolutionary timescale sets on a value less than those found for the activity indicators. Without any ancillary data available as photometry, covering the same timespan 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 it was observed before in previous studies of different systems (e.g. afferetal2016). Also, the hyper-parameter, corresponding to the mean amplitude of the stellar signal, is fully compatible with the expected value derived by suarezmascarenoetal2017b from the mean activity level of the star, , that is m s.
|Jump parameter||prior||Best-fit value|
|[m]||Modified scale invariant||1.68|
|[BJD-2 450 000]||6995.86|
|[m]||Modified scale invariant||2.5|
|[BJD-2 450 000]||10730|
|Derived quantitiesDerived quantities from the posterior distributions. We used the following equations (assuming ): (; , where is the gravitational constant.|
|¡0.13 ( percentile)|
|¡0.40 ( percentile)|
4.4 Comparison with CARMENES results
In a very recent work, trifonovetal2018 have published 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 showed no evidence of the presence of Gl15A b, thus the authors concluded it to be an artifact due to the stellar noise. The authors also identified a long-term trend in the combined Keck+CARMENES dataset which they proposed as due to a distant low-mass companion. As discussed in Sec. 4.3, our HARPS-N 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 HARPS-N+Keck datasets makes a decisive case for the existence of the long-period low-mass giant Gl15A c based on a 5-yr time-span of HARPS-N observations, partly overlapping with the Keck time-series. The orbital elements for Gl15A c derived in this work have larger formal uncertainties than those reported in trifonovetal2018 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 HARPS-N observations, that overlaps with most of the CARMENES data, compared to our two-Keplerian plus correlated stellar noise model. This season is clearly dominated by activity-related variations in the optical HARPS-N spectra, while the internal errors for HARPS-N are typically 3 times smaller than those of CARMENES (0.6 m s vs. 1.8 m s). We nonetheless tested the effect of combining the HARPS-N 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 HARPS-N dataset alone.
We also tested our two planet + stellar noise model described in Sec. 4.3 on the complete dataset obtained combining the HIRES, HARPS-N 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 semi-amplitude and minimum mass for the inner planet are m s and M respectively. These values are well within 1- of the results from the analysis of only the HIRES and HARPS-N datasets: no significant drop in the amplitude of Gl15A b signal is observed. It’s also worth noticing that the residual jitter found in the analysis is small ( m s) 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 Tab. 3 ( m s, d and d), meaning that the activity signal description derived from the first two datasets is robust also in modeling the CARMENES data. Given the intrinsically higher quality of the HARPS-N data, that drive the GP regression modeling 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 Sec. 5.2.
trifonovetal2018 also stated, as an argument in favour of the non-Keplerian interpretation of Gl15A b, that the d signal disappeared when analysing the last two years of HIRES observations, subsequent to the time series used by howardetal2014. Studying the same time span of data, we observed that, when ignoring the outlier described in Sec. 2, the days signal is still clearly visible in the GLS periodogram.
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 probability to be and probability to be . 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 investigate 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, , and the inclination, , with the same period but opposite phase. It was originally studied by lidov1962 and kozai1962 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 modeling from astrometry and RV data
The first obstacle in the dynamical analysis of the system was the poorly constrained orbital parameters of the companion. lippincott1972 estimated a period of yr from yr of astrometric measurements, which cover less than of the orbit.
Dealing with a similar case of an eccentric planet hosted by a wide binary system, hausermarcy1999 developed a technique for constraining long-period binary orbital parameters, combining astrometric and radial velocity measurements. The method is based on the fact that, being newtonian mechanics deterministic, knowing the instantaneous full position and velocity vector, , 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 hausermarcy1999, 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 and which we need to calculate the velocity component in the plane of the sky, we fitted the data with a second order polynomial:
From these we can easily derive and as:
Adopting the parallax value from Table 1, mas, we can derive the Cartesian position and velocity using their Equations (1)-(4):
To obtain the third component of the velocity vector, , we use the combined Doppler information of the two stars: from the same HARPS-N spectra we used to obtain the RV time series for the planet detection, collected as illustrated in Sect. 2, we extract the absolute radial velocities with the DRS pipeline. We use 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 HARPS-N epochs, and subtract from the absolute RV the planetary and stellar signals, as derived in Sect. 4.3. For Gl15B we use the 5 spectra we took on January 2017, as described in Section 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 fit the two RV time series with first-order polynomials. From the difference of the two linear fits we obtain the relative line-of-sight velocity of Gl15B with respect to Gl15A. We can then select an epoch, and compute the values of for that time.
The selected epoch to derive the binary orbit is BJD , that is January 1st 2017, which is well in between 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 line-of-sight separation , which cannot be measured, but can be constrained imposing the condition of a bound orbit for the binary system, which 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:
where , and, of course, .
From this we get a range of acceptable values AU. From every value of is possible to compute all the orbital parameters for the binary system (see hausermarcy1999, Eq. (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 constrain, and this is still insufficient for any meaningful dynamical analysis.
The procedure by hausermarcy1999 provides the orbital solution for every single value of , 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 line-of-sight separation .
To do this we perform 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 apply a prior on , to have the best possible constraint on it. We select the eccentricity distribution from tokovininkiyaeva2016, who studied a sample of 477 wide binaries within 67 pc from the Sun with median projected separation of AU, very close to that of the Gl15AB system. The forementioned prior is:
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 duchenekraus2013 for low-mass binary stars, that is a log-normal distribution, with AU and .
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 almost unbound, so it can be ignored, both because, as we said, 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 intended to study. The latter can be said also on the peak on the left, which correspond to a period of yr, and a semi-major axis of AU, again too distant from the planetary system to have a significant influence.
We thus use the central peak of Fig. 12 to derive the orbital parameters best solutions and error bars. To do this we fit a truncated normal distribution in the range AU, and use the mean value to calculate the corresponding best solution orbital parameters as described before. The upper and lower errors on the orbital parameters are calculated by taking and and deriving the corresponding values of . The orbital parameters solutions and errors are listed in Table 4.
Assuming the best fit orbital parameters listed in Table 4, we can 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 d , which is fully compatible with the result from our MCMC listed in Table 3.
5.2 Lidov-Kozai Interaction modeling
The results of the previous section show that the most likely orbit of the stellar companion Gl15B has a high eccentricity, . 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 literature-coined acronym EKL). This mechanism applies to hierarchical triple-body 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 three body systems.
The EKL mechanism is very sensitive to the mutual orientations of the orbits of the perturber, i.e. the companion star, and of the influenced body, i.e. the planet. We have derived of Gl15B, but we do not know either 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 Sec. 4.3, we calculate the fraction of time spent by Gl15A c below and ( and ), which can be considered a proxy of the probabilities ( and , respectively) to observe the system with eccentricities under those thresholds (e.g. andersonlai2017). The same can be done to study the interaction of Gl15B with Gl15A b, in which case the and probability levels correspond to and .
We perform some preliminary tests to verify the various mechanisms which could be involved in the dynamical evolution of the system: we prove 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 focus our efforts to study the orbit of the newly discovered Gl15A c; we also check 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 and the ones of Gl15A c with . All the EKL integrations were performed with a timescale of 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 . The other unknown is the longitude of the ascending node of the outer planet . 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:
= arccos(sini_c cosΩ_c sini_B cosΩ_B +
sini_c sinΩ_c sini_B sinΩ_B + cosi_c cosi_B), and thus can be ignored as long as is conveniently sampled, which can be obtained by changing . For our analysis we fix and vary the planet orbital inclination in the interval .
The results are shown in Figure 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 consider 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 possibly producing planet-planet scattering events, that is when:
where is the maximum eccentricity reached during the EKL oscillations. In these cases we also consider .
As we can see in both the panels of Fig. 13, there are regions in the space that are unstable regardless of the initial eccentricity of the planet : between and and around . Only the solid black line, corresponding to , behaves differently, showing the system to be stable between and and unstable at larger inclinations. We can also see that the Lidov-Kozai interaction is weak for and strengthen as decreases. The top panel of Fig. 13 shows that for the resulting eccentricity ranges are compatible with the observed values. The constraints for the threshold are somewhat looser, but pointing in the same direction.
Since as previously said we do not know the initial eccentricity of the system, 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 just discussed, with initial inclinations around and between and leading to an unstable system, and inclinations higher than producing a good agreement with the observed orbital configuration.
6 Discussion and conclusions
We present in this paper the fifth planet detected by the HADES programme conducted with HARPS-N at TNG. This long period planet was found orbiting the planet-host M1 star Gl15A, from the analysis of high precision, high resolution 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 long-period 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 M and M for the inner and outer planet respectively. The mass we find for Gl15A b is much smaller than the value found by howardetal2014, which was almost double due to the higher signals amplitude. The smaller value we find, can be easily explained by the additional information brought by the high-precision HARPS-N RVs, along with the new calibration of the archival HIRES data published by butleretal2017. The combined dataset is almost twice as large that the one analysed by howardetal2014, stretched on a significantly longer timespan, with better sampling and precision. This, together with the simultaneous modeling 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 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 howardetal2014 ( at a level of confidence), thus confirming their conclusion that the planet’s motion is best described by a circular orbit.
With its period of yr, Gl15A c is the longest-period sub-Jovian planet detected up to date with the RV method555https://exoplanetarchive.ipac.caltech.edu/ - 28/09/2017, the second being HD 10180 h (lovisetal2011b) with a period of yr and a minimum mass of M (kanegelino2014). With the confirmed presence of two widely spaced planetary mass companions, Gl15A is now the multi-planet system closest to our Sun, at a distance of only pc.
We also compared the results of our analysis of the HIRES and HARPS-N RV data with the recent results by trifonovetal2018 based on CARMENES high-cadence monitoring of the target. trifonovetal2018 found no evidence for the presence of planet , while we showed in Sec. 4.4 that the signal can be clearly detected when combining the HIRES, HARPS-N and CARMENES data, without any loss in significance. Things notwithstanding, however, the reason why CARMENES does not see the signal of Gl15A b is not fully clear, but the higher quality, and much longer timespan, of the HARPS-N data, combined with our modeling of the stellar activity quasi-periodic signal could be a possible explanations for the non-detection based on the CARMENES data alone. Another factor contributing to this non-detection could be the sampling of the CARMENES RV data, which appear not to be homogeneously distributed over the 11.44 d period (see Fig. 10 from trifonovetal2018).
The CARMENES visual arm contains a spectral region extending all the way to m, i.e. a significantly redder spectral range than that covered by HIRES and HARPS-N. As the amplitude of activity induced RV variations is known to be chromatic (reinersetal2010), the non-detection of the 11.44 d signal in the CARMENES time series could be an indication of a wavelength dependent-amplitude of the signal, that would clearly indicate its stellar origin. However, the CARMENES visual arm spectral range still significantly overlaps with that of HIRES and HARPS-N and thus it must be affected by stellar activity in a similar way. Based on the higher RV precision of HARPS-N, allowing a detailed modeling of quasi-periodic 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 HARPS-N time series, the Keplerian origin of the signal still seems the most straightforward explanation for the observed data. It would be however interesting to carry out a systematic study on both CARMENES and HARPS-N time series of this target, adopting the same strategy outlined by fengetal2017 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 for Sun-like stars (duqmay1991) and somewhat lower for M dwarfs (bergforsetal2012). Many young binaries possess either circumstellar or circumbinary disks (e.g. moninetal2007), and the existence of stable planetary orbits in binary systems was postulated well in advance of the first exoplanets discoveries (dvorak1982).
Early studies proposed different mass-period relations for planets around binaries and single stars (zuckermazeh2002) but in the following years the evidence of such diversity decreased (e.g. desiderabarbieri2007; eggenberger2010), until most recently ngoetal2017 claimed not to be any difference of planetary properties between the two kind of systems. On the other hand, recent works like moutouetal2017 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, , 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, and , 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 semi-amplitude 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 long-period planetary-mass companion. Additional RV observations in the years to come will, however, be very helpful to better constrain the 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 of the full orbit. However, based on the torres1999 formalism curvature effects in the stellar motion should typically amount to as yr, thus they should be easily revealed by Gaia, that for such a bright star as Gl15A will be able to deliver end-of-mission proper motion accuracies as yr (e.g. gaiaetal2016).
Yet even with the orbital solution 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/2007-2013) under Grant Agreement No. 313014 (ETAEARTH). The HARPS-N 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 TNG-INAF 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 AYA2014-56359-P. J.I.G.H. also acknowledges financial support from the Spanish MINECO under the 2013 Ramón Cajal program MINECO RYC-2013-14875. 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 HARPS-N 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 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 errors are due to photon noise through error propagation.