© ESO 2021

1 Introduction

To date, more than 4000 extrasolar planets have been confirmed, and more than 800 of these were discovered with Doppler spectroscopy, also known as the radial velocity (RV) method¹. While most of the discovered planets orbit main-sequence (MS) stars, the number of detections around evolved stars has also risen quickly: since the first discovery of a planet around a giant star in 2001 (Frink et al. 2002), 112 exoplanets orbiting evolved intermediate-mass stars have been published².

Giant stars allow us to extend planet surveys to higher stellar masses, while also still being sensitive in the lower mass regime: whereas intermediate-mass MS stars show only few and rotationally broadened absorption lines, their evolved counterparts, such as K- and G-giant stars, have many sharp spectral lines and are therefore perfect targets for RV measurements. Furthermore, the discoveries of planets around these stars help to improve our understanding of the evolution of planetary systems once the stars evolve into giants (Villaver & Livio 2009; Reffert et al. 2015).

Despite the overall large number of extrasolar planets discovered to date, only a comparably small fraction of 3−4% has been detected in stellar binaries³. In the case of giant stars, this number is especially small, with only five known cases: 11 Com (whose companion falls into the brown dwarf regime, Liu et al. 2008), γ Leo (Han et al. 2010), 91 Aqr (Mitchell et al. 2013), 8 UMi (Lee et al. 2015), and HD 59686 (Ortiz et al. 2016), and only one of these (HD 59686) is a spectroscopic, that is, rather close binary. To a large extent, this small fraction of discovered planets in binary star systems can be explained by the fact that most exoplanet surveys focus on single stars, which is unfortunate since binary systems harboring planets serve as good special cases to constrain the theory of planet formation and evolution.

In this work we study the K giant star ε Cyg, which has been observed spectroscopically for more than 100 yr and is known to undergo large RV changes, hinting at the existence of a close stellar companion to the primary component that is not directly visible. We use our own RV measurements of the star both from the Lick and Stellar Observations Network Group (SONG) telescopes to derive a precise orbit of the spectroscopic stellar companion and to investigate whether additional short-period RV variations are caused by an S-type (i.e., circumstellar) planet in the system.

In Sect. 2 we present the known properties of the ε Cyg system and its primary component. Section 3 describes the RV data sets used in this analysis. Next, in Sect. 4, we perform the RV modeling. In Sect. 5 we then analyze the stability of the system using an N-body code. In Sect. 6 we investigate alternative explanations for the short-period signal. Section 7 summarizes and concludes our analysis.

2 Stellar properties

ε Cyg A is a bright (m_V = 2.48 mag) K0 III giant star. Its HIPPARCOS parallax is 44.86 ± 0.12 mas (van Leeuwen 2007), which puts it at a distance of 22.29 ± 0.06 pc. It has been observed spectroscopically since the beginning of the 20th century (see e.g., Campbell & Moore 1906; Kustner 1908; McMillan et al. 1992), and has long been known to host a spectroscopic binary companion, for which Gray (2015) derived an orbital period around 55.1 yr without putting an error on that value. The spectroscopic companion has never been observed visually though, and the spectra do not show any evidence of a second set of spectral lines nor any temporal changes of line asymmetry (Gray 1982). McMillan et al. (1992) also detected short-period RV variations on a time scale on the order of a few 100 d with much smaller amplitude, which they did not analyze in detail, but simply compared them qualitatively to other examples of K-giants with fast variations. Further RV measurements of the ε Cyg system were performed from 1999–2010 by Gray (2015), who present a spectroscopic analysis of the binary system. As the formal uncertainties of these measurements are too large to clearly identify the short-period signal, Gray (2015) only reference the data taken by McMillan et al. (1992) and propose that the variations might be caused by stellar activity, modulated by the rotation of the star.

In addition to the spectroscopic companion, ε Cyg A has an optical companion (ε Cyg B) at a separation of 71′′ (Montes et al. 2018), whose Gaia DR2 parallax and distance are 4.35 ± 0.02 mas and 229.9 ± 1.3 pc (Gaia Collaboration 2018), respectively, and which is therefore not gravitationally associated. A second optical companion C has a separation of 78′′ (Montes et al. 2018), and its Gaia DR2 parallax is 45.51 ± 0.03 mas, putting it at a distance of 22.97 ± 0.01 pc and therefore very close to ε Cyg A (Gaia Collaboration 2018). According to Montes et al. (2018), it is an M4V+ star, and as its proper motion also resembles that of the primary star (ε Cyg C: + 354.62 ± 0.04 mas/yr in RA, + 329.18 ± 0.05 mas yr⁻¹ in declination; Gaia Collaboration 2018; ε Cyg A: + 355.66 ± 0.08 mas yr⁻¹ in RA, + 330.60 ± 0.09 mas yr⁻¹ in declination; van Leeuwen 2007), it can be expected to be a very wide physical companion. We do not pay any further attention to it within this work though, as its separation from the primary at the distance of the system is at least 22.29 pc ⋅ 78 mas = 1739 AU, and any influence therefore can be neglected.

We have two independent sources for the stellar properties of ε Cyg A: the first one comes from Stock et al. (2018), who used a Bayesian interpolation scheme for evolutionary tracks. This method provides a probability for the star to either fall onto the red giant branch (RGB) or onto the horizontal branch (HB). It also delivers nonsymmetrical probability density functions (PDFs), which allows us to compute asymmetric 1σ confidence intervals for each of the stellar parameters. The mode values of the PDFs are used as the most probable values. According to this method, ε Cyg A is most probably an HB star (P = 99.5%), with a stellar radius of and a mass of . This mass estimate is smaller than the numbers that are usually adopted in the literature, which are closer to 2 M_⊙ (see e.g., Gray 2015). It is possible that previous authors assumed ε Cyg A to be an RGB star and thus overestimated its stellar mass. This is plausible considering the (very unlikely) RGB solution by Stock et al. (2018), which yields a mass of , somewhat closer to the past estimates⁴.

In order to get a second independent measurement for the stellar properties, we also performed an asteroseismic campaign on ε Cyg A with SONG, and derived its mass from the frequency of maximum power ν_max of the star, using the relation as in Stello et al. (2017). We used the same approach as in Arentoft et al. (2019), who analyzed a similar SONG data set for the red giant ε Tau, to derive ν_max and its uncertainty. In short, following the methods described in Mosser & Appourchaux (2009) and Stello et al. (2017), we applied a Gaussian fit combined with a linear (background) trend to the oscillation signal in the power spectrum to determine ν_max. The uncertainty was found by performing the same fit to slightly modified versions of the power spectrum, in each of which a single oscillation mode had been subtracted, and using the variations in the determined ν_max -values to estimate the uncertainty; see Arentoft et al. (2019) for details. The analysis resulted in a ν_max = 32.16 ± 0.81 μHz, which gives a mass estimate of M_⋆ = 1.103 ± 0.042 M_⊙; the 1σ-error bars of this result overlap with the ones of the HB mass estimate from Stock et al. (2018), and we can rule out the possibility of ε Cyg A being an RGB star with high certainty. Given the robustness and high fidelity of the asteroseismic method, we adopt that mass estimate for our analysis. Table 1 lists all the stellar parameters.

3 Observations

From 1999 to 2012, we monitored a sample of 373 G- and K-giant stars and measured their RV variations, using the Hamilton Échelle Spectrometer at the Lick observatory in California, USA. The spectrograph was fed by the 0.6 m Coudé Auxiliary Telescope (CAT) and covers a wavelength range of 3755–9590 Å at a resolving power of R ~ 60 000. High RV precision is achieved using the iodine cell method as described by Butler et al. (1996). Our observations led to several planet detections (e.g., Frink et al. 2002; Reffert et al. 2006; Trifonov et al. 2014), among them a massive circumprimary planet in a close and eccentric binary system (Ortiz et al. 2016). The measurements also proved valuable in gaining a statistical understanding of the properties of extrasolar planets around giant stars (Reffert et al. 2015).

The K-giant star ε Cyg (HIP 102488) is a member of our sample, and we collected a total of 109 RV measurements at Lick observatory. They cover a time span of more than eleven years (from June 2000 until November 2011), and their median measurement precision is 4.8 m s⁻¹.

In addition to the Lick data set, we obtained a total of 5 272 RV measurements of ε Cyg with the 1 m robotic SONG telescope on the island of Tenerife, Spain (Andersen et al. 2014; Fredslund Andersen et al. 2019). A set of 5 063 of these individual measurements were part of an asteroseismic campaign to better constrain the stellar properties of the K giant (see Sect. 2); these observations were therefore performed at very high cadence, with about 100 to 150 measurements per night, for a little more than a month from June 20 until July 27, 2018. From these measurements we can get an estimate of the short-term RV jitter, that is, the expected stochastic astrophysical stellar RV variation. Data from a few nights of this campaign are shown in Fig. 1; the stellar oscillations are clearly visible, with peak-to-peak variations of ~ 40−50 m s⁻¹ on timescales of several hours, in agreement with the value for ν_max found above. The observed oscillations are a combination of a number of individual oscillation modes, which together give rise to a short-term variation in the time series with a standard deviation of ~ 11 m s⁻¹. This scatter agrees well with other stars with similar maximum-power frequencies ν_max as ε Cyg (compare Yu et al. 2018). Following Kjeldsen et al. (2008), we isolated the p-mode signal and determined the amplitude per radial mode to be 1.45 ± 0.06 m s⁻¹. That value is approximately 9 times larger than the amplitude of the oscillation signal seen in the Sun, and accordingly the standard deviation and peak-to-peak variations of our asteroseismic measurements are roughly one order of magnitude larger than the solar values (see Kjeldsen et al. 2008; Fredslund Andersen et al. 2019).

For our orbital analysis of the system this large number of data points is unsuitable, as it puts too much weight on the SONG measurements and requires too much computational power when performing orbital fits. For this reason, we computed the median RV for each night with more than one observation. Our updated SONG data set then consists of 228 RV measurements, which fall between April 2015 and December 2018. They cover exactly the last periastron passage of the eccentric stellar companion to ε Cyg, which occurred around February 2017, and thus enable us to determine the binary orbit with very high precision. Just as the Lick measurements, the high-resolution spectra taken by SONG are calibrated by the means of an iodine cell (Grundahl et al. 2017), and the measurements reach a similar precision, with a median of 2.3 m s⁻¹.

In addition to the long-period signal induced by the close stellar companion, both the Lick and SONG data sets show RV variations with periods just below 300 d (see Sect. 4.2), which could be caused by a planetary companion orbiting the primary star in an S-type configuration. To gain a better coverage of the orbits of the binary and the possible planetary companion, we complement our own two data sets with the older measurements from McMillan et al. (1992). The McMillan data set consists of 213 individual RV measurements of ε Cyg, collected between May 1987 and May 1992 with an interferometer at the 0.9 m-telescope of the Steward Observatory at Kitt Peak, USA. Calibrations were done with emission lines of a Fe-Ar hollow cathode lamp, yielding a precision of 12 m s⁻¹ for all individual measurements.

Finally, we use six RV measurements from Kustner (1908) to cross-check our results for the orbital solution of the long-period stellar companion to ε Cyg. The observations by Kustner were carried out in Bonn, Germany, with a 30 cm refractor and a spectrograph constructed out of three prisms, and RVs were computed by comparing the positions of the stellar absorption lines to a reference iron spectrum. The observations of ε Cyg cover the time between July 1904 and September 1906. Their errors are large in comparison to modern values, and vary between 1000 and 2100 m s⁻¹, so we do not include these measurements in our fits. The merit of this data set however lies in the fact that it seems to record the periastron passage of the stellar companion two full periods back, and we can overplot and check our results.

We decided to refrain from using the RVs collected by Gray (2015) from 2001 until 2010 as they cover mostly the same time as our own Lick data and have much larger measurement errors, which is why they do not help to constrain our models. All in all, our combined data set for the fitting procedures consists of 550 RV measurements from McMillan, Lick and SONG.

Fig. 1 Uncorrected SONG RV measurements from the last week of our asteroseismic campaign on ε Cyg. Oscillations with periods shorter than a day are clearly visible.

4 Analysis of the RV data

4.1 Determining the orbit of the close stellar companion

We use a Keplerian model in combination with a Levenberg-Marquardt (LM) minimization scheme to fit the combined RV measurements of ε Cyg and determine orbital parameters for the spectroscopic companion, which we hereafter also refer to as the close stellar companion; the main host star is denoted as in the literature, ε Cyg A. The Keplerian fit incorporates eight free parameters: RV semi-amplitude K, orbital period P, eccentricity e, argument of periastron ω, and mean anomaly M₀ of the orbit, as well as zero-point offsets for the McMillan, Lick and SONG data sets, respectively. In order to account for the RV jitter of the star, we quadratically added a fixed value of 20 m s⁻¹ to all individual measurement errors (see Sect. 4.3). Error estimation was performed with a Markov-chain Monte Carlo (MCMC) approach, using the emcee python package (Foreman-Mackey et al. 2013): we constructed 32 walkers (four times the number of free parameters) with initial values drawn randomly from a small Gaussian ball around the best-fit solution from above, and let the sampler run for altogether 6000 steps (where we discarded the first 1000 steps). This method produced 36 649 unique samples of the parameters; the 1σ intervals of these parameter sets around the mean values of the posterior distribution serve as our uncertainties.

The LM fitting scheme delivers a best-fit model with an orbital period of 19 619 d at , which lies within the errors from the mean period of the MCMC posterior distribution of (see column “Single-Keplerian” of Table 2). Our results therefore are roughly one-and-a-half years less than the orbital period of about 55.1 yr derived by Gray (2015), who used 53 RV measurements taken between 2001 and 2011 and complemented them with the McMillan data set as well as data by Griffin (1994).

When plotting our best-fit solution with a period of 19 619 d over the measurements taken by Kustner (1908), one clearly notices that the last RV from that data set falls far away from the curve, by a value of about 7400 m s⁻¹, which corresponds to more than three times the formal error of 2100 m s⁻¹ (compare Fig. 2). If we fit all four data sets combined, the best Keplerian model delivers an orbital period of 20 175 d and a of 1.95. In this model, the former outlier now also agrees with the fit curve within its 1σ-uncertainty and falls exactly at the sharp peak of the curve; however, the fit leads to systematic (albeit small) deviations from the RVs of the three modern data sets by McMillan, Lick and SONG. As our further analysis focuses on these measurements, we decide to adopt the model with the shorter orbital period of 19 619 d for now. We do not exclude the possibility that the result for the orbital period of the close stellar companion might still change on the level of a few percent with future observations coming in, but given our current data it seems most probable that the last Kustner measurement is actually an outlier.

Using the primary mass of M_⋆ = 1.103 ± 0.042 M_⊙ and solving the posteriorresults of the mass equation f(m) numerically for fixed inclinations i = 90°, we derive a minimum mass of 278 ± 7 M_jup = 0.265 ± 0.007 M_⊙ for the close stellar companion. Our estimate places it therefore well above the brown dwarf regime. As the spectroscopic companion has never been imaged directly to our knowledge, and its absorption lines are not visible in the spectra, its luminosity must be much lower than that of the primary star ε Cyg A. If we additionally assume that both stars have been formed around the same time, as is expected for binaries, this opens up two possibilities for the nature of the companion: the first is a white dwarf, which would mean that it was originally the primary, more massive component in the system and therefore evolved more quickly.

The second possibility is that the spectroscopic companion to ε Cyg A is a main-sequence star with a mass lower than that of the primary, therefore evolving more slowly and being much less luminous. This places it anywhere between an early M-dwarf to an early G-type star, and constrains the possible orbital inclination to values larger than ~ 14°, to keep its mass smaller than the mass of the primary, ~1.103 M_⊙. Therefore the apparent brightness in the V band should be between approximately 6 and 12 mag.

Griffin (1994) speculated about the possibility of imaging the companion directly (either at visual wavelengths or the infrared) and predicted that its angular separation on the sky could get nearly as large as 2′′ during apastron, but he had no certain solution for the orbit yet and used a semi-major axis of approximately 20 AU and eccentricity of 0.9 for his calculations. Our orbit solution puts the apastron at 30.5 AU from the primary. At the distance of the system, 22.1 pc, this would translate to a maximum projected separation of 1.38′′ in the case of an optimal position of the system on the sky. However, we were also able to derive a solid estimate of the argument of periastron of the orbit, which is at 275° and therefore places the periastron of the orbit nearly exactly in-between the observer and the primary component, and the apastron behind the primary ε Cyg A as seen from Earth. Hence the angular separation at apastron will only become large for a low inclination of the orbit. For an inclination higher than 80°, the maximum separation of the two components would be reached around the co-vertices of the orbital ellipse (semi-minor axis: b ≈ 5.8 AU, projected separation: 0.26′′).

Modern direct imaging instruments offer the capability to resolve the two stars: as an example we used the exposure time calculator of the ESO instrument SPHERE-IRDIS⁵ to calculate the achievable contrast for ε Cyg within an exposure time of 1800 s with median weather conditions, which resulted in a magnitude difference of ~ 8.4 mag at a separation of 0.1′′, and ~ 14 mag at the co-vertex separation of 0.25′′. With the magnitude difference of the two stars probably not being higher than 12 mag in the visual, it might be an interesting option to try and image the spectroscopic companion directly. According to our calculations, it will reach the co-vertex around JD 2459797, that is, in August 2022, which is also the best time within the year to observe ε Cyg directly from Cerro Paranal Observatory.

Fig. 2 Top plots: RV measurements of ε Cyg by Kustner (1908) (left), and McMillan et al. (1992), Lick observatory and the SONG telescope (right), along with two orbital models of the binary companion: Fit 1 (black line) is the best-fit solution for the three modern data sets of McMillan, Lick and SONG, with an orbital period of 19 618.6 d; Fit 2 (gray line) was derived from fitting all four data sets and yields a period of 20 175 d. Middle plots: residuals of the RV measurements after removing the orbital solution of Fit 1. The last measurement of Kustner then is a clear outlier. Bottom plots: residuals of the RV measurements after removing the orbital solution of Fit 2. Now all the Kustner RVs agree with the model within their 1σ-uncertainties.

4.2 GLS periodogram of the RV residuals

We aim to find relevant short-period variations in the RV measurements of ε Cyg by calculating generalized Lomb-Scargle (GLS) periodograms as described in Zechmeister & Kürster (2009). As the RV data are largely dominated by the signal of the stellar companion with a semi-amplitude of K ≈ 4.6 km s⁻¹, we first subtracted the best-fit orbital solution of the binary from the measurements and then calculated the periodograms of the residuals. Figure 3 shows the significance of variations with periods between 100 and 500 d, for the three individual data sets by McMillan, Lick, and SONG (lower three panels), as well as for the combined measurements (uppermost panel). False-alarm probabilities(FAPs) of 5, 1, and 0.1% are depicted by the gray dash-dotted, dashed, and dotted lines, respectively.

All three data sets show a highly significant peak at periods slightly shorter than 300 d that greatly exceeds the FAP of 0.1%, but the exact positions of those peaks vary: in the McMillan data the highest peak lies at 277.9 ± 2.8 d, the Lick data set shows its strongest signal at 287.3 ± 1.7 d, and the SONG data set at 291.2 ± 3.0 d. The uncertainty estimation used by Zechmeister & Kürster (2009) is based on the curvature of the peaks in the periodograms (as described in Ivezić et al. 2014), and in simulations, using fixed periods with varying white noise amplitudes and the sampling of the original data sets, we found a good agreement with these errors.

The period of the most dominant signal changes by more than 13 d over the course of 30 yr from the first until the last measurement. All the peak periods lie at least 1σ away from each other; in the case of the McMillan and Lick periods the separation is even larger than 3σ. The change of the period over time therefore seems to be quite significant, and we examine that behavior more closely in Sect. 4.4.

The GLS periodogram of the combined data sets shows a forest of peaks around periods of 300 d, with the strongest one at 291 d. The side peaks can at least partly be attributed to the sampling of the data: Between the individual data sets there are gaps of several years without any data, which can produce such patterns (see e.g., VanderPlas 2018).

Fig. 3 GLS periodograms of the residuals of the RV measurements after removal of the binary signal. The periodogram of the combined data sets (top panel) shows two peaks with similar heights at ~ 282.5 d and ~ 290.8 d. The gray dash-dotted, dashed and dotted lines represent false-alarm probability (FAP) levels of 5, 1 and 0.1%, respectively.

4.3 Fitting a double-keplerian model to the data

We first investigate the hypothesis that this short-period signal is caused by a possible additional companion of planetary nature in the system orbiting the main component. Therefore we use a double-Keplerian model to constrain the orbits of the known stellar component and the possible second companion simultaneously. The fit then incorporates 13 free parameters: RV semi-amplitude K, orbital period P, eccentricity e, argument of periastron ω, and mean anomaly M₀ for each of the companions, as well as zero-point offsets for the McMillan, Lick and SONG data sets, respectively. Again we estimated the errors through the MCMC method, similar to the single Keplerian model, but now with 52 walkers. The 1σ-intervals of the 62 535 samples around the mean values of the posterior distribution serve as our uncertainties.

The best-fit results from the model and the mean orbital parameters from the MCMC posterior distribution are shown in Table 2 (columns “Double-Keplerian”), and they agree within the errors. The parameters of the stellar companion only change slightly as compared to the single-Keplerian model. The best-fit orbital period of the putative planet is 291.2 d, and it has a low eccentricity of ~ 0.173. Using the primary mass of M_⋆ = 1.103 M_⊙ puts the minimum mass of the planetat about 1 M_jup, and its semi-major axis at 0.89 AU.

Without taking any jitter into account, the of the best fit is 37.49. When we fit our Lick data set alone, letting all parameters vary, a jitter value of 17.4 m s⁻¹ brings the of that fit down to unity. This value is larger than the jitter estimate from the asteroseismic measurements (see Sect. 3), which is no surprise as the Lick data set covers a much longer baseline. Stellar jitter stems from a number of stellar phenomena, and oscillations are just one of those – others are granulation, stellar rotation, and magnetic cycles (Dumusque 2016; Dumusque et al. 2017). These jitter sources act on considerably different time scales (few hours to several years), which might well explain the discrepancy between the jitter estimates for the whole Lick data set (covering 11 yr) and the asteroseismic data (~1 month). For consistency, similarly as for the Lick data we also fitted only the SONG measurements, but left out the asteroseismic data, which resulted in a jitter estimate of 15.0 m s⁻¹ to achieve . This value falls in-between the jitter estimates from the longer-baseline Lick data and the shorter-baseline asteroseismic measurements.

We can also compare these results to the color-dependent distribution of jitter estimates found for the stars in our own K-giant sample (see Frink et al. 2001; Trifonov et al. 2014): For the color of ε Cyg, (B − V) = (1.04 ± 0.01) mag, the measured values scatter from about 10 to 30 m s⁻¹, so all of the jitter results from above lie within the distribution. In the end we decided to use a fixed jitter value of 20 m s⁻¹ for all our Keplerian and dynamical fits in this analysis. This corresponds roughly to the mean of the above-mentioned jitter distribution for stars of similar color(B − V), and allows to account for additional jitter sources that only come into play over the 32-yr long baseline of the three data sets McMillan, Lick and SONG combined. The double-Keplerian fit to all these measurements then delivers a of 1.21, which is a clear improvement to the single-Keplerian model with from Sect. 4.1. The combined rms of all data sets is 23.85 m s⁻¹, which corresponds roughly to our adopted jitter value.

Figure 4 shows the RVs of the three data sets along with the best Keplerian fit to the data phased to the best-fit period P = 291.2 d (top) and the residuals of the RVs from the fit (bottom). While many of the data points lie close to the best-fit solution, there are a number of McMillan and Lick measurements that fall far away from the curve. They form an additional “branch”; they are systematically shifted in phase and offset in vertical direction with respect to the fit. The amplitude of the best-fit curve is a little too small to perfectly match the distribution of the SONG RVs, but overall the SONG data set is much better represented by the model than the other two. The fact that the SONG data dominate the fit is not surprising as this is the largest data set with the smallest measurement uncertainties. However, the outliers especially in the McMillan and Lick data sets show that the double-Keplerian model does not provide a satisfactory fit to the combined data. A possible explanation for that would be a temporally changing RV signal; as the Keplerian model assumes undisturbed orbits, changes of the Keplerian elements are not incorporated.

Fig. 4 Top panel: RVs of the three modern data sets McMillan, Lick and SONG, phase-folded by the best-fit period of the putative planet P = 291.2 d, are shown in symbols along with the RV curve of the best fit as a black line. Bottom panel: residuals of the individual data points.

4.4 Investigating the temporal evolution of the Keplerian signal

To further investigate the problems of the Keplerian model, we examine the three data sets individually by fitting a double-Keplerian model to each one of them separately and comparing the resulting parameters. By this we are hoping to find out whethera changing RV signal is the cause for the bad fit of the combined model. In a second step, we split each data set into shorter sections and fit each of these separately, in order to follow changes of the Keplerian elements over even shorter times. It is important to assure a large number of data points that is more or less evenly sampled as well as a time coverage that is longer than the period of the signal for each of the sections. The McMillan data set covers 5 yr with 213 measurements, but the measurement uncertainties are comparatively large (12 m s⁻¹), and the observations were mostly taken in chunks within a few days, leaving gaps of several days or weeks in between. This leads to a potentially poor phase coverage of the period, which is why we split the McMillan data set into only two sections of 106 and 107 data points, respectively. The Lick observations cover the longest time span of the three data sets used here, but also contain the smallest number of measurements. Therefore we also split them into two sections of 55 and 54 data points, respectively. The SONG data set in contrast consists of a large number of measurements with small uncertainties, which allowed us to split it into three sections: The first consists of 46 data points and covers the time right before the very quick RV change during periastron passage of the close stellar companion occurs; the second is made up of the 128 data points taken during the close approach of the two stars; and the third consists of 54 data points and focuses on the time after periastron passage.

Each single fit to a complete data set incorporates 11 parameters, two times five Keplerian elements for the stellar and the planetary companion, and one RV zero point. When fitting the shorter sections of the data sets however, the long-period orbit of the spectroscopic companion is not constrained sufficiently, due to the smaller number of data points and shorter time spans; therefore, inthe models for the individual sections, we left the Keplerian parameters of the stellar companion fixed at the best-fit results from the whole respective data set. This leaves 6 free parameters for each fit to a section: 5 Keplerian elements for the planetary companion and one RV zero point. Uncertainties on the free parameters were again computed from an MCMC analysis for each data set and each section, in analogy to the single and double-Keplerian models to the combined data set.

Each model to a subset of the full data, either individual data set or even shorter section, results in a much better fit than the complete double-Keplerian model, with only very few random outliers. The for the McMillan, Lick and SONG data sets are 0.96, 0.77, and 0.47, respectively. Figure 5 shows the evolution of the orbital period P and the semi-amplitude K of the short-period signal over time: the solid lines depict best-fit results of the individual data sets, while the dashed lines denote the results from fitting the sections. In both cases distinct changes of the parameters over time become apparent: The period clearly follows the pattern already observed in the GLS periodograms of the three data sets, changing from 278 d in McMillan over 286 d in Lick to 296 d in the SONG data. The fits to the shorter sections follow a similar trend, but show more scatter. Some of that variance might be explained by the fewer measurements within each section and the shorter time spans covered by the data. Especially in the case of the SONG data set, the first and the last section only contain around 50 measurements, and none of the three sections covers two full periods of the short-period signal.

The semi-amplitude K of the signal in contrast stays mostly constant for McMillan and Lick around a value of 26− 27 m s⁻¹, before it drastically increases to 53 m s⁻¹ in the SONG data, effectively doubling. Again, the fits to the sections show some more scatter but generally follow the same trend. Interestingly, of the three SONG sections the second one, which covers the time around periastron, results in a much larger K than the other two. This might hint at an actual physical process that increases the semi-amplitude exactly during periastron passage, but could also be simply explained by the fitting problems described above. According to our calculations, for the edge-on configuration the gravitational redshift z_grav only changes by 0.6–0.7 m s⁻¹ over the eccentric orbit of the stellar companion, so it fails as an explanation for the varying short-period semi-amplitude.

Fig. 5 Orbital period P (top) and semi-amplitude K (bottom) of the presumed planet, derived from double-Keplerian fits to the individual data sets (solid lines) and to shorter sections of the data (dashed lines). The gray-shaded areas denote the time spans from the first until the last measurement within each section.

4.5 Fitting a dynamical model to the RV data

One possibility to explain the changes of the Keplerian elements over time would be the gravitational interaction of the planet in question with the close stellar companion: the double-Keplerian model assumes undisturbed orbits around the center of mass and neglects any other forces that might be present. However, for two bodies in such close proximity as the planetary and stellar companion in our model, and with such high masses, this assumption certainly does not hold true anymore, and the orbits of the two bodies would be subject to considerable changes. This is even more true because our total observational time span covers more than 30 yr by now, over which any evolution of the orbits can be traced. Furthermore, the SONG measurements were taken around the periastron passage of the close stellar companion, where all three bodies come closest to each other and the biggest changes can be expected.

To take this effect into account, we attempt to better fit the data using a fully dynamical model, as done before for HD 59686 in Trifonov et al. (2018). It consists of a modified version of the Bulirsch-Stoer N-body integrator in the SWIFT package, that outputs the RV data of the primary component in the system, and a Levenberg-Marquardt minimization algorithm. As for the double-Keplerian model, the free parameters comprise the RV semi-amplitude K, orbital period P, eccentricity e, argument of periastron ω, and mean anomaly M₀ for each of the companions, as well as zero-point offsets for the McMillan, Lick and SONG data sets, respectively.

Additionally, in a dynamical model the data can potentially constrain the inclinations of the orbits relative to the sky plane i_1,2 and the longitudes of the ascending nodes Ω_1,2 of the close stellar companion (subscript 1) and the putative planet (subscript 2). However, it is not feasible to keep both parameters free in our case as there are many local minima in the log(likelihood) plane, and we do not have any independent knowledge of the true orientation of the orbits that we could use to find suitable starting parameters. Therefore, we applied our model only to strictly coplanar prograde orbits, where i₁ = i₂ = 90° and Ω₁ = Ω₂ = 0°, and strictly coplanar retrograde orbits, that is, i₁ = i₂ = 90° and Ω₁ = 0°, Ω₂ = 180°.

Our best dynamical models however show no improvements to the double-Keplerian fits, with a of 1.42 and 1.23 for the prograde and retrograde configuration, respectively. Both models are still left with many systematic outliers from the best-fit curve, so it seems as if the changes of the Keplerian parameters over time cannot be fully accounted for by the expansion to a dynamicalapproach. When integrating the two models further over time, with a modified version of the Wisdom-Holman N-body integrator in the SWIFT package, we find stability time of only 85 and 454 yr for the prograde and retrograde configuration, respectively, after which the planetary eccentricity becomes larger than 1 and the planet is kicked out of the system. The long-term survival of the putative planet therefore seems highly questionable, and we present a more in-depth analysis of the dynamical stability in Sect. 5.

Figure 6 shows the GLS periodograms of the residuals of the dynamical models and the RVs, for both the prograde and retrograde configuration and for all individual data sets as well as the combined data. In all periodograms there are still peaks left that exceed the FAP level of 0.1% and mostly lie at periods close to the fitted orbital periods just below 300 d, which illustrate the above-mentioned systematics in the residuals. It is evident that the dynamical models do not fully describe the signals present in the RV data.

4.6 The co-orbital scenario

In this section we investigate the source of the significant peak of the RV residuals after subtracting the possible planet signal in the context of the potential presence of a co-orbital body. Co-orbital configurations have not yet been detected outside of the Solar System although several candidates have recently been published (of particular interest is the case of TOI-178, Leleu et al. 2019; but see also Hippke & Angerhausen 2015; Janson 2013 or the TROY project⁶ by Lillo-Box et al. 2018a,b). Interestingly, co-orbital planet pairs are stable under a very relaxed condition developed in Laughlin & Chambers (2002), stating that such configurations would remain long-term stable as long as the total mass of the planet and its co-orbital companion is smaller than 3.8% of the mass of the star.

The motivation of this analysis for our system comes from the theoretical study of Leleu et al. (2015), who derived the implications on the RV periodogram of a co-orbital planet. In such case, the long-term libration of the co-orbital motion would introduce two signals in the periodogram at n ± ν, where n is the frequency of the orbital period and ν corresponds to the libration frequency. However, the sampling and precision of the data can hide these two signals in the periodogram. Leleu et al. (2015) propose a technique to enhance their detection through the so-called demodulation technique: the RV residuals are convolved with a sinusoidal function of same frequency as the carrier (the Keplerian frequency of the planet), reproducing the peak at the libration period. The libration period depends on the Trojan, planet and star masses so a minimum (m_t = 0) and maximum (m_t = m_p) libration period can be estimated. We applied this technique and obtained the periodogram of the demodulated RVs, which is shown in Fig. 7. The result displays a peak in the expected period range of possible libration frequencies, which opens up the possibility of a librating co-orbital body to be the cause of these RV residual peaks. A more in-depth analysis is needed to properly investigate this possibility, which is beyond the scope of this paper. It is also necessary to point out that the issue of instability of a single circumprimary planet, which we examine in depth in Sect. 5, potentially also poses a problem for the co-orbital scenario; further investigations into this direction therefore require a thorough dynamical study.

Fig. 6 GLS periodograms of the residuals of the best dynamical models for the prograde (left) and retrograde configuration (right). The top-most plot shows the periodograms of the residuals of the combined data, the bottom three of the individual data sets. FAPs of 5, 1 and 0.1% are indicated by the gray dash-dotted, dashed and dotted lines, respectively.

Fig. 7 Lomb-Scargle periodogram of the radial velocity residuals (blue) and the corresponding demodulated radial velocities (red). The red vertical dotted lines correspond to the range of allowed libration periods for this star-planet system based on a minimum (mt = 0) and a maximum (mt = mp) mass for theco-orbital body.

5 Dynamical stability analysis

5.1 Theoretical considerations about the orbital configuration

An S-type planet around ε Cyg A would be subjected to substantial gravitational forces from the close stellar companion, which could potentially alter its orbit considerably over time, as is explained above in Sect. 4.5. Even though the semi-major axis of the orbit of the stellar companion is a₁ = 15.8 AU, due to its high eccentricity the periastron distance is only q₁ = 1.11 AU, which is just 0.22 AU larger than the proposed semi-major axis of the planet, a₂ = 0.89 AU. According to Hamilton & Burns (1992), the Hill radius of the stellar companion can be approximated as (1)

with m₁ = 0.265 M_⊙ being the minimum mass of the stellar companion, m₀ = 1.103 M_⊙ the mass of the primary, and a₁ and e₁ the semi-major axis and the eccentricity of the orbit of the stellar companion. For the best coplanar double-Keplerian model the orbital path of the putative planet thus would pass through the Hill sphere of the close stellar companion during its periastron passage, making its long-term survival very questionable. For low inclinations and thus higher masses of the stellar companion, the problembecomes even worse as the Hill sphere becomes larger. However, despite the small distance of the orbits during periastron, the two companions could in principle always be much farther apart than the Hill radius of the stellar companion, as they do not necessarily have to have the same true anomalies at the time of closest approach.

We know of a number of extreme multi-companion systems today that are dynamically stable because their orbits are locked in very specific configurations, for instance, secular alignments of their periastra and/or mean-motion resonances (MMR), such as HD 59686 (Trifonov et al. 2018) or HD 82943 (Tan et al. 2013). In order to gain a better understanding of the possible orbital configurations of the ε Cyg system, we ran a comprehensive dynamical analysis: using a modified version of the Wisdom-Holman N-body integrator in the SWIFT package, we tested the temporal evolution of many different orbital configurations by varying the starting parameters in a broad range around our best-fit results. First we tested both prograde and retrograde coplanar orbits, with a special focus on configurations that are locked in secular apsidal alignment or mean-motion resonance. In a second step we extended our analysis to mutually inclined orbits, for which we concentrated on regions of the parameter space near the fixed point of the Kozai-Lidov mechanism as these would offer the greatest chances for long-term stability (for more details see e.g., Kozai 1962; Lidov 1962; Lithwick & Naoz 2011).

Fig. 8 Survival times for test particles placed into the ε Cyg stellar binary, in prograde (left) and retrograde coplanar configuration (right). The eccentricity of the orbit of the stellar companion and the initial semi-major axes of the test particles were varied in order to find the stability boundary. The position of the best-fit solution from Sect. 4.3 is plotted by the white diamond; its size is at least by a factor of 10 larger than the uncertainties on the Keplerian parameters.

5.2 Stability analysis of coplanar orbits

We set up a system of two bodies, using the derived masses of the primary and secondary stellar components as well as the semi-major axis of the secondary from Sects. 2 and 4.3, and inserted a variety of massless test particles around the best-fit orbital solution of the planet. As starting parameters, we chose the periastron of both the test particles and the stellar companion to be at 0° (so Δω = 0°), and the mean anomalies of the test particles to be 0°, while that of the stellar companion was set to 180°. The semi-major axes of the test particles were varied within a range of 0.06 AU around the best-fit solution of the planet, and their eccentricities varied between 0 and 0.8. In order to understand the influence of the high binary eccentricity on the inner orbits, these simulations were repeated for eccentricities of the stellar companion ranging between 0.8 and the best-fit value of 0.93. The whole analysis was done both for prograde and retrograde orbital configurations, and a system was considered stable if the test particle stayed within a critical distance interval {a_c,1, a_c,2} from the primary component for the whole integration time of 10⁴ orbits of the stellar companion (~5.4 × 10⁵ yr). The inner boundary a_c,1 = 0.22 AU was chosen to be the distance at which one orbit of a test particle is only resolved by 10 integration time steps (the time step of the simulation being 3.4 d), the outer boundary a_c,2 corresponds to the semi-major axis of the stellar companion.

As was expected from our considerations in Sect. 5.1, the system generally shows a very high degree of chaotic behavior. The most stable coplanar configurations were achieved for retrograde MMR or secularly aligned orbits with test particle eccentricities around 0.4. A 1:n MMR is established for systems where at least one of the MMR angles (2)

librates around a constant, with λ_1,2 being the mean longitudes of the stellar companion and the test particle, respectively. For prograde orbit, n is positive and m = 1, …, n. For retrograde orbit, n is negative and m = 1, …, |n| + 2. For secular alignment the angle between the pericenters of the two orbits, (3)

librates around zero. For all our simulated systems we checked whether these requirements were fulfilled; we were thus able to identify the MMR and secularly aligned configurations. Nevertheless, even for those configurations general long-term stability only occurred for eccentricities of the stellar companion being smaller than 0.84 (for prograde orbits) or smaller than 0.9 (for retrograde ones). At the best-fit binary eccentricity of 0.93 the instability timescale is on the order of only 10 binary periods (≅500 yr) in the retrograde case and on the order of just 1 binary period in the prograde case. Altering the semi-major axes of the test particles within the range defined above has barely any effect on the stability.

We redid the simulations for 1000 test particles with initial semi-major axes evenly distributed between 0.65 and 1.1 AU and fixed eccentricities at 0.1, while varying the binary eccentricity between 0.8 and 0.93 in steps of 0.01, to illustrate our general findings for the pro- and retrograde case in Fig. 8. Here it becomes clear that orbiting closer-in to the primary component improves the stability of the test particles, but the position of the best fit (indicated by the white diamond) is far away from the stable region even in the retrograde configuration. With the eccentricity of the orbit of the stellar companion being our best-constrained parameter (e_AB = 0.9295 ± 0.0003 in the single-Keplerian model), there is no doubt that the putative planet falls into a highly unstable region of the parameter space.

5.3 Stability analysis of mutually inclined orbits

From our considerations in Sect. 4.1 we know that the inclination of the binary orbit must be larger than 14°, but apart from that we do not have any constraints on the orbital inclinations. Therefore, in order to test the stability of mutually inclined orbits, we decided to focus on the region of possible Kozai librations within the parameter space, which we deem most promising to guarantee long-term stability. We ran simulations for retrograde orbits only (these being more stable than prograde orbits), varying the mutual inclination Δi between 90° and 180° in steps of 5°, with ω₂ = 90° (this corresponds to the stable Kozai libration regime) and ΔΩ = 0° (to avoid close encounters). A smaller time step of 0.35 d was chosen, which puts the inner boundary (for at least 10 time steps per orbit) at a_c,1 = 0.055 AU and therefore allows the possibility of highly eccentric orbits of the test particles, as are expected to occur in Kozai oscillations. The outer boundary a_c,2 again is equal to the semi-major axis of the stellar companion. This time we did not test different binary eccentricities, but only used the best-fit value of 0.93 in our simulations.

The most stable solutions were found for mutual inclinations 135° < Δi < 140°, but even for those configurations the instability timescale is on the order of 10 binary periods (≅ 500 yr), after which the test particles are ejected from the system; only a small percentage of systems survive longer than 100 binary periods, but fully stable regions do not appear in the parameter space. We also did not find any signs of Kozai oscillating orbits; if they exist, they must be very short-lived and are unlikely to provide long-term stable solutions.

6 Possible alternative explanations for the short-period RV variations

6.1 Hierarchical triple

In a series of publications (Morais & Correia 2008, 2011, 2012) it was shown that RV variations of a star can be caused by two bodies orbiting each other while orbiting that star, as is the case in hierarchical triples. In Morais & Correia (2012) the authors develop a secular theory for such systems and derive an expression for the precession rate of the orbit of the companion pair around the main component. Furthermore they suggest this effect as an alternative explanation to the planet hypothesis for the ν Octantis system (Ramm et al. 2009, 2016), which is known to be a close single-line spectroscopic binary where the two stellar components orbit each other with a period of 1050 d. In addition to the long-period RV signal caused by the stellar companion, Ramm et al. (2009) observed short-period RV variations with a period of 417 d, which they interpreted as a possible planet with a minimum mass of 2.5 M_jup. Stable orbital configurations have been published by Ramm et al. (2016), but previously the stability of the proposed planet seemed questionable. Morais & Correia (2012) argue that the short-period RV variations could also be caused by a very close pair orbiting the main component, meaning that the stellar companion ν Octantis B would itself be made up of two bodies. They also show that the precession rate of the long-period orbit as measured from the RVs is in agreement with their hypothesis.

Here we aim to test the hierarchical triple configuration for the ε Cyg system. Unfortunately, our RV data do not even cover one full phase of the long-period orbit, which makes any measurement of the precession rate very unreliable. Therefore we chose a different approach: Using a modified version of the Wisdom-Holman N-body integrator in the SWIFT package we modeled hierarchical systems and tried to reproduce the RV variations observed in our data. As the integrator was originally created to model the Solar System, it only allows to set up orbits of bodies around the central star of the modeled system, not of two outer bodies orbiting each other whilst orbiting together around the central star. By changing from barycentric to Jacobi coordinates however, we are able to choose one of the two components of the companion pair (denoted as m₁) as central object. In this transformed setup, the other component of the companion pair (m₂) as well as the main component (m₀) are then orbiting m₁.

We used a brute-force technique to try to find systems that produce RVs similar to what we observe. The long-period orbit of the companion binary is effectively given by our solution in Sect. 4.1; as described above, these (Jacobian) orbital elements were now ascribed to the main component m₀. We then replaced the single companion by a pair of two bodies m₁ and m₂; m₁ was chosen as the central star in the system and we systematically varied the Keplerian elements of the orbit of m₂ about m₁ as well as the masses of the two bodies. After setting up the system, we let it run for 200 yr with a time step of 0.5 days and recorded the RVs of the main component.

In these simulated RVs the long-period orbit of the m₁ -m₂ pair around the main component m₀ was in general clearly visible and closely resembled our data (see top plot of Fig. 9). We then fitted a Keplerian model to that RV variation and subtracted it from the RVs. By computing a GLS periodogram of the residuals, we searched for short-period variations. While in most simulated systems the binarity of the companion was too small to produce any considerable signature, or one of the components was lost within a very short time, some simulated configurations did show peaks at the periods in question. These were all systems with a comparably large separation a₁ between the two bodies m₁ and m₂, and all of them showed signs of long-term instability. However, the short-period variations were always limited to the time around periastron passage of the m₁ -m₂ binary around the main component.

Figure 9 shows a characteristic example, for a system with a₁ = 0.3 AU, m₁ = 0.16 M_⊙ and m₂ = 0.12 M_⊙: after subtracting the long-period orbit, the residuals show variations with an amplitude of 30 m s⁻¹ near periastron passage (middle plot), with a quasi-periodic nature (lower plot). During that time, all three bodies are very close to each other (the periastron distance of the m₁ -m₂ system to m₀ is approximately 1.2 AU), so the varying quadrupole (and higher) moments of the companion pair induce a large motion on the main component. But this signal quickly decays to amplitudes below 1 m s⁻¹ further away from the time of periastron passage. In contrast, we clearly observe RV variations with amplitudes around 25 m s⁻¹ in the McMillan and Lick data of ε Cyg far from periastron passage.

We conclude that while hierarchical triples may for some systems be a valid alternative explanation for the observed RV signals, this is not a viable solution for the observed RVs of ε Cyg. In contrast to the case discussed in Morais & Correia (2012), the long-period orbit of ε Cyg is far too eccentric to be compatible with large-amplitude short-period RV variations over the whole orbit induced by a companion pair.

Fig. 9 Top: simulated RVs of the main component of a hierarchical triple, with the orbit of the companion pair around the primary chosen according to our solution for the close stellar companion to ε Cyg A. Middle: residuals of the simulated RVs after subtracting the high-amplitude, long-period RV signal modeled by a Keplerian orbit. Bottom: zoom in on the time around periastron passage of the eccentric long-period orbit.

6.2 Stellar spots

Many stars exhibit stellar activity, that is, stellar spots, flares etc., which can mimic a planet in RV observations by blocking or enhancing light from parts of the surface of the stellar photosphere. This effect can become especially problematic when searching for planets around M dwarfs, which are often quite active and therefore can show activity-induced RV signals of several m s⁻¹ (see e.g., Reiners et al. 2010; Barnes et al. 2011; Tal-Or et al. 2018). By checking time-series of activity indicators, such as photometry, bisector spans of the absorption lines, or the depths of the Hα and Ca II lines of these stars, and comparing them to the RV time series one can distinguish between activity-induced and planetary signals. In addition, the so-called chromaticity, that is, wavelength-dependence, of RV signatures caused by stellar activity can serve as a tracer (Reiners et al. 2010; Zechmeister et al. 2018).

Evolved stars in contrast are known to be less affected by large-scale activity. Our sample stars were chosen specifically because they were comparably quiet stars in the HIPPARCOS photometry measurements. The top panel of Fig. 10 displays the HIPPARCOS data of ε Cyg, which have been averaged whenever there was more than one data point within a day. The time series largely overlaps with the second half of the McMillan RVs, which are displayed in the plot underneath. Apart from some outliers, most of the HIPPARCOS measurements lie very close around the mean of 2.643 mag with an rms scatter of 3.7 mmag. This translates to a relative change in flux of about 0.34%.

Nevertheless, when taking a periodogram of the HIPPARCOS photometry, there is a peak around a period of 300 d (see bottom panel of Fig. 10), which matches the observed period of the RV signal (indicated by the red line) quite well. This might hint at a signal induced by spots in the photometry data and would mean that the rotation period of the star would be approximately 300 d or an integer multiple of that. To our knowledge, there are two modern measurements of the projected rotational speed of ε Cyg: Massarotti et al. (2008) determined it to be v_rot sini_⋆ = 1.2 km s⁻¹, while the analysis of Gray (2015) resulted in v_rot sini_⋆ = 1.0 ± 0.2 km s⁻¹ (where i_⋆ denotes the inclination of the stellar spin axis). Assuming i_⋆ = 90° these values correspond to rotation periods of 461 and 554 ± 111 d, so in-between 1 and 2 times the period that is visible in our data. Taking the uncertainty in the measurement of Gray (2015) into account as well as the fact that for smaller values of the unknown inclination i_⋆ the rotation periods become shorter, it seems possible that the peak in the periodogram of the photometry actually traces the rotation of the star.

To further investigate this we searched the literature for additional photometric time series, but only found two more data sets: the BRITE nano-satellites BRITE-Toronto (BTr) and UniBRITE (UBr; Weiss et al. 2014) observed ε Cyg in the years 2014 and 2015, which coincides with the beginning of our SONG measurements (see Fig. 11, top). Unfortunately, the two data sets only cover 50 and 146 d, respectively, which is much shorter than one full phase of the short-period signal and therefore does not help to shed more light on the question whether that signal has a direct counterpart in photometry (see Fig. 11, top). The rms scatter of the data sets BTr and UBr are 2.0 and 1.1 mmag, respectively, which is even lower than in the HIPPARCOS data. Interestingly however the periodograms of the measurements reveal some variability on even shorter periods, with a strong peak around 40 d in both data sets (see Fig. 11, bottom). The HIPPARCOS periodogram seems to have a rough analog with a peak around 60 d, although in that case it might also be caused by the sampling frequency (see Fig. 10, bottom).

Generally, it would be surprising if the observed short-period RV variations of ε Cyg were caused by spots, as it is an evolved star with a radius of about 11 R_⊙. In order to cause RV variations on the order of the observed semi-amplitude of 20− 30 m s⁻¹, stellar spots would have to be extremely large, which in turn should become clearly visible in photometry as well. To study the effect qualitatively, we used the online tool of the star spot simulator SOAP 2.0 (Dumusque et al. 2014) to find a combination of spots that could possibly mimic the observed RV signal. Setting the rotation period of the simulated ε Cyg A to the observed modulation of 291 d, we found that a single spot would need to have a radius of 0.2 R_⋆ to induce a RV semi-amplitude of about 25 m s⁻¹, and the resultant flux change would be around 2.7%, much more than the scatter of 0.34% in the HIPPARCOS photometry. Using more than one single spot to create the same RV amplitude generally allows to decrease the individual spot sizes a little, but the photometric variability remains the same or even increases. Four spots with radii of 0.12 R_⋆ each at very close longitudes and latitudes, for example, keep the RV semi-amplitude roughly constant but lead to a flux variability of 3.5%. Similarly, placing spots on opposite sides of the star, or assuming an inclination i < 90°, only increases the predicted photometric variability.

We conclude that even though there are hints of a signature of the rotation of the star in the HIPPARCOS periodogram, the modeling suggests that spots are an unlikely cause of the observed RV variations.

Fig. 10 Top: HIPPARCOS photometry for ε Cyg, taken around the same time as the last RVs of the McMillan data set. The displayed photometry data are median values for all cases of multiple measurements within one day; the RVs are the residuals of the McMillan measurements after subtracting the long-period signal induced by the binary companion (see Sect. 4.1), and the dashed line denotes the best-fit model. Bottom: GLS periodogram of the HIPPARCOS photometry data and of the window function of the observations. The best-fit Keplerian orbital period of the signal of 291 d is marked with a red line.

Fig. 11 Top: BRITE-Toronto (BTr) and UniBRITE (UBr) photometry for ε Cyg, taken contemporaneously with the first RVs of the SONG data set. The displayed photometry data are median values for all cases of multiple measurements within one day, and both data sets have been normalized by their respective median values; the RVs are the residuals of the SONG measurements after subtracting the long-period signal induced by the binary companion (see Sect. 4.1), and the dashed line denotes the best-fit model. Bottom: GLS periodogram of the Brite photometry data sets. The best-fit Keplerian orbital period of the signal of 291 d is marked with a red line.

6.3 Oscillatory convective modes

In many evolved stars with very high luminosities, photometric oscillations are known with typical periods between about 400 to 1500 d, but the causes of these oscillations are not fully understood (Wood et al. 1999, 2004; Hinkle et al. 2002). Those stars are classified as long secondary period (LSP) variables because their primary oscillation periods are much shorter, with the period ratios between the secondary and primary oscillations ranging from about 5 to 13 (Wood et al. 1999). While the oscillations in most LSP stars have only been detected photometrically, Hinkle et al. (2002) and Wood et al. (2004) also obtained RVs for some of these stars and found the long variations present in most cases, with amplitudes of a few km s⁻¹ and periods consistent with the photometric variations. Saio et al. (2015) examine the observed periodic variations of the LSP stars and suggest so-called oscillatory convective modes, that is, nonadiabatic g⁻ modes in the deep convective envelopes of the stars, as a possible explanation.

However, all of the known LSP stars have luminosities of L_⋆ > 300 L_⊙ (Saio et al. 2015), whereas ε Cyg has a much lower luminosity of L_⋆ = 57 L_⊙. It is therefore not clear whether the mechanism described by Saio et al. (2015) might be responsible for the observed RV variations inthis case, but it is certainly possible that LSP oscillations also exist in stars with lower luminosities, and that the distribution of known LSP stars is constrained by observational biases: most of them have been found and characterized in the OGLE survey (see e.g., Udalski et al. 1997; Soszyński et al. 2009), which is aimed at discovering microlensing events and focuses on very distant stars in fields toward the Galactic Center and the Large Magellanic Cloud; stars similar to ε Cyg therefore are probably just not bright enough to be observed with sufficient S/N in order to detect the photometric variations.

Two other examples of evolved stars whose RV variations might well be caused by this mechanism are γ Draconis (Hatzes et al. 2018) and Aldebaran (Reichert et al. 2019); both show signals with periods around 700 and 600 d, respectively, which could easily be confused with Keplerian signals if one just examined parts of the data. Similar to ε Cyg, in both cases only the very long time span of the collected RV measurements allowed the authors to discover the varying nature of the signal and therefore reject the planet hypothesis. As the luminosities of γ Draconis and Aldebaran are L_γDra = 510 ± 51 L_⊙ (Hatzes et al. 2018) and (Reichert et al. 2019), they fall comfortably inside the distribution of known LSP stars and therefore oscillatory convective modes seem to be a plausible explanation for the observed RV signals. With only these two examples of similar (non-)detection histories, it could be that ε Cyg is another example of that category, even though its stellar parameters might suggest otherwise.

6.4 Potential stellar oscillations through the heartbeat phenomenon

A different mechanism that is known to produce oscillations in some eccentric binary systems is the so-called heartbeat phenomenon, where a stellar companion passes close by the primary component and excites tides within the latter. With the Kepler mission alone, more than 150 of these systems have been observed (Kirk et al. 2016). They are characterized by a sudden increase in brightness of the primary star during periastron passage of the companion, due to the deformation of the star caused by the increased gravitational force, which is called the equilibrium tide. Furthermore, sometimes one or even several periodic luminosity-signals can be observed over the whole orbit, the so-called dynamical tides. These oscillate in characteristic stellar oscillation modes, which are excited through the regular close passage of the stellar companion (for a full theoretical description see e.g., Fuller 2017).

All of the known binary systems with heartbeat variations are however much more compact than ε Cyg, with orbital periods typically shorter than 1 yr (Kirk et al. 2016). Still, in an attempt to investigate this possibility further, we revisit two publications on heartbeat systems that have been monitored in RVs to better constrain their orbits: Shporer et al. (2016) analyze a sample of 19 main-sequence stars in heartbeat binaries, with orbital periods between 8− 90 d, while Beck et al. (2014) examine 18 RGB heartbeat stars whose stellar companions orbit them with periods between 20− 438 d. In analogy to Fig. 18 from Beck et al. (2014), Fig. 12 displays important orbital and stellar parameters of all systems from both samples. The plot shows a clear positive correlation between primary stellar radius and binary orbital period for the RGBstars, while the main-sequence stars seem to be scattered more or less uniformly at the small-radius, short-period end of the distribution. Generally, for a given primary stellar radius heartbeats may occur in systems with longer orbital periods P if the eccentricity e is also larger, as the periastron distance scales as q ∝ P^2∕3(1 − e). ε Cyg however falls far away from all other systems at a much longer orbital period; even though its eccentricity is also larger than any other inthe combined sample and its stellar radius rather large, it is questionable from this figure whether ε Cyg could be a heartbeat system.

Shporer et al. (2016) further investigate the ratio of the tidal forcing due to the stellar companion to the surface gravity of the primary, using the relation (4)

where R_1,2 and M_1,2 are the radius and mass of the primary and secondary component, respectively, G is the gravitational constant and q is the periastron distance of the two stars. To compute the masses of the secondary components, the authors used the mass functions from their RV analysis and assumed sin³i = 0.6495; this value corresponds to the median of the sin³i distribution for the unknown inclinations and it was applied instead of the mean due to the very asymmetric shape of the distribution (see Shporer et al. 2016, for more information). We redid these calculations for both samples of Shporer et al. (2016) and Beck et al. (2014), using the stellar and orbital parameters from their works, as well as for the ε Cyg system (scaling its companion mass M₂ by the same inclination as used for the other systems). In Fig. 13 we plot the primary stellar radius (top) and tidal forcing ratio as calculated in Eq. (4) (bottom) over the periastron distance in units of the primary stellar radius (q∕R_*) for each system; the stellar radii are also indicated by the sizes of the symbols in the plots.

Plotted against q∕R_*, the stars of the two samples now fall much closer to each other, with most of them having periastron distances between 4 and 12 R_*. The RGB heartbeat stars populate thearea at small values of q∕R_*, whereas themain-sequence stars extend the distribution to larger values of q∕R_*. One main-sequence heartbeat system even has a q ≈ 24 R_* (this one being KIC 10334122), which is comparable to the periastron distance of ε Cyg, 22 R_*. Due to its larger radius however, ε Cyg falls above the distribution of the samples of RGB and main-sequence stars.

Finally, as is to be expected from Eq. 4, the ratio of the tidal forcing to the surface gravity of the primary shows a clear correlation with q∕R_*. Most of the main-sequence and RGB stars fall between values of 4 × 10⁻³ (when q∕R_* is small) and 4 × 10⁻⁴ (when q∕R_* is larger), and the rms scatter of the two samples combined is 1.55 × 10⁻³; the main-sequence system at q ≈ 24 R_* has a tidal forcing ratio around 7 × 10⁻⁵. Calculating the same quantity for ε Cyg gives a value of 2.3 × 10⁻⁵, which is not far from the main-sequence outlier. Given that ε Cyg lies so close to a known heartbeat system in these quantities makes it appear plausible that its primary star ε Cyg A might also undergo tidally induced stellar oscillations. Additionally, it is important to keep in mind that the tidal forcing ratio scales with the companion mass M₂, and that ε Cyg might therefore be even closer to the other stars in the case of a low inclination of the system. For example, for the lower limit on the inclination of i = 14°, the tidal forcing would be 1 × 10⁻⁴, which is greater than that of the main-sequence outlier and nearly within the rms scatter from the mean of the RGB and main-sequence distribution.

We also emphasize the fact that the variability of ε Cyg has been observed in RVs, whereas the systems from Beck et al. (2014) and Shporer et al. (2016) have all been discovered through their photometric variations in Kepler observations. Even though both these studies also monitored their systems in RVs to constrain their orbits, they do not mention any detection of RV variations corresponding to the heartbeat oscillations. This could have several reasons: they might not have picked up the signals in RVs due to sparse sampling or large errors of their measurements. Also some of their heartbeat systems barely seem to show any noticeable dynamical tides and mostly just equilibrium tides during periastron passage of the stellar companion, so there might not even be an RV signal to detect during most parts of the orbit. For ε Cyg in contrast we have a large number of RV measurements over most parts of the binary orbit, which clearly record a short-period signal that changes in period and semi-amplitude especially around periastron passage of the stellar companion. Given the differing detection method as compared to the other heartbeat stars, it is possible that ε Cyg experiences a similar phenomenon even though it does not follow all correlations exactly.

Fig. 12 Primary stellar radii over binary orbital periods for a number of heartbeat systems, with the orbital eccentricity color-coded. Triangles mark the main-sequence stars from Shporer et al. (2016), circles the RGB stars from Beck et al. (2014). For comparison ε Cyg has been added, denoted by a cross. The plot has been created in analogy to Beck et al. (2014).

Fig. 13 Primary stellar radii (top) and tidal forcing ratio (bottom) over periastron distances in units of stellar radii (q∕R*) for heartbeat systems, consisting of the main-sequence star sample from Shporer et al. (2016) (green) and the RGB stars from Beck et al. (2014) (blue). ε Cyg has been added for comparison (orange). Symbol sizes are proportional to the stellar radii. The tidal forcing ratio is calculated as in Eq. (4), and the large errors of the main-sequence stars stem from the high uncertainties on the radii R1 and the masses M1 of the primary components. The plot has been created in analogy to Shporer et al. (2016).

7 Summary and conclusions

Short-period RV variability of the spectroscopic binary ε Cyg had already been noticed in past publications (McMillan et al. 1992; Gray 2015) and had been attributed to intrinsic stellar variations rather than an extrinsic source. We observe a similar RV signal in our measurements of the star at Lick and with SONG and further investigated the possibility of a planetary origin.

In a first step we determined the long-period binary orbit with very high precision by fitting a Keplerian model to the three RV data sets of McMillan et al. (1992), Lick and SONG. The results indicate the possibility of directly imaging the secondary component in the years 2021/2022. A direct observation would constrain the orbital inclination and the mass of the spectroscopic companion, possibly revealing whether it is a white dwarf or a main-sequence star, which could give insight into the past evolution of the binary system.

Next we computed GLS periodograms of the residuals after subtracting the Keplerian model from the RV data; they show clear signs of additional periodic signals at periods just below 300 d for all three data sets. We fitted the combined RV data with a double-Keplerian model to explore the possibility of an S-type planetary companion as cause for the short-period RV variations. The best fit leaves a considerable number of systematic outliers that look like a phase shift of the signal at some point in time or a change of orbital period. By modeling each data set, and smaller sections of the data sets individually, we revealed that the planetary period gradually increases over time, whereas the RV semi-amplitude mostly stays constant until it increases by a factor of 2 at periastron passage of the stellar companion. This indicates that a single circumprimary planetary companion is an unlikely explanation for the observed RV signal, but the residuals of the double-Keplerian model could point toward a scenario with a second low-mass, co-orbital companion in the system. Using a fully dynamical model for the planetary companion, treating it as a test particle in the stellar binary system, does not help to solve the issue, as it cannot reproduce the systematic changes of the Keplerian elements, neither for the prograde nor for the retrograde case. Additionally, our stability analysis shows highly unstable behavior in large regions around the best-fit orbital configurations, thus strengthening the arguments against a planet.

Following Morais & Correia (2008, 2011, 2012), we explored whether the RV variations could be due to an additional component not orbiting the primary, but the secondary star, making the system a hierarchical triple. Our models clearly show that in our case such a setup cannot reproduce the data. We therefore investigated possible intrinsic stellar origins of the signal, such as stellar spots. Data from HIPPARCOS suggest that ε Cyg A is a photometrically relatively quiet star, even though a periodogram of the measurements reveals a possible weak signal around 300 d, which is close to the period observed in RVs and could correspond to the rotation period of the star. If the RV variations were due to stellar spots, we would however expect the photometric variation to be much larger than the observed scatter.

Two other examples of ambiguous RV variations of giant stars, γ Draconis (Hatzes et al. 2018) and Aldebaran (Reichert et al. 2019), might belong to the class of LSP variables described by, for example, Hinkle et al. (2002), Wood et al. (1999, 2004), and Saio et al. (2015). As ε Cyg A is much smaller and less luminous than the known LSP stars, this mechanism can offer a satisfactory explanation in this case only if the distribution of LSP stars known today is incomplete due to observational biases.

A more likely, even though also not completely convincing solution is the possibility of ε Cyg being a heartbeat system, with the close stellar companion inducing tides in the primary component during each periastron passage that excite long-standing oscillations with the observed period. However, there are some distinct differences between the known heartbeat systems and ε Cyg, with the periastron distance of the latter being much larger. Also the question remains whether there are oscillation modes at periods around 300 d in HB stars that could easily be excited through this mechanism.

Generally, we showed that none of the phenomena discussed above offer a fully satisfactory explanation for the RV variability by themselves at the moment. Nevertheless, it is possible that a combination of some of the mechanisms could explain the observed RVs while also eliminating some of the problems encountered in the analysis. For example, the oscillatory convective modes might be present in ε Cyg A and, although not strong enough to cause recognizable RV variations by themselves in this type of star, could be excited regularly by the tidal interaction with the stellar companion. On the one hand, this would explain why ε Cyg A is so far away from the distribution of known LSP stars and still shows oscillations; on the other hand, the fact that these oscillation modes are present and probably easily excited offers a solution to the problem of the comparably large separation of the stellar binary when compared to other heartbeat systems. Similarly, Saio et al. (2018) show the presence of tidally excited Rossby waves (so-called r modes) in some heartbeat systems, which typically appear at periods slightly longer than the rotation period. If the observed RV signal of ε Cyg was caused by this mechanism, the star’s rotation period therefore would have to be shorter than ~ 280 d, and the inclination of its spin axis i_⋆ smaller than ~ 30° (see Sect. 6.2, given that the v_rot measurement by Gray 2015 is correct). As the work by Saio et al. (2018) is focused on upper MS stars, it obviously remains questionable whether r modes can also be present in HB stars like ε Cyg A.

Another possibility is that there actually is a planet present around ε Cyg A, and its orbital period corresponds to the rotational period of the star, meaning the system is tidally locked. The observed RV variations could then be a sum of the Keplerian signal of the star and the rotational modulation, therefore solving the mismatch between the low photometric variations and the large RV amplitude for a rotational signal alone. Additionally, a tidally coupled system might improve the stability of a planetary companion in an orbit around the best-fit parameters derived in Sect. 4.3. In order to investigate this, it is necessary to pair an N-body code with a state-of-the-art algorithm for the computation of tidal forcing, which is beyond the scope of this work.

For the sake of completeness, we also mention the recent work by Maciejewski et al. (2020), who show that the slightly eccentric RV signature of WASP-12, which has been explained by a hot Jupiter on an unusual eccentric orbit, could also be explained by a circular orbit of the planet and the signature of tides in the host star. Arras et al. (2012) delivered the theoretical foundation, proving that tides in stars, induced by massive planets on tight orbits, can manifest themselves in RV signals with periods of half the orbital period and an amplitude of a few m s⁻¹; due to a phase shift between the planetary and tidal RV signals the combined signature can then mimic an orbit with nonzero eccentricity and longitude of periastron equal to 270°. This raises the question whether a similar effect could be at work in ε Cyg, where tides raised in the primary ε Cyg A during the close periastron passage of the spectroscopic companion might contribute to the recorded RV signature; with the eccentricity of the long-period orbit being constrained mostly by the data around periastron passage, the true orbital eccentricity might then be lower than our models suggest. This in turn could solve the stability issues of the putativeplanet discussed in Sect. 5 if the eccentricity actually would drop below a value of 0.9. It is questionablethough whether the RV contribution of tides would be large enough, and a more thorough analysis is needed to assess this possibility.

Generally however, we must conclude that the cause of the short-period RV variations of ε Cyg remains a mystery for now. This analysis confirms what has been shown in other recent publications about evolved stars: With more years of continued RV measurements we are now detecting more stars with signals that mimic planets but are more likely caused by other phenomena. Only thanks to the long time spans of the observations and through careful analysis are we able to identify these false positives. This should be a warning, as a considerable number of the published planets around giant stars, especially the early discoveries, had been declared on the basis of very sparse data sets, and some of them might have to be revoked once more measurements are available.

Acknowledgements

Part of this work was supported by the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg, IMPRS-HD, Germany. We would like to thank the staff at the Lick Observatory for their support during the years of this project. We are also thankful to the CAT observers who assisted with this project, including Saskia Hekker, Simon Albrecht, David Bauer, Christoph Bergmann, Stanley Browne, Kelsey Clubb, Dennis Kügler, Christian Schwab, Julian Stürmer, Kirsten Vincke, and Dominika Wylezalek. This work also includes observations made with the Hertzsprung SONG telescope operated at the Spanish Observatorio del Teide on the island of Tenerife by the Aarhus and Copenhagen Universities and by the Instituto de Astrofísica de Canarias. Part of the work has been performed with the RV analysis tool Exo-Striker (Trifonov 2019), and many of the simulations have been carried out on the machine of the Novel Astronomical Instrumentation based on photonic light Reformatting (NAIR) group at the Landessternwarte Heidelberg. This research has been funded by the Spanish State Research Agency (AEI) Projects No.ESP2017-87676-C5-1-R and No. MDM-2017-0737 Unidad de Excelencia “María de Maeztu”- Centro de Astrobiología (INTA-CSIC). Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant agreement no.: DNRF106). S.A. acknowledges the support from the Danish Council for Independent Research through the DFF Sapere Aude Starting Grant No. 4181-00487B. K.H.W. and M.H.L. are supported in part by Hong Kong RGC grant HKU 17305015. V.A. was supported by a research grant (00028173) from VILLUM FONDEN. The research made use of the SIMBAD database operated at CDS, Strasbourg, France, as well as of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. We are grateful to Bogdan Ghiorghiu’s contribution in the analysis of the RV measurements. We thank A. Leleu and A. Correia for private discussions on the co-orbital demodulation technique. We also thank H. Kjeldsen for his feedback concerning the jitter analysis, and R. Kuschnig for making the BRITE data available to us. Finally, we thank the anonymous referee for his very valuable review, especially his questions concerning the jitter analysis.

References

All Tables

All Figures

	Fig. 1 Uncorrected SONG RV measurements from the last week of our asteroseismic campaign on ε Cyg. Oscillations with periods shorter than a day are clearly visible.
In the text

Fig. 2 Top plots: RV measurements of ε Cyg by Kustner (1908) (left), and McMillan et al. (1992), Lick observatory and the SONG telescope (right), along with two orbital models of the binary companion: Fit 1 (black line) is the best-fit solution for the three modern data sets of McMillan, Lick and SONG, with an orbital period of 19 618.6 d; Fit 2 (gray line) was derived from fitting all four data sets and yields a period of 20 175 d. Middle plots: residuals of the RV measurements after removing the orbital solution of Fit 1. The last measurement of Kustner then is a clear outlier. Bottom plots: residuals of the RV measurements after removing the orbital solution of Fit 2. Now all the Kustner RVs agree with the model within their 1σ-uncertainties.

In the text

	Fig. 3 GLS periodograms of the residuals of the RV measurements after removal of the binary signal. The periodogram of the combined data sets (top panel) shows two peaks with similar heights at ~ 282.5 d and ~ 290.8 d. The gray dash-dotted, dashed and dotted lines represent false-alarm probability (FAP) levels of 5, 1 and 0.1%, respectively.
In the text

	Fig. 4 Top panel: RVs of the three modern data sets McMillan, Lick and SONG, phase-folded by the best-fit period of the putative planet P = 291.2 d, are shown in symbols along with the RV curve of the best fit as a black line. Bottom panel: residuals of the individual data points.
In the text

	Fig. 5 Orbital period P (top) and semi-amplitude K (bottom) of the presumed planet, derived from double-Keplerian fits to the individual data sets (solid lines) and to shorter sections of the data (dashed lines). The gray-shaded areas denote the time spans from the first until the last measurement within each section.
In the text

	Fig. 6 GLS periodograms of the residuals of the best dynamical models for the prograde (left) and retrograde configuration (right). The top-most plot shows the periodograms of the residuals of the combined data, the bottom three of the individual data sets. FAPs of 5, 1 and 0.1% are indicated by the gray dash-dotted, dashed and dotted lines, respectively.
In the text

	Fig. 7 Lomb-Scargle periodogram of the radial velocity residuals (blue) and the corresponding demodulated radial velocities (red). The red vertical dotted lines correspond to the range of allowed libration periods for this star-planet system based on a minimum (mt = 0) and a maximum (mt = mp) mass for theco-orbital body.
In the text

Fig. 8 Survival times for test particles placed into the ε Cyg stellar binary, in prograde (left) and retrograde coplanar configuration (right). The eccentricity of the orbit of the stellar companion and the initial semi-major axes of the test particles were varied in order to find the stability boundary. The position of the best-fit solution from Sect. 4.3 is plotted by the white diamond; its size is at least by a factor of 10 larger than the uncertainties on the Keplerian parameters.

In the text

Fig. 9 Top: simulated RVs of the main component of a hierarchical triple, with the orbit of the companion pair around the primary chosen according to our solution for the close stellar companion to ε Cyg A. Middle: residuals of the simulated RVs after subtracting the high-amplitude, long-period RV signal modeled by a Keplerian orbit. Bottom: zoom in on the time around periastron passage of the eccentric long-period orbit.

In the text

Fig. 10 Top: HIPPARCOS photometry for ε Cyg, taken around the same time as the last RVs of the McMillan data set. The displayed photometry data are median values for all cases of multiple measurements within one day; the RVs are the residuals of the McMillan measurements after subtracting the long-period signal induced by the binary companion (see Sect. 4.1), and the dashed line denotes the best-fit model. Bottom: GLS periodogram of the HIPPARCOS photometry data and of the window function of the observations. The best-fit Keplerian orbital period of the signal of 291 d is marked with a red line.

In the text

Fig. 11 Top: BRITE-Toronto (BTr) and UniBRITE (UBr) photometry for ε Cyg, taken contemporaneously with the first RVs of the SONG data set. The displayed photometry data are median values for all cases of multiple measurements within one day, and both data sets have been normalized by their respective median values; the RVs are the residuals of the SONG measurements after subtracting the long-period signal induced by the binary companion (see Sect. 4.1), and the dashed line denotes the best-fit model. Bottom: GLS periodogram of the Brite photometry data sets. The best-fit Keplerian orbital period of the signal of 291 d is marked with a red line.

In the text

	Fig. 12 Primary stellar radii over binary orbital periods for a number of heartbeat systems, with the orbital eccentricity color-coded. Triangles mark the main-sequence stars from Shporer et al. (2016), circles the RGB stars from Beck et al. (2014). For comparison ε Cyg has been added, denoted by a cross. The plot has been created in analogy to Beck et al. (2014).
In the text

Fig. 13 Primary stellar radii (top) and tidal forcing ratio (bottom) over periastron distances in units of stellar radii (q∕R*) for heartbeat systems, consisting of the main-sequence star sample from Shporer et al. (2016) (green) and the RGB stars from Beck et al. (2014) (blue). ε Cyg has been added for comparison (orange). Symbol sizes are proportional to the stellar radii. The tidal forcing ratio is calculated as in Eq. (4), and the large errors of the main-sequence stars stem from the high uncertainties on the radii R1 and the masses M1 of the primary components. The plot has been created in analogy to Shporer et al. (2016).

In the text