# Dynamics and habitability in circumstellar planetary systems of known binary stars

Ákos Bazsó, Elke Pilat-Lohinger, Siegfried Eggl, Barbara Funk, David

Bancelin, Gioia Rau

► To cite this version:

Ákos Bazsó, Elke Pilat-Lohinger, Siegfried Eggl, Barbara Funk, David Bancelin, et al. Dynamics and habitability in circumstellar planetary systems of known binary stars. Monthly Notices of the Royal Astronomical Society, 2017, 466 (2), pp.1555-1566. 10.1093/mnras/stw3095. hal-02262500

HAL Id: hal-02262500

https://hal.science/hal-02262500v1

Submitted on 26 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Dynamics and habitability in circumstellar planetary systems of known binary stars

Ákos Bazsó,<fnref n="1" />* Elke Pilat-Lohinger,<fnref n="1" />* Siegfried Eggl,<fnref n="2" /> Barbara Funk,<fnref n="1" /> David Bancelin<fnref n="1" /> and Gioia Rau<fnref n="1" />

Accepted 2016 November 24. Received 2016 November 23; in original form 2016 September 28

## ABSTRACT

We present a survey on exoplanetary systems of binary stars with stellar separations less than 100 au. For a sample of 11 binaries that harbour detected circumstellar giant planets, we investigate the frequency of systems with secular resonances (SR) affecting the habitable zone (HZ). SR are connected to dynamically unstable or chaotic regions by enforcing highly eccentric motion. We apply a semi-analytical method to determine the locations of linear SR, which is based on finding the apsidal precession frequencies of the massive bodies. For configurations where the giant planet is located exterior to the HZ, we find that there is always an SR interior to its orbit; the exact location of the SR strongly depends on the system's architecture. In systems with the giant planet interior to the HZ, no SR can occur in the Newtonian framework. Taking into account the general relativistic precession of the perihelion, which increases the precession frequencies, planets with $a < 0.1$ au can cause SR in the HZ. We find two cases where the SR is located inside the HZ and some more where it is close to the HZ. Generally, giant planets interior to the HZ are more favourable than exterior planets to avoid SR in the HZ. Around the location of the SR weaker mean-motion resonances are excited and resonance overlap is possible. Existing analytical models are not as accurate as the semi-analytical method in locating the SR and deviate by ~0.1 au or more.

**Key words:** celestial mechanics – planets and satellites: dynamical evolution and stability – binaries: general.

## 1 INTRODUCTION

Extrasolar planets (exoplanets) orbiting other stars exist in a large variety of different configurations. Some exoplanets orbit their host star with periods of only a few days, e.g. the Hot Jupiters (see Beaugé & Nesvorný 2012; Naoz, Farr & Rasio 2012). Other exoplanets were observed on wide orbits at large distances from their host star (e.g. HR 8799; Marois et al. 2008); however, their formation mechanism is still poorly understood (Dodson-Robinson et al. 2009; Vorobyov 2013).

An important point is that exoplanets are not exclusive companions to single stars, but also binary and multiple stars host planetary objects (Desidera & Barbieri 2007; Roell et al. 2012). A considerable fraction of stars in the solar neighbourhood are members of binary and multiple star systems. Observational surveys derived a percentage of multiple stars of about 45 per cent (Duquennoy

& Mayor 1991; Raghavan et al. 2010; Duchêne & Kraus 2013; Tokovinin 2014). Of course, not all of those systems are hosting planets. Recent estimates on the multiplicity of exoplanet host stars yield the percentage of double stars to be 10–15 per cent, while multiple systems make up for about 2 per cent (Raghavan et al. 2006; Mugrauer & Neuhäuser 2009; Roell et al. 2012).

Although exoplanets in binary and multiple star systems have been found, it is not clear whether or not these environments are more hostile for the presence of planets than single stars (Boss 2006; Bromley & Kenyon 2015; Jang-Condell 2015). The planet formation process in binaries is basically proceeding in the same way as for single stars (Thebault & Haghighipour 2015); however, the presence of a companion star can truncate the proto-planetary disc (Kley & Nelson 2010) and influence planet formation (Müller & Kley 2012), possibly limiting the occurrence rate of planets (see Winn & Fabrycky 2015). Wang et al. (2014) found that planets in multiple star systems with stellar separations between 10 and 1000 astronomical units (au) occur about a factor of 2 less frequently compared to single stars. On the other hand, exoplanets may still form in 'tight' binaries with separations of ≲30 au. Jang-Condell

<footnote n="1">¹Institute for Astrophysics, University of Vienna, Türkenschanzstraße 17, A-1180 Vienna, Austria</footnote>

<footnote n="2">²IMCCE, Observatoire de Paris, 77 Avenue Denfert-Rochereau, F-75014 Paris, France</footnote>

<footnote>* E-mail: bazso@astro.univie.ac.at (AB); elke.pilat-lohinger@univie.ac.at (EP-L)</footnote>
(2015) calculated the likelihood that giant planets can form from the truncated disc in such tight binaries. The results indicate that the discs still contain enough material to form giant planets by core accretion, while terrestrial planets can form in even more depleted discs.

Detected exoplanets in binary star systems can be divided into three dynamical classes, of which S-type (Rabl & Dvorak 1988) and P-type (Dvorak, Froeschle & Froeschle 1989) are the most important cases. S-type planets are circumstellar planets; they orbit one component of the binary star system, just as a planet would orbit a single star. P-type planets are circumbinary planets that orbit both stars in some distance. In both cases, long-term stable planetary orbits may occupy only certain regions. Rabl & Dvorak (1988) established an empirical formula for the size of the circumstellar region where stable planetary motion can occur in an eccentric binary system for equal mass stars. Holman & Wiegert (1999) generalized those results to both S-type and P-type binaries with variable mass-ratios and gave formulas for the edge of the stable region as function of binary separation, eccentricity and mass-ratio. Pilat-Lohinger & Dvorak (2002) also included the influence of a non-zero initial planetary eccentricity on the orbital stability region, which can shrink considerably for large mass-ratios and eccentric planets relative to the circular case.

Albeit the orbital stability of exoplanets is an intriguing topic in itself, it is not the only prerequisite for habitability. Substantial changes in insolation that are due to the variations in planetary orbits have to be taken into account as well (Williams & Pollard 2002). This tight coupling between orbital dynamics and habitability is demonstrated by Forgan (2016), who uses a coupled model for spin-orbit-radiative perturbations to investigate Milankovitch cycles of P- and S-type planets in binary stars. The stellar companion induces temperature oscillations on the planet on short and long time-scales via the changes in the obliquity and precession.

Since orbital dynamics is also important to assess habitability, we need to consider the effect of both mean-motion resonances (MMR, involving orbital frequencies) and secular resonances (SR, related to apsidal frequencies) on the planet's orbit. SR act on time-scales comparatively long to a planet's orbital period. In connection with dynamically unstable or chaotic regions in planetary systems, we are primarily focusing on the locations of SR (see Murray & Dermott 1999; Levison & Agnor 2003; Beaugé, Ferraz-Mello & Michtchenko 2012), because SR are usually connected to highly eccentric motion and so they can have a strong influence on a planet's habitability.

In this work, we investigate a sample of binary star systems with a known giant planet. Our goal is to check how frequently do SR affect the orbits of bodies in the HZ and what are the conditions for this to happen. We employ a semi-analytical method (Pilat-Lohinger, Bazsó & Funk 2016) that allows us to quickly identify SR and exclude unstable regions from a stability analysis. A possible application of the method is for newly discovered exoplanetary systems of binary stars, where it would help to exclude those systems that cannot host additional planets in the HZ. This, in turn, makes the method especially interesting as a guide to select such binary star systems for in-depth observations that would allow planets in the HZ from a dynamical point of view.

Section 2 presents the binary star systems selected for our investigation and explains how we determined the HZ for those systems. In Section 3, we give a short summary of the semi-analytical method introduced by Pilat-Lohinger et al. (2016) that allows us to calculate accurately the location of SR in binary star systems. We present the results of the survey in section 4 and discuss the semi-analytical

method and compare its performance relative to some purely analytical models in section 5. Finally, in Section 6, we summarize the results and give the conclusions.

## 2 THE INVESTIGATED SYSTEMS

As of 2016 October, in total, 79 binary star systems with 113 planets had been discovered, of which 58 are of S-type and 21 of P-type.<fnref n="1" />

For this survey, we have selected S-type binary stars with (projected) separations $a_B \le 100$ au; such binaries have typical orbital periods below 1000 yr. This limit is based on the observation that planets in 'wide' binaries (with $a_B \ge 300$ au) have essentially the same mass and period distributions as planets orbiting single stars (see Desidera & Barbieri 2007). Another consideration was that the distance ratio planet/secondary star should not become too small, as we are focusing on secular effects and a small ratio would lead to only weak effects (see Section 3).

Given the above constraints on stellar separation, our sample consisted of 15 candidates. We discarded the Kepler-296 (Rowe et al. 2014; Barclay et al. 2015) and Kepler-444 (Campante et al. 2015) systems, because they are compact multiplanetary systems with five planets. We also rejected OGLE-2008-BLG-092L (Poleski et al. 2014) and OGLE-2013-BLG-0341L (Gould et al. 2014) due to insufficient data regarding the orbits of the companion star and the planet.

### 2.1 System properties

In the end, 11 systems remained to be investigated. Table 1 provides an overview of their basic parameters such as stellar masses ($M_A, M_B$) and the planet's minimum mass ($M_P = m \sin i$), as well as the semimajor axes ($a_B, a_P$) and eccentricities ($e_B, e_P$) for the secondary star and the planet, respectively. This table, in connection with Fig. 1, summarizes the main physical and dynamical aspects of the selected systems. The table includes references to the original detection paper and, in several cases, to additional sources that provide important updates to the parameters (e.g. on stellar or planet masses, distances, etc.). For several binaries, the eccentricity of the secondary star is unknown; in those cases, we assume a moderate value of $e_B = 0.2$ based on Duquennoy & Mayor (1991) and Raghavan et al. (2010).

Fig. 1 shows a graphical overview of the binary star systems. The circles for the primary (host) and secondary star indicate their masses (cf. Table 1); the box in the top left corner displays four typical cases for M-, K-, G- and F-type stars. A large dashed line extending from the left marks the host star's region of orbital stability for planetary motion, beyond which it can no longer host S-type planets (Holman & Wiegert 1999; Pilat-Lohinger & Dvorak 2002). As opposed to this, the secondary star might harbour S-type planets in a sub-region outside this zone (not shown). We indicate the eccentricities of the secondary star and the planet with error bars having the colour of the respective object.

Note that for $\alpha$ Centauri, the existence of the planet is disputed. Dumusque et al. (2012) first claimed the discovery of a close-in planet with a few Earth masses. However, Hatzes (2013) cautioned that the radial velocity (RV) signal may be induced by stellar activity rather than a planet. Recently, Rajpaul, Aigrain & Roberts (2016) presented evidence for a false positive detection. We also caution that in the Kepler-420 (KOI-1257) system, the giant planet and

<footnote n="1">¹ http://www.univie.ac.at/adg/schwarz/multiple.html</footnote>
<table><caption>**Table 1.** List of investigated binary star systems ordered by decreasing stellar separation $a_B$. The planet hosting star is indicated by the capital letter A or B following the system's designation. Stellar masses $M_A$ and $M_B$ are in units of the solar mass, while planet masses $M_P$ are given in Jupiter masses.</caption><thead><tr><th>Name</th><th><i>M</i><sub>A</sub><br>(M<sub>☉</sub>)</th><th><i>M</i><sub>B</sub><br>(M<sub>☉</sub>)</th><th><i>M</i><sub>P</sub><br>(<i>M</i><sub>J</sub>)</th><th><i>a</i><sub>B</sub><br>(au)</th><th><i>a</i><sub>P</sub><br>(au)</th><th><i>e</i><sub>B</sub></th><th><i>e</i><sub>P</sub></th><th>Reference</th></tr></thead><tbody><tr><td>94 Cet Ab</td><td>1.34</td><td>0.90</td><td>1.68</td><td>100</td><td>1.42</td><td>0.26</td><td>0.3</td><td>Mayor et al. (2004); Roell et al. (2012)</td></tr><tr><td>HD 177830 Ab</td><td>1.37</td><td>0.23</td><td>1.49</td><td>97</td><td>1.22</td><td>0.20<sup>a</sup></td><td>0.001</td><td>Vogt et al. (2000); Roberts et al. (2015)</td></tr><tr><td>HD 177830 Ac</td><td>...</td><td>...</td><td>0.15</td><td>...</td><td>0.51</td><td>...</td><td>0.349</td><td></td></tr><tr><td>GJ 3021 Ab</td><td>0.90</td><td>0.13</td><td>3.37</td><td>68</td><td>0.49</td><td>0.20<sup>a</sup></td><td>0.51</td><td>Naef et al. (2001); Chauvin et al. (2007)</td></tr><tr><td>τ Boo Ab</td><td>1.30</td><td>0.40</td><td>4.13</td><td>45</td><td>0.05</td><td>0.20<sup>a</sup></td><td>0.02</td><td>Butler et al. (1997); Rodler, Lopez-Morales &amp; Ribas (2012)</td></tr><tr><td>HD 126614 Ab</td><td>1.15</td><td>0.32</td><td>0.38</td><td>36.2</td><td>2.35</td><td>0.50<sup>b</sup></td><td>0.41</td><td>Howard et al. (2010)</td></tr><tr><td>α Cen Bb</td><td>1.11</td><td>0.93</td><td>0.004</td><td>23.4</td><td>0.04</td><td>0.52</td><td>0.0</td><td>Dumusque et al. (2012); Endl et al. (2015)</td></tr><tr><td>HD 41004 Ab</td><td>0.70</td><td>0.40</td><td>2.54</td><td>23</td><td>1.64</td><td>0.20<sup>a</sup></td><td>0.2</td><td>Zucker et al. (2004)</td></tr><tr><td>HD 196885 Ab</td><td>1.33</td><td>0.45</td><td>2.98</td><td>21</td><td>2.60</td><td>0.42</td><td>0.48</td><td>Correia et al. (2008); Chauvin et al. (2011)</td></tr><tr><td>γ Cep Ab</td><td>1.40</td><td>0.41</td><td>1.85</td><td>20.2</td><td>2.05</td><td>0.41</td><td>0.05</td><td>Hatzes et al. (2003); Endl et al. (2011)</td></tr><tr><td>Gliese 86 Ab</td><td>0.83</td><td>0.49</td><td>4.01</td><td>19</td><td>0.11</td><td>0.40</td><td>0.05</td><td>Queloz et al. (2000); Fuhrmann et al. (2014)</td></tr><tr><td>Kepler-420 Ab</td><td>0.99</td><td>0.70</td><td>1.45</td><td>5.3</td><td>0.38</td><td>0.31</td><td>0.772</td><td>Santerne et al. (2014)</td></tr></tbody></table>

Notes. <sup>a</sup>assumed eccentricity; <sup>b</sup>upper limit for stable planetary motion.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Overview of the binary star systems. The sizes of the circles representing the primary and secondary star indicate their masses according to the scale in the top left corner; the planet is not to scale. Error bars show the minimum and maximum distances of the respective objects, which are related to the eccentricity. The large dashed line marks the region of orbital stability for S-type motion about the host star in the respective system.</figcaption></figure>

secondary star have such high eccentricities (Santerne et al. 2014) that the dynamically stable zone (DSZ) around the host star is most probably devoid of additional planets.

## 2.2 Habitable zone

The classical HZ was defined for a single Earth-like planet on an almost circular orbit around a single star (Kasting, Whitmire &

<table><caption>**Table 2.** Extension of the binary star habitable zone (AHZ) and the DSZ with inner and outer borders, outer limit of the stability region and location of SR in the binary star systems. Dots indicate non-existent values.</caption><thead><tr><th rowspan="2">System</th><th colspan="2">AHZ (au)</th><th colspan="2">DSZ (au)</th><th rowspan="2">Stability limit (au)</th><th rowspan="2">Resonance location (au)</th></tr><tr><th>Inner</th><th>Outer</th><th>Inner</th><th>Outer</th></tr></thead><tbody><tr><td>94 Cet</td><td>1.84</td><td>3.11</td><td>2.62</td><td>3.11</td><td>27.90</td><td>...</td></tr><tr><td>HD 177830</td><td>1.63</td><td>3.72</td><td>2.20</td><td>3.72</td><td>30.02</td><td>2.02, 3.78</td></tr><tr><td>GJ 3021</td><td>0.81</td><td>1.40</td><td>1.24</td><td>1.40</td><td>14.96</td><td>...</td></tr><tr><td>τ Boo</td><td>1.67</td><td>2.84</td><td>1.67</td><td>2.84</td><td>13.30</td><td>0.45</td></tr><tr><td>HD 126614</td><td>1.09</td><td>1.90</td><td>...</td><td>...</td><td>4.63</td><td>0.91</td></tr><tr><td>α Cen</td><td>0.72</td><td>1.26</td><td>0.72</td><td>1.26</td><td>2.45</td><td>1.16</td></tr><tr><td>HD 41004</td><td>0.79</td><td>1.40</td><td>0.79</td><td>0.94</td><td>5.20</td><td>0.37</td></tr><tr><td>HD 196885</td><td>2.06</td><td>3.53</td><td>...</td><td>...</td><td>3.84</td><td>1.11</td></tr><tr><td>γ Cep</td><td>2.51</td><td>5.90</td><td>...</td><td>...</td><td>3.86</td><td>0.78</td></tr><tr><td>Gliese 86</td><td>0.63</td><td>1.10</td><td>0.63</td><td>1.10</td><td>2.34</td><td>...</td></tr><tr><td>Kepler-420</td><td>0.95</td><td>1.67</td><td>...</td><td>...</td><td>1.02</td><td>...</td></tr></tbody></table>

Reynolds 1993; Kopparapu et al. 2013). Various attempts have been made to extend the concept of the HZ to binary stars. The major difference is to take into account the insolation on to the planet from both stars as a function of their luminosity, orbital parameters, as well as the changes in the planet's orbit.

Here we follow Eggl et al. (2012) to calculate the so-called averaged habitable zones (AHZ), where the effective insolation on to the planet is, on average, within the prescribed maximum and minimum habitable insolation limits: $S_{in} \ge \langle S_{eff} \rangle_{i} \ge S_{out}$. The AHZ represents an optimistic estimate of the HZ; besides the AHZ, the extended HZ and permanent HZ pose more restrictive conditions on the effective insolation. In the following, we always mean AHZ when using the term HZ. Alternative methods for calculating HZ in binary star systems were proposed by Kaltenegger & Haghighipour (2013) and Cuntz (2014). In contrast to Eggl et al. (2012), these methods neglect variations of the planet's orbit and would yield slightly different results.

Table 2 summarizes the extent of the binary star HZ for each system (in astronomical units) by specifying the inner and outer border of the AHZ. Besides the HZ, we define the DSZ to identify those parts of the HZ where terrestrial planets can remain on long-term stable orbits subject to the gravitational perturbations of the giant planet and secondary star. This table also indicates the border of the region of orbital stability for S-type motion. We have determined the boundary for regular S-type motion around the host star for
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Binary star HZ for the planet hosting star (hatched boxes). The blue area is the intersection of the HZ with the dynamically stable region (see the text for details). The two red-framed regions are beyond the stability limit. Black dots with horizontal error bars show the position and eccentricity of the planet, a red cross indicates the location of an (linear) SR.</figcaption></figure>

each individual system by using the Fast Lyapunov Indicator (FLI; Froeschlé, Lega & Gonczi 1997) as presented in Pilat-Lohinger & Dvorak (2002). These values are somewhat more restrictive than what the formula of Holman & Wiegert (1999) would predict, as they also account for the planet's eccentricity. The last column shows the location of the linear secular resonance with the detected planet (see Section 4); missing entries for four systems mean that there is no such resonance.

In Fig. 2, we plot the binary star HZ as hatched boxes. The intersection of the HZ with the region of stable orbital motion (DSZ) is plotted as the solid filled area, which we derived from numerical computations with test particles. The black dots show the location of the detected planets, where the error bars (extending from periastron to apoastron) visualize the planet's eccentricity. Interactions of a planet in the HZ with the giant planet can lead to dynamical instability in parts of the HZ by effects such as MMR and SR. In several cases, low-order MMRs with the giant planet occur at various locations inside the HZ; these MMR are not included in Fig. 2 but collected in Table 3. A cross symbol marks the location of linear SR (see Section 4). For the systems γ Cep and Kepler-420, parts of the HZ lie outside the S-type stable region (cf. Table 2); this is indicated by a red border around the corresponding part of the HZ. For HD 177830, the two planets cause SR at least at two locations, these are marked in Fig. 2 and given in Table 2.

Table 3 provides a non-exhaustive list of MMR with the giant planet(s) that happen to fall inside the HZ or nearby the SR. The

<table><caption>**Table 3.** Locations of MMR inside the HZ of the investigated systems or nearby the SR. We indicate MMR by the test planet to giant planet period ratio.</caption><thead><tr><th>System</th><th>MMR</th><th>Resonance location (au)</th></tr></thead><tbody><tr><td rowspan="3">94 Cet</td><td>2:3</td><td>1.86</td></tr><tr><td>1:2</td><td>2.25</td></tr><tr><td>2:5</td><td>2.61</td></tr><tr><td rowspan="4">HD 177830 b</td><td>2:3</td><td>1.60</td></tr><tr><td>1:2</td><td>1.94</td></tr><tr><td>8:17</td><td>2.02</td></tr><tr><td>2:11</td><td>3.81</td></tr><tr><td>HD 177830 c</td><td>1:8</td><td>2.05</td></tr><tr><td rowspan="3">GJ 3021</td><td>1:2</td><td>0.78</td></tr><tr><td>1:3</td><td>1.02</td></tr><tr><td>1:4</td><td>1.23</td></tr><tr><td rowspan="3">HD 126614</td><td>4:1</td><td>0.93</td></tr><tr><td>3:1</td><td>1.13</td></tr><tr><td>2:1</td><td>1.48</td></tr><tr><td rowspan="4">HD 41004</td><td>9:1</td><td>0.38</td></tr><tr><td>2:1</td><td>1.03</td></tr><tr><td>3:2</td><td>1.25</td></tr><tr><td>4:1</td><td>0.81</td></tr><tr><td rowspan="3">γ Cep</td><td>3:4</td><td>2.48</td></tr><tr><td>2:3</td><td>2.69</td></tr><tr><td>1:2</td><td>3.25</td></tr></tbody></table>

numbers in the second column indicate the orbital period ratios; when the latter number is larger, there is an exterior resonance ($a_p < a$), while it is smaller for interior resonances ($a < a_p$). We note that some of the MMR are located close to the border of the HZ or coincide with the transition from unstable to regular motion. For HD 177830, HD 126614, HD 41004 and γ Cephei, we included those MMR that are nearest to the location of the SR (cf. Fig. 2 and Table 2).

We note that HD 177830 A and γ Cephei A are post-main-sequence stars (subgiants) and the companion of Gliese 86 is a white dwarf (WD; Fuhrmann et al. 2014). In these cases, we cannot calculate the HZ as for ordinary main-sequence stars. For Gliese 86, a resort is to ignore the contribution of the WD companion, and to approximate the current HZ by the host star's HZ alone. This is justifiable by the low luminosities of cool WDs, e.g. $L/L_⊙ < 10^{-3}$ for $T_{eff} < 10^4$ K (see Liebert, Bergeron & Holberg 2005). For the other two evolved stars, we use the model of Ramirez & Kaltenegger (2016) to calculate the HZ of post-main-sequence stars. Their model provides the effective stellar flux $S_{eff}$ for the inner and outer boundary of the HZ for a grid of stellar masses between 0.5 and $1.9 M_⊙$, following the stellar evolution from the beginning of the red giant branch until the end of the asymptotic giant branch. We determine the effective stellar flux values for HD 177830 and γ Cep by systematic interpolation between the masses of stars in the model to fit the current masses, which yields the values in Table 2 when combining the results with the published luminosities (Jofré et al. 2015).

In addition, we tried to recover the main-sequence parameters for the two stars in order to obtain an estimate of their main-sequence HZ. A main-sequence mass for γ Cep and HD 177830 can be calculated through a Bayesian estimation of stellar parameters, i.e. calculating the most likely values of stellar mass, age and radius, knowing the effective temperature, the metallicity, the V-photometry and the distance (parallax); see da Silva et al. (2006) for details. The assumed Bayesian priors for the calculations are an initial mass function for single stars from Chabrier (2001) (lognormal), and a
constant star formation rate in the age interval from 0.1 to 12.0 Gyr. For the Bayesian analysis, we adopted the values from Soubiran et al. (2016) for $T_{\text{eff}}$, [Fe/H] and V-magnitude, and we used the parallax from van Leeuwen (2007). Given the wide range of possible metallicities and effective temperatures for $\gamma$ Cep, as well as the varying mass estimates (see Torres 2007), we can only set an upper limit on the main-sequence mass of $M = 1.62 \pm 0.10 M_{\odot}$. For HD 177830, we find $M = 1.41 \pm 0.03 M_{\odot}$, which is consistent with the mass given by Roberts et al. (2015) to within the error bars. We use the Kopparapu et al. (2013) model to calculate HZ boundaries for runaway/maximum greenhouse, where we assume both stars to be of spectral type F0 and having an effective temperature of 7200 K (the upper limit in the model). This results in HZ distances of 1.86–3.20 au for $\gamma$ Cep and 1.77–3.05 au for HD 177830. However, given the many assumptions about metallicities, masses and luminosities, these ranges should be used with caution.

## 3 METHODS

To investigate the secular dynamics of the binary star systems, we use a semi-analytical method that is presented and discussed in more detail in Pilat-Lohinger et al. (2016); here we give just a quick summary: we treat the binary star–giant planet–test planet four-body problem as two coupled three-body problems. First, for the pair giant planet–secondary star, we derive from a single N-body numerical integration the secular precession frequency of the giant planet’s pericentre with a characteristic frequency determination. This precession frequency is approximated as being fixed. Next, for the test planet–giant planet pair, we can apply the Laplace–Lagrange perturbation theory to find the locations of SR.

### 3.1 Laplace–Lagrange theory

The Laplace–Lagrange secular perturbation theory (LL) was developed to serve as a first approximation to the dynamical behaviour of the Solar system (see Murray & Dermott 1999). It is a first-order theory in the masses and includes terms in the eccentricity and inclination of degree $e^2$ and $(\sin i)^2$, while higher order terms are neglected. Nevertheless, this approximation already allows us to determine the location of the most important linear SR.

Let us consider an exoplanetary system consisting of a test planet (with mass $m \ll m_j$) moving under the influence of the host star ($m_0$), a giant planet ($m_1$) and the secondary star ($m_2$). All bodies are assumed to be point mass objects and located in a common plane.

To set up the equations of motion, usually the variables ($h = e \sin \varphi, k = e \cos \varphi$) are defined. The general solution to the equations of motion is then of the form

$$h(t) = e_{\text{free}} \sin(gt + \phi) - \sum_{j=1}^{2} \frac{v_j}{g - g_j} \sin(g_j t + \phi_j), \\ k(t) = e_{\text{free}} \cos(gt + \phi) - \sum_{j=1}^{2} \frac{v_j}{g - g_j} \cos(g_j t + \phi_j). \quad (1)$$

Here $a, e, n, \varphi$ are the semi-major axis, eccentricity, mean motion and longitude of pericentre of the test particle, respectively, and variables with subscript $j$ are associated with the perturbers ($j = 1$ the giant planet, $j = 2$ the secondary star). The constants $e_{\text{free}}$ and $\phi$ are given by the initial conditions for the particle, while $g_j, \phi_j$ and $v_j$ follow from the secular dynamics of the massive perturbers.

The quantity $g$ is the test particle’s proper secular frequency calculated from

$$g = \frac{n}{4} \sum_{j=1}^{2} \frac{m_j}{m_0} \alpha_j \bar{\alpha}_j b_{3/2}^{(1)}(\alpha_j). \quad (2)$$

This frequency $g$ depends on the semimajor axes of the particle and of the perturbers via the ratio $\alpha_j = a/a_j$ (in the case $a < a_j$) or $\alpha_j = a_j/a$ (when $a > a_j$), the Laplace coefficient $b_{3/2}^{(1)}(\alpha_j)$, as well as on the perturber’s mass $m_j$. Equations (1) and (2) can easily be extended to an arbitrary number $N$ of perturbing bodies. From equation (2), we can see that $g \propto O(\alpha^3)$ (using $b_n^{(k)} \propto O(\alpha^k)$), so the secular frequency quickly diminishes with decreasing $\alpha$, i.e. when the perturber is located farther away. For this reason, the main contribution to $g$ stems from the giant planet that is much closer than the secondary star, despite the difference in the masses. Extending this argument to well-separated binaries (with $a_B > 100$ au), it becomes clear that a distant secondary star would have only a negligible influence on the test planet’s secular frequency. In that case, the system can be approximated as a single star with a giant planet having a fixed forced secular frequency $g_1$.

Wherever $g \approx g_j$ in equation (1), the particle’s secular frequency equals one of the perturber’s frequencies, which then gives rise to a secular resonance. In these cases, the particle suffers from large secular perturbations due to the perturber $j$, where the particle’s forced eccentricity can reach values close to $e = 1$. Then, it can cross the orbit of the giant planet or be ejected from the system by a close approach to the host star.

### 3.2 Numerical part

In equation (1), we need to know the frequencies $g_j$ accurately in order to determine the location where $g - g_j \approx 0$, i.e. the location of a linear secular resonance. We chose to determine the fundamental frequencies numerically throughout this study. This is necessary as the eccentricities of the binary stars are large for all selected system ($e_B \ge 0.2$, see Table 1), such that a low-order analytical theory (like the LL theory) has only a limited accuracy.

We use the Lie-series method (Hanslmeier & Dvorak 1984; Eggl & Dvorak 2010; Bancelin, Hestroffer & Thuillot 2012) and the Mercury integrator’s Radau method (Everhart 1974; Chambers 1999) as complementary tools to integrate the equations of motions of the selected binary star–planet system. Next, we perform a frequency analysis on the results of our integrations to extract the fundamental frequencies; this step is accomplished by using the Fast-Fourier Transform (FFT) library FFTW<fnref n="2" /> of Frigo & Johnson (2005). As an independent check we use the tool SigSpec of Reegen (2007) that performs a Discrete Fourier Transform (DFT) and gives the probability that a peak is not due to noise in the time series. Whenever possible, we remove short-period contributions (associated with the orbital frequencies of the giant planet and the secondary star) by a low-pass filter. From the frequency analysis, we obtain the giant planet’s secular frequency $g_P$ (or $g_1$), which is effectively the time-averaged frequency $\langle g(t) \rangle_T$ over the integration time $T$. Finally, to determine the locations of SR for the test particles, we identify the intersection points of equation (2) with the values $g_j$.

Fig. 3 shows an example for the processing of the numerical integrations. The upper panel displays the time evolution of $k(t)$ during the first $10^7$ yr from a 30 Myr integration of the system HD 41004. In the bottom part, we see the DFT (black) and FFT (grey) spectrum of

<footnote n="2">² http://fftw.org/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of DFT and FFT results. Top: time evolution of variable $k(t)$ (from equation 1) for the system HD 41004. Bottom: Fourier amplitudes for the signal from above. The peaks identified as $g_1, g_2$ denote the secular frequencies of the planet and the star, respectively.</figcaption></figure>

The full signal, where frequencies above 330 arcsec yr$^{-1}$ have been filtered out. The two curves are shifted vertically relative to each other, which means that their amplitudes were scaled differently and are not directly comparable; however, the horizontal frequency axis is the same for both. Both techniques agree excellently in the identification of the peaks, where $g_1$ and $g_2$ are the system's main secular frequencies attributable to the giant planet and the secondary star. The third largest peak $g_1 + \gamma_2$ is a linear combination of the planet's secular frequency with a short-period term from the star's signal, but we have not attempted to identify its origin. In the top panel, the effect of $g_2$ is clearly visible as the ~4 Myr oscillation, while $g_1$ is superposed by a longer period.

### 3.3 Relativistic effects

From Table 1, we can observe that there are several planets with $a_p \lesssim 0.5$ au. For such close-in planets the general relativistic precession of the perihelion (GRP) becomes important, like in the case of Mercury in our Solar system (Laskar 2008). We determine this additional precession frequency for the longitude of pericentre in the simplified form taken from Brasser et al. (2009):

$$ \left\langle \frac{d\sigma}{dt} \right\rangle = \frac{3GM}{c^2} \frac{n_p}{a_p(1 - e_p^2)}. \qquad (3) $$

As a check of the formula, we performed direct numerical test calculations of the relativistic equations of motion, solving them either in the full Einstein–Infeld–Hoffmann form (Beutler 2005; Brumberg 2007) or in a more approximate and simplified variant (Bancelin et al. 2012). Results from these calculations agree well with equation (3), so that we use it instead of a more time-consuming numerical integration. We do not include tidal effects between the close-in planets and the host star (see Fabrycky & Tremaine 2007), because the orbits are already quite circular, such that a further circularization is not expected.

Fig. 4 shows the absolute and relative magnitude of the GRP. The upper panel plots contour levels for equation (3) in the $a_p - e_p$ parameter space. The exoplanets are placed according to their parameters in Table 1. In the lower panel of this figure, we plot the ratio of the GRP to the secular precession frequency; for comparison, Mercury is also included. Whenever the GRP is of the same order of magnitude or larger than the secular frequency of the giant

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top: contour levels indicating the order of magnitude of the giant planet's GRP of perihelion (in arcseconds per year) for different binary star systems. Bottom: ratio of the relativistic precession of perihelion relative to the secular precession frequency induced by the binary companion. In the shaded area, the relativistic precession frequency is at least as important as the secular one.</figcaption></figure>

planet, which is indicated by the dashed line and shaded area, we include the general relativistic effect. This figure demonstrates that the GRP can be orders of magnitude larger than the precession frequency induced by purely Newtonian interaction between the giant planet and secondary star. In case of $\tau$ Boo, $\alpha$ Cen, Gliese 86 and GJ 3021 – lying above the dashed line – we have to include this additional contribution.

## 4 RESULTS

It becomes evident from Table 1 that HD 177830 is the only binary with two planets, in all other cases, only a single planet has been detected so far. This leads us to the following questions: can there be yet undetected additional (Earth-like) planets? Would their RV signals be detectable with current detectors?

### 4.1 Observational aspects for finding additional planets

The lack of known Earth-mass planets in the selected binary star systems is a direct consequence of their detection via the RV method (except for Kepler-420). Current RV planet search programs have typical measurement precisions between 0.8 and 5.0 m s$^{-1}$ (Fischer et al. 2016), which favours the detection of massive and/or close-in planets. The discovered exoplanets from Table 1 produce RV signals of 5–500 m s$^{-1}$ with rms fit residuals of 4–20 m s$^{-1}$ (see Butler et al. 2006, Table 3). In contrast, Earth-like and super-Earth-like planets with masses of 1 and 10 M$_\oplus$ at a typical distance of 1 au would produce signals of 0.09 and 0.9 m s$^{-1}$, respectively.
We estimate the maximum RV signal of a terrestrial planet using the approach of Eggl, Haghighipour & Pilat-Lohinger (2013). This method assumes a planet on an initially circular orbit. Perturbations from both the secondary star and the giant planet cause its eccentricity to vary secularly (i.e. over long periods of time) between $e = 0$ and $e = e_{max}$, where the RV signal will be maximal at the latter value:

$$V_r^{\max} = M_P \sin i \sqrt{\frac{G}{a(M_A + M_P)}} \sqrt{\frac{1 + e_{\max}}{1 - e_{\max}}}. \quad (4)$$

Eggl et al. (2012) used an analytic formula due to Georgakarakos (2003) to calculate the value of $e_{max}$, depending on the masses and orbital parameters of the perturbers. Combining these methods, we calculate that planets with $1 M_⊕$ in the HZ of the respective host star never exceed an RV semi-amplitude of $1 \text{ m s}^{-1}$, while planets of $10 M_⊕$ can reach $2 \text{ m s}^{-1}$ only in the innermost part of the HZ.

Consequently, additional terrestrial planets can be hidden in the RV signal residuals and could be detected in the future with more precise measurements. We are aware that very tight binaries can pose a problem as the gas giant might have interfered with the formation of terrestrial planets (Quintana et al. 2007). Nevertheless, Raymond, Quinn & Lunine (2005), Fogg & Nelson (2007) and Ogihara, Kobayashi & Inutsuka (2014) presented promising results that show the feasibility of terrestrial planet formation after migration of a giant planet to orbital distances $≤0.5$ au. Although these simulations covered only discs of single stars, the accretion of a few terrestrial mass objects outside the hot Jupiter’s orbit should be viable even for truncated discs of binaries (Haghighipour & Raymond 2007).

Another possibility is the existence of Neptune- to Jupiter-mass planets at larger orbital distances. Based on the stability limits (cf. Table 2 and Fig. 1), such planets could be imagined for the wider binaries, but excluding HD 126614, HD 196885 and Kepler-420. However, Jupiter-mass planets up to a reasonable distance of about 5 au (for reasons of the long orbital periods) would always produce RV signals of more than $10 \text{ m s}^{-1}$, and hence this kind of planets should have already been detectable so far. Neptune-mass planets (with about $20 M_⊕$) are more difficult to detect and would produce signals of a few metres per second. Anyway, the presence of an additional perturber increases the number of fundamental frequencies and leads to more possible combinations that can become resonant (also the periods become shorter, see the discussion on HD 177830 in Section 4.3). Therefore, the main focus of this survey is on a single giant planet in binary star systems; work on a similar study for two planets in binaries is in progress.

## 4.2 Testing massless versus massive planets

It is common practice to use the test particle approximation for the dynamics of low-mass objects. This approximation is suitable for small mass-ratios, like for asteroids to Jupiter or a terrestrial mass planet to a star. Here we demonstrate that is also applicable in the case of terrestrial planets in binary star systems interacting with a giant planet.

Fig. 5 compares the restricted to the full four-body problem for HD 41004. The black points indicate the initial semimajor axes of massless test particles that start on initially circular orbits in an interval centred on the SR (at $a = 0.37$ au, see Table 2). In the course of the numerical integration, they reach a maximum eccentricity defined as

$$e_{\max} = \max_{t \le T} [e(t)]. \quad (5)$$

<figure><img src="image_7.png" /><figcaption>Figure 5. Maximum eccentricities of massless (black points) and massive (red diamonds) test planets in the system HD 41004, distributed in an interval around the location of the secular resonance.</figcaption></figure>

The curve shows the effect of secular forcing of the eccentricity amplitude when approaching the SR. Particles located further away from the resonance remain on low to moderately eccentric orbits, while those in the vicinity of the resonance are ejected from the system ($e_{max} \ge 1$) in relatively short time. The spike at $a = 0.56$ au is due to the 5:1 MMR with the giant planet and the neighbouring smaller spikes are caused by high-order MMR. The red diamonds represent the maximum eccentricity of Earth-mass planets. Each point was calculated in a separate run to avoid mutual interactions between these planets, which could tamper with the results. It is visible that there is an excellent agreement between the massless and massive planet eccentricities.

## 4.3 Interior versus exterior giant planets

The planetary systems in Table 1 can roughly be distinguished into two classes depending on the location of the giant planet. When the giant planet is exterior to the HZ, like Jupiter in the Solar system, we refer to these systems as ‘exterior cases’. Otherwise, we refer to giant planets located interior to the HZ as ‘interior cases’. Regarding Fig. 2, we have to clarify that we consider HD 196885 and $\gamma$ Cephei as being exterior cases, despite both planets being located exactly inside the (main-sequence) HZ.

Fig. 6 displays the proper secular period depending on the test planet’s semimajor axis. The intervals are chosen in such a way that the test particles never approach the giant planet closer than three Hill radii or alternatively never closer than the pericentre and apocentre distance of the giant planet (for column a and for the other two columns, respectively), whichever is larger. The systems in column (a) all belong to the exterior planet type and, in each case, there is a linear SR interior to the giant planet’s orbit. Contrary to that, for the interior cases in columns (b) and (c), there is no linear SR, except for HD 177830. In the right-hand column, it is visible that the precession periods can decrease appreciably when including the relativistic precession of the perihelion, such that an SR with the close-in planet may occur.

The results suggest that there is a qualitative difference in the dynamics depending on the location of the giant planet. In all of the exterior cases, a single giant planet causes a linear SR at some place interior to its orbit. As opposed to this, a single interior giant planet does not cause a linear SR, at least for $a_P \gtrsim 0.1$ au when the GRP is insufficient to increase the precession frequency by the required amount.

This difference has important implications for the HZ and the possible presence of habitable planets. A hot/warm Jupiter in a
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Plot of the proper secular periods of test planets as a function of their semimajor axis. An intersection of the black curve with the horizontal line indicates the occurrence of a linear secular resonance. The left-hand column (a) collects systems with exterior giant planets, while columns (b) and (c) show systems with a giant planet that is interior to the orbit of the test planet. Column (c) presents cases where the close-in giant planet is affected by the GRP (dashed line without GR, dot-dashed line including GR).</figcaption></figure>

binary star system would not cause an SR in the HZ. Planets in the HZ of Sun-like stars in such systems would only be subject to MMR, but the period ratio to the close-in planet only allows for high-order MMR that have a weak effect. The configuration of giant planets exterior to the orbit of terrestrial planets, based on the architecture of the Solar system, might not be the optimal model for habitable planets. It is well known that Jupiter and Saturn are causing the $v_5$ and $v_6$ SR in the Solar system and the latter plays an important role for perturbing the main belt of asteroids. In a similar way, the accretion of planet-mass objects could be obstructed by exterior giant planets.

Indeed, Bancelin et al. (2015) and Bancelin, Pilat-Lohinger & Bazsó (2016) showed that SR can exist inside the main-belt analogues of certain binary star–giant planet configurations. When the

SR lies in the asteroid belt, it can overlap with MMRs and be responsible for a high rate of asteroids crossing the HZ and therefore for an enhanced impact rate and water transport to the HZ (Bancelin, in preparation). In addition, they compared the water delivery to planets (or planetary embryos) orbiting in the HZ, according to the location of the SR – either inside the HZ or inside the asteroid belt. They concluded the following.

(i) A planet moving on an eccentric orbit (when the SR lies inside the HZ) would increase the probability to collide with incoming wet asteroids. However, the relatively high-impact velocity will strongly reduce the quantity of water finally arriving at the planet’s surface.
(ii) A planet moving on a nearly circular orbit (when the SR lies inside the asteroid belt) can boost its water content due to a higher flux of incoming wet asteroids at relatively low-impact velocities.

But what about multiplanet systems? Do these objections also apply to them? HD 177830 is an analogy to the Solar system, just with the giant planets interior to the HZ. In that case, there are two giant planets whose mutual gravitational interaction increases the apsidal precession frequencies relative to those that would be caused by the secondary star alone. For multiple giant planets, SR are a more common phenomenon and this leads to a large variety of possible planetary system architectures (see Levison & Agnør 2003).

The explanation for the different behaviour of external and internal cases lies in the behaviour of the analytical curve described by equation (2). The secular frequency *g* is the superposition of two functions of *α* that are strictly monotonic. In the exterior case, the periods for the test planets are large nearby the host star for small values of *α*, i.e. the torque is small and the exchange of angular momentum with the giant planet is inefficient. The periods decrease monotonically and approach zero as *a* → *a*<sub>P</sub> (because of the properties of the Laplace coefficient *b*<sub>3/2</sub><sup>(1)</sup>(*α*) when *α* → 1). At some point in between, the curve must necessarily cross the value of *g*<sub>1</sub>, which gives rise to the SR. For the interior case, on the one hand, the periods are first low in the vicinity of the giant planet, but increase monotonically with increasing distance to the planet (i.e. *α* = *a*<sub>P</sub>/a is decreasing). On the other hand, the periods are steadily decreasing under the influence of the secondary star. From the combination of these two components, the secular periods rise to some peak value before they start decreasing again (see Fig. 6). For the systems investigated, the secular periods of the giant planet and the peak value for test planets are differing by more than a factor of '0. Consequently, interior giant planets are unable to cause a linear SR, except for when the GRP is large enough.

## 4.4 Effect of SR

Now let us investigate two typical cases for the effect of the SR for exterior cases. For a more detailed analysis of the dynamical stability, see e.g. Pilat-Lohinger & Funk (2010) and Funk, Pilat-Lohinger & Eggl (2015).

Fig. 7 shows the maximum eccentricities of test particles reached over an integration time of $10^6$ yr in panel (a). All test particles start on initially circular orbits, but those beyond $a \approx 0.9$ au quickly reach a value of $e \ge 1$ and are ejected from the system by a close encounter with the giant planet. At two other locations, we observe a local increase of the maximum eccentricity. These locations correspond to the non-linear<fnref n="3" /> SR $2g - g_P = 0$ and $3g - g_P = 0$, while nearby there are also the 6:1 and 8:1 MMR at $a \approx 0.79$ and 0.65 au, respectively. In the bottom panel (b), we show the intersection points of the test particle's proper secular period curve (dashed black) with multiples of the giant planet's period ($P_{GP}$; the horizontal lines). One can observe that the intersection points correspond to the peaks in maximum eccentricity. This demonstrates that the semi-analytical method is also useful and capable to detect non-linear SR. Note that for this simulation we have reduced the giant planet's initial eccentricity to $e_P = 0.3$, as otherwise the planet's perihelion distance ($q \approx 1.3$ au) and the resonance width of the 3:1 MMR (at $\sim 1.25$ au) would destabilize the majority of test particles in this region.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Effect of the secular resonance on test particles in the HD 196885 system. Over a time interval of $T = 10^6$ yr, the maximum eccentricity rises to unity (top) at the locations where the test particle's proper period (bottom; dashed black curve) is equal to multiples of the giant planet's secular period (horizontal lines).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Effect of the secular resonance on test particles in the $\gamma$ Cephei system for two values of the giant planet's eccentricity (black and grey curves). Top: the maximum eccentricity as a function of the initial semimajor axis. Bottom: normalized FLI of the same region.</figcaption></figure>

Fig. 8 shows the variation of maximum eccentricity (top) and FLI (bottom) as a function of the initial semimajor axis of test particles. In the first case, the giant planet's eccentricity is $e_P = 0.05$ (black curve), which is actually close to the dynamically relaxed value of $e_F \approx 0.056$ for the forced eccentricity obtained by Andrade-Ines et al. (2016) from an analytical model. Increasing the giant planet's eccentricity moderately to $e_P = 0.1$ has a strong effect on the test particles (grey curve). The maximum eccentricities now reach values sufficiently high for ejection from the system. The FLI is a chaos indicator that is sensitive to MMR (Froeschlé et al. 1997); a low FLI value means regular motion, while a high value indicates chaotic motion. Using the FLI, the two peaks of the black curve in panel (b) are clearly identified as corresponding to the 5:1 ($a = 0.71$ au) and 9:2 ($a = 0.75$ au) MMR, respectively.

The main mechanism responsible for the strong increase of the maximum eccentricity is the forced eccentricity component, which

<footnote n="3">³ We follow Knezevic et al. (1991, section 2) for the definition of 'linear' and 'non-linear' SR.</footnote>
depends on the perturber's eccentricity. A minor increase in the giant planet's eccentricity from 0.05 to 0.1 (cf. the black and grey curves in panel a) is sufficient to drive the forced eccentricity to values near unity. This means that already a moderately eccentric giant planet is able to perturb terrestrial planets quite effectively. The overlap of (relatively high-order) MMR and SR (Wisdom 1980) clearly leads to chaotic behaviour and effectively empties the region around the SR. Resonance overlap has also been studied by Mudryk & Wu (2006) for binary systems, but in the case of the secondary star ejecting the exoplanet. They found that secular forcing of the companion shifts the centres of MMR, thus leading to extended resonance overlap regions. We can observe this effect in panel (b) for the 5:1 MMR that is shifted from 0.71 au (for $e_P = 0.05$) to 0.70 au (for $e_P = 0.1$), consistent with the calculations of Mudryk & Wu. This effect is less pronounced for the weaker 9:2 MMR and would only be visible with a higher resolution in the grid of initial conditions.

# 5 DISCUSSION

The semi-analytical method presented in Section 3 relies on a numerical integration of the equations of motion to extract the fundamental frequencies of the perturbers. These numbers are then used in the analytical model to calculate the location of SR. An analytical approach involving simple expressions for the giant planet's secular frequency would facilitate the study of a wide variety of systems, including newly detected ones, without resorting to numerical integrations.

## 5.1 Purely analytical methods

Here we test various analytical methods to calculate the fundamental frequencies $g_j$.

The LL model can also be used to calculate the secular frequencies of an arbitrary number N of mutually interacting massive bodies (see Murray & Dermott 1999, section 7). The resulting frequencies will be reliable only if the assumptions of low eccentricity and inclination are fulfilled.

In the Heppenheimer (1978) model (hereafter HEP), a restricted three-body problem is assumed, where the planet's mass is negligible relative to the masses of the two stars. The disturbing function is limited to terms of $O(e_P^2)$ in the planet's eccentricity, but it can handle an arbitrary eccentricity $e_B$ of the binary (see Andrade-Ines et al. 2016). The expression for the forced secular frequency of the giant planet is

$$g_P = \frac{3}{4} \left( \frac{m_B}{m_A} \right) \left( \frac{a_P}{a_B} \right)^3 n_P (1 - e_B^2)^{-3/2}, \quad (6)$$

where $a_P$ ($n_P$) is the giant planet's semi-major axis (mean motion), $m_A$ and $m_B$ are the stellar masses, $a_B$ is the secondary star's distance from the primary and $e_B$ its eccentricity.

Giuppone et al. (2011) constructed a secular model (hereafter GIU) building on the Heppenheimer model, but extending it to second order in the masses. Instead of using the full analytical expressions, they fitted a simple empirical formula for the secular frequency to their model, which is given by

$$g_P = g_0 \left[ 1 + 32 \left( \frac{m_B}{m_A} \right) \left( \frac{a_P}{a_B} \right)^2 (1 - e_B^2)^{-5} \right], \quad (7)$$

where $g_0$ is given by the expression in equation (6). Strictly speaking, their formula was fitted to the $\gamma$ Cephei system and was not intended to be used elsewhere.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Magnitude of relative errors for the different analytical methods used to determine the fundamental secular frequency of the giant planet. The reference value for the frequency is derived from a numerical integration. Systems are ordered by increasing semimajor axis ratio $a_P/a_B$.</figcaption></figure>

In a series of papers, Georgakarakos (2002, 2003, 2006, 2009) worked out a model (GEO) for hierarchical triple systems, i.e. a close binary star (or star–planet pair) with a much more distant third object. He derived closed formulas for the eccentricity evolution of the inner and outer object, as well as the secular frequency going to second order in the masses. In this model, the eccentricity of the inner object should remain bounded below ~0.2, otherwise the neglected higher order terms in the eccentricity deteriorate the accuracy. The expression for the secular frequency can be written as<fnref n="4" />

$$g_P = g_0 \left[ 1 + \frac{25}{8} \frac{m_B}{\sqrt{m_A(m_A + m_B)}} \left( \frac{a_P}{a_B} \right)^{3/2} \frac{3 + 2e_B^2}{(1 - e_B^2)^{3/2}} \right], \quad (8)$$

where again $g_0$ is given by equation (6). Note that equation (8) contains the complete analytical dependence on masses and semimajor axis ratio to the order of approximation, in contrast to equation (7) that assumes a certain form and is different in important details.

The latter three models all assume a restricted three-body problem ($m_P = 0$). This choice mainly affects and limits the GEO model that is also able to handle a massive planet (see details in Georgakarakos 2003), but is, of course, less accurate when the planet's mass is neglected. We use the restricted problem primarily to permit for a fair comparison between the models.

Using these analytical methods (LL, HEP, GIU, GEO), we calculate the secular frequency of the giant planet and compare the errors relative to the reference value obtained from a numerical simulation (without relativity) in Fig. 9. In the HD 177830 system, there are four massive bodies instead of three as presumed by the HEP, GIU and GEO models; owing to this inherent inability, we do not show HD 177830 in that figure. The LL method shows the largest overall errors, which is due to its construction of neglecting higher order terms in the eccentricity. The one exception to this trend is the Kepler-420 case; here the LL performs best among the four methods. One reason for the mismatch in the analytical frequency for this particularly tight binary system is the combination of the semimajor axis ratio $\alpha = a_P/a_B$ (or orbital period ratio) with the mass of the secondary star. Comparing the other methods, we see that they agree fairly well in most cases. For the largest separation binary, 94 Ceti, all errors are below 10 per cent.

Table 4 compares the performance of the analytical methods for the exterior giant planet cases from Fig. 6. The semi-analytical

<footnote n="4">⁴ Georgakarakos, private communication.</footnote>
<table><caption>**Table 4.** Comparing the location of the secular resonance obtained by combinations of different analytical methods with the reference semi-analytical method. For each pair, the first item indicates the model to calculate the giant planet's frequency, the second item the model for the test planet frequency. Abbreviations: LL = Laplace–Lagrange, HEP = Heppenheimer, GIU = Giuppone, GEO = Georgakarakos.</caption><thead><tr><th rowspan="2">Method</th><th colspan="4">Resonance location (au) for system</th></tr><tr><th>HD 41004</th><th>HD 126614</th><th>HD 196885</th><th>γ Cep</th></tr></thead><tbody><tr><td>Semi-analytical</td><td>0.37</td><td>0.91</td><td>1.11</td><td>0.78</td></tr><tr><td>LL/LL</td><td>0.30</td><td>0.70</td><td>0.87</td><td>0.62</td></tr><tr><td>HEP/HEP</td><td>0.32</td><td>0.98</td><td>1.06</td><td>0.93</td></tr><tr><td>HEP/LL</td><td>0.32</td><td>0.88</td><td>0.99</td><td>0.71</td></tr><tr><td>GIU/GIU</td><td>0.35</td><td>1.08</td><td>1.29</td><td>1.06</td></tr><tr><td>GIU/LL</td><td>0.35</td><td>0.94</td><td>1.17</td><td>0.79</td></tr><tr><td>GEO/GEO</td><td>0.34</td><td>1.02</td><td>1.17</td><td>0.99</td></tr><tr><td>GEO/LL</td><td>0.34</td><td>0.91</td><td>1.07</td><td>0.75</td></tr></tbody></table>

method serves as a reference (cf. locations in Table 2), where we used a combination of a numerical integration of the respective system for the giant planet frequency and the LL model for the test planet frequency. All other entries are pairs of methods, where the first item identifies the model for the giant planet's frequency and the second item the model for the test planet frequencies. Comparing mixed pairs involving the LL method, we see that usually they are closer to the numerical data than pairs of purely analytical methods. This behaviour can be understood when recalling that most importantly, the LL method takes into account the contributions from the giant planet and secondary star (restricted four-body problem), while the HEP, GIU and GEO methods are limited to a restricted three-body problem (giant planet only). Another issue might be the convergence rate of the Legendre expansions; see Andrade-Ines et al. (2016) for a discussion.

## 5.2 Additional improvements

A crucial ingredient is the giant planet's secular frequency, $g_P$; differences in this frequency (cf. Fig. 9) transform into deviations $\Delta a$ of several 0.1 au in the location of the resonance. To improve on the analytical estimate of this frequency, one needs the second-order models in the masses including the higher order terms in the planet's eccentricity.

This path has been followed recently by Andrade-Ines et al. (2016). The authors constructed a second-order model (in the masses) and determined the limits of its validity. The model complexity increased considerably as compared to the first-order model, but they could show that for some of their investigated systems, only the second-order model was capable to correctly describe the dynamics. However, for a few extreme cases, neither the first- nor the second-order model was sufficient; in those cases, a direct numerical integration remains the only viable option. Libert & Henrard (2005) used a similar approach with a high-order expansion of the disturbing function in the eccentricities to cope with moderately to highly eccentric extrasolar systems. Libert & Sansottera (2013) extended the previous work by also moving to a second-order theory with respect to the masses.

Another approach is to stay with a first-order model in the masses, but to use a high-order expansion of the perturbing potential. Georgakarakos, Dobbs-Dixon & Way (2016) derived an expansion up to order 11 in the semimajor axis ratio $\alpha$ valid for a terrestrial planet–giant planet pair. In this way, he obtained refined estimates for the secular frequency and maximum eccentricity of the smaller planet.

So far we have considered only coplanar systems, but the orbital inclination is another important parameter. Pilat-Lohinger et al.

(2016) demonstrate for the HD 41004 system that the secular resonance can be traced numerically for inclinations up to $i \sim 40^{\circ}$. Williams & Faulkner (1981) presented a method to map the SR surfaces in the proper $(a, e, i)$ space. They applied this method to the Solar system and investigated the main linear SR (with Jupiter and Saturn) in the asteroid main belt.

## 6 SUMMARY

We investigated 11 binary star systems with detected planets in circumstellar (S-type) motion and we calculated the location and extent of the HZ for these systems. In most cases, the presence of yet undetected Earth-mass planets cannot be ruled out. Using the semi-analytical method of Pilat-Lohinger et al. (2016), we showed that in two systems (HD 177830 and $\alpha$ Cen), a linear SR is located directly in the HZ. For three other binaries (HD 126614, HD 41004 and HD 196885), we found an SR close to the HZ. Our analysis revealed that giant planets located exterior to the HZ generally cause a linear SR, but not necessarily in the HZ. A moderate eccentricity of the giant planet ($e_P \ge 0.1$) is already sufficient to excite the forced eccentricities of terrestrial planets to values of $e \sim 1$ through the overlapping of the SR with nearby MMRs, where the secular forcing enhances the effect of the MMR. Therefore, it is important to determine the location of the SR accurately, which underlines the usefulness of such studies. In the case of giant planets interior to the HZ, we can exclude SR. However, for close-in giant planets ($a_P \le 0.1$ au), the GRP of the pericentre can still lead to an SR; this is shown for two systems ($\tau$ Boo and $\alpha$ Cen). We compared the accuracy of our semi-analytical method to three analytical models that can be used instead of a numerical integration to calculate the secular frequency of the giant planet. These analytical models are able to determine the resonance location to about $\Delta a \sim 0.1$ au, provided that the LL model is used for the test planet frequencies.

## ACKNOWLEDGEMENTS

The authors acknowledge support by FWF projects S11608-N16 and P22603-N16. Our special thanks go to Nikolaos Georgakarakos for valuable discussions and helpful comments about the application of the analytical models. We also thank our referee for constructive comments that helped to improve the manuscript. For this investigation, we made use of the Extrasolar Planet Encyclopaedia<fnref n="5" /> maintained by J. Schneider (Schneider et al. 2011) and the Catalogue of Exoplanets in Binary Star Systems<fnref n="6" /> maintained by R. Schwarz (Schwarz et al. 2016).

## REFERENCES

* Andrade-Ines E., Beaugé C., Michtchenko T., Robutel P., 2016, Celest. Mech. Dyn. Astron., 124, 405

* Bancelin D., Hestroffer D., Thuillot W., 2012, Celest. Mech. Dyn. Astron., 112, 221

* Bancelin D., Pilat-Lohinger E., Eggl S., Maindl T. I., Schäfer C., Speith R., Dvorak R., 2015, A&A, 581, A46

* Bancelin D., Pilat-Lohinger E., Bazsó Á., 2016, A&A, 591, A120

* Barclay T., Quintana E. V., Adams F. C., Ciardi D. R., Huber D., Foreman-Mackey D., Montet B. T., Caldwell D., 2015, ApJ, 809, 7

* Beaugé C., Nesvorný D., 2012, ApJ, 751, 119

* Beaugé C., Ferraz-Mello S., Michtchenko T. A., 2012, Res. Astron. Astrophys., 12, 1044

<footnote n="5">⁵ http://exoplanet.eu</footnote>

<footnote n="6">⁶ http://www.univie.ac.at/adg/schwarz/multiple.html</footnote>
Beutler G., 2005, Astronomy and Astrophysics Library, Methods of celestial mechanics. Vol. I: Physical, Mathematical, and Numerical Principles. Springer, Berlin

Boss A. P., 2006, ApJ, 641, 1148

Brasser R., Morbidelli A., Gomes R., Tsiganis K., Levison H. F., 2009, A&A, 507, 1053

Bromley B. C., Kenyon S. J., 2015, ApJ, 806, 98

Brumberg V., 2007, Celest. Mech. Dyn. Astron., 99, 245

Butler R. P., Marcy G. W., Williams E., Hauser H., Shirts P., 1997, ApJ, 474, L115

Butler R. P. et al., 2006, ApJ, 646, 505

Campante T. L. et al., 2015, ApJ, 799, 170

Chabrier G., 2001, ApJ, 554, 1274

Chambers J. E., 1999, MNRAS, 304, 793

Chauvin G., Lagrange A.-M., Udry S., Mayor M., 2007, A&A, 475, 723

Chauvin G., Beust H., Lagrange A.-M., Eggenberger A., 2011, A&A, 528, A8

Correia A. C. M. et al., 2008, A&A, 479, 271

Cuntz M., 2014, ApJ, 780, 14

da Silva L. et al., 2006, A&A, 458, 609

Desidera S., Barbieri M., 2007, A&A, 462, 345

Dodson-Robinson S. E., Veras D., Ford E. B., Beichman C. A., 2009, ApJ, 707, 79

Duchêne G., Kraus A., 2013, ARA&A, 51, 269

Dumusque X. et al., 2012, Nature, 491, 207

Duquennoy A., Mayor M., 1991, A&A, 248, 485

Dvorak R., Froeschle C., Froeschle C., 1989, A&A, 226, 335

Eggl S., Dvorak R., 2010, in Souchay J., Dvorak R., eds, Lecture Notes in Physics, Vol. 790, Dynamics of Small Solar system Bodies and Exoplanets. Springer Verlag, Berlin, p. 431

Eggl S., Pilat-Lohinger E., Georgakarakos N., Gyergyovits M., Funk B., 2012, ApJ, 752, 74

Eggl S., Haghighipour N., Pilat-Lohinger E., 2013, ApJ, 764, 130

Endl M., Cochran W. D., Hatzes A. P., Wittenmyer R. A., 2011, in Schuh S., Drechsel H., Heber U., eds, AIP Conf. Series Vol. 1331, Planetary Systems Beyond the Main Sequence. Am. Inst. Phys., New York, p. 88

Endl M. et al., 2015, Int. J. Astrobiology, 14, 305

Everhart E., 1974, Celest. Mech., 10, 35

Fabrycky D., Tremaine S., 2007, ApJ, 669, 1298

Fischer D. A. et al., 2016, PASP, 128, 066001

Fogg M. J., Nelson R. P., 2007, A&A, 472, 1003

Forgan D., 2016, MNRAS, 463, 2768

Frigo M., Johnson S. G., 2005, Proc. IEEE, 93, 216

Froeschlé C., Lega E., Gonczi R., 1997, Celest. Mech. Dyn. Astron., 67, 41

Fuhrmann K., Chini R., Buda L.-S., Pozo Nuñez F., 2014, ApJ, 785, 68

Funk B., Pilat-Lohinger E., Eggl S., 2015, MNRAS, 448, 3797

Georgakarakos N., 2002, MNRAS, 337, 559

Georgakarakos N., 2003, MNRAS, 345, 340

Georgakarakos N., 2006, MNRAS, 366, 566

Georgakarakos N., 2009, MNRAS, 392, 1253

Georgakarakos N., Dobbs-Dixon I., Way M. J., 2016, MNRAS, 461, 1512

Giuppone C. A., Leiva A. M., Correa-Otto J., Beaugé C., 2011, A&A, 530, A103 (GIU)

Gould A. et al., 2014, Science, 345, 46

Haghighipour N., Raymond S. N., 2007, ApJ, 666, 436

Hanslmeier A., Dvorak R., 1984, A&A, 132, 203

Hatzes A. P., 2013, ApJ, 770, 133

Hatzes A. P., Cochran W. D., Endl M., McArthur B., Paulson D. B., Walker G. A. H., Campbell B., Yang S., 2003, ApJ, 599, 1383

Heppenheimer T. A., 1978, A&A, 65, 421 (HEP)

Holman M. J., Wiegert P. A., 1999, AJ, 117, 621

Howard A. W. et al., 2010, ApJ, 721, 1467

Jang-Condell H., 2015, ApJ, 799, 147

Jofré E., Petrucci R., Saffe C., Saker L., de la Villarmois E. A., Chavero C., Gómez M., Mauas P. J. D., 2015, A&A, 574, A50

Kaltenegger L., Haghighipour N., 2013, ApJ, 777, 165

Kasting J. F., Whitmire D. P., Reynolds R. T., 1993, Icarus, 101, 108

Kley W., Nelson R. P., 2010, in Haghighipour N., ed., Astrophysics and Space Science Library Vol. 366, Planets in Binary Star Systems. Springer, Netherlands, Dordrecht, p. 135

Knezevic Z., Milani A., Farinella P., Froeschle C., Froeschle C., 1991, Icarus, 93, 316

Kopparapu R. K. et al., 2013, ApJ, 765, 131

Laskar J., 2008, Icarus, 196, 1

Levison H. F., Agnor C., 2003, AJ, 125, 2692

Libert A.-S., Henrard J., 2005, Celest. Mech. Dyn. Astron., 93, 187

Libert A.-S., Sansottera M., 2013, Celest. Mech. Dyn. Astron., 117, 149

Liebert J., Bergeron P., Holberg J. B., 2005, ApJS, 156, 47

Marois C., Macintosh B., Barman T., Zuckerman B., Song I., Patience J., Lafrenière D., Doyon R., 2008, Science, 322, 1348

Mayor M., Udry S., Naef D., Pepe F., Queloz D., Santos N. C., Burnet M., 2004, A&A, 415, 391

Naoz S., Farr W. M., Rasio F. A., 2012, ApJ, 754, L36

Ogihara M., Kobayashi H., Inutsuka S.-i., 2014, ApJ, 787, 172

Müller T. W. A., Kley W., 2012, A&A, 539, A18

Murray C. D., Dermott S. F., 1999, Solar System Dynamics. Cambridge Univ. Press, Cambridge

Naef D., Mayor M., Pepe F., Queloz D., Santos N. C., Udry S., Burnet M., 2001, A&A, 375, 205

Naoz S., Farr W. M., Rasio F. A., 2012, ApJ, 754, L36

Ogihara M., Kobayashi H., Inutsuka S.-i., 2014, ApJ, 787, 172

Pilat-Lohinger E., Dvorak R., 2002, Celest. Mech. Dyn. Astron., 82, 143

Pilat-Lohinger E., Funk B., 2010, in Souchay J., Dvorak R., eds, Lecture Notes in Physics, Vol. 790, Dynamics of Small Solar System Bodies and Exoplanets. Springer Verlag, Berlin, p. 481

Pilat-Lohinger E., Bazsó Á., Funk B., 2016, AJ, 152, 139

Poleski R. et al., 2014, ApJ, 795, 42

Queloz D. et al., 2000, A&A, 354, 99

Quintana E. V., Adams F. C., Lissauer J. J., Chambers J. E., 2007, ApJ, 660, 807

Rabl G., Dvorak R., 1988, A&A, 191, 385

Raghavan D., Henry T. J., Mason B. D., Subasavage J. P., Jao W.-C., Beaulieu T. D., Hambly N. C., 2006, ApJ, 646, 523

Raghavan D. et al., 2010, ApJS, 190, 1

Rajpaul V., Aigrain S., Roberts S., 2016, MNRAS, 456, L6

Ramirez R. M., Kaltenegger L., 2016, ApJ, 823, 6

Raymond S. N., Quinn T., Lunine J. I., 2005, Icarus, 177, 256

Reegen P., 2007, A&A, 467, 1353

Roberts L. C., Jr, et al., 2015, AJ, 150, 103

Rodler F., Lopez-Morales M., Ribas I., 2012, ApJ, 753, L25

Roell T., Neuhäuser R., Seifahrt A., Mugrauer M., 2012, A&A, 542, A92

Rowe J. F. et al., 2014, ApJ, 784, 45

Santerne A. et al., 2014, A&A, 571, A37

Schneider J., Dedieu C., Le Sidaner P., Savalle R., Zolotukhin I., 2011, A&A, 532, A79

Schwarz R., Funk B., Zechner R., Bazsó Á., 2016, MNRAS, 460, 3598

Soubiran C., Le Campion J.-F., Brouillet N., Chemin L., 2016, A&A, 591, A118

Thebault P., Haghighipour N., 2015, in Jin S., Haghighipour N., Ip W.-H., eds, Planetary Exploration and Science: Recent Results and Advances. Springer Geophysics, Springer-Verlag, Berlin Heidelberg. p. 309

Tokovinin A., 2014, AJ, 147, 87

Torres G., 2007, ApJ, 654, 1095

van Leeuwen F., 2007, A&A, 474, 653

Vogt S. S., Marcy G. W., Butler R. P., Apps K., 2000, ApJ, 536, 902

Vorobyov E. I., 2013, A&A, 552, A129

Wang J., Fischer D. A., Xie J.-W., Ciardi D. R., 2014, ApJ, 791, 111

Williams J. G., Faulkner J., 1981, Icarus, 46, 390

Williams D. M., Pollard D., 2002, Int. J. Astrobiology, 1, 61

Winn J. N., Fabrycky D. C., 2015, ARA&A, 53, 409

Wisdom J., 1980, AJ, 85, 1122

Zucker S., Mazeh T., Santos N. C., Udry S., Mayor M., 2004, A&A, 426, 695

This paper has been typeset from a TeX/LATEX file prepared by the author.