Binary Gravitational Perturbations and Their Influence on the Habitability of Circumstellar Planets

In order to assess the habitability of planets in binary star systems, not only astrophysical considerations regarding stellar and atmospheric conditions are needed, but orbital dynamics and the architecture of the system also play an important role. Due to the strong gravitational perturbations caused by the presence of the second star, the study of planetary orbits in double star systems requires special attention. In this context, we show the important role of the main gravitational perturbations (resonances) and review our recently developed methods which allow a quick determination of locations of secular resonances (SRs) in binary stars for circumstellar planetary motion where a giant planet has to move exterior to the habitable zone (HZ). These methods provide the basis for our online-tool ShaDoS which allows a quick check of circumstellar HZs regarding secular perturbations. It is important to know the locations of SRs since they can push a dynamically quiet HZ into a high-eccentricity state which will change the conditions for habitability significantly. Applications of SHaDoS to the wide binary star HD106515 AB and the tight system HD41004 AB reveal a quiet HZ for both systems. However, the study of these systems indicates only for the tight binary star a possible change of the HZ's dynamical state if the orbital parameters change due to new observational data.


INTRODUCTION
Since the discovery of 51 Peg b (Mayor and Queloz, 1995) thousands of exoplanets have been detected so far [see the Extrasolar Planet Encyclopedia 1 Schneider et al. (2011) or the NASA Exoplanet archive 2 ]. However, of the nearly 3,500 planetary systems none harbors a habitable world similar to planet Earth. This fact raises the question if Earth-like habitability requires a certain planetary system architecture where e.g., giant planets are orbiting the host-star exterior to the habitable zone (HZ). It can be expected that such a configuration has a less violent terrestrial planet formation process than a planetary system where the giant planet is closer to the host-star than the HZ. The giant planet formed probably at a larger distance to the star and migrated inwards by crossing the HZ which could cause problems for the terrestrial planet formation in the HZ.
Thus, in this study we considered only systems where the giant planet is exterior to the HZ and investigate the occurrence of strong perturbations in the HZ. Studies of planetary motion in binary stars are of special interest as a high fraction of stars in the solar neighborhood are members of binary and multiple star systems. Observational surveys by Duquennoy et al. (1991) and Raghavan et al. (2010) established that in the solar neighborhood (up to 25 parsec) about 40-45 % of all Sun-like stars (spectral types F6-K3) are members of binary and multiple star systems. Tokovinin (2014) derived a fraction of 33 % of binary stars from a sample of about 4800 F-/G-type mainsequence stars within 67 parsec of the Sun. For more details see Duchêne and Kraus (2013).
The Catalog of Exoplanets in Binary Star Systems 3 lists about 100 planetary systems which indicates that exoplanets are not restricted to single-stars. But it is doubtful whether these environments are more hostile for the existence of planets or not (see e.g., Boss, 2006;Bromley and Kenyon, 2015;Jang-Condell, 2015). In a study based on Kepler data (Armstrong et al., 2014) the occurrence rate of co-planar circumbinary planets has been found to be similar to that for single stars. Circumbinary motion is also known as P-type motion (Dvorak, 1984(Dvorak, , 1986 where a planet orbits both stars. However, the Catalog of Exoplanets in Binary Star Systems indicate that most of the exoplanets in binary star systems orbit only one star which is known as circumstellar or S-type motion (Dvorak, 1984(Dvorak, , 1986. Theoretical stability studies of S-and P-type motion have been carried out decades before the detection of exoplanets (see e.g., Harrington, 1977). Studies by Rabl and Dvorak (1988) and Dvorak et al. (1989) published the stability boundaries of S-and P-type motion in equal-mass binary stars as a function of the binary's eccentricity. This study has been repeated and extended to a wider range of binary star configurations by Holman and Wiegert (1999). In their study they provide expressions for the stability boundaries which depend on the mass-ratio and the binary's eccentricity for both, circumstellar and circumbinary configurations. In connection to this study, Pilat-Lohinger and Dvorak (2002) analyzed the stability of S-type motion using the Fast Lyapunov indicator (FLI) 4 . In addition, Pilat-Lohinger et al. (2003) investigated the stability of inclined P-type planetary orbits.
The detections of exoplanets in tight binary stars like γ Cephei b (Hatzes et al., 2003), Gliese 86 b (Santos et al., 2000), and HD 41004 Ab (Zucker et al., 2004) have led to a growing interest in understanding planetary formation in such stellar systems. There are binary specific problems for planet formation due to the gravitational interaction of the secondary star resulting in e.g., a truncated protoplanetary disk (see e.g., Artymowicz and Lubow, 1994;Savonije et al., 1994) which influences the formation and evolution of planets throughout several stages of the planet-forming process. For more details about the problems of planet formation in binary stars we refer the reader to Marzari and Thebault (2019) and references therein.
In case of circumstellar motion mainly the outer edge of the disk is influenced (Kley and Nelson, 2010;Müller and Kley, 2012) while for circumbinary disks mainly the inner edge is affected (Rafikov, 2013).
During the planet-formation stage where planetesimals (kmsized bodies) collide and merge to Moon-sized embryos dynamical perturbations play an important role, as planetesimal accretion requires low encounter velocities. Thébault et al. (2006) showed that planetesimal accretion cannot occur at distances beyond 1 au from the host-star in tight binaries with separation of 20 au. This result raises questions about the formation of detected planets in tight binary stars (Thébault et al., 2004(Thébault et al., , 2008(Thébault et al., , 2009Thebault, 2011). Moreover, Fragner et al. (2011) have found an increase of impact velocities when studying the dynamical behavior of planetesimals taking into the self gravity of the gas disk. Some years later, Gyergyovits et al. (2014) studied the full interaction of more than 2,000 embryos with a gas disk in a γ Cephei-like configuration. They concluded that the growth from embryos to planets within a dynamically evolving gas disk is strongly altered by the dynamical evolution of the disk, which leads to a decreased probability for planet formation at least in the inner parts of the gas disk.
Apart from these problems of planetary embryo formation, the late planet-formation stage where embryo-sized bodies grow to planets can be easily simulated by N-body calculations (Haghighipour and Raymond, 2007;Quintana et al., 2007;Pilat-Lohinger et al., 2018). A detailed discussion on planet formation in binary star systems can be found in Thebault and Haghighipour (2015).
A weak point in all these numerical studies is certainly the treatment of collisions where usually the so-called perfect merging of two bodies is assumed without taking into account any fragmentation or mass-loss of the bodies due to the collision. However, Bancelin et al. (2017) studied the collision parameters (impact angle and velocity) in tight binary stars and showed that the perfect merging approach can significantly overestimate the water content and mass of the formed planets. Thus, this simplified approach has to be replaced by a more realistic one that includes results of detailed collision simulations.
The habitability of planets in binary stars is certainly an interesting issue, especially due to the fact that most of the stars in the solar neighborhood build such stellar systems. Considering the different spectral types, especially F to M type stars are of interest for habitability studies, since the life-times of these stars on the main sequence are sufficiently long (see Kasting et al., 1993) to permit the evolution of life on a terrestrial-like planet in the HZ. For a summary of astrophysical conditions for planetary habitability we refer the reader to Güdel et al. (2014).
In a binary star system, an important requirement for habitability is certainly the dynamical stability of the HZ. Moreover, the influence of the secondary star will increase the planet's eccentricity which also affect the insolation on a planet in the HZ where the perturbations depend on the distance and the eccentricity of the secondary star. Eggl et al. (2012) has found a correlation of both, eccentricity and insolation and introduced various classes of HZs in binary star systems: (i) permanent, (ii) extended, and (iii) averaged HZ. An application of this HZ classification to binary star systems in the solar neighborhood is shown in Eggl et al. (2013).
In recent years, different approaches for the determination of HZs in binary stars have been published e.g., Kaltenegger and Haghighipour, 2013;Cuntz, 2014). For a detailed review on this topic we refer the reader to Eggl et al. (2020).
In this study, we focus on the gravitational influence of the secondary star on planetary motion in the HZ. We consider binary star systems with a gas giant in an orbit exterior to the HZ and a terrestrial planet in the HZ. The interplay of these bodies causes perturbations like mean motion resonances (MMRs) and secular resonances (SRs) whose locations strongly depend on the architecture of the system. Such resonances can also influence the habitability of planets in the HZ. Thus, it is important to determine the locations of resonances, which is easy for MMRs but not for SRs. However, for the SRs new methods have been developed that allow a quick determination of the SR location (Pilat-Lohinger et al., 2016;Bazsó et al., 2017;Bazsó and Pilat-Lohinger, 2020).
In this paper, we will briefly discuss these methods and show some applications to real binary systems that host a giant planet. The structure of this paper is the following: First we define the dynamical model and describe the gravitational perturbation. Then we introduce the semi analytical method (SAM) to determine the location of an SR. Subsequently, we describe the basics of the recently published combined analytical method (CAM) and its application for the online tool named SHaDoS. Finally, we apply SHaDoS to some real binary-starplanet systems and discuss the results.

Dynamical Model
In our study, we focus on circumstellar planetary motion where (i) a terrestrial planet is located in the primary stars' HZ, (ii) a gas giant is moving in an orbit exterior to the HZ, and (iii) the secondary star is a far away perturber. First of all, the planets have to move on stable circumstellar orbits around the primary star. The stability of the planetary motion depends strongly on the distance of the two stars, their eccentricity and their masses. To ensure that the outer planet i.e., the giant planet moves inside the stable region of the primary star, it is advisable to apply the investigations by Pilat-Lohinger and Dvorak (2002) and Holman and Wiegert (1999). However, the situation for two planets is more complex due to arising resonances among the planets which may overlap with those of the secondary star. Therefore, Marzari and Gallina (2016) investigated such systems using a Frequency Map Analysis method to study the influence of the secondary star on the two planets orbiting the primary star. In their study, they obtained a semi-empirical equation that defines the minimum semi-major axis of the secondary star for which the two-planet system is stable.
Nevertheless, the presence of a distant secondary star will perturb the giant planet, which might pass the perturbations to the terrestrial planet in the HZ. In our study, we are primarily interested in the orbital behavior of the terrestrial planet in the HZ and how perturbations of the giant planet and the secondary star may affect its motion. In case these perturbations cause high eccentricity motion in the HZ then the planet's habitability would be affected, especially when it leaves the HZ periodically.

Gravitational Perturbations
It is well-known that in N-body systems consisting of more than two massive bodies variations of the orbital parameters occur due to gravitational interactions between these bodies. Thus, in binary star systems that host a planet such variations occur. These mutual gravitational interaction between celestial bodies can lead to resonances which appear when the ratio of two frequencies f 1 and f 2 can be expressed as a rational number: |f 1 /f 2 | = p/q, p, q ∈ N. Depending on the frequencies that are involved, a typical time-scale can be associated to any resonance. As resonances may lead to significant changes in the orbital motion with large variations in eccentricity, the planet's habitability could be affected. Usually, the HZ is a quite narrow region thus, the planet might leave the HZ periodically due to its eccentricity motion. Therefore, it is important to determine the locations of resonances, where mean motion and secular resonances are considered in our study. Mean Motion Resonances (MMRs) occur when n ′ /n ∼ j 1 /j 2 where n and n ′ are the mean motions 5 of e.g., the two planets and j 1 , j 2 are integers. The location of an MMR can be easily derived from the third Kepler law: where a res is the resonant semi-major axis and M, m and m ′ are the masses of the host-star and the two planets, respectively. MMRs can be the source of both stability and chaos, which depends sensitively upon the orbital parameters. Secular Resonances (SRs) arise when one of the precession frequencies of a test-planet (either of the argument of perihelion ω or of the node ) is equal (or a linear combination) of the proper modes of the giant planet. The strongest secular perturbations occur at locations when where q, p are integers, g TP and g GP are the proper frequencies of the test-planet and the giant planet, respectively. For p = q = 1, the proper frequencies of the test-planet and the giant planet are equal which causes the strongest effects. The location of an SR with the giant planet can thus be determined by applying Equation (2) to test-planets of a certain area. Therefore, the proper frequencies of the orbits of all test-planets and the giant planet have to be determined which usually requires long-term computations of the dynamical system where the integration time depends on the distance of the two stars. For a tight binary star like γ Cephei or HD41004 AB with a stellar separation of 20 au, calculations over some 10 6 years are needed. However, the computation time increases with the stellar distance a B (as discussed in Bazsó et al., 2017). Thus, for wide binaries such a study can be quite time-consuming. In addition, the time series of the orbital computations must be analyzed using e.g., a Fast Fourier Transformation to determine the main frequency in the Fourier spectrum of an orbit which represents the proper frequency. Fourier spectra of strongly perturbed orbits are littered with frequencies which makes it is impossible to determine their proper frequencies and is also an indication for a perturbation due to a resonance. To overcome the time-consuming procedure of identifying the proper frequency via Fourier analysis of each orbit, a new semi-analytical approach has been introduced by Pilat-Lohinger et al. (2016) and Bazsó et al. (2017).

QUICK METHODS TO DETERMINE LOCATIONS OF AN SR
For the development of our semi-analytical method, the tight binary star HD41004 AB has been used. In this binary star system a giant planet (of 2.5M Jup ) has been discovered (Zucker et al., 2004)

Semi-Analytical Approach
In the study by Pilat-Lohinger et al. (2016) the numerical effort for the determination of the proper frequencies could be reduced by applying the Laplace-Lagrange secular perturbation theory (see e.g., Murray and Dermott, 1999) to obtain an analytical solution for test-planets (with negligible mass) in the area between 0.2 and 0.6 au in the HD41004 AB system. The secular frequency of test-planets (g TP ) was deduced using the following secular linear approximation (see e.g., Murray and Dermott, 1999): where α i = a/a i for i = 1, 2 and a 1 , a 2 , a are the semimajor axes of the giant planet, the secondary and the test-planet, respectively. m 1 and m 2 are the masses of the giant planet and the secondary and m 0 is the mass of the primary star. b (1) 3/2 is a Laplace coefficient. Considering a planar system, the secular frequency s which is related to the node , has not been taken into account. Moreover, this theory is restricted to low inclination and low eccentricity motion. Pilat-Lohinger et al. (2016) showed that the analytical values of proper frequencies for the test-planets were in good agreement with the numerical solution and save computation time for the determination of the location of SR.
Thus, to satisfy Equation (2) only the numerical integration of the giant planet's orbit in the binary star has to be carried out for a computation time which covers the secular period of the binary star. From this numerical simulation, the proper frequency of the giant planet has to be determined via Fourier analysis. Therefore, this semi-analytical method (SAM) represents a quick approach to define the position of SRs for circumstellar planetary motion in binary stars.

The Location of the SR
The location of the SR is defined by the semi-major axis of the test-planet which proper frequency equals the one of the giant planet. Figure 1 illustrates the application of SAM to different binary-star-planet configurations to figure out the influence of the secondary star on the location of the SR. For these plots the HD41004 AB configuration has been used, where the distance between the two stars was varied between 20, 30, and 40 au (in the left panel) and the mass of the secondary star was changed from an M-star (of 0.42 solar mass) to either a K-(of 0.7 solar mass) or G-type star (of 1 solar mass). The eccentricities of the binary and of the giant planet were set to 0.2.
The left panel of Figure 1 shows the influence of the distance of the two stars on the location of the SR. The diagonal curve represents the analytically derived proper periods of test-planets orbiting the host-star in the region between 0.1 and 1.5 au. Whereas, this curve does not vary significantly for the different stellar separations, strong changes are visible for the proper period of the giant plant which is represented by the horizontal line in the corresponding color of the stellar distance. The intersection of an horizontal line with the same colored curve defines the location of the SR in the system. Different intersection points show clearly that the SR is shifted to smaller semi-major axes-thus toward the host-star-with increasing the distance of the two stars.
In the same way, the right panel of Figure 1 indicates the influence of the mass of the secondary on the location of the SR. For low-mass M-type stars (0.4 M Sun , red line/curve) it can be seen that the SR is closer to the host-star than for a more massive secondary, like a K-type (0.7 M Sun , green line/curve) or a G-type star (1 M Sun , blue line/curve). Furthermore, the study by Pilat-Lohinger et al. (2016) showed that a change in the binary's eccentricity has stronger effects on the location of the SR than a change in the giant planet's eccentricity. Note that eccentricity motion widens the area which is affected by the SR. To determine the width of the SR, Pilat-Lohinger et al. (2016) calculated the SR location by using the periand apo-center positions instead of a GP and found good results for the binary star HD41004 AB.
An application of SAM to binary star systems with stellar separations up to 100 au and a discovered exoplanet in a circumstellar orbit has been published by Bazsó et al. (2017). The aim of this study was to check the HZ of the host-stars regarding secular perturbations that could alter the conditions for habitability. In this context, SAM has been used also for systems where the giant planet orbits the host-star interior to the HZ (= interior case) and not only for systems where the giant planet is exterior to the HZ (= exterior case) like in the test-system HD41004 AB. In a first approach SAM did not work for the numerous interior cases because no crossing of the giant planet's proper frequency with one of the test-planet's proper frequency was found. Only by adding the effect of general relativity for the interior cases, the method worked well for these systems Bazsó et al. (see 2017). Even if the application of SAM saves a lot of computation time when studying numerous binary-star-planet configurations regarding dynamical conditions for habitability. However, the method cannot be used for an internet-tool.
Therefore, further improvements are needed to speed up the procedure which is only feasible with an analytical approach for the giant planet's proper frequency.

Analytical Model for Giant Planet's Proper Frequency
Assuming that the dynamical evolution of the giant planet's orbit is dominated by the secular interaction with the perturber (the secondary star), a simple approach to obtain analytical values for the giant planet's proper frequency has been provided by Heppenheimer (1978): where g H is a first order approximation (in masses) for the proper frequency of the giant planet, n GP is its mean motion, m A and m B are the masses of the primary and the secondary star, a GP and a B are the semi-major axes of the giant planet and the secondary star and e B is the eccentricity of the binary. The Heppenheimer model has been developed for the restricted three body problem, where the planet's mass is negligible compared to the stellar masses. A suitable improvement 6 of the Heppenheimer model is provided in 6 Such a correction to the secular precession frequency can be found in previous attempts, e.g., Giuppone et al. (2011). the study by Andrade-Ines and Eggl (2017) who established a simple and more accurate secular model where the giant planet acts like a mass-less body relative to the two stars. This model can be used also for tight binary stars with stellar separations of only two to three times the distance of the giant planet to the host-star.
The new approach by Andrade-Ines and Eggl (2017) is defined as: where δ g is an empirical correction term. The model expresses the secular precession frequency g GP as a polynomial function of (i) the mass-ratio of the binary star, (ii) the semi-major axis ratio a GP /a B , and (iii) the perturber's eccentricity e B . The expression of δ g has usually less than 20 terms. For the exact values we refer to Equation (24) in Andrade-Ines and . Equation (5) is applicable to a wide range of binary stars. Regarding the semi-major axes ratio, it covers the range of currently observed exoplanets in binary systems. Only the range of binary's eccentricities does not cover all systems where an exoplanet has been discovered.
Applying Equation (5) for the solution of g GP in Equation (2) we obtain a fully analytical approach which is called Combined Analytical Method (CAM) in Bazsó and Pilat-Lohinger (2020). The CAM is a promising method that allows a very quick determination of the location of secular perturbations for circumstellar planetary motion in binary star systems.
Our main interest is to define the dynamical state of circumstellar HZs and to classify perturbed (pHZ) and quiet HZs (qHz). In this regard, CAM might help to exclude binary-starplanet systems (with pHZs) from the observational search of habitable planets. Moreover, CAM is fast enough to be applied for an internet-tool.

ONLINE-TOOL: SHADOS
SHaDoS is an acronym for Secular perturbations in Habitable zones of Double Stars and is accessible at https://www.univie. ac.at/adg/shados/index.html. It implements the CAM model and solves Equation (2) for any given set of orbital parameters. Figure 2 presents the flow chart of this online-tool which shows an object-oriented approach where a four-step input process yields the desired results in 2D diagrams. In each step the tool needs input from the user about the parameters of (i) the hoststar, (ii) the giant planet, (iii) the secondary star, and (iv) the parameter space of the resulting plot.
Step 1: The required input parameters of the host star are the luminosity, effective temperature, and mass. According to these parameters, the tool will determine the borders of the host-star's HZ 7 using the approach by Kopparapu et al. (2014) where the effective temperatures are limited to the range 2600 ≤ T eff ≤ 7200 K. For some stellar types, namely F, G, K, and M-type main-sequence stars the standard input data is provided by the tool which can be modified via the user defined input. Based on this input the HZ will be calculated. 7 Note that the single star HZ corresponds to the averaged HZ in binary stars.
Step 2: The required parameters of the giant planet are the mass (in Jupiter-masses), distance to the host star (semi-major axis in au) which must be exterior to the HZ, and orbital eccentricity which is restricted to elliptic orbits.
Step 3: The required parameters of the secondary star are the mass (in solar mass), distance to the host-star (in au), and the eccentricity which is restricted to values 0 ≤ e B ≤ 0.6 due to limitations in the analytical model.
Step 4: Following parameters for the 2D-plots can be chosen by the user: there are five options for the x-axis, which are planet mass, secondary star mass, planet distance, secondary star distance, and secondary star eccentricity for which the user has to specify the minimum and maximum values and the steps in-between (for the grid). And for the y-axis there are three options which are planet distance, secondary star distance, and secondary star eccentricity.
With the input of all parameters the online-tool performs the calculations and provides a 2D plot depending on the selection in Step 4 which can be saved if desired. This plot usually shows a gray-shaded area which defines the parameters that might perturb the host-star's HZ.
Due to the constraints of the methods used in SHaDoS and the requirement that the giant planet has to be exterior to the HZ, the application to real systems is still very limited. In addition, for most binary star systems the eccentricity is not known which also has strong influence on dynamical studies. Even if nowadays the number of suitable circumstellar systems is rare, we expect an increase of such systems in the near future, especially when the PLATO 2.0 mission starts the observations.

APPLICATION OF SHADOS
To demonstrate the functionality of SHaDoS, we have selected two binary system from the Binary Catalog of Exoplanets which fulfill the requirements for an application of the online-tool. These are the wide binary star HD106515 AB and the tight binary HD41004 AB. The latter has also been used for the development of SAM. For the first system, we show all possible combinations provided in Step 4 while for the other for HD41004 AB we present only the most interesting plots.

The Wide Binary Star HD106515 AB
HD106515 AB consists of two G-type stars, with masses of 0.91 and 0.88 solar-masses for the primary and the secondary star, respectively. The distance of the two stars is 345 au and the binary's eccentricity is 0.42. Desidera et al. (2012) detected a quite massive planet (9.61 Jupiter-masses) orbiting the primary star at a distance of 4.59 au in an eccentric orbit (e = 0.572). Since the HZ of a G-type star is between 0.95 and 1.7 au, the giant planet is exterior to the HZ as required for the application of SHaDoS.
The first 2D-plot of SHaDoS (Figure 3) shows the secondary's semi-major axis vs. the giant planet's semi-major axis, where the gray-shaded area indicates all combinations of semi-major axes in the chosen range of the x-and y-axes which would cause a secular perturbation in the HZ of the primary star. Considering the observational data of the HD106515 AbB system, Figure 3 shows clearly, that one of the two celestial bodies must be moved quite a lot, either closer to the primary (i.e., the secondary to 50 au ) or farther away (i.e., the planet to 18 − 22 au) to find an SR in the primary's HZ. Since the orbital parameters of exoplanets show large uncertainties and some of them might change when new observational data is available, the plots from the application of SHaDoS indicate also to what extend such changes of the orbital parameters might influence the dynamical state of the HZ (e.g., if a qHZ turns into a pHZ or vice versa). For the HD106515 AbB system such changes must be significantly large which seems unlikely to happen thus HD106515 A has a quiet HZ.
More possible parameter combinations-according to the input in Step 4-are shown in Figures 4-6. Each panel of these figures shows a gray-shaded area from which we get the information whether the host-star's HZ is perturbed or not. For the observed configuration of HD106515 AbB the resulting 2Dplots of SHaDoS clearly indicate a qHZ of HD106515 A. Secular perturbations in the HZ would occur only if at least one of the parameters shown in the various panels is changed to values in gray area.

The Tight Binary Star HD41004 AB
HD41004 AB consists of a K-and an M-type main-sequence star as primary and secondary, respectively. The masses are 0.7 solar mass for the K star and 0.4 solar mass for the M star. Actually, both stars have a sub-stellar companion, but in our study, the close-in brown dwarf orbiting the M-type star has been ignored.
In the vicinity of the primary (between 1.3 and 1.64 au) a gas giant (of about 2.5 Jupiter-mass) has been detected (Zucker et al., 2004). It orbits the K star exterior to the HZ. Thus, the system allows an application of SHaDoS. Different values for the eccentricities of the binary and the giant planet have FIGURE 3 | Result of the online-tool SHaDoS for the binary system HD106515 AB. The gray stripe shows which combinations may lead to an SR in the primary star's HZ for the defined semi-major axis range of the secondary star (x-axis) and the giant planet (y-axis). Figure 3 but for the variation of the stellar mass (x-axis) vs. the variation of the semi-major axis either of the giant planet (Upper) or of the secondary star (Lower) in the y-axes.

FIGURE 4 | Same as
been published but these are quite uncertain. Thus, in the application of SHaDoS, the eccentricities were set to 0.2 for both bodies. Figure 7 shows three 2D-plots of the online-tool SHaDoS, where the top panel displays the semi-major axis of the secondary star on the x-axis and that of the giant planet on the y-axis. This plot indicates that for the published values of the semi-major axes (a GP = 1.64 au and a B = 23 au) there are secular perturbations in the vicinity of the HZ of HD41004 A but the HZ itself is not perturbed. The middle panel displays the secondary star's semimajor axis vs. binary eccentricity which indicates that in case of high eccentricity motion of the binary (0.4 < e B < 0.9) there would be a high probability that the HZ of the K-type star would be perturbed by an SR. The bottom panel indicates the same but for different masses of the giant planet on the x-axis.
Since the eccentricities of the HD41004 AbB system are not well-determined, it is quite probable that the HZ of the K-type star might change into a pHZ if the orbital parameters vary due to new observational data.

DISCUSSION AND CONCLUSION
In this study, we presented recently developed methods to determine the location of secular perturbations for circumstellar (or S-type) planetary motion in binary stars with a special emphasis to planets in the cirumstellar HZ in such systems. We did not take into account circumbinary (or P-type) planetary FIGURE 6 | Same as Figure 3 but for the variation of the giant planet's mass (x-axis) vs. a variation of the binary's eccentricity on the y-axis (Upper). The (Lower) shows a variation of the secondary star's semi-major axis (x-axis) vs. the binary's eccentricity (y-axis). motion as we assume stronger stellar perturbations for a planet in the HZ due to interactions of the stars (e.g., colliding stellar winds; Johnstone et al., 2015) especially in the early phase of a system which could severely affect the conditions for planetary habitability.
The appearance of secular resonances (SR) in circumstellar HZs was investigated for certain binary-star-planet configurations where the giant planet orbits the host-star exterior to the HZ. Considering planar systems, SRs are caused by the precession of the giant planet's perihelion which results from gravitational interaction with the secondary star. Moreover, limitations of the adopted methods have restricted our study to eccentricities ≤ 0.6 for the planet and the binary star.
In a first step we replaced the fully numerical approach by a semi-analytical method (SAM), which uses the Laplace-Lagrange theory (Murray and Dermott, 1999) to define the proper frequencies/periods of all test-planets moving in a certain area of the host-star. Pilat-Lohinger et al. (2016) showed that FIGURE 7 | Results for HD41004 AbB. The gray stripe in each panel indicates conditions that cause perturbations in the primary's HZ due to an SR. The top panel shows the result for the semi-major axis of the secondary star (x-axis) and the giant planet (y-axis). In the middle panel the y-axis shows the binary's eccentricity. In the bottom panel, the x-axis shows the giant planet's mass. these analytically determined values of proper frequencies are in good agreement with the numerical results. Thus, the application of SAM removed a huge amount of time-consuming numerical computations. Only one computation for the determination of the giant planet's proper frequency remains.
However, a general parameter study to distinguish between systems with either qHZ or pHZ needs a parameter space of at least five variables (i.e., mass-ratio of the binary star, mass of the giant planet, semi-major axes of the secondary star and the giant planet, and the eccentricity of the binary) where a change of each parameter will modify the location of a resonance makes therefore, the application of SAM is still too costly in terms of time. Only a fully analytical approach could realize such a study in a reasonable time.
We strove to replace the numerical part for the determination of the giant planet's proper frequency by an analytical method. In this context, the approach by Andrade-Ines and  has been found to be appropriate for our studies. By joining the Laplace-Lagrange method and the Andrade-Ines-Eggl approach we obtained a combined analytical method (CAM) which can easily perform all computations of the 5-dimension parameter space very quickly.
Instead of doing a huge amount of calculations in order to compile a catalog which states the dynamical behavior of the circumstellar HZ in various binary-star-planet configurations, we decided to develop an online-tool which applies CAM. The online-tool is named SHaDoS (=Secular perturbations in Habitable zones of Double Stars) 8 and is easy to use. After a 4-Step input procedure the user gets a 2D-plot which indicates the parameters which would cause a secular resonance in the circumstellar HZ of the host-star by a gray-shaded area. To demonstrate SHaDoS we showed applications to a wide binary star (HD106515 AB) and a tight binary star (HD41004 AB) and analyzed the primaries' HZ regarding secular perturbations. Using the observational data, we have found a qHZ for both systems. 8 https://www.univie.ac.at/adg/shados/index.html Moreover, the plots of the wide binary star indicate that the probability of a change in the dynamical behavior of the HZ is extremely low-only if new observations should yield significantly different orbital parameters. While in case of the tight binary star changes of the orbital parameters could lead to a pHZ.
Thus, the presented methods can be considered as an important contribution to the habitability research of exoplanets as an SR could affect the motion of a planet in the HZ by increasing the eccentricity. Such variations of eccentricity due to SRs usually happen on a long time-scale and could lead to a misestimation of planetary habitability in such systems. Detected planets in the HZ on nearly circular orbits could move onto highly eccentric orbits after some thousands to millions of years. This would certainly influence the conditions of habitability of these planets.
Even if the number of real systems for which SHaDoS can be applied is still very low, we expect that the number of such systems will increase significantly in the near future when large telescopes like ELT or the Rubin Observatory start observing and especially, when the space mission PLATO 2.0 is launched and starts to explore the sky.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.