The Formation of Habitable Planets in the Four-Planet System HD 141399

TL;DR

This study uses numerical simulations to explore the potential formation of habitable terrestrial planets in the HD 141399 system's habitable zone, revealing that such planets could have formed there through natural processes.

Contribution

It demonstrates the feasibility of terrestrial planet formation in HD 141399's habitable zone using detailed N-body simulations including water content and gravitational interactions.

Findings

Terrestrial planets could form in the habitable zone of HD 141399.

Planets tend to form near mean motion resonances.

A belt-like structure analogous to the Solar System's asteroid belt appears.

Abstract

The presented work investigates the possible formation of terrestrial planets in the habitable zone (HZ) of the exoplanetary system HD 141399. In this system the HZ is located approximately between the planets c (a = 0.7 au) and d (a = 2.1 au). Extensive numerical integrations of the equations of motion in the pure Newtonian framework of small bodies with different initial conditions in the HZ are performed. Our investigations included several steps starting with 500 massless bodies distributed between planets c and d in order to model the development of the disk of small bodies. It turns out that after some 10^6 years a belt-like structure analogue to the main belt inside Jupiter in our Solar System appears. We then proceed with giving the small bodies masses (Moon-mass) and take into account the gravitational interaction between these planetesimal-like objects. The growing of the objects - with certain percentage of water - due to collisions is computed in order to look for the formation of terrestrial planets. We observe that planets form in regions connected to mean motion resonances (MMR). So far there is no observational evidence of terrestrial planets in the system of HD 141399 but from our results we can conclude that the formation of terrestrial planets - even with an appropriate amount of water necessary for being habitable - in the HZ would have been possible.

Figures

Tables

Table 1. Table 1: A few examples of confirmed planetary systems with at least one gas giant at an orbit with a > > 2.0 au (check https://www.exoplanet.eu/ for updates and more detailed information). The system of 47 UMa is topic of another paper within our group (Cuntz et al., 2018 ) .

System	$m_{\mathrm{star}}$	nr of	nr with	m $>$ 0.8 $m_{\mathrm{Jup}}$
	($m_{\mathrm{Sun}}$)	Planets	m $>$ 0.8 $m_{\mathrm{Jup}}$	a $>$ 2.0 au
HD 125612	1.1	3	2	1
HD 204313	1.045	3	2	2
47 UMa	1.03	3	2	1
GJ 676 A	0.71	4	2	1
mu Ara	1.08	4	2	1
ups And	1.27	4	3	2
HR 8799	1.56	4	4	4
55 Cnc	0.95	5	2	1
HD 131499	1.07	4	2	1

Table 2. Table 2: Properties of the planets in the system HD 141399. The 3:2 MMR is called ”Hilda region” and resides in the HZ of HD 141399.

Planets	Period (d)	e	a (au)	mass ($m_{\mathrm{Jupiter}}$)
HD 141399 b	94	0.04	0.4225	0.46
HD 141399 c	202	0.05	0.7023	1 36
3:2 MMR	732		1.63
HD 141399 d	1070	0.06	2.1348	1.22
HD 141399 e	3717	0.00	4.8968	0.69

Table 3. Table 3: The proper frequencies of four planets.The f 𝑓 f , g 𝑔 g and s 𝑠 s represent the frequencies of mean motion, periastron precession and nodal precession, respectively. The values are all given in a unit of 2 π / 2\pi/ yr. The sign ± plus-or-minus \pm indicates the precession direction, ‘ + + ’ for prograde and ‘ − - ’ for retrograde. The nodal precession of Planet d is so slow that we cannot detect it.

	Planet b	Planet c	Planet d	Planet e
$f$	$+3.866383$	$+1.811227$	$+3.461740\times {10}^{-1}$	$+9.438578\times {10}^{-2}$
$g$	$+3.751551\times {10}^{-3}$	$+3.406273\times {10}^{-4}$	$+6.832154\times {10}^{-5}$	$+2.066966\times {10}^{-5}$
$s$	$-2.276786\times {10}^{-3}$	$-1.240933\times {10}^{-4}$	$0.000000$	$-3.178531\times {10}^{-5}$

Full text

1 Introduction

The ’race’ to find a second Earth led to the discovery of many such candidates with different techniques from space and from ground. Nevertheless, although we have knowledge of systems hosting terrestrial111using the definition of a terrestrial planet meaning a planet smaller than 5 Earth-masses with a rocky surface planets without gas giants – for example Trappist 1 2227 terrestrial planets moving partly in the HZ, see http://exoplanet.eu/ (e.g., Gillon et al., 2016) – our Solar System is an example of hosting big planets outside of the habitable zone (HZ). This zone was proposed by Kasting et al. (1993) and later revised by e.g., Kopparapu et al. (2013) and describes the zone where water can exist in liquid form on a terrestrial planet orbiting the star.

Since the first discovery of a big planet in 1995 (Mayor & Queloz, 1995) – a so called hot Jupiter – many gas giants have been found around stars of different spectral type.

The discoveries of many hot Jupiters over the past 20 years and the implications of their existence for the process of planet formation make it difficult to believe that terrestrial planets may co-exist in the HZ of those systems (Dawson & Johnson, 2018).

Systems with only one gas giant are well represented in the long list of extrasolar planets333see http://exoplanet.eu/, but in our work we concentrate on systems with more than one gas giant of which several ones are known (see Table 1). Of special interest to us seems the recently discovered four-planet system HD 141399 (Vogt et al., 2014).

Vogt et al. (2014) carefully describe the method of determining the orbital elements of the planets and check their stability which seems to be guaranteed because of the almost circular orbits of the four planets.

The aim of our research is to find whether terrestrial planets could be have been formed in the HZ around the late K-type star.

The structure of this research paper is as follows: In Section 2 we describe the system of HD 141399 in detail. The method of investigation is presented in Section 3. The results of the integration of massless bodies in the HZ between planet c and d are shown in Section 4, where we present a thorough analysis of the outcomes and show examples of places where terrestrial planets can form on stable orbits. The results were gained using extensive numerical simulations with integration times of up to 10 Myr. In Sub - section 4.2 we take into account that the observations were made using radial velocity measurements which yield only a minimum mass for the planets and concetrate on the formation of bodies in the system with assumed larger masses for the giant planets.

Section 5 shows the results of computations with massive bodies which formed terrestrial planets through collisions. In Section 6 we describe the growth process from planetesimals to protoplanets in more detail. Finally in Sections 7 and 8 we show the main results and draw the conclusions with focus on possible terrestrial planets in the HZ.

In Appendix A.1 we investigate the dynamics around the Lagrangian Point $L_{4}$ of planet d in detail and describe the acting secular resonances in the system in Appendix A.2. We note that the Trojan region of planet d is just outside of the HZ. In the presented article we included the computations for possibly habitable planets formed in this region taking into account the eccentricity of planet d. On its orbit it will sometimes be inside and sometimes just outside the defined HZ.

2 The System of HD 141399

It is important to have an overview of possibly stable terrestrial planets in known multiple planetary systems with gas giants (Agnew et al., 2018). Nevertheless, we think that each system should be treated separately because of the different architecture. The dynamical stability of additional – possibly terrestrial – planets in a planetary systems with gas giants is only one side of the coin, the other side is to know how these planets could have formed. The present knowledge is that in our Solar System the gas giants formed prior to the terrestrial planets in the Solar environment of gas and dust. Different scenarios are possible (e.g., Grand Tack as described by Morbidelli et al., 2011). We know that in our Solar System no stable terrestrial planet could exist between Jupiter and Saturn – they all formed inside the Jovian zone.

We looked into the literature to find an interesting system quite similar to our own and found that HD 141399 would be a good example (Vogt et al., 2014). The HZ in HD 141399 is located between two gas giants (see Table 2), namely between planet c (a = 0.7 au) and planet d (a = 2.1 au). Calculation of the boundary values for the habitable zone shows that it lies between 0.8 au and 2.0  au from the star (Selsis et al., 2007).

HD 141399 is a K0V star with m = 1.07 $\mathrm{m_{Sun}}$, r = 1.46 $\mathrm{r_{Sun}}$, and $\mathrm{T_{eff}}$ = 5600 K. The system was presented by Vogt et al. (2014) who used radial velocity data sets from Keck-HIRES and the Lick Observatory’s Automated Planet Finder Telescope and Levy Spectrometer on Mt. Hamilton. 91 observations from over 10.5 yrs were analyzed and the parameters of four planets shown in Table 2 have been found. Interestingly the planets HD 141399 b and HD 141399 c are close to the 2:1 mean motion resonance (MMR) which might give a clue on the history of the system. Vogt et al. (2014) argue, that a possible inward migration would point towards a capture in resonance, which cannot be the case, given the values derived from the observations. Conversely, Batygin & Morbidelli (2013) show, that a divergence away from exact resonance behavior is a natural outcome of dissipative evolution of resonant planetary pairs. Vogt et al. (2014) carefully describe the method of determining the orbital elements of the planets and checked their stability, which seems to be guaranteed because of the almost circular orbits of the four planets.

3 The Method

The method used for the investigation of the formation of possible terrestrial planets between HD 141399 c and HD 141399 d was the Lie-ingration method which has already been used extensively by our group (e.g., Dvorak, 1986; Delva, 1984; Hanslmeier & Dvorak, 1984; Lichtenegger, 1984). It is well adapted for the integration of the equations of motion for planetary systems hosting massless and massive bodies and treats close encounters and collisions with a high precision due to its automatically chosen step size. The method has been compared to other methods in detail by Eggl & Dvorak (2010).

We ask the question where and how these terrestrial planets could have formed in HD 141399. For being habitable the planets need to have water and therefore we included water content of the planetesimals in our formation computations.

Making use of the Lie integration method four main runs were undertaken:

• •

First we distributed 500 massless bodies in the region between 0.7 au and 2.1 au (the location of planet c and planet d) with initial conditions randomly distributed for the orbital elements: semi-major axis a ($0.75\,au<\mathrm{a}<1.75\,au$), eccentricity $\mathrm{e}<0.13$, inclination i ($0.5^{\circ}<\mathrm{i}<1.5^{\circ}$), $\Omega=154^{\circ}$, and $\omega=0.5^{\circ}$. The mean anomaly was randomly distributed. The whole system was integrated for 1 Myr for 8 differently chosen initial conditions. See Figure 1.

• •

In the next step we replaced the massless bodies with massive ones ($0.92\,\mathrm{m_{Moon}}<m<1.23\,\mathrm{m_{Moon}}$) with a total mass of approximately 5 $\mathrm{m_{Earth}}$, which was for example also used in Quintana & Lissauer (2014)444several authors used initially different values of up to 10 Earth-masses (e.g., Raymond et al., 2006).. The orbital elements were initially chosen randomly like in the first run. These computations were performed for a total set of 50 different initial conditions.

• •

Again the initial conditions were chosen randomly but we took into account that the estimated masses by the observers are only spectroscopically determined. Thus three additional computations were performed for massless bodies but with more massive four giant planets (the most probable masses given by the most probable inclinations, Voigt, 2012).

• •

In a forth run we checked the outer region between planets d and e by populating it with massless bodies (see Figure 3).

The water content of the planetesimals was calculated to be initially between 5% to 10% according to the mass of the body and the distance to the central star.

We note that the number of bodies diminished with integration time because of

• •

collisions amongst them resulting in formation of bigger planets (for the computations with massive planetesimals – perfect merging was assumed555the parameters of each collision (impact velocities, impact angle) have been stored for possible follow up studies with SPH (Smoothed Particle Hydrodynamics).)

• •

collisions with one of the gas giants

• •

collisions with the host star

• •

escape from the system after a close encounter with one of the giant planets (as escape criterion for a body we fixed: e ¿ 0.99).

4 The ’Belt’ of Massless Bodies

The reason for the choice of a close investigation of the region between planets c and planet d (we call it the belt) is the following: we are interested in the formation of possible terrestrial planets in the HZ of HD 141399.

In a recent work of Agnew et al. (2018) amongst other multi-planetary systems with gas giants the dynamical stability of HD 141399 is proven, using best-fit planetary and stellar parameters from the NASA Exoplanet Archive. They show that a potential terrestrial planet (m $\sim$ 1 $\mathrm{m_{Earth}}$) can inhabit a stable orbit inside the HZ of this particular system for the integration period of $10^{7}$ yrs.

In order to have a first overview of the dynamics of orbits inside this belt we compute 8 ’packages’ of 500 massless bodies each. All the orbits were integrated for 1 Myr.

4.1 MMR Induced Structures in the Belt

The structure visible in Figures 2 and 4 looks very similar to the Kirkwood gaps in the main belt of asteroids in our Solar System. In Figure 2 we marked the MMRs.

When we compare the results of an integration of massless bodies between the planets d and e in Figure 3 with the region between Jupiter and Saturn in our Solar System there is a big difference evident. In the Solar System there are no stable families of asteroids with small eccentricities between $5.2\,au<a<9.52\,au$. In HD 141399 we found accumulations of bodies in certain resonances. Even in the equilateral Lagrange points of planet d and planet e many planetesimals are captured and form Trojan families.

The extension of computing bodies outside of planet d has been undertaken for the sake of completeness of the dynamical study of HD 141599. It becomes relevant as one recognizes the Trojans around planet d (see Figure 3). This was one of the reasons to explore the Trojan region of planet d in more detail (see Appendix). Note that around Saturn there is no a Trojan cloud (Dvorak et al., 2008, 2014). Nevertheless, it would be worth to do a more detailed investigation.

4.2 Larger Masses

The following Section deals with a check of stability and formation of families of massless bodies in the system of HD 141399 assuming larger masses for the four gas giants. We performed these computations with three different initial conditions (compare Section 3 and Figure 4).

For these computations we assumed the most probable inclination of the orbital plane of the giant planets – the statistical mean value. We assume equally distributed inclinations ($\overline{\sin^{3}(i)}=0.56$, Voigt, 2012) which leads to an increasing factor of 0.59 and gives $\mathrm{m_{planet\,b}}$ =  8.035, $\mathrm{m_{planet\,c}}$ = 2.378, $\mathrm{m_{planet\,d}}$ = 2.127 and $\mathrm{m_{planet\,e}}$ = 1.15 (masses in $\mathrm{M_{Jupiter}}$). Because of the small eccentricities, the gas giants are stable for Gyrs in this configuration, too. For the distribution of massless bodies in the HZ one can observe that gaps and groups appear as in our previous computations (see Figure 4). A comparison of Figures 2 and 4 shows similarities of the gaps and groups (areas of accumulation of bodies) in the distribution of bodies after an integration time of 1 Myr. These gaps and groups result from the acting of resonances (MMRs) with the giant planets.

5 Some Examples of New Terrestrial Planets

As an example of the formation process taking place in the system of HD 141399 we picked out one simulation for discussion. Semi-major axis versus mass of the bodies of this particular simulation are shown in Figure 5. One can observe that after 0.1 Myr most of the original Moon-sized objects are gone either due to collisions (with one of the gas giants or the host star) or they were thrown out of the system after close encounter(s) with one of the gas giants. A close encounter can lead to large eccentricity (e $>$ 0.7) and subsequent ejection from the system. In the presented simulation two massive planets were formed – one with the mass of Mars and one with approximately 0.5 Earth-masses. They stay in stable orbits between planet c and planet d for 1 Myr. An extension of the integration time to 10 Myr in this particular case shows that they obtain stable orbits over the whole integration time.

In Figure 6 mass versus semi-major axis of the remaining bodies after 0.05 Myr, 0.1 Myr, and 1 Myr are shown. We remark that after 0.05 Myr (red points) there are still many – now protoplanets – present in the system. Additionally still some of the originally Moon-sized bodies remain. We observe the formation of a planet of 20 Moon-masses at $\sim$3.8 au, another one with 16 Moon-masses at $\sim$4.4 au, and one with 10 Moon-masses at $\sim$4.3 au within the first 0.05 Myr.

Due to their eccentricities, which allow close encounters with planet d and subsequent collision or ejection (not checked in detail here) the number of remaining bodies after 0.1 Myr (green points) is small.

Note that Figure 5 and Figure 6 show results of the same simulation while Figure 7 depicts results from another run.

6 The Growing from Planetesimals to Protoplanets

There are two ways of looking at the the growth of planets:

• a)

follow one particular body during its evolution (monitoring each collision in order to watch the growth to the final protoplanet/ planet)

• b)

follow the development of the whole belt of bodies (checking the number of remaining bodies and the growth after several time steps)

We follow the growth of a planetesimal to the moment of reaching the mass of a terrestrial planet. A special case is shown in Figure 7 where after 9 collisions (perfect merging was assumed) a terrestrial planet is formed. The body considered has a mass of 25 $\mathrm{m_{Moon}}$ (blue lines in Figure 7). The second most massive one has only 20 $\mathrm{m_{Moon}}$ (green lines in Figure 7).

The growth process can be divided into different steps: In the beginning small bodies collide and form larger and larger objects until they reach a mass of several Moon-masses. The collision of two bigger bodies leads to the formation of a terrestrial planet like object. Yet, there are still many small bodies around which are ’added’ to the big body via collisions.

6.1 Example of the Evolution of the Belt

We show one out of 50 integrations and follow the 500 Moon-sized bodies over time in order to depict the evolution of the belt. We show the initial conditions and snapshots after 1 kyr, 10 kyr, and 1 Myr in Figure 8.

One can observe that the eccentricities grow quite fast when we compare in Figure 8 the initial distribution to the distribution semi-major axis versus eccentricity after 1000 years (left and right upper panels in Figure 8). Bodies of the size of several times the Moon are formed during such a short time. After 10 kyrs we can see that already some Mars-sized protoplanets were formed; like in the former plot one recognizes also that the majority of bodies is now confined primarily between 0.9 au and 1.5 au (lower left panel in Figure 8).These larger bodies disappear during the continued integration of the equations of motion up to 1 Myrs; we can explain it by encounters with the gas giants which throw them far away from the central star to large semi-major axes and eccentricities and may also lead to escapes. But even collisions with the four large planets let them ’disappear’. The left over after 1 Myr consists now of only two Mars-sized planets and two Moon-like ones (lower right panel in Figure 8).

7 Main Results for the Formation of Terrestrial Planets

The results from the integrations are summarized in the following plots which show (note that the water content will be discussed in Section 7.1)

• •

Figure 9: We show the resulting planets formed out of 50 x 500 planetesimals distributed between 0.7 au and 2.1 au (between planet c and d). Only two of the planets with mass comparable to the mass of Earth were formed via collisions and consecutive merging (these planets gained 0.9 Earth-masses). In addition many Mars-sized planets are formed. One can observe that most of the more massive planets inhabit low eccentric orbits. The biggest bodies have eccentricities e $<$ 0.1, which assures that they are in relatively stable orbits. Nevertheless, there are some bodies with mass of $\sim$ 0.5 Earth-masses with eccentricities up to $e\sim 0.3$ which might experience close encounters or even collision with one of the gas giants. Eventually they might fall into the star or even be ejected from the system.

• •

Figure 10: The distribution of bodies at the end of the integration (0.5 Myr) shows that most of the bodies remain between the planets c and d. An accumulation of objects can be observed at certain semi-major axes. The majority of the formed planets with up to 0.9 Earth-masses reside in the area $0.8\,\textrm{au}<a<1.6\,\textrm{au}$. Some objects accumulate around the gas giant d. They have been captured into the 1:1 MMR as Trojan planets. Furthermore there appears to be a more or less empty region inside of planet d. Clearly, empty regions analogue to the ’Kirkwood gaps’ in our Solar System are seen. The groups where many larger planets are formed, are observable at 1.15 au (by far the biggest one), 1.35 au, and 1.5 au.

• •

Figure 11: The distribution of bodies at the end of the integration shows that most of them stay in the plane of the gas giants. In particular the more massive ones show only small values for inclination (i $<$ $10^{\circ}$). In some simulations bodies were scattered to retrograde orbits during the formation process.

• •

Figure 12: The plot shows, that some of the bodies captured as Trojans of planet d have orbits with inclinations of up to $80^{\circ}$. Some are even on retrograde orbits. Nevertheless, the majority of bodies stays in the plane of the gas giants with $i<20^{\circ}$.

• •

Figure 13: Statistical summary of all the remaining bodies after the end of the integration (included are all 50 runs). One can see that most of the planetesimals remained small. Nevertheless, some (to be exact: 2) managed to gain a mass of up to 0.9 Earth-masses.

7.1 The Role of Water Content

At the beginning of the integration each of the approximately Moon-sized bodies in the belt between the planets c and d is assigned with a water mass fraction of 5 – 10 % depending on the distance from the star. We assume that after each collision the bodies perfectly merge. Figure 9 and 10 show the summary of all 50 runs. Depicted are the mass and water content of the bodies after the end of the integration time. One can see that, obviously, the more massive ones collected more water (in units of Earth-masses). Nevertheless, the water content in percentages is the same for small and massive bodies.

SPH (smooth particle hydrodynamics) computations of collisions between bodies show that they suffer from essential loss of water depending on the impact velocity, the impact angle etc. (Burger et al., 2018; Maindl et al., 2013, 2014; Schäfer et al., 2016). All these parameters were not yet included in our simulations. As we assume that no mass is lost during a collision we clearly overestimate the water mass fraction present in each planet at the end of the simulation.

The more realistic models of collisions with the aid of SPH-codes need time consuming simulations for every single collision and should be treated in combination with a n-body code.

Nevertheless, our results show that the planets formed in between the planets c and d could foster life as they remain on stable low eccentric orbits in the HZ around HD 141399 with possibly water on their surfaces.

8 Conclusions

The goal of the presented work was to investigate the probability of the existence and formation of terrestrial planets in an extrasolar planetary system with several gas giants. We found a very interesting example: The system HD 141399, where a K0V star hosts 4 giant planets in distances to the central star between 0.4 au and 5.2 au. The small eccentricities of the planets guarantee stable orbits up to Gyrs. This is true, even when it is taken into account that the determined masses via spectral analysis from observations with large telescopes are only minimum values.

In order to answer the question about existence and formation of terrestrial planets in the system of HD 141399 we assumed a leftover belt of hundreds of small bodies between the two planets c and d orbiting at 0.7 au respectively 2.1 au. Furthermore we assume that the formation of the four gas giants has been finished (see Section 3). The region between the planets c and d includes the habitable zone (HZ) where water on a rocky planet can be present in liquid form. This might give the planet the chance to develop and sustain life on its surface.

The work was done in several steps. First of all we studied the evolution of orbits of massless bodies in the HZ which yielded a structure similar to the ’Kirkwood Gaps’ found in the main belt of asteroids in our Solar System – all associated to mean motion resonances (MMRs) with the gas giants. These simulations have been done with the primarily determined masses of the four gas giants. In addition the integrations were repeated with masses of the four giants estimated for the statistical mean value of inclination. There is only a small difference in the final distribution of the remaining massless bodies inside the belt (comparison of Figure 2 and 4).

As a second step we integrated Moon-sized bodies for 500 different initial conditions between the planets c and d with randomly chosen orbital elements, small initial eccentricities ($e\sim 0.1$), and small initial inclinations ($i\sim 1^{\circ}$). We provide numerical integrations of the equations of motion in the Newtonian framework up to the formation of terrestrial planets. We carefully studied the collisions between these bodies under the assumption of complete merging. Being aware of this rough estimation we started SPH (Smooth Particle Hydrodynamics) computations which will provide more realistic results for the collisions in a future project. We performed 50 different runs with different initial conditions for the 500 planetary embryos. In several plots we show the respective results – e.g after how many collisions such a terrestrial planet is formed, how fast the formation happens in different stages of the evolution, how many planets form, and in which distance to the host star most of the planets are formed.

We found that there is primarily one region around 1.15 au where many terrestrial planets formed within approximately 1 Myr (most of the massless bodies survived in this region in previous computations (see Figure 10). Also in two other regions around 1.35 au and 1.55 au planets formed. The latter one corresponds to the so-called Hilda region (named after the Hilda asteroids in the main belt of asteroids in our Solar System) in the 2:3 MMR with planet d. Those Hildas are less numerous and less massive.

Concerning formation computations with larger masses for the gas giants we found that it is far more difficult to form such planets, because relatively soon the objects undergo perturbations by the giant planets leading to large eccentricities and subsequent close encounters and/or ejections. Our test computations with a system with larger masses of the primary bodies has shown that already after some thousand years only half of the original planetesimals survived. Although this number is comparable to the computation with the nominal masses (see Figure 8) there is a big difference insofar that for the latter test 95 % of the original bodies left the central region because of close encounters respectively collisions with the gas giants and also the central star and not even one bigger body was formed within this time (see Section 4.2).

The special case of a MMR with one of the giant planets – the capture of a small body in 1:1 resonance – a so called Trojan, was investigated separately. One can observe a few aggregated bodies up to the size of Mars in the Trojan region around planet d. In the results of our integrations of fictitious massless Trojans for 1 Myr we found stable regions around planet d up to an inclination of the Trojans of about $80^{\circ}$ (compare Figure 12).

In our computations we also included the water content of the bodies and followed the collisional process of accumulation of more and more mass up to the final formation of a habitable terrestrial planet. In conclusion we can say that terrestrial like planets may exist in the extrasolar planetary system of HD 141399 and it might even harbor life.

Acknowledgements

This research is supported by the Austrian Science Fund (FWF) trough grant S11603-N16 (B.L. and R.D. and L.Z.). The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC). This work has been supported by the National Natural Science Foundation of China (NSFC, Grants No. 11473016, and No. 11333002) Special thanks to T. I. Maindl for critically reading the manuscript and his help with Gnuplot.

Appendix A Trojans Around Planet d

Inspired by the presented results on the accumulation of bodies around planet d (shown in Figure 3) and the importance of possible Trojan planets in extrasolar planetary systems (Dvorak et al., 2008; Schwarz et al., 2007; Dvorak et al., 2004; Laughlin & Chambers, 2002) we checked the stability of the Lagrange point $L_{4}$: for a grid of initial conditions 100 by 60 (semi-major axis versus initial inclination) we integrated the respective orbits of fictitious massless Trojans for 1 Myr (Figure 14). There is a big stable region around a = 2.1 au visible between 1.95 au $<$ a $<$ 2.175 au (for the plane case). The horseshoe-like contour of the stable region extends up to an inclination of $\mathrm{i}=57^{\circ}$ . The stability criterion was the libration in the angle $\lambda$ around the equilibrium point $60^{\circ}$ ahead in the planets orbit on one side (not shown here) and the largest eccentricity achieved during the integration. The almost black region shows that the libration point itself is very stable; more to the edge the eccentricity reaches values up to 0.3 before they escape (yellow region). All the structures visible are comparable to those found in other studies of Trojans in the Solar system (e.g., Dvorak et al., 2012).

A.1 Dynamical Map

Beside the maximum eccentricity map, we also construct the dynamical map of the Trojan region using the spectral number (SN) as the indicator of regularity of orbits. Concretely, the SN is defined to be the number of peaks over a given threshold in the power spectrum of a specific variable that can represent some basic properties of the orbital motion. The method of SN, which could reflect the long term stability within a relatively short integration time, was introduced by Michtchenko & Ferraz-Mello (1995) and has been successfully applied (e.g., Michtchenko et al., 2002; Zhou et al., 2009, 2011).

In this study, we calculate the power spectrum of the resonant angle of each test particle in the Trojan region, $\sigma=\lambda-\lambda_{d}$ where $\lambda$ and $\lambda_{d}$ are the mean longitudes of the test particle and the planet HD141399 d. After some test runs, we set our integration time to be $4.6\times 10^{6}$ yr to cover at least two complete periods of all the secular variation terms in the system.

With the same initial conditions as the one we set for obtaining the maximum eccentricity ($e_{\max}$) map, we perform two runs of calculations and two dynamical maps are obtained, as shown in Figure 15. The SN in logarithm is indicated by colours. The red and blue imply the most irregular (chaotic) and most regular (stable) motions, respectively, while the intermediates are shown in the colour bar. The empty area (white in colour) are those orbits that cannot survive our integration time. In the left panel, all planets in this system are nearly coplanar with small inclination of $0.1^{\circ}$, and in the right panel we show the dynamical map for the non-coplanar system, in which the planets in the HD141399 system from inside out are assumed to have the same inclinations and ascending nodes as Jovian planets in our Solar System, i.e. Jupiter, Saturn, Uranus and Neptune, respectively. We will call the latter one the Jovian-analogue system hereinafter.

The SN dynamical map reveals richer details than the maximum eccentricity map. As shown in Figure 14, the $\mathrm{e_{\max}}$ does not change with the inclination but has a monotonic increase from the resonance centre ($a=2.09$ au) outwards to both boundaries along the semi-major axis. Contrarily, in Figure 15 evident arc-structures can be seen, indicating that the orbital stability depends sensitively both on a and i. Introducing chaos into the Trojan region, these distinct structures shrink and shred the stable region. Particularly, the left and right branches of the broad orange arcs corresponding to the most irregular motion are connected to each other at $(a_{0},i_{0})\sim(2.083\textrm{ au },38^{\circ})$, making the stable region at high inclination disconnected from the one at low inclination.

To our experience, the arc structures in the left panel of Figure 15 arise due to some secondary resonances, where the frequency of $\sigma$ is commensurable to some secular frequencies of the system. In addition, two V-shaped thin stripes can be clearly seen, which most probably is related to the nodal secular resonances. Since the system is nearly coplanar (planets’ inclinations are all $0.1^{\circ}$), these structures may be enhanced if the planets’ inclinations are higher.

In the right panel, we show the dynamical map in the Jovian-analogue system, in which the inclinations of the planets are taken as the same as the Jovian planets (higher than $0.1^{\circ}$). Apparently, the V-shaped stripes are enhanced as expected. Compared to the left panel, we find that the irregularity has developed considerably. Although the Jovian inclination is still quite small in fact, the stable region in high inclination $\mathrm{i}_{0}>25^{\circ}$ almost disappears now. But in the low-inclination region, the stability increases. As a result of the disappearance of those cyan-yellow arc structures corresponding to irregular motion, the blue stable region extends both to lower inclination and to a wider range in semi-major axis. It is also worth to note that the resonance centre in this Jovian-analogue system shifts a little from the near-coplanar system, because here all the initial conditions are osculating orbital elements.

A.2 Resonance Map

Generally, a resonance happens when two (or multiple) frequencies in a dynamical system match each other or they are commensurable by some simple integers. So, it is possible to find the resonances if the frequencies of the system are known. Following the similar method and procedure as in Zhou et al. (2009, 2018), we numerically determine the frequencies of four planets in this system, which are listed in Table 3.

From Table 3, we notice that the mean motions of Planet c and Planet d are very close to a ratio of 5:1. In fact, in our numerical simulations, we find that this near resonance plays a significant role in the dynamics of Trojans around Planet d, as we will see below. The frequency of this near resonance, defined as the frequency of the angle $(5\lambda_{d}-\lambda_{c})$ and denoted by $f_{5:1}$, is calculated as $f_{5:1}=8.0\times 10^{-2}2\pi/$yr.

To find out the resonances that sculpt the features in the dynamical map, we calculate the proper frequencies of the test Trojans on the $(\mathrm{a}_{0},\mathrm{i}_{0})$-plane as well. For each test Trojan, the frequencies of mean motion $f$, apsidal precession $g$ and nodal precession $s$ are calculated and finally the empirical expressions of these frequencies as functions of the initial conditions $\mathrm{a}_{0},\mathrm{i}_{0}$ are obtained. With all the frequencies in hand, we can determine the resonances on the $(\mathrm{a}_{0},\mathrm{i}_{0})$-plane by carefully checking the relationships among these frequencies. The most significant resonances are found and their locations are plotted on the $(\mathrm{a}_{0},\mathrm{i}_{0})$-plane.

The most evident structure on the dynamical map is the family of grand ‘arcs’ crossing the stable region. We find that these arcs are basically related to the secondary resonances between the libration of $\sigma$ (critical angle of the 1:1 MMR) and the near resonance angle $(5\lambda_{d}-\lambda_{c})$. As illustrated in Figure 16, the line $3f_{\sigma}-f_{5:1}=0$ just crosses the unstable arc, and some other lines labelled by C1 to C4 that stand for the following equations, emplace the similar structures nearby

[TABLE]

The secular resonance that has the strongest dynamical effects is the $\nu_{2}$ resonance. As explained above, the $\nu_{2}$ resonance indicates the situation where the Trojan’s periastron processes in the same rate as Planet c’s periastron, $\dot{\varpi}=\dot{\varpi}_{2}$. The $\nu_{2}$ drives the Trojan’s eccentricity to vary in a large range and destablizes its orbit, thus carves the upper boundary of the stable region. The Kozai mechanism, occurring just above the $\nu_{2}$ curve in Figure 16, may also contribute to the instability in this neighbourhood.

Other conspicuous resonances include the nodal secular resonances $\nu_{12},\nu_{13}$ and $\nu_{14}$. Very clearly, the $\nu_{13},\nu_{14}$ curves just sit exactly inside the V-shaped stripes. The effects of the lowest order secular resonances, particularly the $\nu_{12}$ and $\nu_{14}$ are greatly enhanced when the inclinations of the planets increase to the Jovian-analogue system (right panel in Figure 15), while those high-order secular resonances, e.g. the one related to $2s+g_{2}-3g_{3}=0$, fade away in a relative sense.