© ESO, 2017

1. Introduction

Exoplanet searches by radial velocity (RV) or transit methods have revealed the presence of giant planets (GPs) orbiting very close to their stars. Core accretion (CA), the commonly accepted formation mechanism for GPs in the solar system, is expected to occur in a region of a few au in width around the star beyond the iceline (see Kennedy & Kenyon 2008). To explain the existence of the hot Jupiters, the migration of the planet within the protoplanetary disk at an early stage of evolution has been proposed (see e.g. Lin et al. 1996).

Recently, several GPs on very remote orbits (dozens to hundreds of au from their stars) have been imaged. These GPs cannot have formed in situ by CA as this process would take too much time at these distances. They may have formed via another mechanism such as gravitational instability (GI) within the protoplanetary disk (Cameron 1978; Boss 1997) or within collapsing clouds, i.e. like binary systems. They may alternatively have formed closer in and have then been kicked away by some dynamical mechanism such as planet-planet scattering within the planetary system or via passing-by stars. A recent example of such a scenario was extensively investigated by Rodet et al. (2017) for the HD 106906 system, which hosts a planet orbiting at more than 650 au (see Bailey et al. 2014) and a debris disk (see Lagrange et al. 2016); we recently showed, using HARPS and VLTI/PIONIER data, that the central star is in fact a binary (Lagrange et al. 2016). HD 106906 is located in the Lower Centaurus-Crux LCC (Sco-Cen) association, the nearest (≃100–200 pc) OB association and star forming region (de Zeeuw et al. 1999).

We recently acquired HARPS data of the young (16 ± 1Myr, Pecaut et al. 2012), 1.82M_Sun, A2-type star star HD 131399 (HIP 72940), another member of Sco-Cen which hosts a remote giant planet recently discovered by the SPHERE extreme adaptive optics imager (Wagner et al. 2016); the planet’s semi-major axis and orbital parameters are not well constrained, but its projected separation is greater than 80 au. This star forms a hierarchical triple system with the tight binary HD 131399BC, located about 350 au away. Interestingly, there is so far no evidence of dust around HD 131399, conversely to HD 106906. Wagner et al. (2016) have suggested different scenarios to explain the current planet orbit: “Scenario A) the planet formed on a short orbit around star A, and subsequently underwent a planet-planet scattering event that ejected it onto its current orbit; Scenario B) The planet formed as a circumbinary planet around components B and C and underwent a scattering event via interactions with another planet or with the binary itself, and was furthermore captured by star A; Scenario C) The planet formed around either component before the triple system arrived in its present configuration. The stellar orbits could have evolved subsequently due to interactions with the natal disks or secular effects”. Noticeably, Wagner et al. (2016) and Veras et al. (2017) have shown that the planet orbit may not be stable over the system lifetime, due to the high ratio between the planet-HD 131399 A separation and the HD 131399BC–HD 131399 A separation.

We report the results of our HARPS observations. The data and results are described in Sect. 2. In Sect. 3 we list some consequences of these results on the HD 131399 system.

2. Observations and analysis

2.1. Observations

We acquired 116 high resolution (R ≃ 110 000), high S / N spectra with HARPS, between July 2016 and January 2017. Each spectrum contains 72 spectral orders, and covers the 380–690 nm wavelength range. The S / N of the spectra at 550 nm varies between 70 and ≥200; the average S / N is 153. For each pointing, at least two consecutive spectra of a few minutes in length were recorded. In addition, two continuous 1h30 sequences of individual spectra were obtained to check for possible short-term variations.

The HD 131399 spectrum is characteristic of a low vsin(i) early-type star¹. The vsin(i) measured on the cross-correlation functions (CCF) is about 19 km s^-1. No short-term variations are detected on timescales of minutes to hours, but huge RV variations are observed, with an amplitude of more than 60 km s^-1, on timescales of days. No correlation between the RV and bisector velocity span (hereafter BVS) measured on the CCFs are observed (see Fig.1). Finally, the Lomb-Scargle periodogram and clean periodograms of the RV indicate a peak at about 10 days (Fig. 1), which is not present in the BVS periodogram.

The RVs measured with our SAFIR software (Galland et al. 2005) are provided in Appendix A and shown in Fig. 1.

Fig. 1 1st row: HARPS RV time series of HD 131399 (the error bars are lower than 30 m/s). 2nd row: Bisector Velocity Span (BVS) as a function of RV. The y-scale is different from that in the first panel. 3rd row: Clean Periodogram of the RV. 4th row: Clean Periodogram of the BVS.

2.2. Radial velocity analysis

The absence of correlation between RV and BVS indicates that the RV variations are not due to spots. Spots can also be ruled out given the star’s early spectral type and the extremely high amplitude of the RV variations. Pulsations with periods as long as 10 days are not known today among A-type main sequence stars. We therefore attribute the RV variations to a close companion.

The RV data are very well fitted by a low mass star companion (M sini = 0.43 M_Sun) orbiting with a period of 9.9 days (hence with a semi-major axis a = 0.12 au) and an eccentricity e = 0.49. The amplitude of the residuals is about 160 m s^-1, compared to the 60 km s^-1 amplitude of the input RV data. The rms of the residuals is 30.6 m s^-1, compared with the 19 km s^-1 rms of the input RV data. Figure 2 shows the results of the fit made using the yorbit software (Ségransan et al. 2011), and Table 1 summarizes the main parameters of the orbit.

Fig. 2 Phase folded fit of HD 131399A RV variations.

3. Discussion

3.1. Companion properties

The true mass M_p of the companion cannot be straightforwardly derived from the RV orbital fit. The mass derived from the Yorbit fit, M_psini = 0.43 M_Sun (Table 1) is obtained assuming M_p ≪ M_∗, where M_∗ is the mass of the central star. This assumption is valid for exoplanets, but may be questionable here. Classical Keplerian formalism gives the RV signal induced by a companion onto a star as (1)where K is the amplitude, ω is the argument of periastron, e is the eccentricity, and ν is the true anomaly on the orbit. The yorbit fits provide K, ω, e, the orbital period P, and the time of periastron passage t_p. The amplitude K is related to the mass M_p as (2)Whenever M_p ≪ M_∗, the (M_p + M_∗)^{2 / 3} term is usually approximated as , which allows M_psini to be derived from the fit, provided M_∗ is known. We must consider here the full formula, which changes the determination of M_p for a given sini. Figure 3 illustrates this correction, showing the relationship between sini and M_p using the standard M_psini determination (blue) and the full formula (red). We see that for any given inclination, the use of the standard approximation leads to an underestimation of M_p. Noticeably, the minimum M_p value obtained for i = 90° (edge-on configuration) is now 0.51 M_Sun instead of 0.43 M_Sun.

Fig. 3 Relation between the mass Mp of the companion of HD 131399A and its orbital inclination with respect to the plane of the sky. Blue curve: assuming Mp ≪ M∗; red curve: general case (see text).

Determining precisely the maximum mass of the companion is not possible with the present data. To get a rough upper limit, we estimate the mass that would lead to a detectable spectrum of the secondary in the best conditions, i.e. at maximum velocities of both components and for a low vsin(i) companion. We only consider G-M companions as spectra of earlier type companions would be detectable. The maximum velocity of the secondary is estimated using the relation V_B/V_A = M_A/M_B (which is only true for circular orbits, and approximate for an eccentric orbit), where V and M refer to the maximum velocities and masses of A and B. Given the low vsin(i) of HD 131399A, both spectra should be largely resolved (no blend) at maximum velocities, so the CCF signal of the companion should only be compared to the noise. We then approximate that the companion is detectable if C/d ≥ S/N, where d is the depth of the secondary CCF, C is the contrast between the primary and the secondary, and S/N is the signal-to-noise ratio. To get an order of idea of the depth of the companion CCF, we considered the CCFs of the G-M stars of our Harps survey of members of young and close-by associations that have vsin(i) ≤ 5 km s^-1. A large scatter is found, from 0.2 to 0.5, for the CCF depths. In a favourable case where the depth is 0.5, we find that companions earlier than K0 should give detectable spectroscopic signatures. Hence the companion spectral type is earlier than K0 and its mass lower than about 0.6–0.8 M_⊙.

HD 131399A has a vsini of 19 km s^-1, which is lower than the mean vsini of early A-type stars. From the study of 191 stars with masses in the range 1.6–2.0 M_Sun, Zorec & Royer (2012) get a distribution of equatorial velocities peaking at about 150 km s^-1. Assuming this value for HD 131399A would imply an inclination of HD 131399 rotation axis i ≲ 10°. In addition, given the Msin(i) of the companion and its maximum mass, we conclude that the companion orbit is very inclined with respect to the equatorial plane of the primary. This conclusion is also supported by the eccentric nature of the orbit. Tight binaries are known to undergo tidal circularization during the pre-main sequence phase up to periods of 8–10 days (Melo et al. 2001). The HD 131399A pair falls close to this threshold and thus might not be circularized yet. However, due to the high eccentricity of HD 131399A companion, tides should be active at periastron. In this case, even if it is not fully circularized today, the orbit should have moderate eccentricity (≲0.1). The observed high eccentricity is an indication of strong ongoing or past gravitational perturbations onto this system that could also have caused a misalignment with the star’s equatorial plane.

These perturbations could arise or have arisen from the distant pair HD 131399BC. Wagner et al. (2016) report for the wide A-BC orbit an inclination between 45° and 65° with respect to the plane of the sky, i.e. also with respect to the HD 131399A equatorial plane if we assume a pole-on configuration. This could trigger a misalignment of the stellar companion, possibly via Kozai-Lidov (KL) cycles (Kozai 1962), which are known to cause large amplitude eccentricity and inclination variations. We note that we do not expect KL cycles to be currently active on the orbit. KL cycles in triple systems are characterized by an oscillation period that basically scales as within a factor of order unity (see Beust & Dutrey 2006; Krymolowski & Mazeh 1999; and Ford et al. 2000), where P_out is the orbital period of the outer perturber, P_in is the inner orbital period, M_tot is the total mass, and M_in is the mass of the central body. Taking P_in = 9.9 days, P_out = 3556 yr, M_in = 1.82 M_Sun, M_tot = 1.82 + 0.96 + 0.6 M_Sun, this gives a KL oscillation period of ~1 Gyr. This very slow precession induced to the central binary is likely to be overriden by the orbital advance due to General Relativity if not by tides, as shown by Fabrycky & Tremaine (2007) and Beust et al. (2012). The orbit could nevertheless have been wider in the past and been subject to KL cycles from the outer pair, thanks to a larger P_in. In this case, tides started to act at periastron in peak eccentricity phases of the KL cycles. This caused the orbit to shrink to its present-day state. This process, known as Kozai migration, was described by Wu & Murray (2003), Fabrycky & Tremaine (2007), and Beust et al. (2012). It naturally generates short-period but misaligned orbits, and could possibly have been at work in the HD 131399 system.

3.2. Link with planetary systems

HD 131399 thus appears to be a hierarchical multiple system made of two pairs of stars. While the HD 131399BC pair is a binary separated by a few au, HD 131399A is a very tight binary, with a separation of 0.12 au. The 17 confirmed imaged planets orbiting stars and with masses ≤16 M_Jup (to take into account possible metallicity effects on the limit mass for deuterium burning, Spiegel et al. 2011) orbit further than 20 au, except β Pic b (9 au), HR8799 e (14 au), and 51 Eri b (14 au). These remote GPs cannot a priori have formed situ by classical CA (Kennedy & Kenyon 2008), even though pebble accretion could help form CA planets at larger separations (tens of au) from their host stars than previously expected (Bitsch et al. 2015). These GPs are then regarded as good candidates for GI formation. Yet, among these 17 systems, we note that 7 are members of binary systems, with separations between fractions of au (HD 106906, HD 131399) and hundreds of au (2000 au for 51 Eri, Feigelson et al. 2006). Of course, the actual number of binaries could be higher (hence greater than 50%).

3.3. Impact on the planet dynamical history

The planet HD 131399A b is thus a rare example of a circumbinary planet imbedded in a quadruple stellar system. Two major dynamical issues can be raised concerning the planet. The first is its stability, which is not clear (Wagner et al. 2016) given its not-so-small separation compared to the separation between the two binaries. Even though the planet’s current orbit (eccentricity, semi-major axis, etc.) is not well characterized and given the very short separation between both stars, it is very probable that the binary nature of HD 131399A does not play a role in the current stability of the planet. We also note that considering the probable misalignment of the inner binary orbit with respect to the central star’s equatorial plane, no hypothesis about the current orbital plane of the planet can be proposed.

The second issue is the past dynamical history of the planet and its present-day orbital evolution. If we assume that it did not form in situ, then it might either have formed orbiting the HD 131399BC pair and then been captured by HD 131399A, or it may have formed by conventional CA as a circumbinary planet orbiting HD 131399A, but closer than today, and been subject to dynamical ejection to reach its current remote orbit. In the latter scenario, the binary nature of HD 131399A could have played an active role. Rodet et al. (2017) recently investigated a model to explain the formation of the distant planet (~700 au) imaged around HD 106906. In this scenario, the recently discovered binary nature of HD 106906 (Lagrange et al. 2016) plays a crucial role. The planet first forms at a standard distance in the circumbinary disk. Then it is subject to inward migration due to viscous interaction with the disk, and gets captured in a mean-motion resonance with the binary. The resonant interaction causes an increase in the planet’s eccentricity that drives it too close to the binary at periastron, so that it finally gets ejected by a close encounter. A similar scenario could easily be considered for HD 131399Ab. Rodet et al. (2017) has shown that a key point for this scenario is the eccentricity of the HD 106906 binary. The resonant eccentricity increase of the planet that furthermore leads to a close encounter is a direct consequence of the binary’s eccentricity. With a circular binary, the planet, once captured in resonance with it, would remain there forever. The HD 131399A pair’s high eccentricity obviously reinforces the similarity between the two cases. The main difference is the presence of the outer pair HD 131399BC. The weakest point in the Rodet et al. (2017) scenario is the stabilization of the planet on its current wide bound orbit. Most of the time, the close encounter with the central binary actually leads to a definitive ejection of the planet. Rodet et al. (2017) invoke stabilization via interaction with passing-by stars, and reach the conclusion that it is unlikely although not to be excluded. In the HD 131399 case, the outer pair could easily have played a similar role to the hypothetical stellar fly-bys in the HD 106906 system; the main difference is that the outer pair is known and bound to the system, thus permanently exerting gravitational perturbations onto the planet.

If the central binary itself could have undergone Kozai migration triggered by the outer pair (see above), this would increase the likelihood of the proposed scenario. Resonant capture of the planet could indeed be due to the migration of the planet or to the Kozai-induced migration of the binary. Moreover, the planet eccentricity could have been enhanced by KL cycles induced by the outer pair with no need of resonant capture. Obviously, a dedicated deeper dynamical study is required to investigate the proposed scenario and its variants.

The present-day dynamical interaction of the planet with the outer pair also needs to be investigated. As noted above, no relevant hypothesis about the coplanarity of the various orbits can be made. If the current orbit of HD 131399Ab is significantly inclined with respect to the wide orbit between the two binaries, then the planet could still be subject to KL cycles due to the outer pair. In that case, the planet would have regular high eccentricity phases with repeated periastron passages down to a few au. These episodes could easily have caused the destruction of any remaining circumstellar (CS) disk around HD 131399A and thus explain the current absence of CS material. Assuming the orbital parameters given by Wagner et al. (2016), the period of these KL oscillations should be ~40 000 yr. This is short enough for them to have carved the remaining CS disk.

4. Concluding remarks

HD 131399 is a fascinating complex system with many interrelated dynamical issues. The binary nature of the central star that we evidence in this paper may have played an important role in the generation of the system’s current configuration. Various dynamical studies need to be conducted to investigate its history and its present-day dynamics. Also, many still unknown parameters need to be better constrained to allow a more focused study. A direct detection using interferometers of the HD 131399 A companion, if possible, would help constrain its mass and orbit. Longer term astrometric follow-up of the orbital motion of the planet and of the outer pair will help refine the tentative orbital determination of Wagner et al. (2016) and probe the hypothesis of ongoing KL cycles for the planet.

Note added in proof. Shortly after this letter was accepted for publication, Nielsen et al. (2017) questioned the reality of the HD 131399 b planet and claimed the source is instead a background star. Monitoring of the astrometry, as well as a more detailed spectroscopy, will be decisive in order to disentangle the two hypotheses. In any case, the results on the binarity of the central star presented in the present paper are not affected.

Acknowledgments

We acknowledge support from the French CNRS and the support from the Agence Nationale de la Recherche (ANR-14-CE33-0018), the OSUG LABEX. A.M.L. thanks P. Rubini for his help on the developments of our RV measurement tool SAFIR, and D. Ségransan for providing us with Yorbit.

References

Appendix A: Additional tables

All Tables

All Figures

	Fig. 1 1st row: HARPS RV time series of HD 131399 (the error bars are lower than 30 m/s). 2nd row: Bisector Velocity Span (BVS) as a function of RV. The y-scale is different from that in the first panel. 3rd row: Clean Periodogram of the RV. 4th row: Clean Periodogram of the BVS.
In the text

	Fig. 2 Phase folded fit of HD 131399A RV variations.
In the text

	Fig. 3 Relation between the mass Mp of the companion of HD 131399A and its orbital inclination with respect to the plane of the sky. Blue curve: assuming Mp ≪ M∗; red curve: general case (see text).
In the text