New Dynamical State and Habitability of the HD 45364 Planetary System

Over the past decade, there has been a dramatic increase in the number of exoplanets detected, largely attributable to the Kepler (Borucki et al. 2010) and Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) missions. The huge inventory of exoplanets we see today shows a vast diversity in terms of planetary intrinsic properties (Konacki & Wolszczan 2003; Kipping & Bakos 2011; Barclay et al. 2013; Masuda 2014; Gaudi et al. 2017), orbital characteristics (Jones et al. 2006; Smith et al. 2018; Zhang et al. 2021), and system architectures (Muirhead et al. 2012; Schwamb et al. 2013; Cabrera et al. 2014; Gillon et al. 2017; Feng et al. 2017; Bohn et al. 2020). From the over 5000 confirmed exoplanets within over 3000 planetary systems so far, there are over 800 systems that host multiple planets, many of which host planets in mean motion resonance (MMR) chains. MMR can occur when planets within a system have orbital periods of near integer ratios with each other. Planets that are within an MMR chain often exert significant gravitational influence on each other and exchange angular momentum periodically (Raymond et al. 2008; Petrovich et al. 2013; Goldreich & Schlichting 2014). Of particular interest are giant planets exhibiting MMR orbital configurations at locations within the snow line. According to standard planetary formation theories, giant planets are believed to form at distances far away from the star where the temperature is low enough for condensation of volatile compounds, such as water ice (Sasselov & Lecar 2000). Therefore, the giant planets within an MMR chain at locations within the snow line are likely to have initially formed beyond the snow line, then gradually migrated inwards through convergent migration while being embedded within the protoplanetary disk. Such interacting planets undergoing migration thus provide important constraints on planetary formation, planet migration history, system evolution, and architecture (Correia et al. 2009; Rosenthal et al. 2019). Additionally, the presence of giant planets within the system could greatly influence the habitability of terrestrial planets (Georgakarakos et al. 2018; Sánchez et al. 2018; Kane & Blunt 2019; Kane et al. 2020; Bailey & Fabrycky 2022). Giant planets passing through the system's habitable zone (HZ) during migration provide important clues to the viability of forming additional terrestrial planets in the HZ around the host star and, if so, the dynamical stability of these low-mass planets. Systems such as these with giant planets interacting with each other in MMR, particularly in the HZ of the star, offer rare opportunities to simultaneously study the intriguing dynamical history of the planets as well as the habitability of the system.

The star HD 45364 was discovered to host two gas giant planets by Correia et al. (2009, C09 hereafter) using radial velocity (RV) data from the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph (Pepe et al. 2000). The two planets were found to be trapped in a 3:2 MMR orbiting inside the HZ of the system, and it was the first exoplanetary system found to host planets exhibiting such an orbital configuration. Orbital parameters published by C09 predicted strong planet–planet interactions, and the interesting dynamics of the system prompted several studies into the possible planetary formation and migration route that could explain the observed MMR. Notably, Rein et al. (2010, R10 hereafter) proposed a formation scenario where the outer planet underwent a very rapid inward type III migration (Masset & Papaloizou 2003) that allowed it to quickly pass through the very stable 2:1 resonance and finally settle in a 3:2 resonance with the inner planet. This scenario was verified by R10 through both N-body and hydrodynamical simulations under the assumptions that the surface density of the protoplanetary disk was about five times higher than the minimum solar nebula (MMSN; Hayashi 1981) at 1 au and that there was no active mass accretion during the process of migration. That is, the two planets achieved their final orbital states with giant planet masses before the start of the migration. Instead of the unconventional type III migration model that R10 proposed, Correa-Otto et al. (2013, C13 hereafter) suggested that a more traditional model with the combination of type I and type II migration is possible to explain the observed 3:2 MMR. In this scenario, both planets were assumed to be in their planet embryo states initially with masses of 0.1 M_⊕ and started the migration under the type I regime within the disk with surface density equivalent to MMSN at 5 au. The two planets were allowed to accrete mass while migrating inward under the two-stage giant planet formation model by Alibert et al. (2005), where the migration effects are included during the formation. The planet pair was able to pass through 2:1 and other MMRs due to their small masses and minimal interaction with each other. This scenario by C13 suggests that the two planets were able to reach near the 3:2 MMR at the end of the first stage of planetary formation before reaching their critical masses and the onset of runaway growth. Both planets would then be massive enough to continue migration toward the observed orbital configuration slowly under type II migration.

One of the intriguing aspects about the different proposed migration scenarios of the HD 45364 planet pair is that they predicted different orbital configuration outcomes. The outer planet type III migration to a 3:2 resonance capture with the inner planet modeled in R10 resulted in much-lower-eccentricity, near-circular orbits for both planets and, as a result, different libration patterns from that reported in C09. With the planet pair undergoing type I and II migration in C13, two types of simulations were conducted, K3 and K100, with different ratios (3 and 100) of the e-folding times of semimajor axis decay and eccentricity damping for both planets assumed for their type II migration stage. Interestingly, the K100 simulation reproduced a similar result to that in R10 while K3 simulation was consistent with the configuration originally derived in C09. It is worth noting that, at the time of these publications, RV data were limited to an observation baseline of ∼4 yr; perhaps too short to detect the longer-term effect of the system dynamics and distinguish between the different proposed models. However, the previous works have stressed that the effort to resolve the orbital configuration and migration model degeneracies could greatly benefit from further RV observations.

In this work, we provide a clearer picture of the orbital configuration of the two planets as well as the system architecture using the latest RV data set. This data set, as will be described in Section 3, doubles the amount of RV data and quadruples the observation baseline compared to the previously published data. Such a data expansion allows us to revisit this fascinating system after over ten years to refine the orbital parameters for both planets and attempt to detect the predicted strong planet–planet interaction through both Keplerian and dynamical modeling of the RV data. We also provide constraints on the system's orbital inclination using only RV information by fitting the data with dynamical models at different inclinations. An RV survey completeness analysis is conducted, and we explore the potential of additional undetected planets within the system. The new orbital configuration has particular implications for the dynamical state of the system, which we also study in detail. We are especially interested in the possibility of detecting additional low-mass terrestrial planets in the HZ of the system. As mentioned earlier, the presence of two giant planets in the HZ as a result of inward migration may have a huge impact on the formation of other, smaller planets in the HZ. The confirmation or exclusion of possible additional low-mass planets in the system could provide a clue to the habitability prospects of this system and will be useful for target selection processes related to future missions that search for potentially habitable terrestrial planets.

The paper is organized as follows: In Section 2, we describe the HD 45364 system architecture and the extent of the HZ. We present the revised RV solutions of the system with the new RV data collected using both Keplerian and dynamical models in Section 3, including orbital inclination constraints and RV completeness calculations. Section 4 discusses an analysis of the system's new dynamical state near MMR, and the planet–planet interaction within that resonance regime. Prospects for terrestrial bodies within the HZ of the system, including planets and exomoons, are investigated in Section 5. We finally discuss the implications of the results from this work and provide concluding remarks in Section 6.

HD 45364 is a nearby star with a V-band magnitude of 8.08 (Correia et al. 2009), and is located at a distance of 34.35 ± 0.04 pc (Gaia Collaboration et al. 2018). The star has a spectral type K0V, an effective temperature of 5466${}_{-32}^{+59}$ K, and a luminosity of 0.637 ± 0.002 L_⊙ (Gaia Collaboration et al. 2018). The mass and radius were estimated to be 0.82 ± 0.05 M_⊙ and ${0.89}_{-0.02}^{+0.01}$ R_⊙, respectively (Correia et al. 2009; Gaia Collaboration et al. 2018). The system was reported by C09 to host two giant planets through the HARPS RV search with minimum masses of 0.1872 M_J and 0.6579 M_J for the (b) and (c) planets, respectively. The orbits of the (b) and (c) planets have semimajor axes of 0.6813 and 0.8972 au, respectively, and eccentricities of 0.168 ± 0.019 and 0.097 ± 0.012, respectively. Error estimates unfortunately were not provided by C09 for both minimum masses and semimajor axes. These reported orbits are summarized in Table 2 and depicted in a top-down view of the system, shown in the left panel of Figure 1, showing the proximity of the two planetary orbits that potentially render the system unstable. However, dynamical simulations performed by C09 found that the two planets, with orbital periods 226.93 ± 0.37 and 342.85 ± 0.28 days, reside within a 3:2 MMR such that the planets never pass within 0.37 au of each other. The system configuration was thus determined to be stable over Gyr timescales, making it the first detected exoplanetary system to exhibit such a MMR and an interesting case study of planetary formation and orbital dynamical scenarios.

Zoom In Zoom Out Reset image size

In addition to the MMR configuration, the locations of the two planets are intriguing to the study of planetary habitability in the HZ (Kopparapu et al. 2013, 2014; Kane et al. 2016a), and also for the study of potentially habitable exomoons (Heller et al. 2014; Hill et al. 2018). We calculated the boundaries of the conservative HZ (CHZ) and optimistic HZ (OHZ) in the system following the definitions described by Kopparapu et al. (2013, 2014). We adopted the aforementioned stellar luminosity and effective temperature values. The inner and outer boundaries for the CHZ were determined to be 0.77 and 1.38 au, respectively. For the OHZ, the inner and outer boundaries were found to be 0.61 and 1.45 au, respectively. The extent of these HZ boundaries is depicted in Figure 1, along with the C09 planetary orbits described above. The left panel of Figure 1 demonstrates that the orbit of planet (c) lies completely within the bounds of the CHZ, while planet (b) spends 86% of its orbital period within the CHZ and the remaining time either in the inner OHZ region or interior to the HZ entirely. Although giant planets are likely poor considerations for habitable conditions, additional terrestrial exoplanets within the HZ or exomoons orbiting the known giant planets could potentially have feasible environments for life (Heller et al. 2014). The right panel of Figure 1 shows a depiction of the revised orbital solution for both planets from this work, and is described in detail in Section 3.

The original discovery of the two-planet system by C09 utilized 58 RVs acquired using the HARPS spectrograph (Pepe et al. 2000) on the ESO 3.6 m telescope at La Silla, Chile in the span of 1583 days from 2003 December to 2008 April. The best-fit orbital solution at that time (see Table 2 and Figure 1) places the two planets in a 3:2 MMR and slightly eccentric orbits with the potential for close encounters. Since then, HD 45364 had been continuously monitored by the HARPS team for an extended period of time until 2017 September, extending the RV baseline to over 5000 days with 114 data points in total. The full HARPS data set for this target was later rereduced with the new Spectrum Radial Velocity Analyser pipeline (Zechmeister et al. 2018) and had several systematics corrected for including nightly zero-point variations and intra-night drifts. The newly reduced HARPS data achieved slightly higher precision and was published as part of the HARPS RV database (Trifonov et al. 2020). The fiber upgrade to the HARPS instrument in 2015 May (Lo Curto et al. 2015) modified the instrument profile and created a vertical offset between the pre- and post-upgraded RV data (Trifonov et al. 2020). Because of the upgrade, we treated the pre- and post-upgraded HARPS RV as two separate instruments during the fitting process. In addition to the HARPS data set, we obtained 7 RVs taken using High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on Keck I at Maunakea, Hawaii from 2009 December to 2021 September (Table 1). The data reduction follows the similar procedure as in Howard et al. (2010). In total, 121 RV data points were collected for this system with a baseline of about 6500 days. The improved RV precision as well as the extended observational baseline, 5000 days longer than the original data set, allows the opportunity to revisit and provide a more accurate orbital solution to the two planets, as suggested by the original paper.

3.1. Keplerian Model

Given the extended baseline and additional observations obtained since the original publication, we conducted a Fourier analysis for the entire data set to check for possible additional periodicities or linear trends that may be indicative of additional companions in the system. We used RVSearch (Rosenthal et al. 2021), an iterative planet searching tool that uses the change in the Bayesian Information Criterion (ΔBIC) between models as a measure of the goodness of the fit. We set a period grid from 2 to 10,000 days and allowed linear trends and curvatures to be incorporated in the period search. Signals were considered significant if they peak above the 0.1% false-alarm-probability (FAP) level. Two signals with periods similar to those of the originally reported planets were present and were picked up by the periodogram as expected, but with no long-term linear or quadratic trend detected.

We then sampled model posteriors using the RV modeling toolkit RadVel (Fulton et al. 2018). We used five parameters: orbital period P, time of inferior conjunction T_c , $\sqrt{e}$ cosω, $\sqrt{e}$ sinω, where e and ω are orbital eccentricity and argument of periastron, respectively, and RV semiamplitude K as the fitting basis to better sample the near-circular orbits and to minimize the bias toward the higher orbital eccentricities during the fitting process (Lucy & Sweeney 1971). A Markov Chain Monte Carlo (MCMC) search of the parameters' posterior space was carried out using parameters returned from RVSearch as the initial guess. All parameters including instrumental offset and jitter terms were allowed to explore freely with no additional constraint other than forcing the positive RV semiamplitude (K) values, eccentricities (e) to be between 0 and 1, and jitter terms between 0 and 10 m s⁻¹. The chain converged relatively quickly, and Figure 2 shows the best-fit model along with the data used. The model returned approximately similar result as reported by the original paper except, in this case, the orbital eccentricities for both planets were favored to be near circular, as supported by a model comparison, in contrary to the mild eccentricities of e ∼ 0.17 and e ∼ 0.1 for planets (b) and (c), respectively, that were originally reported in C09. In addition, the derived mass of the outer planet is slightly smaller than the previously reported value. As a consequence of more circular orbits, the fitted arguments of periastron (ω) and time of periastron (T_p ) values were not as well constrained as those reported in C09. The fit produces a log-likelihood (ln ${ \mathcal L }$) = −240.67, rms = 1.76 m s⁻¹, and ${\chi }_{\nu }^{2}$ = 1.15. We present the Keplerian model fit results including the solution with 16%, 50%, and 84% quantiles as well as that for the maximum likelihood in Table 2. We also show the results from C09 in the same table for comparison.

Zoom In Zoom Out Reset image size

Table 2. Orbital Parameters of the Two Planets in HD 45364 from C09 and This Work

Source	Planet	P (days)	e	ω (deg)	Tp (BJD)	K (m s−1)	a (au)	Mp (MJ)
C09	b	226.93 ± 0.37	0.168 ± 0.019	162.6 ± 6.3	2,453,464 ± 4	7.22 ± 0.14	0.6813	0.187
	c	342.85 ± 0.28	0.097 ± 0.012	7.4 ± 4.3	2,453,407 ± 4	21.92 ± 0.43	0.8972	0.658
This Work (K) Quantiles	b	${227.93}_{-0.17}^{+0.16}$	${0.035}_{-0.024}^{+0.033}$	${232}_{-80}^{+63}$	$2,453,{471}_{-35}^{+53}$	${7.27}_{-0.24}^{+0.25}$	0.684 ± 0.014	0.191 ± 0.010
	c	343.84 ± 0.18	${0.025}_{-0.015}^{+0.016}$	${123}_{-31}^{+36}$	$2,453,{200}_{-25}^{+180}$	18.13 ± 0.25	${0.899}_{-0.019}^{+0.018}$	${0.548}_{-0.024}^{+0.023}$
This Work (K) Max Likelihood	b	227.95	0.051	232	2,453,462	7.32	0.662	0.180
	c	343.93	0.033	120	2,453,186	18.17	0.871	0.513
This Work (D) Quantiles	b	${225.79}_{-0.76}^{+0.81}$	0.067 ± 0.016	${92}_{-25}^{+22}$	$2,453,{374}_{-16}^{+15}$	7.23 ± 0.24	${0.6793}_{-0.0015}^{+0.0016}$	${0.1893}_{-0.0063}^{+0.0062}$
	c	${345.43}_{-0.57}^{+0.54}$	${0.019}_{-0.010}^{+0.011}$	${244}_{-68}^{+54}$	$2,453,{306}_{-65}^{+52}$	${18.15}_{-0.24}^{+0.25}$	0.9020 ± 0.0010	${0.5490}_{-0.0074}^{+0.0075}$
This Work (D) a Max Likelihood	b	225.34	0.070	92	2,453,375	7.26	0.6784	0.1897
	c	345.76	0.010	276	2,453,336	18.17	0.9026	0.5497

Note. Letters K and D in the Source column denote Keplerian and dynamical models from this work, respectively. ω values reported here are those of the planets, not those of the star as usually reported in RV discoveries. T_p values for C09 were derived using λ values reported in Table 2 of C09. For the parameters from the dynamical model in this work, e and ω were derived from fitted esinϖ and ecosϖ (ω was used instead of ϖ for the Keplerian model since RV data contain no information regarding Ω; Ω was assumed to be 0 for the dynamical fit); T_p , K, and a were derived from λ, M_p , and P, respectively. All models in this table assume edge-on orbital inclinations.

^a The best-fit model we employ for the rest of this work.

Download table as:  ASCIITypeset image

3.2. Dynamical Model

The pure Keplerian model assumes planets orbiting around their host stars are unperturbed by other external forces other than that from the central star. However, there were indications from previous works that, due to the close proximity of the two planets with high minimum masses in MMR, strong planet–planet interactions are expected (Correia et al. 2009). RVs of planets undergoing strong interactions may gradually deviate from the true Keplerian model, and such deviations can be detected in the RV data set if the observational baseline is long enough. Parameters derived using a Keplerian model would thus become an inaccurate representation of the system of interest. To account for the potential presence of planet–planet interaction between the (b) and (c) planets, a full dynamical model is needed to integrate the planets along their orbits and obtain the induced motion on the host star at each step.

We made use of the code RVStab (Rosenthal et al. 2019) as the backbone of our analysis to dynamically fit the RV data by integrating the entire system with the inclusion of planet–planet interaction using an N-body integrator REBOUND (Rein & Liu 2012), and to carry out a posterior parameter search with MCMC exploration. We used IAS15 (Rein & Spiegel 2015), a nonsymplectic integrator with adaptive time stepping, and chose the coordinate system to be with reference to the central star. We assumed coplanarity for the system and tested the case where the system inclination is edge on (i = 90°) first. The dynamical model for the system was set up using a set of five parameters of the orbital period P, planetary mass M_p , planet's mean longitude λ (referenced to the first data point in our data set), Lagrange orbital elements esinϖ, and ecosϖ for each planetary body with the addition of three jitter terms for the three instruments. ϖ is the longitude of periastron, which is defined as ϖ = ω+Ω with Ω being longitude of ascending node, which is assumed to be zero in our model. We used the results from the Keplerian fit (Table 2) and introduced a small Gaussian deviation from the parameters as an initial guess to initialize the MCMC walkers. At each step, the set of orbital parameters was fitted to the RV data, and a ln ${ \mathcal L }$ value was calculated for each instrument using Equation (1) below and was summed for all instruments as a measure of the goodness of the fit with the parameters at the current step:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L } & = & -\displaystyle \sum _{i}^{N}\displaystyle \frac{{({\upsilon }_{o,i}-{\upsilon }_{\mathrm{off}}-{\upsilon }_{m,i})}^{2}}{2({\sigma }_{i}^{2}+{\sigma }_{\mathrm{jit}}^{2})}\\ & & -\displaystyle \sum _{i}^{N}\mathrm{ln}(\sqrt{2\pi ({\sigma }_{i}^{2}+{\sigma }_{\mathrm{jit}}^{2})})+\mathrm{ln}(\sqrt{2\pi {\sigma }_{z}^{2}})\end{array}\end{eqnarray} \tag{ 1 }$$

where υ_o,i , υ_off, and υ_m,i are the ith observed velocity, instrumental velocity offset, and model velocity, respectively. σ_i is the ith measurement error associated with υ_o,i ; σ_jit is the estimated jitter for the instrument; and ${\sigma }_{z}^{2}$ is defined as follows:

$$\begin{eqnarray}&&{\sigma }_{z}^{2}={\left(\sum _{i}^{N}\displaystyle \frac{1}{{\sigma }_{i}^{2}+{\sigma }_{\mathrm{jit}}^{2}}\right)}^{-1}.\end{eqnarray} \tag{ 2 }$$

The walkers were able to explore freely under uniform priors between the specified lower and upper bounds of fitted parameters: between 100 and 500 days for P; 0.0157 and 5 Jupiter masses for M_p ; 0° and 360° for λ; −1 and 1 for esinϖ and ecosϖ; and between 0 and 10 for the jitters. The chains were considered converged once the length of the chain was at least 50 times the autocorrelation time, and the change in autocorrelation time between steps was smaller than 1% for five consecutive steps. The chain was able to converge successfully, and the fit is shown in Figure 3. The orbital solution with 16%, 50%, and 84% quantiles as well as maximum likelihood values is reported in Table 2. The orbital parameters from the dynamical fit are largely consistent with the one from the Keplerian fit, where we see much more circular orbits for both the (b) and (c) planets compared to the solution from C09. The mass of the outer planet appears to be smaller as well. However, we do notice the difference between the fitted orbital periods of the two planets from the dynamical model and the Keplerian model, possibly indicating the presence of planet–planet interaction. The difference in orbital periods derived from the Keplerian and dynamical models is 3σ significant for both planets if quoting the period uncertainty values from the dynamical model. If using the uncertainties from the Keplerian model, the period difference amounts to 13σ and 9σ for planets (b) and (c), respectively. The dynamical model with edge-on inclination yielded ln ${ \mathcal L }$ = −232.49, rms = 1.75 m s⁻¹, and ${\chi }_{\nu }^{2}$ = 1.12. The dynamical model appears to be a slightly better fit to the RV data than the Keplerian model from Section 3.1 based on ln ${ \mathcal L }$.

Zoom In Zoom Out Reset image size

In order to detect the presence of strong planet–planet interaction claimed by previous works, we computed the time series RV difference (ΔRV) between the two models, and the result is shown is Figure 4. It was predicted, based on the limited data at the time of discovery in Figure 2 of C09, that the interaction between the two planets would start to become detectable in the RV data around 2015, and the ΔRV between Keplerian and dynamical models would then gradually increase to an amplitude of around 15 m s⁻¹. However, as can be seen in the bottom panel of Figure 4, the ΔRV between the two models is fluctuating on a similar level over time with low variations from sub-1 to 2.5 m s⁻¹. This value is consistent with the estimated stellar jitter of ∼2.34 m s⁻¹ that we calculated for the host star HD 45364 following the methodology in Isaacson & Fischer (2010). Even at the time of our last data point, which is taken by HIRES near the end of 2021 September, ΔRV does not become any more significant. Therefore, we conclude that, contrary to the previous prediction, interaction between the two planets is not as strong as expected and is in fact negligible within the sensitivity of the current RV data set. The nondetection of planet–planet interaction is possibly due to the more circular orbits and a lower mass for the outer planet derived from the new RV data in this work.

Zoom In Zoom Out Reset image size

3.3. Inclination Constraint

Typical RVs are insensitive to orbital inclination information so that systems observed with RV data have only minimum masses M_p sini, instead of true masses of planets reported. However, when planets are massive enough and are orbiting close to each other, the potential dynamical interaction between the planets could provide a rare opportunity for the system inclination to be constrained. This can be done by introducing orbital inclination as an additional variable and running the dynamical fit under different inclination assumptions. Since varying the assumed inclination changes the true planetary-mass values, planets would start to strongly interact with each other at some inclination angles. The interaction at those angles would become so much, such that the system would be in unstable states or the resulting RV contribution from the planets would deviate from the observed RV data points, allowing us to rule out inclination cases when these situations happen. In some situations where planets are interacting strongly, the orbital inclination of the system could be identified with high certainty using RV data alone. In others where no planet–planet interaction can be detected, the inclination can still be constrained to within a certain range.

In this case, we first varied the system's inclination, and the masses of two planets accordingly, from edge on to face on to test the range of inclinations where the system would become unstable. We assumed coplanarity for the planets and used the orbital solutions that assumed an edge-on case derived from the dynamical work as presented in Section 3.2 and Table 2. Once again, we employed REBOUND to carry out the dynamical integration using a symplectic integrator WHFast (Rein & Tamayo 2015) with an integration time step of ∼5.5 days, 1/40 of the (b) planet's orbital period, and half of the recommended step size from Duncan et al. (1998) to ensure proper orbital sampling. The system was integrated for 10⁷ yr for each inclination angle, and the system stability breaks down at angles smaller than 9°.

Next, we carried out similar dynamical fits to the RV data as presented earlier, except this time the dynamical models at different inclination values up to the instability angle are tested against the RV data, rather than just the edge-on case. Inclination angles were varied every 5°, and the ln ${ \mathcal L }$ value was recorded for each model. In the end, BIC values were calculated for all models as a comparison of goodness of fit. The result is shown in Figure 5. Inclination angles leading to system instability are shaded in red, and the angles with ΔBIC > 2 compared to the best-fit model are shaded in orange and are disfavored. Although there is no distinct peak that would allow us to pinpoint the actual inclination value, the inclination of the system could still be constrained to ≥40° given the models with inclinations between 40° and 90° are indistinguishable according to ΔBIC. It now makes sense that we adopt the edge-on case as the orbital solution for this system since it shares similar statistical significance as the other models, but without taking the inclination as one of the free parameters. We therefore employ the maximum likelihood result from the dynamical model as the best-fit orbital parameters for the rest of this work.

Zoom In Zoom Out Reset image size

3.4. Stellar Activity

3.5. RV Search Completeness

In addition to planetary discoveries, knowing the detection limit or sensitivity of the data is equally important because they entail what other bodies could be detected but were not by the data, or what other objects could potentially be lurking in the system below the detection threshold (Wittenmyer et al. 2020; Li et al. 2021; Rosenthal et al. 2021). Such studies not only extend the knowledge of exoplanet population and system architecture but also serve as a guideline to target the selection of future space missions regarding the search for habitable terrestrial planets in the HZ for example, where the likelihood of the presence of such planets could determine the observing priorities of candidate systems. Here, we carried out an injection–recovery test within the M_p sini versus a parameter space using RVSearch to determine the search completeness of our RV data. We injected 3000 fictitious planets into the current data set and ran the iterative planet search within RVSearch in the same way as described in Section 3.1 to see whether the injected planet could be recovered. The fictitious planet was given an orbit with parameters drawn from log-uniform distributions for the planet's orbital period and minimum mass. The eccentricity was drawn from a Beta distribution following the result from Kipping (2013). The result of the injection–recovery test was then used to compute the RV search completeness contour, which is shown in Figure 6. The two known giant planets that have been the centerpiece of this paper are clearly above the detection threshold of our RV data set. It is interesting to point out from this figure that the data cannot rule out the possibility of additional very long orbital period companions since the current data is insensitive to objects, if there are any, at the very large separations from the host star. Thus a dim, low-mass stellar object that is gravitationally bound to the system is still a possible explanation to the observed trend in two of the activity indicators as we mentioned in the last subsection. In addition, the data suggests that additional gas giants with 1 Jupiter mass could be detected with 90% detection probability under the current data sensitivity for separations up to ∼7 au. However, the nondetection of such planets within our data set indicates the two we have now may be the only giant planets in the system within this separation range.

Zoom In Zoom Out Reset image size

The new RV models, either Keplerian or dynamical, both yield slightly different orbital parameter values than those previously reported, especially for the orbital periods, eccentricity, and planetary masses of the two giant planets. Combined with the fact that we do not see the predicted strong planet–planet interaction, the new picture leads us to wonder if the system could be in a different dynamical state. One way to study the dynamical behavior of the system is by calculating the apsidal trajectory and plotting the behavior in a polar plot of e_b e_c sin Δϖ against e_b e_c cos Δϖ, where Δϖ = ϖ_b − ϖ_c (Figure 7). Typically, if the trajectory encompasses the origin, the system is deemed to be circulating. If not, the system is librating, and the type of libration depends on the location relative to the origin and the shape of the trajectory. These two modes of apsidal behavior are separated by a boundary called "separatrix," for which the system could have a behavior of both libration and circulation if the trajectory comes close to the origin of the polar plot (Barnes & Greenberg 2006a, 2006b). In Figure 7, the polar plot includes the result for 10,000 yr with each plotted point being 1 yr apart. Interestingly, the apsidal trajectory is showing two prominent oval shapes, one to the left of the origin, and the other one slightly to the right but encompassing the origin. These ovals may represent the two modes of apsidal behavior, antialigned libration (left oval) and circulation (right oval). However, the plot shows neither libration nor circulation seems to be in a stable state. Based on the number of data points that lie in between the two modes, the system does not appear to be confined in either one of the modes. This presents a possible but seemingly uncommon scenario for the system's dynamical state in which neither libration nor circulation is dominating the overall behavior, and the system is constantly seeking the balance between the two, thus forming a cone shape in the polar plot where the system is moving back to forth between the two modes. The feature seen in Figure 7 is still present if we extend the simulation duration to 10⁷ yr.

Zoom In Zoom Out Reset image size

The interesting dynamics seen in Figure 7 can be further investigated by calculating the evolution of resonance angles. For an interacting planet pair near a 3:2 MMR, the resonance angles can be calculated as follows:

$$\begin{eqnarray}&&{\phi }_{b}=3{\lambda }_{c}-2{\lambda }_{b}-{\varpi }_{b};\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\phi }}_{c}=3{\lambda }_{c}-2{\lambda }_{b}-{\varpi }_{c}.\end{eqnarray} \tag{ 4 }$$

We ran a short integration of 1000 yr and recorded the resonance angles every 0.1 yr. The evolution of both angles and ${\rm{\Delta }}\varpi$ is shown in Figure 8. Indeed, the top panel of Figure 8 shows that φ_b is librating around 0° whereas φ_c seems to be stuck in between circulation and libration, with the libration pattern around 180°. Upon closer examination, φ_c exhibits several similar but not the same libration and circulation cycles every ∼250 yr. Within each quasi-period, φ_c goes through mixed periods of libration in between each circulation. Such "nodding" behavior where the resonance angle goes through phases of libration and circulation repeatedly over time was studied in detail in Ketchum et al. (2013) and appears to be related to the perturber's orbital eccentricity variation. This is reflected in the eccentricity variation of both planets (b) and (c) in Figure 9 where the duration and time step are the same as those for Figure 8. The (b) and (c) planets go through eccentricity variation of ∼0–0.07 and ∼0–0.03, respectively, and we verified these variations last for simulations of 10⁷ yr. Looking closely at both eccentricities in Figure 9, both planets are going through similar quasi-periodic cycles of eccentricity variation every ∼250 yr, suggesting the orbital shape and orientation change could be the cause of circulation we see in Figures 7 and 8. Since circulations happen when the planet pair's orbits are close to alignment (Δϖ = 0) according to Figure 8, this suggests φ_c goes through circulation when the interaction between the two planets is at a minimum. Due to the eccentricity variation, φ_c cycles back into libration when the orbits are not aligned and planets undergo interaction with each other for a period of time, until the next orbit alignment and circulation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The libration patterns we see here in Figure 8 when the system goes through the libration phase are largely similar to those in previous works (see for example Figure 9 in R10), except that both resonance angles in this case are librating with a much larger amplitude, especially in the case of φ_c where libration becomes so large that sometimes it goes into the circulation regime. Both larger libration amplitudes and two apsidal modes, one after another, indicate that, instead of being within a 3:2 MMR, the planet pair is moving in and out of the resonance repeatedly, in a quasi-resonance state. Such type of motion between libration and circulation near the "separatrix" concerning the HD 45364 system was first suggested in C13 regarding the predicted orbital solution under the type III migration mechanism proposed by R10 and also one of the solutions under the type I, type II combination mechanism by C13, although not directly observed in their models. Given the similarity of the system's dynamical behavior between our model and some of those from R10 and C13, we therefore conclude that the system is indeed in a quasi-resonance dynamical state, and any migration mechanisms predicted near-circular orbit solutions should be preferred as the formation pathway to the configuration that we see today in HD 45364.

The previous orbital configuration with two planets strongly interacting on mildly eccentric orbits may present a challenge to allowing potentially habitable locations within the system. However, the revised orbital solution for the system with near-circular orbits (right panel of Figure 1) and, more importantly, the minimal gravitational interaction between the two giant planets, as shown in previous sections, open up such a possibility. Here, we briefly discuss potential habitability aspects of the system, including the prospect of terrestrial planets within the HZ and moons harbored by the known giant planets.

5.1. Additional Planets

Giant planets observed to orbit within the HZ of the system are believed to have formed farther out beyond the snow line, then gradually migrated inwards through planet–disk interaction. Such migration history may present a significant challenge to the formation of terrestrial planets within the HZ because of the constant gravitational perturbation to the nearby material in the disk, inhibiting the material buildup that would have eventually formed other planets. However, past works have shown that Earth-mass terrestrial planets could indeed survive the giant planet migration and start the formation process within the HZ during or after a giant planet's passage of the HZ (Raymond 2006; Raymond et al. 2006). Here, we conducted a dynamical simulation to test for the viable locations where another Earth-mass terrestrial planet could be present.

We made use of the REBOUND package to carry out the dynamical simulation. For the two known giant planets in the system, we assumed an edge-on orbital configuration derived from our dynamical model presented earlier (Table 2). We then injected an Earth-mass planet in a circular orbit into the existing system at 1000 different locations, one location at a time, evenly spaced within the range of the OHZ from 0.578 to 1.376 au (see Section 2). We used the WHFast integrator with a fixed time step of ∼8 days, roughly 1/20 the orbital period of a test particle orbiting at the inner edge of the OHZ. The system with the hypothetical Earth-mass planet was integrated for 10⁷ yr to test for system stability and the survival of the terrestrial planet, similar to the methodology described by Kane et al. (2015), Kane (2015), Kane et al. (2022). Results were recorded every 100 yr, and the simulation for each test location would stop if any one of the three planets were ejected from the system. After the simulation, we computed the survival rate of the test planet at each test location, defined as the percentage amount of time the test planet would remain in the system throughout the entire integration duration. Any locations with the Earth-mass planet having 100% survival rate are considered dynamically viable locations. The simulation results are shown in Figure 10. As expected, most of the locations within the HZ are prohibitive to the presence of an Earth-mass terrestrial planet due to the influence of the two giant planets within the HZ, as indicated by the red and orange vertical lines in Figure 10. In the vicinity of the two known planets, the injected planet had near-zero survival rates from the inner edge of the HZ up to ∼1.2 au, meaning that the small planet was ejected momentarily after the start of the integration. However, as the test planet was moved farther away from the influence of the two giants and toward the outer edge of the HZ, the dynamically stable area opened up for the terrestrial planet, and it was able to achieve full survival at the end of the simulation at locations near the outer edge of the HZ, some even within the outer boundaries of the CHZ. These results show that, provided the migration of the giant planets allow subsequent formation of terrestrial planets within the system, those planets may retain orbits within the outer edge of the HZ.

Zoom In Zoom Out Reset image size

5.2. Implications for Exomoons

The giant planets orbiting within the HZ of the system are now increasingly scrutinized for their prospects as hosts of potential habitability. These planets themselves may not have the environment suitable for habitability studies, but the potential exomoons orbiting around these giant planets might do thanks to the combined global energy budget from the star, the planet, and tidal energy (Tinney et al. 2011; Heller et al. 2014; Hill et al. 2018). Significant studies have been carried out regarding the orbital stability of moons (Barnes & O'Brien 2002; Kipping 2009; Gong et al. 2013; Kane 2017; Hong et al. 2018; Quarles et al. 2020; Dobos et al. 2021), especially those around planets that have experienced migration (Spalding et al. 2016). For HD 45364, a detailed investigation of the formation and dynamics of planet–moon systems is beyond the scope of this work. We can, however, provide an estimate of the Hill radius of both giant planets and the critical orbits of potential exomoons, provided by the following equations:

$$\begin{eqnarray}&&{r}_{H}\approx a(1-e){\left(\displaystyle \frac{{M}_{p}}{3{M}_{\star }}\right)}^{1/3};\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{a}_{\mathrm{crit}}\approx {{fr}}_{H}.\end{eqnarray} \tag{ 6 }$$

According to Equation (5) above, and using parameters from the best dynamical fit from Table 2, the sizes of the Hill radii for planets (b) and (c) are r_H,b ≈ 0.026 au and r_H,c ≈ 0.053 au, respectively. However, the calculation of Hill radius is only an estimation that does not take into account any other external perturbations (Kipping 2009). Following Barnes & O'Brien (2002) and Kipping (2009), we adopt f = 1/3 and use a_crit from Equation (6) as a conservative outer boundary for the possible locations of exomoons around exoplanets. For the (b) planet, we estimated a_crit,b ≈ 0.009 au, and for the (c) planet, a_crit,c ≈ 0.018 au. Despite the best estimate, this upper bound estimate should not be mistaken as the limit below which an exomoon, if there exists one, is guaranteed to be stable around the planet it orbits. As pointed out by Spalding et al. (2016), the survival of exomoons from type II migration of giant planets depends on the location the planets were migrating from. However, if there are multiple moons orbiting around a giant planet in an MMR configuration, such moon destruction through "evection resonance" may be quenched due to moon–moon interactions (Spalding et al. 2016).

The orbital solution published by C09 in 2009 resulted in two giant planets orbiting in the HZ of the system with mildly eccentric orbits. The planet pair was found to reside in a 3:2 MMR that prevents the planets from experiencing close encounters with each other. Several works have since then attempted to reconstruct the formation pathway to the observed configuration at the time through different planetary migration models. R10 proposed a type III migration model and predicted near-circular orbits for both planets rather than eccentric ones. C13 instead proposed a combination of type I and type II migration scenario and predicted both eccentric and circular orbits are possible under different assumptions. All these works have stressed the importance of future RVs in helping resolve the model degeneracy.

In this work, we revisited the system 13 yr later using a substantially improved RV data set. We presented a new orbital solution for the two planets within the HD 45364 system using new RV data from HARPS and HIRES. The latest RVs allowed us to extend the observational baseline to almost 18 yr, almost four times the baseline used by the original publication by C09 13 yr ago. We conducted a reanalysis of the RV data through both Keplerian and dynamical models. Both models point to near-circular orbits for the two planets instead of having mild eccentricities. In addition, small changes were observed for the orbital periods and planetary masses. As a consequence, the orbits are now more separated, and interactions between the planets may not be as strong as expected. Indeed, by comparing our Keplerian and dynamical models, we did not observe the predicted large amplitude of ΔRV over time between the two models. The ΔRV we determined was consistent with the noise level from both the star and the instruments. We thus conclude that there are negligible planet–planet interactions within the sensitivity of our RV data set. We also attempted to determine the orbital inclination of the system through dynamical fitting. However, due to much weaker interactions between the planet pair, orbital inclination can only be constrained to ≥40°, and an edge-on case is so far preferred.

The orbital dynamics was studied in detail given the change in the orbital configuration of the system and planetary masses. We integrated the system for 10⁷ yr and tracked the evolution of orbital elements for both planets. By tracking the apsidal trajectory and resonance angle evolution, we found that the system was not quite in a 3:2 MMR. Rather, the system exhibits both libration and circulation, and thus appears to be moving in and out of the 3:2 MMR. The libration and circulation pattern is in a quasi-periodic cycle that is consistent with a similar cycle observed in both planets' eccentricity variations, suggesting the observed new dynamical behavior of the system is related to the periodic change in the magnitude of the interaction between the two planets. Therefore, HD 45364, the first exoplanetary system discovered to host a 3:2 MMR, is actually in a quasi-resonance state, rather than being in a true resonance state. The result of our new RV models and dynamical analysis successfully confirmed predictions from previous works. More importantly, our result suggests that any migration models that predict near-circular orbits for both planets should be the preferred migration scenarios for the HD 45364 system. This work demonstrates the importance of continued RV monitoring of a system and how it could impact our understanding of the system's dynamics.

Our orbital solution indicates the (b) planet resides within the inner OHZ, the (c) planet resides within the CHZ, and the orbits of the two planets no longer experience signs of a close encounter (Figure 1). The new system architecture opens up possible habitable locations in the HZ, and we investigated the habitability prospects of the system. An Earth-mass terrestrial planet was injected into the HZ and tested for dynamically viable locations near the two giant planets. The simulation indicated such a planet is indeed possible near the outer edge of the HZ, and detection of such a planet would be challenging because the expected RV amplitude would be well below the current RV sensitivity (e.g., K ≈ 8.5 cm s⁻¹ at 1.35 au). Exomoons around the two giant planets present another possibility for potentially habitable locations within the system, based on the simple estimate of the Hill radius of the two planets. If future observations could detect such objects, the new piece of data could provide further insight into the rich dynamical history that HD 45364 presents. This system therefore holds great value for future space-based missions that search for potentially habitable bodies in nearby exoplanetary systems.

There remains significant further work regarding the HD 45364 system. Transit observations could reveal more about this system. Based on the stellar radius value from Section 2 and orbital parameters of the two planets from Table 2, the transit probability is only 0.61% for the (b) planet and 0.46% for the (c) planet, assuming randomly oriented orbits. Assuming planetary radii of R_b = 0.81 R_J and R_c = 1.22 R_J, derived from the best-fit mass values using the mass–radius relationship of Chen & Kipping (2017), the two planets would have transit depths of 0.87% and 1.98%, respectively. HD 45364 was observed by TESS in sectors 6 and 7 from 2,458,468 to 2,458,516 Barycentric Julian Date (BJD). According to our best dynamical model fit, the inferior conjunction times are T_c,b = $2,458,{453}_{-24}^{+23}$ BJD and T_c,c = $2,458,{487}_{-66}^{+53}$ BJD for planets (b) and (c), respectively. The estimated conjunction times of both planets sit very close to the TESS observation window. We checked TESS data and no transit was detected for impact parameter b < 1. The transit nondetection suggests that either TESS observations had a near miss of the event in time or the system is not edge on. Given the large transit depth of both planets, ground-based follow-up observations could help confirm or rule out the transit scenarios. Future transit observations could potentially disclose more information about this system, such as planetary radius, atmosphere, transit timing variations, and even additional terrestrial planet and exomoon detections, if the system is edge on. In addition, continued RV monitoring by the precision RV facilities could further help refine the orbital solutions, track dynamics of the system, and investigate the source of the mysterious long-term trend seen in some of the activity indicators to determine whether it is of companion or activity origin. HD 45364 is without a doubt one of the most interesting systems in many areas of exoplanet study.

The authors wish to thank Hanno Rein for discussion and suggestions regarding the use of the REBOUND package. The authors also would like to thank Rory Barnes for conversation regarding the system's dynamics. We thank the anonymous referee for the valuable comments that greatly improve the presentation of this work. P.D. acknowledges support from a 51 Pegasi b Postdoctoral Fellowship from the Heising-Simons Foundation. Dynamical simulations in this paper made use of the REBOUND code, which is freely available at http://github.com/hannorein/rebound. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work has also made use of The Habitable Zone Gallery at hzgallery.org (Kane & Gelino 2012). The results reported herein benefited from collaborations and/or information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate.

Software: GLS (Zechmeister & Kürster 2009), RadVel (Fulton et al. 2018), REBOUND (Rein & Liu 2012), RVSearch (Rosenthal et al. 2021), RVStab (Rosenthal et al. 2019).