title: Fear the Shadows of the Giants: On Secular Perturbations in Circumstellar Habitable Zones of Double Stars
authors: Ăkos BazsĂł Elke Pilat-Lohinger
text: 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Ákos Bazsó
  akos.bazso@univie.ac.at
  
  
  
  
      Department of Astrophysics 
  
      University of Vienna 
  
      Türkenschanzstrasse 17 
  
      A-1180 Vienna, Austria
  
  
  
  
      Department of Astrophysics 
  
      University of Vienna 
  
      Türkenschanzstrasse 17 
  
      A-1180 Vienna, Austria
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  After the detection of extrasolar planets in binary and multiple star systems questions arose about their dynamics and habitability.
  
  In this study we investigate a five-dimensional parameter space for wide binary stars with a massive planet beyond the habitable zone (HZ).
  Our aim is to reveal those orbital and physical parameter combinations that expose bodies in the primary star's HZ to secular perturbations.
  
  Building upon an established semi-analytical model, we combine two separate analytical models into a new one to treat the restricted four-body problem.
  We apply this new Combined Analytical Model (CAM) to different synthetic systems and systematically map the occurrence of secular resonances (SR).
  These maps are then visualized as two-dimensional sections of the parameter space.
  
  The CAM model has a median error below 3 % relative to numerical reference simulations.
  We also derive a simplified CAM that performs well for hierarchical systems with moderate to large separations between the bodies.
  Our results show that SR appear in the HZ even for large secondary star distances (up to > 1000 au) if either (i) the planet's distance is larger than Jupiter's, or (ii) its mass is about Saturn's or lower.
  Changes in the secondary star's eccentricity by stellar flybys or galactic tides can push a formerly dynamically quiet HZ to a high-eccentricity state.
  
  Based on these results we provide the easy-to-use online tool `SHaDoS' that evaluates the CAM for a given user input and traces the locations of SR in the HZ in two-dimensional parameter space plots.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
   
  
  
  
  
  Observational surveys revealed that a considerable fraction of stars in the solar neighbourhood are members of binary and multiple star systems.
  In this context, <cit.> and <cit.> established that in the solar neighbourhood 1up to  distances of d ≤ 25 pc about 40-45 % of all Sun-like stars 1with spectral types F6–K3 are members of binary and multiple star systems, independent of whether or not they are hosting extrasolar planets (or short exoplanets).
  <cit.> found for a sample of about 4800 F- and G-type main-sequence stars within 67 pc of the Sun that 33 % of the targets belong to binary star systems.
  Further details on the frequency and main characteristics of multiple star systems can be found in the review of <cit.>.
  
  
  
  
  The detection of exoplanets also in binary and multi-stellar systems showed that planetary companions are not restricted to single stars.
  However, it is uncertain how favourable stellar multiplicity really is for the presence of exoplanets.
  Based on the discovered exoplanet population <cit.>, <cit.>, <cit.>, and <cit.> estimated the stellar multiplicity.
  These studies yield that about 10-15 % of the detected exoplanets reside in binaries, while 1only ∼ 2 % of all 1known exoplanet host stars are members of multiple star systems with three or more components.
  The Catalogue of Exoplanets in Binary and Multiple Star Systems[https://www.univie.ac.at/adg/schwarz/multiple.html] of <cit.> provides an up-to-date compilation of all binary and multiple star systems hosting exoplanets.
  
  
  
  
  Although many exoplanets have been found in both binary and multiple star systems, it is unclear whether or not these environments are more hostile than single star systems for the presence of planets.
  There are different views on the efficiency of planet formation in binary stars and the impact of the stellar companion's perturbations.
  On the one hand, simulations of protoplanetary disks were able to form both terrestrial and giant planets even in tight binary systems, see e.g. <cit.>.
  These studies relate to rather tight systems with stellar separations of < 100 astronomical units (au), such as the poster cases γ Cephei and HD 41004 with ∼ 20 au.
  On the other hand there are indications that the occurrence rate of circumstellar planets is strongly influenced by close binary companions <cit.>.
  For more comprehensive reviews on planet formation in binary star systems refer to <cit.>, <cit.>, and <cit.>.
  
  
  
  
  Apart from the circumstances of planet formation in binaries also the long-term stability of planetary systems is an issue.
  In this study we focus exclusively on circumstellar (S-type) planetary systems of binary stars <cit.>.
  Different kinds of resonances play a key role for the orbital evolution of the planets, 1e.g. mean motion and secular resonances.
  
  
  Mean motion resonances (MMR) 1involve combinations of orbital frequencies and their overlap has already been identified to be crucial in triggering chaos <cit.>.
  The occurrence of circumstellar planets depends sensitively on the interactions with the outer stellar companion and on its orbital parameters <cit.>.
  1
  All these studies have focused on the stability of prograde planets though, but retrograde resonances offer a stabilizing mechanism beyond the established stability borders for the prograde case.
  <cit.> and <cit.> studied retrograde resonances for binary stars in the context of the circular spatial restricted three-body problem.
  These works were based on the study of <cit.> who demonstrated that retrograde MMR can stabilize compact multi-planetary systems.
  In the former studies of binary systems the (retrograde) resonant terms were identified in the expansion of the disturbing function and two main aspects of the stabilizing mechanism were discussed: first, the eccentricity forcing for a prograde planet is larger than for a retrograde one (i.e. different exponents for the eccentricity in the leading term of the expansion), and second the overlap among nearby resonances occurs at larger perturber masses for retrograde planets.
  
  
  
  Secular interactions 1are related to apsidal and nodal precession frequencies and constitute another way how the stellar companion can influence circumstellar planetary systems.
  1
  Similar to MMR, SR affect not only prograde but retrograde orbits as well, and can even induce a switch between these two types of orbital motion.
  <cit.> investigated the secular dynamics of hierarchical three-body systems and derived an octupole order secular approximation to the equations of motion.
  From their model they found a qualitatively different dynamical evolution for an inner body that can become highly eccentric and would then oscillate between prograde and retrograde motion when it is perturbed by an eccentric outer body.
  <cit.> applied a semi-analytical model to investigate SR for prograde and retrograde orbits in a restricted four-body problem.
  They identified the critical angles of the strongest SR that are formed by linear combinations of apsidal and nodal differences of the two adjacent orbits.
  
  In general, even initially circular and coplanar planetary orbits can be strongly altered by secular perturbations from a massive, eccentric and inclined outer companion.
  <cit.> investigated the angular momentum transfer from a stellar companion to weakly or rigidly coupled planetary systems.
  They identified that this cascaded angular momentum transfer through the system leads to a significant eccentricity excitation for the innermost planet for a certain range of outer planet parameters.
  In a similar study, <cit.> found that eccentricities and mutual inclinations of weakly coupled multi-planet systems are more effectively excited by external companions than for their strongly coupled counterparts (i.e. compact planetary systems).
  
  
  Still another type of resonant forcing exerted by the binary companion involves the evection resonance <cit.> that combines the orbital precession frequency of an inner planetary orbit with the orbital motion of the outer secondary star (or another perturber).
  This process can lead to large amplitude oscillations of the planetary eccentricities and mutual inclinations, or even to the disruption of the planetary system.
  Recently, <cit.> demonstrated that the outer stability border of planetary orbits in multi-planet systems of binaries is destabilized by overlap of MMR and SR (both linear and non-linear ones).
  
  
  
  
  In spite of all these different mechanisms to disrupt planetary systems it is not hopeless to find terrestrial planets in binaries, preferably even in the habitable zone (HZ).
  <cit.> examined four different tight binary systems with a known massive planet to explore their potential of being able to host additional habitable terrestrial planets.
  They concluded that two systems are promising candidates and permit dynamically stable orbits in the HZ.
  However, the identification of suitable systems requires a thorough knowledge of the locations of resonances.
  In case of MMR this is trivial, but it is more difficult to find the locations of SR.
  To tackle this task, <cit.> and <cit.> 1described a semi-analytical method that simplifies matters.
  That method determines the orbital precession frequencies for the perturbers and allows to estimate the positions of linear SR for test particles.
  However, the method is valid for the system as it is, i.e. for the instantaneous orbital parameters of the perturbers; if the parameters change with time or by improved observations, then the analysis would need to be repeated.
  
  
  In this work we aim at a more global description of the dynamics in a binary star system with a giant planet.
  Our goal is to reveal parameter combinations that expose bodies (terrestrial planets) in the primary star's HZ to secular perturbations.
  These parameters cover all relevant observables, in particular masses, orbital distances, and eccentricities, such that observational uncertainties are implicitly included into the analysis.
  We investigate this multi-dimensional parameter space by means of a combination of analytical models.
  This allows to quickly identify 1the locations of SR and to 1indicate whether to safely exclude systems with a-priori perturbed HZs 1or to perform a subsequent detailed stability analysis 1with complementary tools.
  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 planets in the HZ in low eccentricity orbits.
  This, in turn, renders the method especially valuable 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.
  
  
  
  
  The outline of the paper is as follows.
  We describe the basic dynamical model along with the analytical methods employed for the parameter search in section <ref>.
  Subsequently, in section <ref>, we first investigate the performance of the analytical models relative to reference data from numerical simulations, and then present the results on the SR location depending on the system parameters.
  Section <ref> contains a discussion of the results and how they compare to other studies.
  Finally, section <ref> highlights the main results.
  Additionally, in Appendix <ref> we introduce an online tool that is based on the models from section <ref> and can be used to create plots similar to those in this article.
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  § THEORY AND MODELS
   
  
  Our main goal is to identify those combinations of physical and orbital parameters for binary star systems with a circumstellar giant planet that lead to secular perturbations inside the host star's habitable zone.
  First, we provide the HZ model and state the HZ location and its extent in section <ref>.
  Next, we introduce the dynamical model for the secular evolution of the system in section <ref>, and define the main relation for the secular frequencies.
  1On the basis of this relation we show how to use the combined analytical model (CAM) that consists of two consecutive steps, each of them given by a simple analytical expression, that are described in sections <ref> and <ref>.
  Finally, we demonstrate some of the properties of the CAM and derive a simplified form for it.
  1Table <ref> contains a list of the variable names used in this work, together with a short description for a quick overview.
  
  
  ll
      Glossary of important variables used in this work.
      Variable     Description
      a_B             semi-major axis of secondary star (component B) 
  
      a_GP            semi-major axis of giant planet 
  
                  inner border of host star's habitable zone 
  
                 outer border of host star's habitable zone 
  
      a_TP            semi-major axis of test planet 
  
      α            semi-major axis ratio, e.g. a_GP / a_B 
  
      b_3/2^(1)     Laplace coefficient 
  
      γ            mass ratio, e.g. m_B / m_A 
  
      δ_g        empirical correction term for the secular frequency g, see eq:andeggcorr 
  
      e_B             eccentricity of secondary star (component B) 
  
      e_GP            eccentricity of giant planet 
  
      g_GP            secular frequency of giant planet 
  
      g_TP            secular frequency of test planet 
  
      m_A             mass of primary star (component A) 
  
      m_B             mass of secondary star (component B) 
  
      m_GP            mass of giant planet 
  
      m_TP            mass of test planet, equals zero by definition 
  
      n_GP            mean motion of giant planet 
  
      n_TP            mean motion of test planet 
  
      
  
  
  
  
  
   §.§ Habitable zone
   
  
  There exist various models for calculating HZs for binary stars and in particular circumstellar HZs, e.g. <cit.>, <cit.>, <cit.>, and its update in <cit.>.
  A common point of all these models is that they extend the classical concept of habitable zones of single stars, see e.g. <cit.> or <cit.>, by taking into account multiple sources of radiation, and to some extent also the dynamical evolution of objects.
  
  In case of circumstellar planets, the HZ location still primarily depends on the luminosity (mass) and effective temperature of the host star.
  Thus, from a practical point of view, it is permissible to approximate the binary star HZ as that of a single star HZ (SSHZ), because in the following we will typically deal with rather large stellar separations (on the order of some hundred au).
  In such cases the secondary's contribution is practically negligible; for stellar separations > 20 au the binary star HZ borders typically shift by < 0.1 au relative to the SSHZ.
  
  For our purposes we fix the host star's spectral type[Without loss of generality; this choice only affects the HZ limits that can be easily calculated for other spectral types (masses).] to G2 V, which gives a solar analogue star with mass 1   M_⊙, luminosity 1   L_⊙, and effective temperature 5780 K.
  With this set of parameters we obtain a HZ that ranges from the inner HZ border (IHZ), = 0.95 au, to the outer HZ border (OHZ), = 1.68 au, according to the model of <cit.>.
  
  
  
  
  
   §.§ Dynamical model
   
  
  Let us assume a binary star system that consists of two stellar components and a (Jupiter-like) giant planet (GP), 1with the restrictions that all bodies are in coplanar and prograde motion.
  We will denote the planet hosting star (component A) the primary star, while the distant perturber (component B) will be named the secondary star.
  The planet is in circumstellar motion (S-type; ) about the host star, and we further assume that it is located exterior to the habitable zone.
  Figure <ref> shows a sketch of this situation.
  In the figure, the horizontal errorbars indicate the range of distances from pericenter to apocenter that each object can cover due to its eccentricity.
  The secondary star must maintain a minimum distance from the primary, in order for the giant planet to remain inside the host star's stability region for all times (see , , and  for details).
  
  
  
  
  A binary star with a circumstellar giant planet forms a hierarchical system.
  In terms of dynamical perturbations, the distant secondary star perturbs the giant planet, which itself would pass the perturbations on to any bodies in the HZ.
  By virtue of the vastly different masses of the bodies entangled in this interaction, an initially small perturbation from the secondary star can result in huge perturbations in the HZ.
  The focus is on secular interactions, where orbital precession frequencies play a major role over extended periods of time, as opposed to short period effects by orbital resonances.
  
  Foremost, the strongest secular perturbations arise 1at locations where the linear combinations of a test planet's orbital precession frequency, g_TP, and that of the giant planet, g_GP, form small numbers
  
      q   g_TP - p   g_GP≈ 0,
  
  1i.e. when the two frequencies are close to form an integer combination g_TP / g_GP≈ p / q (with p, q ∈ℤ).
  This relation forms the basis for all subsequent investigations.
  Similar to MMR, the strongest effects arise when p = q = 1, 1which gives rise to a linear SR.
  1The general case with p ≠ q (p > 0, q > 0) is called a non-linear secular resonance <cit.>.
  1
  Any given SR's strength depends on the specific frequency combination (p, q).
  In case of apsidal frequencies it depends on different powers of the eccentricities,
  according to <cit.> the leading term is ∝ e_TP^q e_GP^p.
  So we can appreciate that for fixed (e_TP, e_GP) the SR with p = q = 1 will always be dominating, but the p = 2, q = 1 and p = 1, q = 2 resonances might show very different effects when e_TP≠ e_GP.
  
  
  1In case of an exact resonance we can recast eq:freqequal into the form q   g_TP = p   g_GP.
  The right-hand side 1of this equation constitutes a three-body problem including the two stars and a giant planet; if the latter is taken to be massless, the dynamical model reduces to a restricted three-body problem.
  Conversely, the left-hand side inevitably forms a restricted four-body problem, which follows from our convention to set test planet masses m_TP = 0.
  This duality makes it difficult to treat eq:freqequal using a single formula or secular model.
  An unified formula would need to be applicable to both the three- and four-body problem and to be valid for a large range of parameters with high accurracy.
  
  In fact, up to now our method of choice was the semi-analytical model (SAM) 1described in <cit.> and <cit.>.
  The SAM starts from the basic equation above and determines the exact location for which q   g_TP - p   g_GP = 0 is fulfilled in terms of a_TP.
  It treats the global restricted four-body problem binary-star – giant planet – test planet as two coupled three-body problems.
  In a first step, the giant planet's precession frequency g_GP in eq:freqequal is determined numerically by frequency analysis from a numerical simulation of the binary-star – giant planet subsystem.
  Then the calculation of g_TP is performed for fixed parameters of the secondary star and giant planet with the determined (and constant) value of g_GP by means of any suitable analytical formula.
  
  The trade-off with an approach like the SAM is that it is less flexible when it comes to scan through a multi-dimensional parameter space.
  However, a combination of analytical models for both frequencies in eq:freqequal would facilitate such studies and allow to investigate a large variety of binary star configurations.
  In <cit.> we have investigated a number of analytical secular models for the giant planet precession frequency, g_GP.
  It turned out that two models are promising candidates for this task, namely the models of <cit.> and <cit.>.
  The latter one provides a framework for the eccentricity evolution (on short and secular time-scales) in the hierarchical full three-body problem, while the former one is only valid for the restricted problem.
  In the following, we describe the transition from the SAM to the new Combined Analytical Model (CAM) that involves a combination of two different analytical models that are connected in order to solve eq:freqequal.
  
  
  
  
  
    §.§.§ Analytical model for giant planet frequency
   
  
  <cit.> established a simple and accurate secular model for the restricted three-body problem, i.e. the giant planet acts like a massless body relative to the two stars.
  Their model is also applicable in case of strongly perturbed three-body systems, i.e. when the secondary star's distance from the primary is just two to three times that of the giant planet.
  The main relation is given by
  
      g_GP = g_H (1 - δ_g),
  
  where g_H is a first order approximation (in masses) to g_GP due to Heppenheimer's formula <cit.>,
  
      g_H = 3/4 n_GP( m_B/m_A) ( a_GP/a_B)^3 (1 - e_B^2)^-3/2,
  
  and n_GP is the giant planet's mean motion.
  The term δ_g is an empirical correction to the frequency, and has the general form
  
      δ_g = ∑_k = 1^M A_k ( m_B/m_A)^p_k( a_GP/a_B)^q_k e_B^r_k ,
  
  with typically M ≤ 20 coefficients A_k and various powers (p_k, q_k, r_k); refer to <cit.> for the actual values of these coefficients and exponents.
  This kind of correction to the secular precession frequency follows the spirit of previous attempts, e.g. <cit.>.
  
  The model expresses the secular precession frequency g_GP as a polynomial function of the basic parameters mass ratio γ = m_B / m_A, semi-major axis ratio α = a_GP / a_B, and perturber eccentricity e_B.
  Formally, <cit.> state the validity ranges of their fit as 0.1 ≤γ≤ 10, 0.01 ≤α≤ 0.4, and 0 ≤ e_B≤ 0.6.
  This interval of mass ratios covers the full range of interest for binary star systems with a G-type primary component, from the lowest mass M-dwarves (for γ = 0.1) up to F- and A-type stars (for γ≈ 2).
  Larger values of γ would lead to an exceedingly fast stellar evolution of the secondary, such that it were no more a main-sequence star.
  The interval for the distance ratio α also covers most – if not all – currently observed exoplanets in binary systems, including cases like γ Cephei (α = 0.1) and HD 196885 (α = 0.12), see <cit.> for a list of tight binary systems.
  One remaining caveat is that for some observed binary star systems the secondary's eccentricity might be larger than e_B = 0.6 (e.g. for τ Bootis)[Another issue stems from the fact that the Washington Double Star catalogue lists less than 1000 systems (from a total of more than 130,000) with well-constrained stellar eccentricity. At larger stellar separations (above 100 au, i.e. in the distance range considered here) there is a deficit of observational information on the eccentricities. Therefore, there is a clear need to estimate perturbing effects on the HZ due to stellar eccentricities.].
  
  
  
  
  
    §.§.§ Analytical model for test planet frequency
   
  
  According to the Laplace-Lagrange model (LLM) the secular precession frequency for a massless test planet (TP) influenced by N massive perturbers is given by
  
      g_TP = n_TP/4∑_i=1^N( m_i/m_A)
                   ( a_TP/a_i)^2
                   b_3/2^(1)( a_TP/a_i)
                  (1 - e_i^2)^-3/2,
  
  where n_TP is the test planet's mean motion, and the function b_3/2^(1) (α_i) is a Laplace coefficient of the single scalar variable 0 ≤α_i = a_TP / a_i < 1 <cit.>.
  For our dynamical model with N = 2 (giant planet and secondary star) it follows that
  
      g_TP = n_TP/4   [
                  ( m_GP/m_A)
                  ( a_TP/a_GP)^2
                  b_3/2^(1)( a_TP/a_GP)
                  (1 - e_GP^2 )^-3/2
              + . 
         . +
                  ( m_B/m_A)
                  ( a_TP/a_B)^2
                  b_3/2^(1)( a_TP/a_B)
                  (1 - e_B^2)^-3/2].
  
  Note that this is an ad-hoc modification of the classical LLM.
  Originally, the classical secular theory is of first-order in the eccentricity and hence is valid only for low values of e_i.
  However, by adding the extra terms (1 - e_i^2)^-3/2 with the eccentricity dependence we can extend its range of validity.
  The motivation is that eq:laplag looks and behaves similar to eq:heppenh with the additional terms added.
  In fact, <cit.> have shown that such a modification acts as an interpolating formula between the classical LLM (in the limit e_i→ 0) and the Heppenheimer model (for a_TP / a_i→ 0).
  
  Equation (<ref>) can handle arbitrary mass ratios and is applicable for large distance ratios up to the limit a_TP / a_i→ 1, which becomes important when the giant planet approaches the outer HZ border (OHZ), but we strictly require that < a_GP.
  Alternatively, the LLM can be replaced by any suitable secular model that can handle at least a four-body problem, e.g. the model of <cit.>.
  
  
  
  
  
    §.§.§ Combined Analytical Model (CAM)
   
  
  After having introduced both analytical sub-models, we now combine equations (<ref>) and (<ref>) to express eq:freqequal in terms of the five parameters (m_GP, m_B, a_GP, a_B, e_B).
  We assume that the eccentricity of the giant planet (e_GP) is small enough to neglect its effect, i.e. e_GP→ 0, which reduces the dynamical system by this variable.
  
  In order to have a more compact notation, we introduce the following abbreviations: γ_1 := m_GP / m_A, γ_2 := m_B / m_A, α_1 := a_TP / a_GP, α_2 := a_GP / a_B, so that we can express a_TP / a_B = (a_TP / a_GP) (a_GP / a_B) = α_1 α_2.
  Furthermore, the eccentricity e_B appears on both sides of the equation in form of the expression ϵ_B := (1 - e_B^2)^-3/2.
  
  With all these definitions we can write eq:freqequal 1for linear[Without loss of generality we focus here on linear secular resonances, because these show the maximum effect in terms of forced eccentricities (cf. Fig. <ref>). In a more general setting, the coefficients (p, q) can be easily reintroduced into eq:freqeqmain to study also nonlinear SR (see Fig. <ref>).] SR (p = q = 1) in terms of (γ_1, γ_2, α_1, α_2, ϵ_B) to obtain
  
      1/4 n_TP(
                  γ_1α_1^2 b_3/2^(1) (α_1) +
                  γ_2( α_1α_2)^2 b_3/2^(1) (α_1α_2) ϵ_B) =
      3/4 n_GPγ_2α_2^3ϵ_B (1 - δ_g).
  
  We will refer to eq:freqeqmain from now on as Combined Analytical Model (CAM).
  
  The expression for δ_g is given in eq:andeggcorr as a function of (γ_2, α_2, e_B), where e_B appears explicitly.
  From now on we skip the subscript g from δ_g, and are going to use simply δ when referring to eq:andeggcorr.
  Note that the frequencies (mean motions) n_TP = n(a_TP) and n_GP = n(m_GP, a_GP) are the only dimensional quantities in this equation, and that they are known functions of the system parameters.
  
  
  
  
  
    §.§.§ Derivation of simplified CAM
   
  
  Naturally, it would be helpful to have at hand some sort of analytical expressions that allow to investigate the qualitative behaviour of the CAM solutions.
  These expressions are not available in general, but they can be found for an approximate form of eq:freqeqmain that we will derive now.
  
  The primary obstacle in eq:freqeqmain are the Laplace coefficients b_3/2^(1).
  We can get rid of them by assuming that the conditions α_1≪ 1 and α_2≪ 1 hold, such that we can replace these functions by their MacLaurin series expansions
  
      b_3/2^(1)(α) ≈ 3 α + 45/8α^3 + 𝒪(α^5)
  
  up to the lowest order in α.
  This approximation is justified because of the hierarchical nature of the dynamical system.
  In case of well separated orbits we expect quite small values of α < 0.1, for which this approximation yields relative errors of roughly 1 - 2 % compared to the values evaluated numerically from the full form of the Laplace coefficient.
  Deviations of this magnitude are tolerable given the inherent fit-errors of the secular model in eq:andeggmain, that itself has median errors of 2 - 3 % relative to a fully numerical solution (compare section <ref>).
  
  Starting from eq:freqeqmain, we replace the Laplace coefficients by their approximations, insert n_TP = (G m_A / a_TP^3)^1/2 and n_GP = (G m_A (1 + γ_1) / a_GP^3)^1/2, and eliminate common factors on both sides of the equation to obtain
  
      γ_1α_1^3 + γ_2α_1^3α_2^3ϵ_B = (1 + γ_1)^1/2γ_2α_1^3/2α_2^3ϵ_B (1 - δ).
  
  This equation, termed the simplified CAM, represents an approximate form of eq:freqeqmain that shares all of its properties in the limit α_1,2→ 0.
  Appendix <ref> provides (approximate) analytical solutions to this simplified CAM with respect to the parameters (γ_1, γ_2, α_1, α_2, ϵ_B).
  
  
  
  
  
   §.§ Properties of the CAM
   
  
  We are now able to investigate the binary-star – giant planet – test planet system as a coupled dynamical problem given by eq:freqeqmain depending on five parameters.
  The goal is to identify combinations of those parameters that would lead to secular perturbations inside the HZ, i.e. for ≤ a_TP≤.
  This knowledge is useful to exclude such systems from searches for potentially habitable planets, and allows for a quick dynamical classification of the HZ (resonant vs non-resonant).
  
  
  
  
  
    §.§.§ Searching for CAM solutions
   
  
  Our approach is to find two sets of solutions of eq:freqeqmain: in the first case for a_TP =, and in the second case by setting a_TP =.
  The two solutions obtained in this way mark the extremal cases where the perturbation just touches the HZ borders.
  For parameter values intermediate to the two identified sets the perturbation will fall somewhere into the HZ.
  However, one must be aware that there are circumstances for which these solutions are still too conservative.
  An example is that the SR falls just slightly outside of the HZ, such that our model would predict that it does not have any influence.
  Resonances have a certain resonance width in real systems though, this width depends on the strength of the perturbation <cit.>.
  So, if the SR were located slightly outside the HZ and had a resonance width that is sufficiently large then it would reach into the HZ and does affect the dynamics there.
  Such situations must be carefully investigated on an individual basis; often numerical simulations will be necessary to decide whether or not a resonance near the HZ has an influence on it <cit.>.
  
  To simplify the search for solutions to the CAM in the five-dimensional parameter space, we fix 4 of the 5 variables and consider eq:freqeqmain to be a one-dimensional function of the remaining variable.
  We then solve the equation via an iterative process, because it depends in a non-trivial way on the variables α_1 and α_2 (via the Laplace coefficients).
  This procedure yields parametric solutions for the selected variable as a function of the other four parameters, and the solutions can be visualized in a two-dimensional parametric space like shown in Fig. <ref> and others in section <ref>.
  
  
  
  
  
    §.§.§ A motivating example
   
  
  
  
  
  We present 1as motivation a first application of the CAM in Fig. <ref>.
  In this example a binary star system is depicted that consists of a G-type primary and an M-type secondary star.
  The main feature is a diagonal stripe (bottom left to top right) that represents the set of parameters in the (a_B, a_GP) parameter space that lead to an SR in the HZ.
  1Three specific examples serve to highlight the strength of this approach.
  In the first case (medium grey or orange colour) a Jupiter mass planet at 5.2 au (1Jupiter's actual semi-major axis in the solar system) would perturb the HZ if the secondary was at a distance between 60-80 au (see filled rectangle in the middle).
  The second case (light grey or yellow colour) places a Jupiter mass planet at Saturn's position at 9.6 au, where this planet induces an SR in the HZ for stellar distances of roughly 170-220 au (lighter filled rectangle to the right).
  1
  Another case shows the same Jupiter-mass planet at a distance of only 3 au (dark grey or purple colour) with perturbations occurring in the HZ for stellar distances between 22 - 36 au (darker filled rectangle to the left).
  
  Note that the band is 1artificially truncated on both ends; for a_B < 20 au, because such a tight binary system would require to place the GP below 2 au 1to maintain its orbital stability for higher secondary eccentricities, but in turn this limits orbital stability in parts of the HZ.
  In practice, the stripe could be traced out to much larger secondary star (and planet) distances, but the more distant the giant planet becomes the weaker its influence on the HZ would be.
  
  
  
  
  
    §.§.§ Effect of a secular perturbation in the HZ
   
  
  
  
  
  Apart from being able to locate and identify secular perturbations some important questions remain:
  How much does it matter? Would potential planets in the HZ even care about such perturbations?
  The answer is likely yes!
  To understand the reason, we need to consider the driver of habitable conditions: the stellar irradiation.
  
  It became evident from the investigations of <cit.>, <cit.>, <cit.>, and others that a static definition of habitable zones in binary star systems is only a rough approximation, and for some applications it is unjustified.
  The gravitational interaction of the (at least) three bodies results in quasi-periodic oscillations of the planet's eccentricity even for initially circular planetary orbits.
  These oscillations in eccentricity are correlated with phases of larger amplitude variations in irradiation.
  Consequently, a planet that may sustain habitable conditions on a low eccentricity orbit could leave the permanently habitable zone (PHZ) for a temporarily increased eccentricity – just to return to a moderate eccentricity some time later.
  It is this effect on the eccentricity that effectively limits the habitability.
  
  In Fig. <ref> we provide the results of calculations for (massless) test planets distributed across the HZ for the 1three cases from Fig. <ref>.
  1These fully numerical simulations were performed with the N-body code Mercury <cit.> using the adaptive time-step Radau method with local error tolerance 10^-12.
  Each figure shows the maximum eccentricities (max-e) of test planets as a function of their initial location.
  The max-e is color coded according to four different scenarios related to the giant planet's orbital eccentricity (from bottom to top): (1) e_GP = 0.01 (square symbols), (2) e_GP = 0.05 (circles), (3) e_GP = 0.1 (triangles), and (4) e_GP = 0.2 (diamond symbols).
  Panel (a) for case 1 refers to the giant planet at 5.2 au, in panel (b) the planet is located at 9.6 au, while in panel (c) the planet is at 3.0 au.
  A triangle-shape overlay (in light gray) delineates the region of the PHZ for the test planets.
  Any test planet that acquires a max-e value outside of this region must be considered uninhabitable (at least temporarily).
  
  In the first example case the perturbing giant planet causes a linear SR near 1.1 - 1.2 au that is marked by the spikes in max-e.
  Even for e_GP = 0.01 the GP forces test planets locally to max-e values of 0.2.
  An increase of e_GP to 0.05 then boosts the eccentricities to values beyond 0.6.
  In these latter cases some small intervals with moderate max-e remain, e.g. from 1.2 - 1.55 au (1.3 - 1.4 au) for e_GP = 0.05 (0.1), respectively.
  For even higher e_GP no test planet can maintain max-e values low enough to stay in the PHZ; some even escape from the system on unbound trajectories.
  
  The second example with a more distant GP in panel (b) exhibits a linear SR at about the same location as before.
  The spike in max-e is not as pronounced as in the other example, but it is still visible even for relatively small e_GP.
  An increase in e_GP first leads to the deterioration of habitability near the edges of the HZ, but increasingly also in the inner regions.
  For e_GP = 0.1 already many test planets achieve too high max-e values and only in an interval from 1.25 - 1.5 au they would remain permanently habitable; for an even more eccentric GP no habitable test planet remain any more.
  
  1
  Case 3 in panel (c) represents a strongly perturbed HZ by a GP at a distance of only 3 au.
  Besides the main linear SR at about 1.1 au also numerous narrow peaks appear that correspond to MMR with the GP, see Table <ref> for a list of MMR.
  The width and max-e amplitude of the SR are strongly modified when increasing the GP's eccentricity, such that for e_GP > 0.1 no habitable test planets remain.
  Note that apparently the 5:1 MMR has the largest effect (max-e > 0.9 already for e_GP≥ 0.05), while the theoretically stronger 3:1 and 5:2 MMR are visible only for e_GP = 0.1 or higher.
  Although the 5:1 MMR is a fourth order MMR, here it interacts with the SR and becomes the dominant MMR locally, apart from the 2:1 MMR that is located slightly beyond the right edge of this figure.
  
  
  
  
  
  These examples show clearly that the GP's eccentricity plays a crucial role for secular perturbations in the HZ.
  1It is obvious that slight variations in e_GP can result in large perturbations that involve the large forced max-e of the TP induced by low to moderately eccentric GP, and the feedback effects on e_GP discussed in section <ref>.
  
  
  
  
  1
  Figure <ref> depicts a binary star system with a K-type secondary (m_B = 0.7   M_⊙) and the bands of various SR.
  Besides the linear SR (middle curve) the figure also shows the (q, p) = (1, 3) (top curve) and (q, p) = (3, 1) (bottom curve) nonlinear SR.
  Note that the vertical width of the curves in a_GP does not indicate their strength (nonlinear SR are rather weak, cf. Fig. <ref>), but rather the feasible range of GP distances that induce these SR in the HZ.
  At small secondary star distances these three curves overlap which means that two (or more) SR may appear simultaneously in the HZ.
  
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  § RESULTS
   
  
  1In a first part (section <ref>) we assess the analytical model's accuracy (i.e. the predicted secular frequencies g_GP and g_TP) relative to numerical reference results, while in section <ref> we apply the CAM to different binary systems and present the effect of varying different system parameters.
  
  
  
  
  
   §.§ Model accuracy
   
  
  1
  As the full CAM consists of two models we check the performance of each of them independently.
  Section <ref> contains the checks related to the CAM and g_GP, while section <ref> contains results pertaining to the g_TP part.
  Additionally, we compare the results of the full CAM (eq:freqeqmain) to those of the simplified CAM (eq:freqeqapprox) in section <ref>.
  
  
  
  
  
  
    §.§.§ Accuracy of full CAM: GP frequency
   
  
  We performed a series of test computations on a two-dimensional grid in the (a_B, e_B) parameter space for four different masses of the secondary star.
  The main goal was to assess the accuracy of the analytical model in eq:andeggmain when its assumptions of zero mass and eccentricity for the GP are not met.
  In fact, we tested several different GP masses and non-zero eccentricities for a total of more than 14,000 test runs.
  1
  The simulations were done in the framework of the full three-body problem, purely Newtonian gravity, and the equations of motion were solved with the Radau method from the Mercury package <cit.>.
  The secular peridos vary over several orders of magnitude (roughly 10^3 - 10^7 years) in the parameter range of interest, so we made sure to cover a simulation time of at least 20 secular periods (estimated from eq:heppenh) and to obtain at least 2000 data points for the time evolution of the heliocentric orbital elements.
  
  Each of these individual test runs was subject to a frequency analysis to obtain a reference value for g_GP.
  1The reference frequency was obtained by the FMFT method <cit.> from the time series of the complex variable z(t) = h(t) + i k(t), where i = √(-1), with h(t) = e(t) sinϖ(t) and k(t) = e(t) cosϖ(t).
  
  
  
  
  Figure <ref> displays a comparison of the CAM and the numerical reference values.
  The contours show the relative error in percent for the case of an M-type secondary, where the giant planet is located at 3 au on an initially circular orbit.
  In this case the agreement is very good, with a maximum deviation of less than 8 %.
  Generally, we observe the largest errors in the top-left corner for tight and highly eccentric secondaries, where orbital instability of the GP becomes an issue (cf. hatched area).
  When increasing the giant planet's mass to 5 M_J also the maximum error increases to ∼ 30 %, while the average error stays below 10 %.
  
  
  
  
  We show a more global comparison in Fig. <ref>, that gives an overview on the statistical distributions of the relative errors for the four secondary star masses.
  For this figure we selected the subset of the simulations with m_GP = 1   M_J and distances a_GP of 3, 5, and 7 au.
  The vertical axis shows that the relative errors are usually quite small with only few exceptions.
  For each of the four cases a boxplot indicates the range of errors in between the first and third quartiles, horizontal lines represent the median values that lie around 2.5 % in all cases.
  Furthermore, the region enclosed by the errorbars contains 95 % of the data points, while black dots mark individual outliers.
  The normalized frequency distributions of data points are shown as shaded regions.
  It is thus clearly visible that the distributions are all unimodal and that in > 95 % of cases the analytical model deviates by less than 10 % from the numerical reference value.
  The remaining outliers stem mainly from those simulations that combine a highly eccentric and tight secondary with a more distant giant planet.
  These edge cases lie close to the instability border (see top left corner in Fig. <ref>) and are characterized by large values of α_2.
  
  
  
  
  
    §.§.§ Accuracy of full CAM: TP frequency
   
  
  1
  A second set of test simulations covers eq:laplag.
  The aim is to assess how well this equation performs relative to numerical reference results in the presence of two massive perturbers.
  All simulations are performed in a restricted four-body problem with the Radau method (see section <ref> for details) and cover a duration of 2 Myr, with output every 100 years in order to resolve short period variations in eccentricity.
  We obtain the reference frequencies again with the FMFT method.
  
  
  
  
  
  1
  Figure <ref> presents the comparisons for two sample cases.
  Panel (a) shows case 1 from Figs. <ref> and <ref>, where we have a_GP = 5.2 au, e_GP = 0.05, and e_B = 0.01.
  The solid grey curve gives the secular frequencies as predicted by eq:laplag (left y-axis).
  Black dots represent numerically determined apsidal precession frequencies for test planets in the HZ range (0.95 - 1.70 au).
  We overplot the max-e curve (with light grey squares, see values on the right y-axis) to aid the identification of relevant features.
  Horizontal dashed lines indicate the GP's proper precession frequency (g_GP) and multiples of its value, whereas vertical dashed lines indicate MMR (see respective labels at the top).
  It is visible that there is an excellent agreement between the analytical and numerical frequencies, except for the interval between the 11:1 and 10:1 MMR.
  The nominal linear SR location is at 1.12 au where the grey curve equals g_GP; the nonlinear SR g_TP - 2   g_GP≈ 0 lies near to 1.65 au.
  According to the numerical data it seems that the exact SR location is slightly shifted to smaller distance (≲ 1.1 au).
  Note that the peak in max-e does not coincide with the location of the SR, moreover it is asymmetric with a sharp edge short of the 9:1 MMR.
  A similar qualitative behaviour has also been demonstrated for secular resonances in the solar system <cit.>.
  The max-e curve is more symmetric for low GP eccentricity (cf. Fig. <ref>), but the asymmetry becomes more pronounced with increasing e_GP.
  
  
  1
  The second panel (b) depicts case 3, with the same parameters as above except for a_GP = 3 au.
  Again the light grey points represent the max-e values for individual test particles, where we now have several narrow peaks due to MMR besides the dominating broad SR.
  The analytical and numerical frequencies are in rather good agreement, except at locations of MMR and starting from about 1.5 au outwards to larger distances.
  A main cause of this disagreement is the influence of the strong 2:1 MMR at 1.89 au (outside of the plot range).
  The analytical model in eq:laplag only accounts for secular terms of the disturbing function, but here also the resonant (short period) terms play an increasingly important role.
  An estimate of the libration width of the 2:1 MMR gives δ a ≈ 0.06 au <cit.>, but this is not wide enough to reach into the HZ, as we can also see from the flat part of the max-e curve at distances larger than 1.5 au.
  The width of the linear SR (at 1.1 au) roughly corresponds to the interval between the 5:1 and 17:4 MMR, compare Tab. <ref> for their exact locations.
  In this region the test planet's proper precession frequency does not change.
  A similar behaviour occurs near the 3:1 MMR that falls into the vicinity of the g_TP - 2 g_GP nonlinear SR.
  Such flat parts (or plateaus) of the frequency curve correspond to crossing of resonance islands, as demonstrated for different MMR in <cit.>.
  Another interesting case is the 3:1 MMR, where the max-e curve shows a `pulse-like' effect.
  <cit.> found a similar behaviour when they analysed the effect of near MMR on the secular precession frequency.
  An explanation for the pulse is that the secular frequency profile is modified by the sign-change of the MMR's small divisor, which – in case of second-order MMR – gives a negative contribution before the resonance and a positive one afterwards.
  This effect on the secular frequency directly relates to the max-e curve where it is visible, too.
  
  
  1
  In general, missing data points in Fig. <ref> indicate a complex dynamical behaviour of the respective particle, where the FMFT was unable to determine a unique dominating frequency.
  Note that in a purely secular setting the d'Alembert rules require additional frequencies (like nodal precession frequencies) to satisfy the resonance relation in eq:freqequal for a nonlinear SR.
  Alternatively, a nonlinear SR may couple with a nearby MMR to fulfill the exact resonance condition, such as for the (p, q) = (2, 1) SR and the 3:1 MMR.
  This coupling gives rise to the max-e peaks visible in Fig. <ref> (c) for e_GP≥ 0.1, but also a local minimum in max-e is associated to the same 3:1 MMR.
  The presence of SR and nearby MMR does not necessarily lead to adverse effects, because MMR can also have a stabilizing effect.
  One particular example was shown in <cit.> for the γ Cephei system, where a stable island around the 3:1 MMR (near to an SR) persists also for non-coplanar motion up to an inclination of 15^∘ relative to the GP.
  
  
  1
  Overall, Figs. <ref> and <ref> clearly show that the linear SR is the dominant perturbation and that higher order SR play a secondary role.
  This finding is backed-up by the results of <cit.> whose max-e plots also show a dominant linear SR that locally couples to various MMR.
  Nonlinear SR become relevant for much higher e_GP when the resonance widths of SR and MMR are wide enough to overlap and cause extended regions of chaotic motion.
  Such an example was presented in <cit.> for the HD 196885 system[Radial velocity observations of the system suggest an exoplanet with minimum mass 3   M_J and eccentricity e_GP∼ 0.5.] where the nonlinear SR 2   g_TP - g_GP and 3   g_TP - g_GP overlap with nearby MMR.
  For another example of nonlinear SR see Fig. <ref>.
  
  
  
  
  
  
    §.§.§ Accuracy of simplified CAM
   
  
  Figures <ref> and <ref> demonstrate that the full CAM provides an accurate description of the secular dynamics.
  The open question is how well the simplified CAM reproduces results of the full equation.
  
  
  
  
  A first test is shown in Fig. <ref> (a).
  The plot uses eq:gamma1sol to express the dependence of the mass ratio γ_1 on the secondary star distance (via α_2).
  It shows two shaded regions that overlap each other partially.
  The darker shade represents the solutions according to the simplified CAM, while the lighter shade are the results from the full CAM.
  Apparently the simplified CAM shows the same qualitative behaviour and follows the full solutions closely, but with 1an offset in vertical direction (note the logarithmic scale of the y-axis).
  This deviation does not show any dependence on the secondary star distance (i.e. α_2) or on the secondary star's eccentricity (not shown).
  We have to bear in mind that eq:gamma1sol involves 1strong approximations with respect to the distance ratios α_1 and α_2.
  Accordingly, the major effect causing the visible deviation is the magnitude of α_1 (fixing a_GP = 5.2 au).
  1The deviation becomes larger with increasing distance ratio α_1, i.e. when the GP is near to the HZ and its perturbations are stronger.
  
  Figure <ref> (b) shows another comparison of the accuracy.
  In this case we use eq:alpha2sol0 to trace the distance ratio α_2 as a function of secondary star mass in the range 0.08 ≤ m_B≤ 2.1   M_⊙.
  We set the giant planet to 5 au and always assign an eccentricity of e_B = 0.25 to the secondary.
  The results demonstrate that the simplified CAM is able to reliably describe the dynamics, albeit with a small deviation.
  It is visible from the figure that the simplified CAM retrieves a narrower region in the parameter space than the full CAM, but the deviation is only minor especially at the lower border.
  Similar to the first test case the parameter α_1 plays a major role for the magnitude of the deviation: the larger α_1 the larger is the deviation.
  
  We can conclude that the simplified CAM is able to 1reliably describe the dynamics, provided that the giant planet does not approach the HZ too closely.
  1Its limitations are clearly connected to the parameter α_1 = a_TP / a_GP that should be α_1 < 1/3 (in particular < 0.1) in order to avoid the effect of MMR that would invalidate a purely secular model.
  The approximate solutions can thus serve to quickly estimate the variation in 1any of the five system parameters, and to scan a large parameter space for a first overview 1on possible secular perturbations in the HZ.
  
  
  
  
  
   §.§ Application to exoplanet systems of binary stars
   
  
  Before continuing to present the model outcomes applied to synthetic systems, we need to be more explicit about the choice of system parameters.
  As mentioned above, we consider only perturbers external to the HZ 1of a G-type main sequence star.
  This puts a constraint on the location of the giant planet 1to orbit beyond about 2 - 3 au depending on its eccentricity.
  Thus, in a first step, we need to get an overview of suitable exoplanet candidates.
  
  
  
  
  
    §.§.§ Selection process
   
  
  
  
  
  We selected from 1various online catalogs[
  For instance, The Extrasolar Planets Encyclopaedia of <cit.> at <https://exoplanet.eu/catalog/>, or the NASA Exoplanet Archive <https://exoplanetarchive.ipac.caltech.edu/index.html>.] all confirmed exoplanets that fulfill the following two criteria: (1) masses between 0.1-10 Jupiter masses (M_J), and (2) semi-major axes between 2 and roughly 10 au, irrespective of planet or stellar multiplicity.
  This gives a total of 220 exoplanets in our sample.
  The choice of parameters is motivated first by the evident incompleteness of planet detections at larger distances, and second by the survey of <cit.>.
  From their combined radial velocity – imaging survey those authors concluded that the frequency of giant planets peaks between 3 and 10 au for the mass range up to 20 M_J.
  
  Figure <ref> shows the scatter plot and respective parameter distributions of these exoplanets.
  The histogram on the right shows that the median (minimum) mass is around 2.5 M_J, with an extended tail reaching up to the selection limit.
  The distribution of semi-major axes in the bottom histogram is more concentrated to lower values, with about 90 % of the sample having values of a_GP≤ 6 au.
  This is a consequence of the detection by mainly radial velocity surveys that favour exoplanets with shorter orbital periods <cit.>.
  
  This selection of giant planet properties has implications on the secondary star parameters.
  We restrict ourselves to masses in the interval 0.4 ≤ m_B≤ 1.3 M_⊙ (roughly spectral types M2V – F5V), and a lower limit on the stellar separation of a_B≥ 20 au.
  The eccentricities can vary between 0 ≤ e_B≤ 0.6 according to the fit limitation in the model of <cit.>.
  
  
  
  
  
    §.§.§ Binary star parameters
   
  
  As we can appreciate from eq:freqeqmain the right-hand side has a non-linear dependence on γ_2, α_2, and  e_B.
  This means that the location of the SR will strongly depend on the secondary star's parameters.
  We will first investigate the qualitative behaviour of the solutions.
  
  In a former study, <cit.> have investigated the system HD 41004 AB as a proxy for binary stars with circumstellar exoplanets.
  They have demonstrated that the positions of SR are connected to the orbital parameters of the secondary star, especially to its distance.
  Accordingly, they have investigated the qualitative behavior of SR location for a shift in a_B (for a_GP = const.).
  Here, our purpose is different: We constrain the SR to lie in the HZ and show how the system parameters (masses and distances of GP and secondary) need be coupled in order to maintain that situation.
  
  From Fig. <ref> we can observe the correlation between a_B and a_GP that results in the diagonal stripe.
  This means that for each value of a_B (for fixed a_TP) there is a unique solution to eq:freqequal that yields a distance a_GP causing an SR in the HZ.
  The figure contains another important detail: When increasing the secondary's distance from the host star (for fixed GP distance along the horizontal lines, from left to right) we can observe that the SR moves from outside in.
  This means that first the outer part of the HZ (upper edge of the stripe) would be affected, then some intermediate region of the HZ, and finally the inner border (lower edge).
  Increasing the distance even more would then remove the SR from the HZ and would place it closer to the host star.
  Such an inverse relation between a_B and the SR location is not immediately obvious from eq:freqeqmain, but we can understand it in terms of α_2.
  The variable α_2 is simply the ratio a_GP / a_B, hence it is inversely proportional to the secondary's distance.
  In case that the other parameters are held fixed (as in Fig. <ref>), an increase of a_B leads to a decrease of g_GP because of the term ∝α_2^3, cf. eq:heppenh.
  Conversely, the frequency g_TP also has to decrease to satisfy eq:freqeqmain.
  The balance is restored through the term ∝γ_1α_1^2 (that is a monotonic function and generally dominates over the second term) via the decrease of α_1, i.e. the resonance location moves away from the giant planet and closer to the host star.
  A physical interpretation of this behaviour is that a more distant perturbing secondary star reduces the rate of angular momentum exchange with the giant planet – its secular frequency g_GP decreases.
  Similarly, test planets have lower orbital precession frequencies the further away they are from the GP (which is their main perturber), such that necessarily the resonance location moves away from the GP.
  
  At the left edge of the Fig.<ref>, the red line denotes the critical distance, , of the GP according to the fits of <cit.> and <cit.>.
  This line separates the dynamically unstable area (hatched, to the left) from the region of orbital stability.
  The indicated example with e_B = 0 is an optimistic best-case scenario, typically the secondary star has a non-zero eccentricity.
  Furthermore, the results of <cit.> suggest that the hatched region would also be more extended if we accounted for the GP's elliptic motion.
  
  One of the main observables of binary star systems is the mass of the secondary star.
  This parameter is of utmost importance to dynamical studies, because it contributes significantly to the size of the host star's region of orbital stability and determines the strength of gravitational perturbations.
  In the following example we show the variation of the SR's location with the secondary star's mass.
  
  
  
  
  Figure <ref> presents two binary stars in the same parameter space like before.
  The G+M (upper curve; m_B = 0.4   M_⊙) and G+F (lower curve; m_B = 1.3   M_⊙) systems relate to edge cases with lowest and highest secondary star masses, respectively.
  All other mass combinations (e.g. G+G, G+K) fall in between these two curves.
  In this example there is a giant planet with one Jupiter mass at 10 au (dashed horizontal line).
  This planet induces an SR in the HZ for secondary star distances of 180 - 240 au in the G+M case, and 260 - 350 au in the G+F case.
  There is a clear trend with secondary mass: the curve becomes flatter with increasing mass, and the absolute width (in terms of a_B) of the SR generating region broadens.
  
  We can thus derive a general picture from Figs. <ref> and <ref>: (1) there is a positive correlation between giant planet and secondary star distance for fixed secondary star mass; and (2) the secondary's mass directly relates to the interval of stellar distances generating SR for fixed giant planet distance.
  This means that a higher mass secondary at a distance of several hundreds of astronomical units might still perturb the GP in such a way as to induce an SR in the HZ.
  
  
  
  
  
    §.§.§ Binary star and planet parameters
   
  
  Until now some planet-related parameters were held constant or ignored altogether.
  Next we will consider also the variation of these planet parameters and investigate their influence on the results.
  
  
  
  
  We have already demonstrated in Fig. <ref> that there is a correlation of stellar mass and resonance location.
  A similar correlation also exists for the planet mass, although with the inverse effect as Fig. <ref> shows.
  This figure displays the same parameter space (a_B, a_GP) like the ones before.
  In this example we have a couple of G-type stars (each with 1   M_⊙) with a planet located at 10 au whose mass is a free parameter.
  The top curve is for a planet with 3 times the mass of Jupiter (M_J), the middle curve for 1 M_J, and the bottom curve corresponds to a Saturn-mass planet.
  As the GP's mass increases (from bottom to top) the curve becomes steeper and the resonance location shifts to the left towards smaller secondary star distances.
  The consequences of this behaviour for different planet masses are twofold: (1) a negative feedback for more massive planets, and (2) a feedback through coupling of planetary mass and eccentricity.
  To clarify both issues we will compare two systems for equal planetary distances but with different masses.
  
  It is intuitively clear that a more massive planet is a more efficient perturber of the HZ area than its lower mass counterpart.
  However, at the same time, a higher mass planet also requires smaller stellar distances 1to force the SR into the HZ, as it is visible from the figure.
  These two facts entail stronger perturbations in the HZ due to larger planet mass and closer secondary star.
  Hence, the negative feedback consists in an increased perturbation strength for a higher mass planet in combination with a closer secondary star.
  
  The second feedback effect concerns the opposite effect.
  With decreasing planet mass the perturbing secondary star has to shift outwards to a distance of several hundred au to place an SR into the HZ.
  Such extended systems have a weaker dynamical coupling than more compact ones, hence one would not expect strong perturbations in the HZ.
  However, the giant planet's eccentricity plays a crucial role and compensates for the weaker coupling to the secondary star, see Fig. <ref> (b).
  In a wide binary system the second star is more susceptible to external perturbations e.g. from passing stars (see section <ref>).
  Such perturbations can act to increase the secondary star's eccentricity.
  The giant planet's forced eccentricity is approximately ∝ e_B / (1 - e_B^2) <cit.>, so it increases with larger e_B.
  Additionally, in former studies we have found that observed exoplanets have higher eccentricities in general than what their forced eccentricity would suggest <cit.>.
  Although the secondary star must be farther away for lower mass planets it can effectively induce a higher planetary eccentricity.
  In this way the increased eccentricity can balance the magnitude of perturbations that otherwise would decrease with smaller planetary mass.
  
  
  
  
  The eccentricity does not only play a role for the planet, but – more importantly – also for the stellar orbit.
  Figure <ref> shows a different cut through the (e_B, a_B) parameter space.
  Three cases present the resonance location as a function of the secondary's eccentricity and the planet's location (here the GP has one Jupiter mass).
  When placing the GP to 3, 5, or 7 au (from bottom to top), respectively, the SR is induced by secondaries in a wide range of distances (y-axis).
  The SR location also depends on the secondary star eccentricity: as e_B grows the SR moves outwards with respect to the circular case.
  We can identify the cause for this variation with the terms ∝ (1 - e_B^2)^-3/2 in eq. (<ref>).
  These two terms are responsible for the weak deviation from the circular value at low e_B, and the increasingly larger difference at high e_B.
  
  A striking characteristic of this figure pertains to the large shifts in a_B for a minor change of a_GP.
  For instance, the SR is displaced from about 50 au (for a_GP = 3 au) to roughly 200 au (for the GP at 7 au).
  These results indicate potentially dramatic changes of the system's secular architecture by slight shifts in a_GP (see section <ref>).
  
  
  
  
  
   §.§ Resonance generating distance of giant planet
   
  
  In Figures <ref>, <ref>, and <ref> we can observe that the distances of GP and secondary star show a positive correlation, but also that this relation is basically nonlinear.
  It is very important to understand the occurrence of SR in the HZ as a function of the GP's semi-major axis, a_GP, for actual applications.
  Hence, we need to understand the relation between a_GP and a_B as a function of the other parameters.
  For this we cannot use the solutions presented in Appendix <ref>, because they express α_2 in terms of α_1 which itself contains a_GP, so those equations only give implicit relations.
  Conversely, in the following we will derive an explicit relation on the basis of eq:freqeqapprox.
  
  We start from the simplified CAM in the limit δ→ 0 to get rid of the additional dependence on α_2.
  Recall the definitions α_1 = a_TP / a_GP and α_2 = a_GP / a_B, and that the constant a_TP is either  or  for the inner or outer HZ border, respectively.
  After inserting the corresponding variables into eq:freqeqapprox, 1we obtain a cubic equation for z = a_GP^3/2.
  The general solution is thus given by
  
      α_2 = 3^-2/3( a_TP/a_B) ( (u + u^-1 - 1)^2/1 + γ_1)^1/3,
  
  with the abbreviations
  
      u = ( -2 c_2^3 - 27 c_0 c_3^2 + 3 √(3)√(4 c_0 c_2^3 c_3^2 + 27 c_0^2 c_3^4))^1/3 / ( 2^1/3 c_2),
  
  and c_0 = - γ_1 a_B^3, c_2 = - γ_2ϵ_B, c_3 = (1 + γ_1)^1/2γ_2ϵ_B a_TP^-3/2.
  
  An important property of solution (<ref>) is that the expression γ_2ϵ_B appears only in this form, i.e. these parameters are never separated from each other.
  Furthermore, we could expand the denominator into a power series (1 + γ_1)^-1/3≈ 1 - γ_1 / 3 + O(γ_1^2) since γ_1 = m_GP / m_A≪ 1.
  Note that the ratio a_TP / a_B also plays a role in u via c_0 c_3.
  
  Equation (<ref>) is useful to calculate the orbital distances of GP that would generate secular perturbations in the HZ for any given binary star system.
  Once the corresponding value of α_2 is known, we can determine the GP location from a_GP = α_2 a_B and construct parametric plots like in Fig. <ref>.
  
  Tables <ref> – <ref> in Appendix <ref> serve to assess the accuracy of eq:agpsol.
  1
  Each table lists the critical values α_2,IHZ and α_2,OHZ (for the inner and outer HZ border, respectively) for a given binary star system.
  These solutions of eq:freqeqmain determine the location of the GP, a_GP, as a function of the secondary star's semi-major axis, a_B, such that an SR affects some part of the HZ.
  The main parameters are the GP's mass (in units of Jupiter's mass), and the secondary star's eccentricity e_B.
  
  The general trends of (1) decreasing α_2 with increasing a_B (i.e. the negative curvature, also visible in Fig. <ref>), and (2) the slight shift in α_2 with increasing e_B, like in Fig. <ref>, are clearly discernible from the data, too.
  
  
  
  
  The numbers in the tables have been used as reference values for comparison with predictions  by employing the model in eq:agpsol.
  Figure <ref> shows the absolute deviations | Δα_2 | between the model and the reference values from the tables.
  For each distance bin in a_B we computed the average error and its standard deviation (1-σ) taken over all different combinations of m_GP and e_B.
  The figure shows that for a_B≥ 100 au the error is less than 10^-3, i.e. the location of the GP that generates an SR can be determined with a typical accuracy of < 0.1 au.
  With smaller secondary star distance the errors are larger though, because this is a direct consequence of neglecting higher order terms in the expansion of the Laplace coefficients when deriving the simplified CAM.
  A closer look at the figures also reveals the trend that the errors grow steadily from G+F to G+M stars (left to right).
  This behaviour can be understood  by the dependence of the SR region on the secondary's mass as mentioned in connection with Fig. <ref>: a more distant secondary must be more massive to generate the same effect as a more proximal lower-mass star.
  Thus, M-type companion stars are naturally closer to the GP, and so the ratio α_2 is larger, which leads to somewhat higher deviations when using eq:agpsol.
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  § DISCUSSION
   
  
  Results from earlier investigations of binary star systems assumed that well-separated binaries (beyond roughly 200 au) can be simply treated as single star systems <cit.>.
  The new results presented in the previous section shed a different light on those assumptions, they demonstrate that secular perturbations do affect the HZ for a wide range of binary system orbital and physical parameters.
  In the following sections we will discuss some of the implications of those resonances on the dynamical evolution of test planets in the HZ.
  
  
  
  
  
   §.§ Effect of stellar fly-bys
   
  
  Passing stars can affect planetary systems severely and can lead to direct or indirect ejections <cit.>.
  A subtle effect relates to the delayed instability of planetary systems of wide binary stars induced by external perturbations, as discussed by <cit.>.
  They explain the apparent overabundance of eccentric (circumstellar) exoplanets in wide binary systems by the action of galactic tides that lead to an increase of the secondary star's eccentricity.
  At a certain point its eccentricity is large enough for the periastron distance to drop below a critical value, below which the star triggers the instability of the exoplanetary system.
  
  <cit.> reiterate on this topic to disentangle the individual contributions of two types of external perturbations: (1) the galactic potential, and (2) stellar encounters via fly-bys.
  They find that stellar passages have a destabilizing effect, because 40 % of their test binary systems are disrupted, and this value is (nearly) independent of the systems's initial eccentricity <cit.>.
  This means that a stellar fly-by, by itself, is effective in modifying the secondary's eccentricity.
  
  These results are verified by <cit.>, who study the dynamical effects of a sequence of stellar encounters between a binary star system and a passing star.
  Their results show that the stellar eccentricity has a maximum variance that reaches up to 0.45 for an initial stellar separation of ≥ 150 au.
  Hence, stellar fly-bys over an extended period of time can significantly modify the companion star's orbital eccentricity.
  
  Returning to Fig. <ref>, we can put these findings about stellar fly-bys into context with the results shown in that figure.
  Let us assume a binary system with an architecture such that initially the HZ is free of SR, i.e. the system would be located somewhere in the white area above one of the coloured zones.
  The combined effect of stellar fly-bys then can lead to a significant change (increase) of the secondary star's eccentricity, perhaps in such a way that it enters one of the coloured bands.
  Consequently, an exoplanetary system with a dynamically calm HZ can transition into one that is affected by strong secular perturbations.
  This could turn a planet that has always had habitable conditions into a potentially uninhabitabe one, depending on how much its eccentricity increases by the perturbations.
  The next section considers how to exactly quantify the effect of perturbations.
  
  
  
  
  
   §.§ Related works
   
  
  The scope of the presented model is not limited to identify and trace SR throughout the HZ of binary star systems.
  1
  In an alternative embodiment of this method we could prescribe any interval of interest, not just the HZ.
  For instance, an application could be to check whether or not an SR occurs at or around the location of a detected or candidate exoplanet in a multi-planet system.
  
  This can be achieved by selecting for a_TP an interval of distances that is different from the HZ range ≤ a_TP≤.
  
  In the present study a major premise is that all bodies are coplanar or nearly so.
  As this is not generally the case, we have to consider the effect of mutual inclinations of the planets and the outer perturber.
  <cit.> have investigated the dynamics and stability of initially coplanar two-planet systems in tight binaries with an inclined stellar companion.
  They identified a narrow zone for the second planet that suppresses the Kozai-Lidov effect on an inner hot-Jupiter and helps to stabilize the system.
  <cit.> improve on these former results by studying basically the same situation, but they derive an analytic stability criterion that takes into account the perturber's separation and eccentricity.
  Their study might be viewed as an extension of the current work to more degrees of freedom by adding an inclination to the system.
  They include the additional perturbations due to the Kozai-Lidov mechanism, however they only consider planetary mass perturbers in lieu of stellar masses like in our work.
  1
  An extension of the current model could follow the lines of <cit.> who studied SR for prograde and retrograde orbits.
  They constructed a semianalytical model for the restricted four-body problem that allows to describe both types of motion with a unified formalism following a suitable change of reference system.
  In this model resonant angles consist of linear combinations of angular differences of apsidal and nodal precession frequencies.
  Although the dynamical model of <cit.> is similar to the one used in this paper (see section <ref>), they studied planetary satellites in very close orbits whereas we focus on rather hierarchical systems with period ratios P_out/P_in > 2.
  
  
  Furthermore, in analogy to the current work, <cit.> have studied the influence of SR on a debris disc for the case of interior planetary perturbers.
  In particular, they mapped how the locations of SR depend on the masses, semi-major axes, and eccentricities of the planets.
  
  The work of <cit.> also contains certain aspects from the current study.
  When they investigate the stability of multiplanet systems in binaries by frequency map analysis, they identify cases when the two planets are not long-term stable although both of them are located well inside the critical stability border.
  They attribute this to the effect of the binary's eccentricity and establish an empirical relationship for (a_B, e_B) to predict a critical value a_B,crit below which the planets become unstable.
  For computational reasons they restricted their parameter range to 10 < a_B < 100 au and to only three different stellar mass ratios.
  The CAM can be regarded as a generalization of their empirical relation that takes into account a much larger parameter space and hence covers a broader variety of use-cases.
  Our CAM model could also offer an explanation for the bifurcated behaviour of eccentricity excitations of exoplanets under binary perturbations that <cit.> describe in their Fig. 14.
  They noted that Kozai-cycles can be either suppressed or fully active in the same parameter space region for a small change in the planet's parameters.
  
  In <cit.> the authors focus on the same topic of a two-planet system subject to the perturbations of an external perturber.
  Their analysis aims to derive analytical expressions for the eccentricities of the inner planetary system as a function of the perturber's orbital elements and mass.
  A main result is that a weak secular coupling between the planets (as is the case in hierarchical systems) leads to a higher susceptibility for eccentricity excitations, especially if the innermost planet is the least massive one, similar to what we have demonstrated in section <ref>.
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  § SUMMARY
   
  
  In this study, we investigate circumstellar planetary systems of binary stars 1(with coplanar and prograde motion) and the conditions that lead to secular perturbations in the habitable zone of the host star.
  Our main aim is to identify zones in the five-dimensional parameter space (orbital and physical parameters of the massive bodies, cf. Fig. <ref>) that generate 1dynamically unfavourable conditions for potentially habitable planets, 1such as elevated orbital eccentricities.
  This contribution helps to identify such systems in large-scale surveys and to exclude them from in-depth observations.
  
  We apply the restricted four-body problem as dynamical model and construct a Combined Analytical Model (CAM) to systematically map the properties of different synthetic systems (see section <ref>).
  Subsequently, we introduce additional simplifications to the CAM (section <ref>) and derive explicit analytical formulas to express the resonance locations in each system parameter as a function of the others (see Appendix <ref>).
  The CAM is assessed in comparison to fully numerical results which shows that it has a median relative error less than 3 % (cf. section <ref>).
  
  Our results are based on the properties of a sample of detected exoplanetary systems from different online databases.
  Secular perturbations in the HZ appear even for hierarchical wide binary systems, where stellar separations are on the order of several hundred astronomical units.
  We give details on the parameter combinations that generate SR in the HZ in the form of two-dimensional sections through the parameter space (section <ref>).
  Additionally, we derive analytical expressions to calculate the set of unfavourable giant planet positions for any given mass, semi-major axis, and eccentricity of the secondary star (in section <ref>).
  
  Based on these results we introduce the online tool SHaDoS (Appendix <ref>) that implements the CAM and provides a simple user interface to reproduce the results presented in this study.
  
  
   
  
  
  
  
  
  
  
  
  AB and EPL are grateful for the continuous support from Austrian Science Fund (FWF) project S11608-N16, part of the National Research Network “Pathways to habitability” (PatH).
  We wish to thank our colleagues D. Bancelin, N. Georgakarakos, and C. Lhotka for inspiring discussion in the course of preparing this manuscript.
  The authors thank an anonymous reviewer whose comments and suggestions helped to improve the manuscript.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § ONLINE TOOL
   
  
  We provide the online tool SHaDoS[This is an acronym for Secular perturbations in Habitable zones of Double Stars.
  The tool is accessible at <https://www.univie.ac.at/adg/shados/index.html>] and allows to reproduce the results presented in this article.
  It implements the CAM model introduced in section <ref> and solves eq:freqequal for any given set of parameters.
  In the following we describe the main features of this tool.
  
  The application implements an object-oriented approach with a step-by-step process for user input.
  There are four consecutive steps tied to the main objects of interest: (1) the host star, (2) the perturbing planet beyond the HZ, (3) the distant secondary star, and (4) the parameter space.
  In each of these substeps the application queries for some system parameter, refer to Fig. <ref> for an overview:
  
      
    * The relevant host star parameters are mass, and spectral type.
      Four different spectral types are preset, namely F, G, K, and M-type main-sequence stars with their respective stellar luminosity and effective temperature.
      Users can manually enter any other stellar characteristics, but effective temperatures are limited to the interval 2600 ≤ T_eff≤ 7200 K.
      The extent of the HZ is then computed based on this input.
  
      
    * Giant planet parameters include mass (in units of Jupiter's mass), distance to the host star (semi-major axis), and orbital eccentricity.
      The eccentricity is restricted to elliptic orbits, and the distance must be larger than the outer HZ border .
  
      
    * The input for the secondary star covers its mass, distance, and eccentricity.
      Eccentricity is restricted to the interval 0 ≤ e_B≤ 0.6 due to limitations in the analytical model (cf. section <ref>).
  
      
    * The parameter space consists of a two-dimensional grid, where a primary parameter must be selected as the independent variable, and a secondary parameter is defined as the sought dependent variable.
      Available parameters are any of those depicted in Fig. <ref> except for the already fixed host star mass m_A.
      Additionally, the minimum and maximum value of the independent variable must be entered, as well as the number of discretization steps.
  
  
  After having finished all of these steps the online tool performs the calculations on-the-fly and finally displays the results.
  There is an option to save the resulting plot as an image in the widely supported PNG format.
  Any comments, suggestions, or bug reports are welcome and should be addressed to the first author.
  
  
  
  
  
  § SOLUTIONS FOR SIMPLIFIED CAM
   
  
  Here we provide analytical solutions of the simplified CAM from eq:freqeqapprox.
  The main goal is to express each of the parameters (γ_1, γ_2, α_1, α_2, ϵ_B) in turn as a function of the other four variables.
  
  
  
  
  
   §.§ Solutions for γ_1
   
  
  There exist two real-valued solutions for γ_1,
  
      γ_1, ± = η[ -2 α_1^3± (1 - δ) ( η (1 - δ) + √(4 α_1^3 (1 - η) + η^2 (1 - δ)^2)) ] / 2 α_1^3,
  
  where η = γ_2α_2^3ϵ_B, and δ is given by eq:andeggcorr.
  The square-root is definitely positive if 1 - η≥ 0, i.e. if γ_2α_2^3ϵ_B≤ 1.
  This puts a coupled constraint on these three parameters, because it holds strictly that 0 ≤α_2 < 1, and ϵ_B≥ 1, while γ_2 can be either smaller or larger than 1.
  Moreover, we can discard γ_1,- < 0 as an unphysical solution, because the mass-ratio cannot be negative.
  The remaining γ_1,+ is the preferred solution to eq:freqeqapprox.
  
  
  
  
  
   §.§ Solution for α_1
   
  
  Apart from the trivial solution α_1 = 0, we find
  
      α_1 = (1 + γ_1)^1/3( η (1 - δ)/γ_1 + η)^2/3.
  
  Note that additionally there exists also a pair of complex conjugated solutions for α_1, but again these are not of physical interest.
  We caution that this functional form is only a first approximation to the `true' solution of the (simplified) CAM.
  The formula above will give satisfactory results only if α_1≪ 1, i.e. typically for α_1 < 0.1, because among all parameters α_1 is most sensitive to the assumptions made about the magnitudes of the distance ratios α_i (i = 1,2).
  In many cases α_1 will be larger than this threshold, compare the histogram of semi-major axes in Fig. <ref>.
  
  
  
  
  
   §.§ Solutions for γ_2
   
  
  While it is straightforward to obtain the solutions for γ_1 and α_1, there is a caveat regarding the other parameters.
  The problem is the explicit polynomial dependence of δ on the parameters γ_2, α_2, and e_B (cf. eq:andeggcorr), which leads to rather complicated equations.
  One way to circumvent this issue is to linearize the equation in γ_2 by the approximation δ→ 0.
  This is equivalent to a fall-back to the Heppenheimer formula eq:heppenh and hence gives an approximation of first-order in masses.
  The unique linearized solution for γ_2 reads
  
      γ_2^(0) = γ_1α_1^3/2/α_2^3ϵ_B( (1 + γ_1)^1/2 - α_1^3/2),
  
  where the superscript 0 denotes the approximative nature of this solution in the limit δ→ 0.
  This solution is defined for any combination of parameters, since (1 + γ_1)^1/2≥ 1 and 0 ≤α_1^3/2 < 1, so the denominator is always positive.
  
  When considering the full equation for γ_2 we only sketch the steps to find the solutions.
  From <cit.> we observe that the correction term is δ(γ_2) ∝ F(γ_2^1/2) + G(γ_2^1) + H(γ_2^2), where the functions F,G,H also depend on α_2 and e_B.
  Using the substitution z = γ_2^1/2 this results in a polynomial equation of degree 6 in z, for which the roots can be found with the Newton-Raphson method with initial guess γ_2^(0).
  
  
  
  
  
   §.§ Solution for α_2
   
  
  On inserting the secular frequency correction term δ from eq:andeggcorr into eq:freqeqapprox we face the problem that it would introduce fractional powers of α_2.
  Hence, we restrict ourselves to the same approximation as before,δ→ 0, such that no additional terms proportional to α_2 arise.
  The resulting equation possesses the solution
  
      α_2^(0) = ( γ_1α_1^3/2/γ_2ϵ_B ((1 + γ_1)^1/2 - α_1^3/2))^1/3.
  
  Note that it resembles eq:gamma2sol, but the roles of γ_2 and α_2 are exchanged.
  A better approximation is found if we kept in δ the term ∝α_2^3/2, which would give a cubic equation in z = α_2^3/2.
  
  
  
  
  
   §.§ Solutions for ϵ_B
   
  
  The function δ in eq:andeggcorr depends on the parameter e_B and involves its powers e_B^0, e_B^2, and e_B^4, but eq:freqeqapprox itself contains ϵ_B = (1 - e_B^2)^-3/2.
  Our goal is to achieve a simpler form of eq:freqeqapprox that should only depend on ϵ_B.
  In fact, the MacLaurin series expansion of ϵ_B with respect to e_B is
  
      (1 - e_B^2)^-3/2 = 1 + 3/2 e_B^2 + 15/4 e_B^4 + 𝒪(e_B^6)
  
  and this expansion recovers the correct powers of e_B in δ.
  There are other evidences that can justify the use of ϵ_B instead of e_B.
  First, the correction δ extends the <cit.> formula, which is a linear function of ϵ_B.
  Second, <cit.> and <cit.> constructed secular models for the three-body problem that extend Heppenheimer's model to higher order in the masses, and both new models contain some power of the expression (1 - e_B^2).
  Following these arguments we set δ = ϵ_B( α_2^3/2Γ_1 (γ_2) + α_2^9/2Γ_2 (γ_2) ), where functions Γ_1 and Γ_2 collect all terms that depend on the parameter γ_2.
  In this way δ becomes a linear function of ϵ_B.
  Then the full solutions are
  
      ϵ_B, ± =
          γ_2α_2^3Δ±√(
                  (γ_2α_2^3Δ)^2
              - 4 (1 + γ_1)^1/2γ_1γ_2α_1^3/2α_2^9/2Γ)/2 (1 + γ_1)^1/2γ_2α_2^9/2Γ,
  
  with the abbreviations Δ = (1 + γ_1)^1/2 - α_1^3/2 (Δ > 0) and Γ = Γ_1 + α_2^3Γ_2.
  Note that ϵ_B≥ 1 in order to be a physically relevant solution.
  
  
  
  
  
  § TABLES
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  apalike
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  